Stiefel

Stiefel{n,k,𝔽} <: AbstractDecoratorManifold{𝔽}

The Stiefel manifold consists of all $n × k$, $n ≥ k$ unitary matrices, i.e.

\[\operatorname{St}(n,k) = \bigl\{ p ∈ 𝔽^{n × k}\ \big|\ p^{\mathrm{H}}p = I_k \bigr\},\]

where $𝔽 ∈ \{ℝ, ℂ\}$, $\cdot^{\mathrm{H}}$ denotes the complex conjugate transpose or Hermitian, and $I_k ∈ ℝ^{k × k}$ denotes the $k × k$ identity matrix.

The tangent space at a point $p ∈ \mathcal M$ is given by

\[T_p \mathcal M = \{ X ∈ 𝔽^{n × k} : p^{\mathrm{H}}X + \overline{X^{\mathrm{H}}p} = 0_k\},\]

where $0_k$ is the $k × k$ zero matrix and $\overline{\cdot}$ the (elementwise) complex conjugate.

This manifold is modeled as an embedded manifold to the Euclidean, i.e. several functions like the inner product and the zero_vector are inherited from the embedding.

The manifold is named after Eduard L. Stiefel (1909–1978).

Constructor

Stiefel(n, k, field = ℝ)

Generate the (real-valued) Stiefel manifold of $n × k$ dimensional orthonormal matrices.

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

When vector_at is nothing, return a random (Gaussian) point x on the Stiefel manifold M by generating a (Gaussian) matrix with standard deviation σ and return the orthogonalized version, i.e. return the Q component of the QR decomposition of the random matrix of size $n×k$.

When vector_at is not nothing, return a (Gaussian) random vector from the tangent space $T_{vector\_at}\mathrm{St}(n,k)$ with mean zero and standard deviation σ by projecting a random Matrix onto the tangent vector at vector_at.

uniform_distribution(M::Stiefel{n,k,ℝ}, p)

Uniform distribution on given (real-valued) Stiefel M. Specifically, this is the normalized Haar and Hausdorff measure on M. Generated points will be of similar type as p.

The implementation is based on Section 2.5.1 in ^{[Chikuse2003]}; see also Theorem 2.2.1(iii) in ^{[Chikuse2003]}.

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

Check whether p is a valid point on the Stiefel M=$\operatorname{St}(n,k)$, i.e. that it has the right AbstractNumbers type and $p^{\mathrm{H}}p$ is (approximately) the identity, where $\cdot^{\mathrm{H}}$ is the complex conjugate transpose. The settings for approximately can be set with kwargs....

check_vector(M::Stiefel, p, X; kwargs...)

Checks whether X is a valid tangent vector at p on the Stiefel M=$\operatorname{St}(n,k)$, i.e. the AbstractNumbers fits and it (approximately) holds that $p^{\mathrm{H}}X + \overline{X^{\mathrm{H}}p} = 0$, where $\cdot^{\mathrm{H}}$ denotes the Hermitian and $\overline{\cdot}$ the (elementwise) complex conjugate. The settings for approximately can be set with kwargs....

default_inverse_retraction_method(M::Stiefel)

Return PolarInverseRetraction as the default inverse retraction for the Stiefel manifold.

default_retraction_method(M::Stiefel)

Return PolarRetraction as the default retraction for the Stiefel manifold.

default_vector_transport_method(M::Stiefel)

Return the DifferentiatedRetractionVectorTransport of the PolarRetraction as the default vector transport method for the Stiefel manifold.

inverse_retract(M::Stiefel, p, q, ::PolarInverseRetraction)

Compute the inverse retraction based on a singular value decomposition for two points p, q on the Stiefel manifold M. This follows the folloing approach: From the Polar retraction we know that

\[\operatorname{retr}_p^{-1}q = qs - t\]

if such a symmetric positive definite $k × k$ matrix exists. Since $qs - t$ is also a tangent vector at $p$ we obtain

\[p^{\mathrm{H}}qs + s(p^{\mathrm{H}}q)^{\mathrm{H}} + 2I_k = 0,\]

which can either be solved by a Lyapunov approach or a continuous-time algebraic Riccati equation.

