Grassmannian manifold

Manifolds.Grassmann β€” Type
Grassmann{T,𝔽} <: AbstractDecoratorManifold{𝔽}

The Grassmann manifold Gr(n,k)\mathrm{Gr}(n,k) consists of all subspaces spanned by kk linear independent vectors 𝔽n𝔽^n, where π”½βˆˆ{R,C}𝔽 ∈ \{ℝ, β„‚\} is either the real- (or complex-) valued vectors. This yields all kk-dimensional subspaces of Rnℝ^n for the real-valued case and all 2k2k-dimensional subspaces of Cnβ„‚^n for the second.

The manifold can be represented as

Gr(n,k):={span⁑(p):pβˆˆπ”½nΓ—k,pHp=Ik},\mathrm{Gr}(n,k) := \bigl\{ \operatorname{span}(p) : p ∈ 𝔽^{nΓ—k}, p^\mathrm{H}p = I_k\},

where β‹…Hβ‹…^{\mathrm{H}} denotes the complex conjugate transpose or Hermitian and IkI_k is the kΓ—kkΓ—k identity matrix. This means, that the columns of pp form an unitary basis of the subspace, that is a point on Gr⁑(n,k)\operatorname{Gr}(n,k), and hence the subspace can actually be represented by a whole equivalence class of representers. Another interpretation is, that

Gr(n,k)=St(n,k)/O⁑(k),\mathrm{Gr}(n,k) = \mathrm{St}(n,k) / \operatorname{O}(k),

i.e the Grassmann manifold is the quotient of the Stiefel manifold and the orthogonal group O⁑(k)\operatorname{O}(k) of orthogonal kΓ—kkΓ—k matrices. Note that it doesn't matter whether we start from the Euclidean or canonical metric on the Stiefel manifold, the resulting quotient metric on Grassmann is the same.

The tangent space at a point (subspace) pp is given by

TpGr(n,k)={Xβˆˆπ”½nΓ—k:XHp+pHX=0k},T_p\mathrm{Gr}(n,k) = \bigl\{ X ∈ 𝔽^{nΓ—k} : X^{\mathrm{H}}p + p^{\mathrm{H}}X = 0_{k} \bigr\},

where 0k0_k is the kΓ—kkΓ—k zero matrix.

Note that a point p∈Gr⁑(n,k)p ∈ \operatorname{Gr}(n,k) might be represented by different matrices (i.e. matrices with unitary column vectors that span the same subspace). Different representations of pp also lead to different representation matrices for the tangent space TpGr(n,k)T_p\mathrm{Gr}(n,k)

For a representation of points as orthogonal projectors. Here

Gr⁑(n,k):={p∈RnΓ—n:p=p˜T,p2=p,rank⁑(p)=k},\operatorname{Gr}(n,k) := \bigl\{ p \in \mathbb R^{nΓ—n} : p = p^˜\mathrm{T}, p^2 = p, \operatorname{rank}(p) = k\},

with tangent space

TpGr(n,k)={X∈RnΓ—n:X=XT and X=pX+Xp},T_p\mathrm{Gr}(n,k) = \bigl\{ X ∈ \mathbb R^{nΓ—n} : X=X^{\mathrm{T}} \text{ and } X = pX+Xp \bigr\},

see also ProjectorPoint and ProjectorTangentVector.

The manifold is named after Hermann G. Graßmann (1809-1877).

A good overview can be found in[BZA20].

Constructor

Grassmann(n, k, field=ℝ, parameter::Symbol=:type)

Generate the Grassmann manifold Gr⁑(n,k)\operatorname{Gr}(n,k), where the real-valued case field=ℝ is the default.

source
Base.convert β€” Method
convert(::Type{ProjectorPoint}, p::AbstractMatrix)

Convert a point p on Stiefel that also represents a point (i.e. subspace) on Grassmann to a projector representation of said subspace, i.e. compute the canonical_project! for

