Orthogonal and Unitary matrices

 OrthogonalMatrices{n} = GeneralUnitaryMatrices{ℝ, n, AbsoluteDeterminantOneMatrixType}

The manifold of (real) orthogonal $n×n$ matrices $\mathrm{O}(n)$. They are precisely the $n×n$ real matrices $M$ that satisfy

\[ M^{T}M=MM^{T}= \mathrm{I}_n,\]

where $M^{T}$ is the transpose of $M$ and $\mathrm{I}_n$ is the $n×n$ identity matrix. Thus, their columns and rows represent $n$ pairwise orthonormal vectors. Such matrices $M$ have the property that $\lVert \det(M) \rVert = 1$.

This is the same manifold as the Stiefel(n, n, ℂ) manifold.

Constructor

OrthogonalMatrices(n)

const UnitaryMatrices{𝔽, T} = GeneralUnitaryMatrices{𝔽, T, AbsoluteDeterminantOneMatrixType}

The manifold $U(n,𝔽)$ of $n×n$ complex matrices (when 𝔽=ℂ) or quaternionic matrices (when 𝔽=ℍ) such that

\[ p^{\mathrm{H}}p = \mathrm{I}_n,\]

where $p^{\mathrm{H}}$ is the conjugate transpose of $p$ and $\mathrm{I}_n$ is the $n×n$ identity matrix. Such matrices p have a property that $\lVert \det(p) \rVert = 1$.

The tangent spaces are given by

\[ T_pU(n) \coloneqq \bigl\{ X \big| pY \text{ where } Y \text{ is skew symmetric, i. e. } Y = -Y^{\mathrm{H}} \bigr\}\]

But note that tangent vectors are represented in the Lie algebra, i.e. just using $Y$ in the representation above. If you prefer the representation as X you can use the Stiefel(n, n, ℂ) manifold.

Constructor

UnitaryMatrices(n, 𝔽::AbstractNumbers=ℂ)

see also OrthogonalMatrices for the real valued case.

rand(::Unitary; vector_at=nothing, σ::Real=1.0)

Generate a random point on the UnitaryMatrices manifold, if vector_at is nothing, by computing the QR decomposition of an $n×n$ matrix.

Generate a tangent vector at vector_at by projecting a normally distributed matrix onto the tangent space.

riemannian_Hessian(M::UnitaryMatrices, p, G, H, X)

The Riemannian Hessian can be computed by adopting Eq. (5.6) [Ngu23], so very similar to the complex Stiefel manifold. The only difference is, that here the tangent vectors are stored in the Lie algebra, i.e. the update direction is actually $pX$ instead of just $X$ (in Stiefel). and that means the inverse has to be applied to the (Euclidean) Hessian to map it into the Lie algebra.

Weingarten(M::UnitaryMatrices, p, X, V)

Compute the Weingarten map $\mathcal W_p$ at p on the Stiefel M with respect to the tangent vector $X \in T_p\mathcal M$ and the normal vector $V \in N_p\mathcal M$.

The formula is due to [AMT13] given by

\[\mathcal W_p(X,V) = -\frac{1}{2}p\bigl(V^{\mathrm{H}}X - X^\mathrm{H}V\bigr).\]

manifold_dimension(M::UnitaryMatrices{n,ℂ}) where {n}

Return the dimension of the manifold unitary matrices.

\[\dim_{\mathrm{U}(n)} = n^2.\]

manifold_dimension(M::UnitaryMatrices{<:Any,ℍ})

Return the dimension of the manifold unitary matrices.

\[\dim_{\mathrm{U}(n, ℍ)} = n(2n+1).\]

const SpecialUnitaryMatrices{T} = GeneralUnitaryMatrices{ℂ, T, DeterminantOneMatrixType}

The manifold $SU(n)$ of $n×n$ complex matrices such that

\[ p^{\mathrm{H}}p = \mathrm{I}_n \text{ and } \det(p) = 1,\]

where $p^{\mathrm{H}}$ is the conjugate transpose of $p$ and $\mathrm{I}_n$ is the $n×n$ identity matrix.

The tangent spaces are given by

\[ T_pU(n) \coloneqq \bigl\{ X \big| pY \text{ where } Y \text{ is skew symmetric and traceless, i. e. } Y = -Y^{\mathrm{H}} \text{ and } \operatorname{tr}(Y) = 0 \bigr\}\]