This implementation follows the Lyapunov approach.

inverse_retract(M::Stiefel, p, q, ::QRInverseRetraction)

Compute the inverse retraction based on a qr decomposition for two points p, q on the Stiefel manifold M and return the resulting tangent vector in X. The computation follows Algorithm 1 in ^{[KanekoFioriTanaka2013]}.

manifold_dimension(M::Stiefel)

Return the dimension of the Stiefel manifold M=$\operatorname{St}(n,k,𝔽)$. The dimension is given by

\[\begin{aligned} \dim \mathrm{St}(n, k, ℝ) &= nk - \frac{1}{2}k(k+1)\\ \dim \mathrm{St}(n, k, ℂ) &= 2nk - k^2\\ \dim \mathrm{St}(n, k, ℍ) &= 4nk - k(2k-1) \end{aligned}\]

representation_size(M::Stiefel)

Returns the representation size of the Stiefel M=$\operatorname{St}(n,k)$, i.e. (n,k), which is the matrix dimensions.

retract(::Stiefel, p, X, ::CayleyRetraction)

Compute the retraction on the Stiefel that is based on the Cayley transform^[Zhu2017]. Using

\[ W_{p,X} = \operatorname{P}_pXp^{\mathrm{H}} - pX^{\mathrm{H}}\operatorname{P_p} \quad\text{where}  \operatorname{P}_p = I - \frac{1}{2}pp^{\mathrm{H}}\]

the formula reads

\[ \operatorname{retr}_pX = \Bigl(I - \frac{1}{2}W_{p,X}\Bigr)^{-1}\Bigl(I + \frac{1}{2}W_{p,X}\Bigr)p.\]

It is implemented as the case $m=1$ of the PadeRetraction.

retract(M::Stiefel, p, X, ::PadeRetraction{m})

Compute the retraction on the Stiefel manifold M based on the Padé approximation of order $m$^{[ZhuDuan2018]}. Let $p_m$ and $q_m$ be defined for any matrix $A ∈ ℝ^{n×x}$ as

\[ p_m(A) = \sum_{k=0}^m \frac{(2m-k)!m!}{(2m)!(m-k)!}\frac{A^k}{k!}\]

and

\[ q_m(A) = \sum_{k=0}^m \frac{(2m-k)!m!}{(2m)!(m-k)!}\frac{(-A)^k}{k!}\]

respectively. Then the Padé approximation (of the matrix exponential $\exp(A)$) reads

\[ r_m(A) = q_m(A)^{-1}p_m(A)\]

Defining further

\[ W_{p,X} = \operatorname{P}_pXp^{\mathrm{H}} - pX^{\mathrm{H}}\operatorname{P_p} \quad\text{where } \operatorname{P}_p = I - \frac{1}{2}pp^{\mathrm{H}}\]

the retraction reads

\[ \operatorname{retr}_pX = r_m(W_{p,X})p\]

retract(M::Stiefel, p, X, ::PolarRetraction)

Compute the SVD-based retraction PolarRetraction on the Stiefel manifold M. With $USV = p + X$ the retraction reads

\[\operatorname{retr}_p X = U\bar{V}^\mathrm{H}.\]

retract(M::Stiefel, p, X, ::QRRetraction)

Compute the QR-based retraction QRRetraction on the Stiefel manifold M. With $QR = p + X$ the retraction reads

\[\operatorname{retr}_p X = QD,\]

where $D$ is a $n × k$ matrix with

\[D = \operatorname{diag}\bigl(\operatorname{sgn}(R_{ii}+0,5)_{i=1}^k \bigr),\]

where $\operatorname{sgn}(p) = \begin{cases} 1 & \text{ for } p > 0,\\ 0 & \text{ for } p = 0,\\ -1& \text{ for } p < 0. \end{cases}$

vector_transport_direction(::Stiefel, p, X, d, ::DifferentiatedRetractionVectorTransport{CayleyRetraction})

Compute the vector transport given by the differentiated retraction of the CayleyRetraction, cf. ^[Zhu2017] Equation (17).

The formula reads

