Stiefel
Manifolds.CaleyRetraction
— TypeCaleyRetraction <: AbstractRetractionMethod
A retraction based on the Caley transform, which is realized by using the PadeRetraction
{1}
.
Manifolds.CaleyVectorTransport
— TypeCaleyVectorTransport <: AbstractVectorTransportMethod
A vector transport that one obtains by differentiating a CaleyRetraction
.
Manifolds.PadeRetraction
— TypePadeRetraction{m} <: AbstractRetractionMethod
A retraction based on the Padé approximation of order $m$
Manifolds.Stiefel
— TypeStiefel{n,k,𝔽} <: AbstractEmbeddedManifold{𝔽,DefaultIsometricEmbeddingType}
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_tangent_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.
Base.exp
— Methodexp(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.
Manifolds.uniform_distribution
— Methoduniform_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].
ManifoldsBase.check_manifold_point
— Methodcheck_manifold_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...
.
ManifoldsBase.check_tangent_vector
— Methodcheck_tangent_vector(M::Stiefel, p, X; check_base_point = true, 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 optional parameter check_base_point
indicates, whether to call check_manifold_point
for p
. The settings for approximately can be set with kwargs...
.
ManifoldsBase.get_basis
— Methodget_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.
ManifoldsBase.inverse_retract
— Methodinverse_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.
ManifoldsBase.inverse_retract
— Methodinverse_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].
ManifoldsBase.manifold_dimension
— Methodmanifold_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}\]
ManifoldsBase.project
— Methodproject(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.
ManifoldsBase.project
— Methodproject(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}$.
ManifoldsBase.representation_size
— Methodrepresentation_size(M::Stiefel)
Returns the representation size of the Stiefel
M
=$\operatorname{St}(n,k)$, i.e. (n,k)
, which is the matrix dimensions.
ManifoldsBase.retract
— Methodretract(::Stiefel, p, X, ::CaleyRetraction)
Compute the retraction on the Stiefel
that is based on the Caley transform[Zhu2016]. 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
.
ManifoldsBase.retract
— Methodretract(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\]
ManifoldsBase.retract
— Methodretract(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}.\]
ManifoldsBase.retract
— Methodretract(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}$
ManifoldsBase.vector_transport_direction
— Methodvector_transport_direction(::Stiefel, p, X, d, ::CaleyVectorTransport)
Compute the vector transport given by the differentiated retraction of the CaleyRetraction
, cf. [Zhu2016] Equation (17).
The formula reads
\[\operatorname{T}_{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 CaleyRetraction
.
Literature
- Chikuse2003
Y. Chikuse: "Statistics on Special Manifolds", Springer New York, 2003, doi: 10.1007/978-0-387-21540-2.
- KanekoFioriTanaka2013
T. Kaneko, S. Fiori, T. Tanaka: "Empirical Arithmetic Averaging over the Compact Stiefel Manifold", IEEE Transactions on Signal Processing, 2013, doi: 10.1109/TSP.2012.2226167.
- Zhu2016
X. Zhu: A Riemannian conjugate gradient method for optimizazion on the Stiefel manifold, Computational Optimization and Applications 67(1), pp. 73–110, 2016. doi 10.1007/s10589-016-9883-4.
- ZhuDuan2018
X. Zhu, C. Duan: On matrix exponentials and their approximations related to optimization on the Stiefel manifold, Optimizazion Letters 13(5), pp. 1069–1083, 2018. doi 10.1007/s11590-018-1341-z.