Grassmannian manifold

Manifolds.GrassmannType
Grassmann{T,𝔽} <: AbstractDecoratorManifold{𝔽}

The Grassmann manifold $\mathrm{Gr}(n,k)$ consists of all subspaces spanned by $k$ linear independent vectors $𝔽^n$, where $𝔽 ∈ \{ℝ, ℂ\}$ is either the real- (or complex-) valued vectors. This yields all $k$-dimensional subspaces of $ℝ^n$ for the real-valued case and all $2k$-dimensional subspaces of $ℂ^n$ for the second.

The manifold can be represented as

\[\mathrm{Gr}(n,k) := \bigl\{ \operatorname{span}(p) : p ∈ 𝔽^{n×k}, p^\mathrm{H}p = I_k\},\]

where $⋅^{\mathrm{H}}$ denotes the complex conjugate transpose or Hermitian and $I_k$ is the $k×k$ identity matrix. This means, that the columns of $p$ form an unitary basis of the subspace, that is a point on $\operatorname{Gr}(n,k)$, and hence the subspace can actually be represented by a whole equivalence class of representers. Another interpretation is, that

\[\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 $\operatorname{O}(k)$ of orthogonal $k×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) $p$ is given by

\[T_p\mathrm{Gr}(n,k) = \bigl\{ X ∈ 𝔽^{n×k} : X^{\mathrm{H}}p + p^{\mathrm{H}}X = 0_{k} \bigr\},\]

where $0_k$ is the $k×k$ zero matrix.

Note that a point $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 $p$ also lead to different representation matrices for the tangent space $T_p\mathrm{Gr}(n,k)$

For a representation of points as orthogonal projectors. Here

\[\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

\[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 ProjectorTVector.

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 $\operatorname{Gr}(n,k)$, where the real-valued case field=ℝ is the default.

source
Base.convertMethod
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

\[ π^{\mathrm{SG}}(p) = pp^{\mathrm{T}}.\]

source
Base.convertMethod
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

\[ π^{\mathrm{SG}}(p) = pp^{\mathrm{T}}.\]

source
ManifoldsBase.change_metricMethod
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_representerMethod
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_methodMethod
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

The Grassmanian represented as points on the Stiefel manifold

Manifolds.StiefelPointType
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.StiefelTVectorType
StiefelTVector <: TVector

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.expMethod
exp(M::Grassmann, p, X)

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

\[z = p V\cos(S)V^\mathrm{H} + U\sin(S)V^\mathrm{H},\]

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

\[\exp_p X = Q.\]

source
Base.logMethod
log(M::Grassmann, p, q)

Compute the logarithmic map on the Grassmann M $= \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_p q = V⋅ \operatorname{atan}(S) ⋅ U^\mathrm{H},\]

where $⋅^{\mathrm{H}}$ denotes the complex conjugate transposed or Hermitian. The matrices $U$ and $V$ are the unitary matrices, and $S$ is the diagonal matrix containing the singular values of the SVD-decomposition

\[USV = (q^\mathrm{H}p)^{-1} ( q^\mathrm{H} - q^\mathrm{H}pp^\mathrm{H}).\]

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

source
Base.randMethod
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 $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_HessianMethod
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$ in their formula. Let $\nabla f(p)$ denote the Euclidean gradient G, $\nabla^2 f(p)[X]$ the Euclidean Hessian H. Then the formula reads

\[ \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.distanceMethod
distance(M::Grassmann, p, q)

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

The distance is given by

\[d_{\mathrm{Gr}(n,k)}(p,q) = \operatorname{norm}(\log_p(q)).\]

source
ManifoldsBase.innerMethod
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

\[g_p(X,Y) = \operatorname{tr}(X^{\mathrm{H}}Y),\]

where $⋅^{\mathrm{H}}$ denotes the complex conjugate transposed or Hermitian.

source
ManifoldsBase.inverse_retractMethod
inverse_retract(M::Grassmann, p, q, ::PolarInverseRetraction)

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

\[\operatorname{retr}_p^{-1}q = q*(p^\mathrm{H}q)^{-1} - p,\]

where $⋅^{\mathrm{H}}$ denotes the complex conjugate transposed or Hermitian.

source
ManifoldsBase.inverse_retractMethod
inverse_retract(M, p, q, ::QRInverseRetraction)

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

\[\operatorname{retr}_p^{-1}q = q(p^\mathrm{H}q)^{-1} - p,\]

where $⋅^{\mathrm{H}}$ denotes the complex conjugate transposed or Hermitian.

source
ManifoldsBase.parallel_transport_directionMethod
parallel_transport_direction(M::Grassmann, p, X, Y)

Compute the parallel transport of $X \in T_p\mathcal M$ along the geodesic starting in direction $\dot γ (0) = Y$.

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

\[\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 $S$ are meant to be elementwise

source
ManifoldsBase.projectMethod
project(M::Grassmann, p)

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

source
ManifoldsBase.projectMethod
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

\[\operatorname{proj_p}(X) = X - pp^{\mathrm{H}}X,\]

where $⋅^{\mathrm{H}}$ denotes the complex conjugate transposed or Hermitian.

source
ManifoldsBase.retractMethod
retract(M::Grassmann, p, X, ::PolarRetraction)

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

\[\operatorname{retr}_p X = UV^\mathrm{H},\]

where $⋅^{\mathrm{H}}$ denotes the complex conjugate transposed or Hermitian.

source
ManifoldsBase.retractMethod
retract(M::Grassmann, p, X, ::QRRetraction )

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

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

where D is a $m×n$ matrix with

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

source
ManifoldsBase.riemann_tensorMethod
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 = (XY^\mathrm{T} - YX^\mathrm{T})Z + Z(Y^\mathrm{T}X - X^\mathrm{T}Y).\]

source
ManifoldsBase.vector_transport_toMethod
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_vectorMethod
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.ProjectorPointType
ProjectorPoint <: AbstractManifoldPoint

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

source
Base.expMethod
exp(M::Grassmann, p::ProjectorPoint, X::ProjectorTVector)

Compute the exponential map on the Grassmann as

\[ \exp_pX = \operatorname{Exp}([X,p])p\operatorname{Exp}(-[X,p]),\]

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

For details, see Proposition 3.2 in [BZA20].

source
Manifolds.horizontal_liftMethod
horizontal_lift(N::Stiefel{n,k}, q, X::ProjectorTVector)

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

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

source
ManifoldsBase.check_pointMethod
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 ∈ \mathbb F^{n×n}$, $\mathbb F ∈ \{\mathbb R, \mathbb C\}$ has to fulfill $p^{\mathrm{T}} = p$, $p^2=p$, and `\operatorname{rank} p = k.

source
ManifoldsBase.check_vectorMethod
check_vector(::Grassmann, p::ProjectorPoint, X::ProjectorTVector; kwargs...)

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

\[Xp + pX = X\]

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

source
ManifoldsBase.parallel_transport_directionMethod
parallel_transport_direction(
    M::Grassmann,
    p::ProjectorPoint,
    X::ProjectorTVector,
    d::ProjectorTVector
)

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

\[\mathcal{P}_{q ← p}(X) = \operatorname{Exp}([d,p])X\operatorname{Exp}(-[d,p]),\]

where $\operatorname{Exp}$ denotes the matrix exponential and $[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.