\[\operatorname{T}_{p,d}(X) = \Bigl(I - \frac{1}{2}W_{p,d}\Bigr)^{-1}W_{p,X}\Bigl(I - \frac{1}{2}W_{p,d}\Bigr)^{-1}p,\]

with

\[ W_{p,X} = \operatorname{P}_pXp^{\mathrm{H}} - pX^{\mathrm{H}}\operatorname{P_p} \quad\text{where } \operatorname{P}_p = I - \frac{1}{2}pp^{\mathrm{H}}\]

Since this is the differentiated retraction as a vector transport, the result will be in the tangent space at $q=\operatorname{retr}_p(d)$ using the CayleyRetraction.

vector_transport_direction(M::Stiefel, p, X, d, DifferentiatedRetractionVectorTransport{PolarRetraction})

Compute the vector transport by computing the push forward of retract(::Stiefel, ::Any, ::Any, ::PolarRetraction) Section 3.5 of ^[Zhu2017]:

\[T_{p,d}^{\text{Pol}}(X) = q*Λ + (I-qq^{\mathrm{T}})X(1+d^\mathrm{T}d)^{-\frac{1}{2}},\]

where $q = \operatorname{retr}^{\mathrm{Pol}}_p(d)$, and $Λ$ is the unique solution of the Sylvester equation

\[ Λ(I+d^\mathrm{T}d)^{\frac{1}{2}} + (I + d^\mathrm{T}d)^{\frac{1}{2}} = q^\mathrm{T}X - X^\mathrm{T}q\]

vector_transport_direction(M::Stiefel, p, X, d, DifferentiatedRetractionVectorTransport{QRRetraction})

Compute the vector transport by computing the push forward of the retract(::Stiefel, ::Any, ::Any, ::QRRetraction), See ^{[AbsilMahonySepulchre2008]}, p. 173, or Section 3.5 of ^[Zhu2017].

\[T_{p,d}^{\text{QR}}(X) = q*\rho_{\mathrm{s}}(q^\mathrm{T}XR^{-1}) + (I-qq^{\mathrm{T}})XR^{-1},\]

where $q = \operatorname{retr}^{\mathrm{QR}}_p(d)$, $R$ is the $R$ factor of the QR decomposition of $p + d$, and

\[\bigl( \rho_{\mathrm{s}}(A) \bigr)_{ij} = \begin{cases} A_{ij}&\text{ if } i > j\\ 0 \text{ if } i = j\\ -A_{ji} \text{ if } i < j.\\ \end{cases}\]

vector_transport_to(M::Stiefel, p, X, q, DifferentiatedRetractionVectorTransport{PolarRetraction})

Compute the vector transport by computing the push forward of the retract(M::Stiefel, ::Any, ::Any, ::PolarRetraction), see Section 4 of ^{[HuangGallivanAbsil2015]} or Section 3.5 of ^[Zhu2017]:

\[T_{q\gets p}^{\text{Pol}}(X) = q*Λ + (I-qq^{\mathrm{T}})X(1+d^\mathrm{T}d)^{-\frac{1}{2}},\]

where $d = \bigl( \operatorname{retr}^{\mathrm{Pol}}_p\bigr)^{-1}(q)$, and $Λ$ is the unique solution of the Sylvester equation

\[ Λ(I+d^\mathrm{T}d)^{\frac{1}{2}} + (I + d^\mathrm{T}d)^{\frac{1}{2}} = q^\mathrm{T}X - X^\mathrm{T}q\]

vector_transport_to(M::Stiefel, p, X, q, DifferentiatedRetractionVectorTransport{QRRetraction})

Compute the vector transport by computing the push forward of the retract(M::Stiefel, ::Any, ::Any, ::QRRetraction), see ^{[AbsilMahonySepulchre2008]}, p. 173, or Section 3.5 of ^[Zhu2017].

\[T_{q \gets p}^{\text{QR}}(X) = q*\rho_{\mathrm{s}}(q^\mathrm{T}XR^{-1}) + (I-qq^{\mathrm{T}})XR^{-1},\]

where $d = \bigl(\operatorname{retr}^{\mathrm{QR}}\bigr)^{-1}_p(q)$, $R$ is the $R$ factor of the QR decomposition of $p+X$, and

