Symmetric, Skew-Symmetric and Triangular Matrices

Among the special arrays implemented in GeometricOptimizers SymmetricMatrix, SkewSymMatrix, UpperTriangular and LowerTriangular are the most common ones and similar implementations can also be found in other libraries; LinearAlgebra.jl has an implementation of a symmetric matrix called Symmetric for example. The versions of these matrices in GeometricOptimizers are however more memory efficient as they only store as many parameters as are necessary, i.e. $n(n+1)/2$ for the symmetric matrix and $n(n-1)/2$ for the other three. In addition, GeometricMachineLearning implements matrix and tensor multiplication for these matrices so that they work in parallel on GPU; see Tensors there. We here give an overview of elementary custom matrices that are implemented in GeometricOptimizers. More involved matrices are the so-called global tangent spaces.

Custom Matrices

GeometricOptimizers has two types of triangular matrices. The first one is UpperTriangular:

\[U = \begin{pmatrix} 0 & a_{12} & \cdots & a_{1n} \\ 0 & \ddots & & a_{2n} \\ \vdots & \ddots & \ddots & \vdots \\ 0 & \cdots & 0 & 0 \end{pmatrix}.\]

And the second one is LowerTriangular:

\[L = \begin{pmatrix} 0 & 0 & \cdots & 0 \\ a_{21} & \ddots & & \vdots \\ \vdots & \ddots & \ddots & \vdots \\ a_{n1} & \cdots & a_{n(n-1)} & 0 \end{pmatrix}.\]