But note that tangent vectors are represented in the Lie algebra, i.e. just using $Y$ in the representation above.

Constructor

SpecialUnitaryMatrices(n; parameter::Symbol = :type)

see also Rotations for the real valued case.

manifold_volume(::SpecialUnitaryMatrices)

Volume of the manifold of complex general unitary matrices of determinant one. The formula reads [BST03]

\[\sqrt{n 2^{n-1}} π^{(n-1)(n+2)/2} \prod_{k=1}^{n-1}\frac{1}{k!}.\]

injectivity_radius(G::SpecialUnitaryMatrices)

Return the injectivity radius for general complex unitary matrix manifolds, where the determinant is $+1$, which is^[1]

\[ \operatorname{inj}_{\mathrm{SU}(n)} = π \sqrt{2}.\]

^[1] > For a derivation of the injectivity radius, see sethaxen.com/blog/2023/02/the-injectivity-radii-of-the-unitary-groups/.

manifold_dimension(M::SpecialUnitaryMatrices)

Return the dimension of the manifold of special unitary matrices.

\[\dim_{\mathrm{SU}(n)} = n^2-1.\]

AbsoluteDeterminantOneMatrixType <: AbstractMatrixType

A type to indicate that we require (orthogonal / unitary) matrices with normed determinant, i.e. that the absolute value of the determinant is 1.

AbstractMatrixType

A plain type to distinguish different types of matrices, for example DeterminantOneMatrixType and AbsoluteDeterminantOneMatrixType.

DeterminantOneMatrixType <: AbstractMatrixType

A type to indicate that we require special (orthogonal / unitary) matrices, i.e. of determinant 1.

GeneralUnitaryMatrices{𝔽, T, S<:AbstractMatrixType} <: AbstractDecoratorManifold

A common parametric type for matrices with a unitary property of size $n×n$ over the field $𝔽$ which additionally have the AbstractMatrixType, e.g. are DeterminantOneMatrixType.

exp(M::Rotations, p, X)
exp(M::OrthogonalMatrices, p, X)
exp(M::UnitaryMatrices, p, X)

Compute the exponential map, that is, since $X$ is represented in the Lie algebra,

\[exp_p(X) = p\mathrm{e}^X\]

For different sizes, like $n=2,3,4$, there are specialized implementations.

The algorithm used is a more numerically stable form of those proposed in [GX02] and [AR13].

log(M::Rotations, p, X)
log(M::OrthogonalMatrices, p, X)
log(M::UnitaryMatrices, p, X)

Compute the logarithmic map, that is, since the resulting $X$ is represented in the Lie algebra,

\[\log_p q = \log(p^{\mathrm{H}}q)\]

which is projected onto the skew symmetric matrices for numerical stability.

log(M::Rotations, p, q)

Compute the logarithmic map on the Rotations manifold M which is given by

\[\log_p q = \log(p^{\mathrm{T}}q)\]

where $\log$ denotes the matrix logarithm. For numerical stability, the result is projected onto the set of skew symmetric matrices.

For antipodal rotations the function returns deterministically one of the tangent vectors that point at q.

cos_angles_4d_rotation_matrix(R)

4D rotations can be described by two orthogonal planes that are unchanged by the action of the rotation (vectors within a plane rotate only within the plane). The cosines of the two angles $α,β$ of rotation about these planes may be obtained from the distinct real parts of the eigenvalues of the rotation matrix. This function computes these more efficiently by solving the system

\[\begin{aligned} \cos α + \cos β &= \frac{1}{2} \operatorname{tr}(R)\\ \cos α \cos β &= \frac{1}{8} \operatorname{tr}(R)^2 - \frac{1}{16} \operatorname{tr}((R - R^T)^2) - 1. \end{aligned}\]

By convention, the returned values are sorted in decreasing order. See also angles_4d_skew_sym_matrix.

manifold_volume(::GeneralUnitaryMatrices{ℝ,<:Any,AbsoluteDeterminantOneMatrixType})

Volume of the manifold of real orthogonal matrices of absolute determinant one. The formula reads [BST03]:

\[\begin{cases} \frac{2^{k}(2\pi)^{k^2}}{\prod_{s=1}^{k-1} (2s)!} & \text{ if } n = 2k \\ \frac{2^{k+1}(2\pi)^{k(k+1)}}{\prod_{s=1}^{k-1} (2s+1)!} & \text{ if } n = 2k+1 \end{cases}\]