\[\bigl( \rho_{\mathrm{s}}(A) \bigr)_{ij} = \begin{cases} A_{ij}&\text{ if } i > j\\ 0 \text{ if } i = j\\ -A_{ji} \text{ if } i < j.\\ \end{cases}\]

vector_transport_to(M::Stiefel, p, X, q, ::ProjectionTransport)

Compute a vector transport by projection, i.e. project X from the tangent space at p by projection it onto the tangent space at q.

exp(M::Stiefel, p, X)

Compute the exponential map on the Stiefel{n,k,𝔽}() manifold M emanating from p in tangent direction X.

\[\exp_p X = \begin{pmatrix} p\\X \end{pmatrix} \operatorname{Exp} \left( \begin{pmatrix} p^{\mathrm{H}}X & - X^{\mathrm{H}}X\\ I_n & p^{\mathrm{H}}X\end{pmatrix} \right) \begin{pmatrix} \exp( -p^{\mathrm{H}}X) \\ 0_n\end{pmatrix},\]

where $\operatorname{Exp}$ denotes matrix exponential, $\cdot^{\mathrm{H}}$ denotes the complex conjugate transpose or Hermitian, and $I_k$ and $0_k$ are the identity matrix and the zero matrix of dimension $k × k$, respectively.

get_basis(M::Stiefel{n,k,ℝ}, p, B::DefaultOrthonormalBasis) where {n,k}

Create the default basis using the parametrization for any $X ∈ T_p\mathcal M$. Set $p_\bot \in ℝ^{n\times(n-k)}$ the matrix such that the $n\times n$ matrix of the common columns $[p\ p_\bot]$ is an ONB. For any skew symmetric matrix $a ∈ ℝ^{k\times k}$ and any $b ∈ ℝ^{(n-k)\times k}$ the matrix

\[X = pa + p_\bot b ∈ T_p\mathcal M\]

and we can use the $\frac{1}{2}k(k-1) + (n-k)k = nk-\frac{1}{2}k(k+1)$ entries of $a$ and $b$ to specify a basis for the tangent space. using unit vectors for constructing both the upper matrix of $a$ to build a skew symmetric matrix and the matrix b, the default basis is constructed.

Since $[p\ p_\bot]$ is an automorphism on $ℝ^{n\times p}$ the elements of $a$ and $b$ are orthonormal coordinates for the tangent space. To be precise exactly one element in the upper trangular entries of $a$ is set to $1$ its symmetric entry to $-1$ and we normalize with the factor $\frac{1}{\sqrt{2}}$ and for $b$ one can just use unit vectors reshaped to a matrix to obtain orthonormal set of parameters.

inverse_retract(M::Stiefel, p, q, method::ProjectionInverseRetraction)

Compute a projection-based inverse retraction.

The inverse retraction is computed by projecting the logarithm map in the embedding to the tangent space at $p$.

project(M::Stiefel,p)