Ο€SG(p)=ppT. Ο€^{\mathrm{SG}}(p) = pp^{\mathrm{T}}.

source
Base.convert β€” Method
convert(::Type{ProjectorPoint}, ::Stiefelpoint)

Convert a point p on Stiefel that also represents a point (i.e. subspace) on Grassmann to a projector representation of said subspace, i.e. compute the canonical_project! for

Ο€SG(p)=ppT. Ο€^{\mathrm{SG}}(p) = pp^{\mathrm{T}}.

source
ManifoldsBase.change_metric β€” Method
change_metric(M::Grassmann, ::EuclideanMetric, p X)

Change X to the corresponding vector with respect to the metric of the Grassmann M, which is just the identity, since the manifold is isometrically embedded.

source
ManifoldsBase.change_representer β€” Method
change_representer(M::Grassmann, ::EuclideanMetric, p, X)

Change X to the corresponding representer of a cotangent vector at p. Since the Grassmann manifold M, is isometrically embedded, this is the identity

source
ManifoldsBase.default_retraction_method β€” Method
default_retraction_method(M::Grassmann)
default_retraction_method(M::Grassmann, ::Type{StiefelPoint})
default_retraction_method(M::Grassmann, ::Type{ProjectorPoint})

Return ExponentialRetraction as the default on the Grassmann manifold for both representations.

source
ManifoldsBase.manifold_dimension β€” Method
manifold_dimension(M::Grassmann)

Return the dimension of the Grassmann(n,k,𝔽) manifold M, i.e.

dim⁑Gr⁑(n,k)=k(nβˆ’k)dim⁑R𝔽,\dim \operatorname{Gr}(n,k) = k(n-k) \dim_ℝ 𝔽,

where dim⁑R𝔽\dim_ℝ 𝔽 is the real_dimension of 𝔽.

source

The Grassmanian represented as points on the Stiefel manifold

Manifolds.StiefelPoint β€” Type
StiefelPoint <: AbstractManifoldPoint

A point on a Stiefel manifold. This point is mainly used for representing points on the Grassmann where this is also the default representation and hence equivalent to using AbstractMatrices thereon. they can also used be used as points on Stiefel.

source
Manifolds.StiefelTangentVector β€” Type
StiefelTangentVector <: AbstractTangentVector

A tangent vector on the Grassmann manifold represented by a tangent vector from the tangent space of a corresponding point from the Stiefel manifold, see StiefelPoint. This is the default representation so is can be used interchangeably with just abstract matrices.

source
Base.exp β€” Method
exp(M::Grassmann, p, X)

Compute the exponential map on the Grassmann M =Gr(n,k)= \mathrm{Gr}(n,k) starting in p with tangent vector (direction) X. Let X=USVX = USV denote the SVD decomposition of XX. Then the exponential map is written using

z=pVcos⁑(S)VH+Usin⁑(S)VH,z = p V\cos(S)V^\mathrm{H} + U\sin(S)V^\mathrm{H},

where β‹…Hβ‹…^{\mathrm{H}} denotes the complex conjugate transposed or Hermitian and the cosine and sine are applied element wise to the diagonal entries of SS. A final QR decomposition z=QRz=QR is performed for numerical stability reasons, yielding the result as

exp⁑pX=Q.\exp_p X = Q.

source
Base.log β€” Method
log(M::Grassmann, p, q)

Compute the logarithmic map on the Grassmann M =M=Gr(n,k)= \mathcal M=\mathrm{Gr}(n,k), i.e. the tangent vector X whose corresponding geodesic starting from p reaches q after time 1 on M. The formula reads

log⁑pq=Vβ‹…atan⁑(S)β‹…UH,\log_p q = Vβ‹… \operatorname{atan}(S) β‹… U^\mathrm{H},

where β‹…Hβ‹…^{\mathrm{H}} denotes the complex conjugate transposed or Hermitian. The matrices UU and VV are the unitary matrices, and SS is the diagonal matrix containing the singular values of the SVD-decomposition