manifold_volume(::GeneralUnitaryMatrices{ℝ,<:Any,DeterminantOneMatrixType})

Volume of the manifold of real orthogonal matrices of determinant one. The formula reads [BST03]:

\[\begin{cases} 2 & \text{ if } n = 0 \\ \frac{2^{k-1/2}(2\pi)^{k^2}}{\prod_{s=1}^{k-1} (2s)!} & \text{ if } n = 2k+2 \\ \frac{2^{k+1/2}(2\pi)^{k(k+1)}}{\prod_{s=1}^{k-1} (2s+1)!} & \text{ if } n = 2k+1 \end{cases}\]

It differs from the paper by a factor of sqrt(2) due to a different choice of normalization.

manifold_volume(::GeneralUnitaryMatrices{ℂ,<:Any,AbsoluteDeterminantOneMatrixType})

Volume of the manifold of complex general unitary matrices of absolute determinant one. The formula reads [BST03]

\[\sqrt{n 2^{n+1}} π^{n(n+1)/2} \prod_{k=1}^{n-1}\frac{1}{k!}.\]

volume_density(M::GeneralUnitaryMatrices{ℝ, TypeParameter{Tuple{2}}}, p, X)

Volume density on O(2)/SO(2) is equal to 1.

volume_density(M::GeneralUnitaryMatrices{ℝ, TypeParameter{Tuple{3}}}, p, X)

Compute the volume density on O(3)/SO(3). The formula reads [FdHDF19]

\[\frac{1-1\cos(\sqrt{2}\lVert X \rVert)}{\lVert X \rVert^2}.\]

volume_density(M::GeneralUnitaryMatrices{<:Any,ℝ}, p, X)

Compute volume density function of a sphere, i.e. determinant of the differential of exponential map exp(M, p, X). It is derived from Eq. (4.1) and Corollary 4.4 in [CLLD22]. See also Theorem 4.1 in [FdHDF19], (note that it uses a different convention).

check_point(M::Rotations, p; kwargs...)

Check whether p is a valid point on the UnitaryMatrices M, i.e. that $p$ has a determinant of absolute value one, i.e. that $p^{\mathrm{H}}p = \mathrm{I}_n$

The tolerance for the last test can be set using the kwargs....

check_point(M::UnitaryMatrices, p; kwargs...)
check_point(M::OrthogonalMatrices, p; kwargs...)
check_point(M::GeneralUnitaryMatrices, p; kwargs...)

Check whether p is a valid point on the UnitaryMatrices or [OrthogonalMatrices] M, i.e. that $p$ has a determinant of absolute value one.

The tolerance for the last test can be set using the kwargs....

check_vector(M::UnitaryMatrices, p, X; kwargs... )
check_vector(M::OrthogonalMatrices, p, X; kwargs... )
check_vector(M::Rotations, p, X; kwargs... )
check_vector(M::GeneralUnitaryMatrices, p, X; kwargs... )

Check whether X is a tangent vector to p on the UnitaryMatrices space M, i.e. after check_point(M,p), X has to be skew symmetric (Hermitian) and orthogonal to p.

The tolerance for the last test can be set using the kwargs....

embed(M::GeneralUnitaryMatrices, p, X)

Embed the tangent vector X at point p in M from its Lie algebra representation (set of skew matrices) into the Riemannian submanifold representation

The formula reads

\[X_{\text{embedded}} = p * X\]

get_coordinates(M::Rotations, p, X)
get_coordinates(M::OrthogonalMatrices, p, X)
get_coordinates(M::UnitaryMatrices, p, X)

Extract the unique tangent vector components $X^i$ at point p on Rotations $\mathrm{SO}(n)$ from the matrix representation X of the tangent vector.

The basis on the Lie algebra $𝔰𝔬(n)$ is chosen such that for $\mathrm{SO}(2)$, $X^1 = θ = X_{21}$ is the angle of rotation, and for $\mathrm{SO}(3)$, $(X^1, X^2, X^3) = (X_{32}, X_{13}, X_{21}) = θ u$ is the angular velocity and axis-angle representation, where $u$ is the unit vector along the axis of rotation.