Projects p from the embedding onto the Stiefel M, i.e. compute q as the polar decomposition of $p$ such that $q^{\mathrm{H}q$ is the identity, where $\cdot^{\mathrm{H}}$ denotes the hermitian, i.e. complex conjugate transposed.

project(M::Stiefel, p, X)

Project X onto the tangent space of p to the Stiefel manifold M. The formula reads

\[\operatorname{proj}_{\mathcal M}(p, X) = X - p \operatorname{Sym}(p^{\mathrm{H}}X),\]

where $\operatorname{Sym}(q)$ is the symmetrization of $q$, e.g. by $\operatorname{Sym}(q) = \frac{q^{\mathrm{H}}+q}{2}$.

retract(M::Stiefel, p, X, method::ProjectionRetraction)

Compute a projection-based retraction.

The retraction is computed by projecting the exponential map in the embedding to M.

ApproximateLogarithmicMap <: ApproximateInverseRetraction

An approximate implementation of the logarithmic map, which is an inverse_retraction. See inverse_retract(::MetricManifold{ℝ,Stiefel{n,k,ℝ},CanonicalMetric}, ::Any, ::Any, ::ApproximateLogarithmicMap) where {n,k} for a use case.

Fields

• max_iterations – maximal number of iterations used in the approximation
• tolerance – a tolerance used as a stopping criterion

CanonicalMetric <: AbstractMetric

The Canonical Metric refers to a metric for the Stiefel manifold, see^{[EdelmanAriasSmith1998]}.

q = exp(M::MetricManifold{ℝ, Stiefel{n,k,ℝ}, CanonicalMetric}, p, X)
exp!(M::MetricManifold{ℝ, Stiefel{n,k,ℝ}, q, CanonicalMetric}, p, X)

Compute the exponential map on the Stiefel(n,k) manifold with respect to the CanonicalMetric.

First, decompose The tangent vector $X$ into its horizontal and vertical component with respect to $p$, i.e.

\[X = pp^{\mathrm{T}}X + (I_n-pp^{\mathrm{T}})X,\]

where $I_n$ is the $n\times n$ identity matrix. We introduce $A=p^{\mathrm{T}}X$ and $QR = (I_n-pp^{\mathrm{T}})X$ the qr decomposition of the vertical component. Then using the matrix exponential $\operatorname{Exp}$ we introduce $B$ and $C$ as

\[\begin{pmatrix} B\\C \end{pmatrix} \coloneqq \operatorname{Exp}\left( \begin{pmatrix} A & -R^{\mathrm{T}}\\ R & 0 \end{pmatrix} \right) \begin{pmatrix}I_k\\0\end{pmatrix}\]

the exponential map reads

\[q = \exp_p X = pC + QB.\]

For more details, see ^{[EdelmanAriasSmith1998][Zimmermann2017]}.

inner(M::MetricManifold{ℝ, Stiefel{n,k,ℝ}, X, CanonicalMetric}, p, X, Y)

Compute the inner product on the Stiefel manifold with respect to the CanonicalMetric. The formula reads

\[g_p(X,Y) = \operatorname{tr}\bigl( X^{\mathrm{T}}(I_n - \frac{1}{2}pp^{\mathrm{T}})Y \bigr).\]

X = inverse_retract(M::MetricManifold{ℝ, Stiefel{n,k,ℝ}, CanonicalMetric}, p, q, a::ApproximateLogarithmicMap)
inverse_retract!(M::MetricManifold{ℝ, Stiefel{n,k,ℝ}, X, CanonicalMetric}, p, q, a::ApproximateLogarithmicMap)

Compute an approximation to the logarithmic map on the Stiefel(n,k) manifold with respect to the CanonicalMetric using a matrix-algebraic based approach to an iterative inversion of the formula of the exp.

The algorithm is derived in^{[Zimmermann2017]} and it uses the max_iterations and the tolerance field from the ApproximateLogarithmicMap.

StiefelSubmersionMetric{T<:Real} <: RiemannianMetric

The submersion (or normal) metric family on the Stiefel manifold.

The family, with a single real parameter $α>-1$, has two special cases:

• $α = -\frac{1}{2}$: EuclideanMetric
• $α = 0$: CanonicalMetric

The family was described in ^{[HüperMarkinaLeite2021]}. This implementation follows the description in ^{[ZimmermanHüper2022]}.

Constructor

StiefelSubmersionMetric(α)

Construct the submersion metric on the Stiefel manifold with the parameter $α$.

q = exp(M::MetricManifold{ℝ, Stiefel{n,k,ℝ}, <:StiefelSubmersionMetric}, p, X)
exp!(M::MetricManifold{ℝ, Stiefel{n,k,ℝ}, q, <:StiefelSubmersionMetric}, p, X)

Compute the exponential map on the Stiefel(n,k) manifold with respect to the StiefelSubmersionMetric.

The exponential map is given by

\[\exp_p X = \operatorname{Exp}\bigl( -\frac{2α+1}{α+1} p p^\mathrm{T} X p^\mathrm{T} + X p^\mathrm{T} - p X^\mathrm{T} \bigr) p \operatorname{Exp}\bigl(\frac{\alpha}{\alpha+1} p^\mathrm{T} X\bigr)\]

This implementation is based on ^{[ZimmermanHüper2022]}.

For $k < \frac{n}{2}$ the exponential is computed more efficiently using StiefelFactorization.

log(M::MetricManifold{ℝ,Stiefel{n,k,ℝ},<:StiefelSubmersionMetric}, p, q; kwargs...)

Compute the logarithmic map on the Stiefel(n,k) manifold with respect to the StiefelSubmersionMetric.

The logarithmic map is computed using ShootingInverseRetraction. For $k ≤ \lfloor\frac{n}{2}\rfloor$, this is sped up using the $k$-shooting method of ^{[ZimmermanHüper2022]}. Keyword arguments are forwarded to ShootingInverseRetraction; see that documentation for details. Their defaults are:

• num_transport_points=4
• tolerance=sqrt(eps())
• max_iterations=1_000

inner(M::MetricManifold{ℝ, Stiefel{n,k,ℝ}, X, <:StiefelSubmersionMetric}, p, X, Y)

Compute the inner product on the Stiefel manifold with respect to the StiefelSubmersionMetric. The formula reads

\[g_p(X,Y) = \operatorname{tr}\bigl( X^{\mathrm{T}}(I_n - \frac{2α+1}{2(α+1)}pp^{\mathrm{T}})Y \bigr),\]

where $α$ is the parameter of the metric.

inverse_retract(
    M::MetricManifold{ℝ,Stiefel{n,k,ℝ},<:StiefelSubmersionMetric},
    p,
    q,
    method::ShootingInverseRetraction,
)

Compute the inverse retraction using ShootingInverseRetraction.

In general the retraction is computed using the generic shooting method.

inverse_retract(
    M::MetricManifold{ℝ,Stiefel{n,k,ℝ},<:StiefelSubmersionMetric},
    p,
    q,
    method::ShootingInverseRetraction{
        ExponentialRetraction,
        ProjectionInverseRetraction,
        <:Union{ProjectionTransport,ScaledVectorTransport{ProjectionTransport}},
    },
)

Compute the inverse retraction using ShootingInverseRetraction more efficiently.

For $k < \frac{n}{2}$ the retraction is computed more efficiently using StiefelFactorization.

StiefelFactorization{UT,XT} <: AbstractManifoldPoint

Represent points (and vectors) on Stiefel(n, k) with $2k × k$ factors.^{[ZimmermanHüper2022]}

Given a point $p ∈ \mathrm{St}(n, k)$ and another matrix $B ∈ ℝ^{n × k}$ for $k ≤ \lfloor\frac{n}{2}\rfloor$ the factorization is

\[\begin{aligned} B &= UZ\\ U &= \begin{bmatrix}p & Q\end{bmatrix} ∈ \mathrm{St}(n, 2k)\\ Z &= \begin{bmatrix}Z_1 \\ Z_2\end{bmatrix}, \quad Z_1,Z_2 ∈ ℝ^{k × k}. \end{aligned}\]

If $B ∈ \mathrm{St}(n, k)$, then $Z ∈ \mathrm{St}(2k, k)$. Note that not every matrix $B$ can be factorized in this way.

For a fixed $U$, if $r ∈ \mathrm{St}(n, k)$ has the factor $Z_r ∈ \mathrm{St}(2k, k)$, then $X_r ∈ T_r \mathrm{St}(n, k)$ has the factor $Z_{X_r} ∈ T_{Z_r} \mathrm{St}(2k, k)$.

$Q$ is determined by choice of a second matrix $A ∈ ℝ^{n × k}$ with the decomposition

\[\begin{aligned} A &= UZ\\ Z_1 &= p^\mathrm{T} A \\ Q Z_2 &= (I - p p^\mathrm{T}) A, \end{aligned}\]

where here $Q Z_2$ is the any decomposition that produces $Q ∈ \mathrm{St}(n, k)$, for which we choose the QR decomposition.

This factorization is useful because it is closed under addition, subtraction, scaling, projection, and the Riemannian exponential and logarithm under the StiefelSubmersionMetric. That is, if all matrices involved are factorized to have the same $U$, then all of these operations and any algorithm that depends only on them can be performed in terms of the $2k × k$ matrices $Z$. For $n ≫ k$, this can be much more efficient than working with the full matrices.

stiefel_factorization(p, x) -> StiefelFactorization

Compute the StiefelFactorization of $x$ relative to the point $p$.