USV=(qHp)βˆ’1(qHβˆ’qHppH).USV = (q^\mathrm{H}p)^{-1} ( q^\mathrm{H} - q^\mathrm{H}pp^\mathrm{H}).

In this formula the atan⁑\operatorname{atan} is meant elementwise.

source
Base.rand β€” Method
rand(M::Grassmann; Οƒ::Real=1.0, vector_at=nothing)

When vector_at is nothing, return a random point p on Grassmann manifold M by generating a random (Gaussian) matrix with standard deviation Οƒ in matching size, which is orthonormal.

When vector_at is not nothing, return a (Gaussian) random vector from the tangent space TpGr(n,k)T_p\mathrm{Gr}(n,k) with mean zero and standard deviation Οƒ by projecting a random Matrix onto the tangent space at vector_at.

source
ManifoldDiff.riemannian_Hessian β€” Method
riemannian_Hessian(M::Grassmann, p, G, H, X)

The Riemannian Hessian can be computed by adopting Eq. (6.6) [Ngu23], where we use for the EuclideanMetric Ξ±0=Ξ±1=1Ξ±_0=Ξ±_1=1 in their formula. Let βˆ‡f(p)\nabla f(p) denote the Euclidean gradient G, βˆ‡2f(p)[X]\nabla^2 f(p)[X] the Euclidean Hessian H. Then the formula reads

Hess⁑f(p)[X]=proj⁑TpM(βˆ‡2f(p)[X]βˆ’XpHβˆ‡f(p)). \operatorname{Hess}f(p)[X] = \operatorname{proj}_{T_p\mathcal M}\Bigl( βˆ‡^2f(p)[X] - X p^{\mathrm{H}}βˆ‡f(p) \Bigr).

Compared to Eq. (5.6) also the metric conversion simplifies to the identity.

source
ManifoldsBase.distance β€” Method
distance(M::Grassmann, p, q)

Compute the Riemannian distance on Grassmann manifold M=Gr(n,k)= \mathrm{Gr}(n,k).

The distance is given by

dGr(n,k)(p,q)=norm⁑(log⁑p(q)).d_{\mathrm{Gr}(n,k)}(p,q) = \operatorname{norm}(\log_p(q)).

source
ManifoldsBase.inner β€” Method
inner(M::Grassmann, p, X, Y)

Compute the inner product for two tangent vectors X, Y from the tangent space of p on the Grassmann manifold M. The formula reads

gp(X,Y)=tr⁑(XHY),g_p(X,Y) = \operatorname{tr}(X^{\mathrm{H}}Y),

where β‹…Hβ‹…^{\mathrm{H}} denotes the complex conjugate transposed or Hermitian.

source
ManifoldsBase.inverse_retract β€” Method
inverse_retract(M::Grassmann, p, q, ::PolarInverseRetraction)

Compute the inverse retraction for the PolarRetraction, on the Grassmann manifold M, i.e.,

retr⁑pβˆ’1q=qβˆ—(pHq)βˆ’1βˆ’p,\operatorname{retr}_p^{-1}q = q*(p^\mathrm{H}q)^{-1} - p,

where β‹…Hβ‹…^{\mathrm{H}} denotes the complex conjugate transposed or Hermitian.

source
ManifoldsBase.inverse_retract β€” Method
inverse_retract(M, p, q, ::QRInverseRetraction)

Compute the inverse retraction for the QRRetraction, on the Grassmann manifold M, i.e.,

retr⁑pβˆ’1q=q(pHq)βˆ’1βˆ’p,\operatorname{retr}_p^{-1}q = q(p^\mathrm{H}q)^{-1} - p,

where β‹…Hβ‹…^{\mathrm{H}} denotes the complex conjugate transposed or Hermitian.

source
ManifoldsBase.parallel_transport_direction β€” Method
parallel_transport_direction(M::Grassmann, p, X, Y)