adjoint swaps between the two: L' is an UpperTriangular and U' is a LowerTriangular. That swap is built around the same storage vector rather than a copy, so parent(L') === parent(L) holds and writing into L' also writes into L. Reusing the storage transposes without conjugating, so that swap is bound to a real element type; a complex one falls through to LinearAlgebra's lazy Adjoint, which conjugates and does not alias.

An instance of SkewSymMatrix can be written as $A = L - L^T$ or $A = U^T - U$:

\[A = \begin{pmatrix} 0 & - a_{21} & \cdots & - a_{n1} \\ a_{21} & \ddots & & \vdots \\ \vdots & \ddots & \ddots & \vdots \\ a_{n1} & \cdots & a_{n(n-1)} & 0 \end{pmatrix}.\]

And lastly a SymmetricMatrix:

\[B = \begin{pmatrix} a_{11} & a_{21} & \cdots & a_{n1} \\ a_{21} & \ddots & & \vdots \\ \vdots & \ddots & \ddots & \vdots \\ a_{n1} & \cdots & a_{n(n-1)} & a_{nn} \end{pmatrix}.\]

Note that any matrix $M\in\mathbb{R}^{n\times{}n}$ can be written

\[M = \frac{1}{2}(M - M^T) + \frac{1}{2}(M + M^T),\]

where the first part of this matrix is skew-symmetric and the second part is symmetric. This is also how the constructors for SkewSymMatrix and SymmetricMatrix are designed. Consider an arbitrary matrix:

M = [1; 2; 3;; 4; 5; 6;; 7; 8; 9]

3×3 Matrix{Int64}:
 1  4  7
 2  5  8
 3  6  9

Calling SkewSymMatrix on $M$ is equivalent to doing $M \to \frac{1}{2}(M - M^T)$:

A = SkewSymMatrix(M)

3×3 SkewSymMatrix{Float64, Vector{Float64}}:
  0.0   1.0  2.0
 -1.0   0.0  1.0
 -2.0  -1.0  0.0

And calling SymmetricMatrix on $M$ is equivalent to doing $M \to \frac{1}{2}(M + M^T)$:

B = SymmetricMatrix(M)

3×3 SymmetricMatrix{Float64, Vector{Float64}}:
 1.0  3.0  5.0
 3.0  5.0  7.0
 5.0  7.0  9.0

We can further confirm the identity above:

M  ≈ A + B

true

Note that for LowerTriangular and UpperTriangular no projection step is involved, which means that if we start with a matrix of type AbstractMatrix{Int64} we will end up with a matrix that is also of type AbstractMatrix{Int64}. The type changes however when we call SkewSymMatrix and SymmetricMatrix:

(typeof(A) <: AbstractMatrix{Int64}, typeof(B) <: AbstractMatrix{Int64})

(false, false)

For the triangular matrices:

U = UpperTriangular(M)
L = LowerTriangular(M)
(typeof(U) <: AbstractMatrix{Int64}, typeof(L) <: AbstractMatrix{Int64})

(true, true)

How are Special Matrices Stored?

The following image demonstrates how a skew-symmetric matrix is stored in GeometricOptimizers:

So what is stored internally is a vector of size $n(n-1)/2$ for the skew-symmetric matrix and the triangular matrices, and a vector of size $n(n+1)/2$ for the symmetric matrix.

Sample Random Matrices

We can sample a random skew-symmetric matrix:

A = rand(SkewSymMatrix, 3)

3×3 SkewSymMatrix{Float64, Vector{Float64}}:
 0.0       -0.521214  -0.586807
 0.521214   0.0       -0.890879
 0.586807   0.890879   0.0

and then access the vector:

A.S

3-element Vector{Float64}:
 0.521213795535383
 0.5868067574533484
 0.8908786980927811

This is equivalent to sampling a vector and then assigning a matrix^[1]:

S = rand(3 * (3 - 1) ÷ 2)
SkewSymMatrix(S, 3)

3×3 SkewSymMatrix{Float64, Vector{Float64}}:
 0.0       -0.521214  -0.586807
 0.521214   0.0       -0.890879
 0.586807   0.890879   0.0

These special matrices are what the layers of GeometricMachineLearning are parametrized by: SympNets, the volume-preserving transformer and the linear symplectic transformer all use one or more of them. That package also batches them over the third axis of a tensor, with mat_tensor_mul and tensor_mat_mul; see Tensors there.

Where a sampled array lands, and what it holds

rand and zeros here come in four shapes, by whether the call names a backend and whether it names an element type. The shape decides both answers, and it decides them the same way for the structured matrices above, for the manifolds and for the horizontal lifts wherever a type offers the shape at all. Not every type offers all four — the second shape is the manifolds' alone, and a shape a type does not offer is a MethodError rather than a different answer:

The third and fourth shapes place on the host without saying so, and that is deliberate. They mirror Base, where zeros(Float32, 3) is a host array and nothing about the call suggests otherwise; these types present as AbstractMatrix, so zeros(SkewSymMatrix{Float32}, n) should read as the Array case does. A host placement also cannot quietly corrupt a device computation: mixing one with a device array throws at the first arithmetic — + and add! both reach the storage arrays and fail there — so the loud failure already gives the guarantee that making these shapes take a backend would buy. copyto! and assign! are the deliberate exception, because they are the transfer: moving a host-built structured matrix onto a device is what they exist for.

The second shape is the only one where the package decides something the caller did not, which is why the choice is a stated rule rather than a literal: Float64 on the host, Float32 on a device. See default_eltype for why each value is what it is, and note that a backend being able to hold a Float64 is not one of the reasons.

The first shape has a rule of its own, in the other direction: an element type the caller names and the backend declares it cannot hold is refused, not narrowed, so rand(MetalBackend(), SkewSymMatrix{Float64}, n) and rand(MetalBackend(), StiefelManifold{Float64}, N, n) are both an ArgumentError. A narrowed result would have a different type from the one asked for, which is exactly what naming the element type rules out. Every allocator of this shape carries the check; KernelAbstractions.supports_float64 is what it asks, so it can only fire where a backend's own package has declared the limitation.

None of this reaches an allocation the package makes for itself. zero, similar, _zero and _similar all take an instance, so the backend and the element type both come from the argument and there is nothing to default — which is what every optimizer cache and every state allocates through. A parameter set on a device stays there.

Why the storage matters here

Because these types keep only their free parameters, they are also what an optimizer has to be able to update — and the generic array methods cannot do it: three of the four have no setindex! for an elementwise operation to broadcast through, similar has to preserve the type rather than widen to a dense Matrix, and $n(n\pm1)/2$ numbers do not reshape back to $n \times n$. See VectorStorageMatrix for the methods that make them usable as optimizer parameters.

The same storage is what a flat parameter vector and a saved file have to hold, for the same reason: $n(n\pm1)/2$ numbers are the whole content of one of these matrices, and the $n^2$ entries of the dense interface are neither the right length nor, for three of the four types, writable at all. NeuralNetworkParameters asks a leaf type for exactly that relation, through its freeparameters/rebuild pair, and loading it alongside this package brings in an extension that answers for all three families here — these matrices, the manifolds, and the horizontal lifts. Flattening, differentiating and saving a parameter set that contains them therefore needs no case per type in the package doing the training.

Library functions

AbstractTriangular, UpperTriangular, LowerTriangular, SkewSymMatrix, SymmetricMatrix and VectorStorageMatrix. Their docstrings are on the reference page, where every docstring in the package is rendered once; the names above link to them.