For $\mathrm{SO}(n)$ where $n ≥ 4$, the additional elements of $X^i$ are $X^{j (j - 3)/2 + k + 1} = X_{jk}$, for $j ∈ [4,n], k ∈ [1,j)$.

get_embedding(M::OrthogonalMatrices)
get_embedding(M::Rotations)
get_embedding(M::UnitaryMatrices)

Return the embedding, i.e. The $\mathbb F^{n×n}$, where $\mathbb F = \mathbb R$ for the first two and $\mathbb F = \mathbb C$ for the unitary matrices.

get_vector(M::OrthogonalMatrices, p, Xⁱ, B::DefaultOrthogonalBasis)
get_vector(M::Rotations, p, Xⁱ, B::DefaultOrthogonalBasis)

Convert the unique tangent vector components Xⁱ at point p on Rotations or OrthogonalMatrices to the matrix representation $X$ of the tangent vector. See get_coordinates for the conventions used.

injectivity_radius(G::GeneraliUnitaryMatrices)

Return the injectivity radius for general unitary matrix manifolds, which is^[1]

\[ \operatorname{inj}_{\mathrm{U}(n)} = π.\]

injectivity_radius(G::SpecialOrthogonal)
injectivity_radius(G::Orthogonal)
injectivity_radius(M::Rotations)
injectivity_radius(M::Rotations, ::ExponentialRetraction)

Return the radius of injectivity on the Rotations manifold M, which is $π\sqrt{2}$. ^[1]

is_flat(M::GeneralUnitaryMatrices)

Return true if GeneralUnitaryMatrices M is SO(2) or U(1) and false otherwise.

manifold_dimension(M::Rotations)
manifold_dimension(M::OrthogonalMatrices)

Return the dimension of the manifold orthogonal matrices and of the manifold of rotations

\[\dim_{\mathrm{O}(n)} = \dim_{\mathrm{SO}(n)} = \frac{n(n-1)}{2}.\]

project(M::OrthogonalMatrices, p, X)
project(M::Rotations, p, X)
project(M::UnitaryMatrices, p, X)

Orthogonally project the tangent vector $X ∈ 𝔽^{n×n}$, $\mathbb F ∈ \{\mathbb R, \mathbb C\}$ to the tangent space of M at p, and change the representer to use the corresponding Lie algebra, i.e. we compute

\[ \operatorname{proj}_p(X) = \frac{p^{\mathrm{H}} X - (p^{\mathrm{H}} X)^{\mathrm{H}}}{2}.\]

 project(G::UnitaryMatrices, p)
 project(G::OrthogonalMatrices, p)

Project the point $p ∈ 𝔽^{n×n}$ to the nearest point in $\mathrm{U}(n,𝔽)=$UnitaryMatrices(n,𝔽) under the Frobenius norm. If $p = U S V^\mathrm{H}$ is the singular value decomposition of $p$, then the projection is

\[ \operatorname{proj}_{\mathrm{U}(n,𝔽)} \colon p ↦ U V^\mathrm{H}.\]

retract(M::Rotations, p, X, ::PolarRetraction)
retract(M::OrthogonalMatrices, p, X, ::PolarRetraction)

Compute the SVD-based retraction on the Rotations and OrthogonalMatrices M from p in direction X (as an element of the Lie group) and is a second-order approximation of the exponential map. Let

\[USV = p + pX\]

be the singular value decomposition, then the formula reads

\[\operatorname{retr}_p X = UV^\mathrm{T}.\]

retract(M::Rotations, p, X, ::QRRetraction)
retract(M::OrthogonalMatrices, p. X, ::QRRetraction)

Compute the QR-based retraction on the Rotations and OrthogonalMatrices M from p in direction X (as an element of the Lie group), which is a first-order approximation of the exponential map.

This is also the default retraction on these manifolds.

riemann_tensor(::GeneralUnitaryMatrices, p, X, Y, Z)

Compute the value of Riemann tensor on the GeneralUnitaryMatrices manifold. The formula reads [Ren11]

\[R(X,Y)Z=\frac{1}{4}[Z, [X, Y]].\]

mean(
    M::Rotations,
    x::AbstractVector,
    [w::AbstractWeights,]
    method = GeodesicInterpolationWithinRadius(π/2/√2);
    kwargs...,
)

Compute the Riemannian mean of x using GeodesicInterpolationWithinRadius.