Compute the parallel transport of X∈TpMX \in T_p\mathcal M along the geodesic starting in direction Ξ³Λ™(0)=Y\dot Ξ³ (0) = Y.

Let Y=USVY = USV denote the SVD decomposition of YY. Then the parallel transport is given by the formula according to Equation (8.5) (p. 171) [AMS08] as

Pp,YX=βˆ’pVsin⁑(S)UTX+Ucos⁑(S)UTX+(Iβˆ’UUT)X\mathcal P_{p,Y} X = -pV \sin(S)U^{\mathrm{T}}X + U\cos(S)U^{\mathrm{T}}X + (I-UU^{\mathrm{T}})X

where the sine and cosine applied to the diagonal matrix SS are meant to be elementwise

source
ManifoldsBase.project β€” Method
project(M::Grassmann, p)

Project p from the embedding onto the Grassmann M, i.e. compute q as the polar decomposition of pp such that qHqq^{\mathrm{H}}q is the identity, where β‹…Hβ‹…^{\mathrm{H}} denotes the Hermitian, i.e. complex conjugate transposed.

source
ManifoldsBase.project β€” Method
project(M::Grassmann, p, X)

Project the n-by-k X onto the tangent space of p on the Grassmann M, which is computed by

projp⁑(X)=Xβˆ’ppHX,\operatorname{proj_p}(X) = X - pp^{\mathrm{H}}X,

where β‹…Hβ‹…^{\mathrm{H}} denotes the complex conjugate transposed or Hermitian.

source
ManifoldsBase.representation_size β€” Method
representation_size(M::Grassmann)

Return the representation size or matrix dimension of a point on the Grassmann M, i.e. (n,k)(n,k) for both the real-valued and the complex value case.

source
ManifoldsBase.retract β€” Method
retract(M::Grassmann, p, X, ::PolarRetraction)

Compute the SVD-based retraction PolarRetraction on the Grassmann M. With USV=p+XUSV = p + X the retraction reads

retr⁑pX=UVH,\operatorname{retr}_p X = UV^\mathrm{H},

where β‹…Hβ‹…^{\mathrm{H}} denotes the complex conjugate transposed or Hermitian.

source
ManifoldsBase.retract β€” Method
retract(M::Grassmann, p, X, ::QRRetraction )

Compute the QR-based retraction QRRetraction on the Grassmann M. With QR=p+XQR = p + X the retraction reads

retr⁑pX=QD,\operatorname{retr}_p X = QD,

where D is a mΓ—nmΓ—n matrix with

D=diag⁑(sgn⁑(Rii+12)i=1n).D = \operatorname{diag}\left( \operatorname{sgn}\left(R_{ii}+\frac{1}{2}\right)_{i=1}^n \right).

source
ManifoldsBase.riemann_tensor β€” Method
riemann_tensor(::Grassmann{<:Any,ℝ}, p, X, Y, Z)

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

R(X,Y)Z=(XYTβˆ’YXT)Z+Z(YTXβˆ’XTY).R(X,Y)Z = (XY^\mathrm{T} - YX^\mathrm{T})Z + Z(Y^\mathrm{T}X - X^\mathrm{T}Y).

source
ManifoldsBase.vector_transport_to β€” Method
vector_transport_to(M::Grassmann, p, X, q, ::ProjectionTransport)

compute the projection based transport on the Grassmann M by interpreting X from the tangent space at p as a point in the embedding and projecting it onto the tangent space at q.

source
ManifoldsBase.zero_vector β€” Method
zero_vector(M::Grassmann, p)

Return the zero tangent vector from the tangent space at p on the Grassmann M, which is given by a zero matrix the same size as p.

source

The Grassmannian represented as projectors

Manifolds.ProjectorPoint β€” Type
ProjectorPoint <: AbstractManifoldPoint

A type to represent points on a manifold Grassmann that are orthogonal projectors, i.e. a matrix p∈Fn,np ∈ \mathbb F^{n,n} projecting onto a kk-dimensional subspace.

source
Base.exp β€” Method
exp(M::Grassmann, p::ProjectorPoint, X::ProjectorTangentVector)

Compute the exponential map on the Grassmann as

exp⁑pX=Exp⁑([X,p])pExp⁑(βˆ’[X,p]), \exp_pX = \operatorname{Exp}([X,p])p\operatorname{Exp}(-[X,p]),

where Exp⁑\operatorname{Exp} denotes the matrix exponential and [A,B]=ABβˆ’BA[A,B] = AB-BA denotes the matrix commutator.

For details, see Proposition 3.2 in [BZA20].

source
Manifolds.horizontal_lift β€” Method
horizontal_lift(N::Stiefel{n,k}, q, X::ProjectorTangentVector)

Compute the horizontal lift of X from the tangent space at p=Ο€(q)p=Ο€(q) on the Grassmann manifold, i.e.

Y=Xq∈TqSt(n,k)Y = Xq ∈ T_q\mathrm{St}(n,k)

source
ManifoldsBase.check_point β€” Method
check_point(::Grassmann, p::ProjectorPoint; kwargs...)

Check whether an orthogonal projector is a point from the Grassmann(n,k) manifold, i.e. the ProjectorPoint p∈FnΓ—np ∈ \mathbb F^{nΓ—n}, F∈{R,C}\mathbb F ∈ \{\mathbb R, \mathbb C\} has to fulfill pT=pp^{\mathrm{T}} = p, p2=pp^2=p, and `\operatorname{rank} p = k.

source
ManifoldsBase.check_vector β€” Method
check_vector(::Grassmann, p::ProjectorPoint, X::ProjectorTangentVector; kwargs...)

Check whether the ProjectorTangentVector X is from the tangent space TpGr⁑(n,k)T_p\operatorname{Gr}(n,k) at the ProjectorPoint p on the Grassmann manifold Gr⁑(n,k)\operatorname{Gr}(n,k). This means that X has to be symmetric and that

Xp+pX=XXp + pX = X

must hold, where the kwargs can be used to check both for symmetrix of XX` and this equality up to a certain tolerance.

source
ManifoldsBase.parallel_transport_direction β€” Method
parallel_transport_direction(
    M::Grassmann,
    p::ProjectorPoint,
    X::ProjectorTangentVector,
    d::ProjectorTangentVector
)

Compute the parallel transport of X from the tangent space at p into direction d, i.e. to q=exp⁑pdq=\exp_pd. The formula is given in Proposition 3.5 of [BZA20] as

Pq←p(X)=Exp⁑([d,p])XExp⁑(βˆ’[d,p]),\mathcal{P}_{q ← p}(X) = \operatorname{Exp}([d,p])X\operatorname{Exp}(-[d,p]),

where Exp⁑\operatorname{Exp} denotes the matrix exponential and [A,B]=ABβˆ’BA[A,B] = AB-BA denotes the matrix commutator.

source

Literature

[AMS08]
P.-A. Absil, R. Mahony and R. Sepulchre. Optimization Algorithms on Matrix Manifolds (Princeton University Press, 2008), available online at press.princeton.edu/chapters/absil/.
[BZA20]
T. Bendokat, R. Zimmermann and P.-A. Absil. A Grassmann Manifold Handbook: Basic Geometry and Computational Aspects, arXiv Preprint (2020), arXiv:2011.13699.
[Ngu23]
D. Nguyen. Operator-Valued Formulas for Riemannian Gradient and Hessian and Families of Tractable Metrics in Riemannian Optimization. Journal of Optimization Theory and Applications 198, 135–164 (2023), arXiv:2009.10159.
[Ren11]
Q. Rentmeesters. A gradient method for geodesic data fitting on some symmetric Riemannian manifolds. In: IEEE Conference on Decision and Control and European Control Conference (2011); pp. 7141–7146.