# Grassmannian manifold

`Manifolds.Grassmann`

β Type`Grassmann{n,k,π½} <: AbstractDecoratorManifold{π½}`

The Grassmann manifold $\operatorname{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

\[\operatorname{Gr}(n,k) := \bigl\{ \operatorname{span}(p) : p β π½^{n Γ k}, p^\mathrm{H}p = I_k\},\]

where $\cdot^{\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

\[\operatorname{Gr}(n,k) = \operatorname{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.

The tangent space at a point (subspace) $x$ 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^{[BendokatZimmermannAbsil2020]}.

**Constructor**

`Grassmann(n,k,field=β)`

Generate the Grassmann manifold $\operatorname{Gr}(n,k)$, where the real-valued case `field = β`

is the default.

`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

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

`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

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

`Manifolds.get_total_space`

β Method`get_total_space(::Grassmann{n,k})`

Return the total space of the `Grassmann`

manifold, which is the corresponding Stiefel manifold, independent of whether the points are represented already in the total space or as `ProjectorPoint`

s.

`ManifoldsBase.injectivity_radius`

β Method```
injectivity_radius(M::Grassmann)
injectivity_radius(M::Grassmann, p)
```

Return the injectivity radius on the `Grassmann`

`M`

, which is $\frac{Ο}{2}$.

`ManifoldsBase.is_flat`

β Method`is_flat(M::Grassmann)`

Return true if `Grassmann`

`M`

is one-dimensional.

`ManifoldsBase.manifold_dimension`

β Method`manifold_dimension(M::Grassmann)`

Return the dimension of the `Grassmann(n,k,π½)`

manifold `M`

, i.e.

\[\dim \operatorname{Gr}(n,k) = k(n-k) \dim_β π½,\]

where $\dim_β π½$ is the `real_dimension`

of `π½`

.

`Statistics.mean`

β Method```
mean(
M::Grassmann,
x::AbstractVector,
[w::AbstractWeights,]
method = GeodesicInterpolationWithinRadius(Ο/4);
kwargs...,
)
```

Compute the Riemannian `mean`

of `x`

using `GeodesicInterpolationWithinRadius`

.

## 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.

`Manifolds.StiefelTVector`

β Type`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.

`Base.exp`

β Method`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 $\cdot^{\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.\]

`Base.log`

β Method`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\cdot \operatorname{atan}(S) \cdot U^\mathrm{H},\]

where $\cdot^{\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.

`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 $T_p\mathrm{Gr}(n,k)$ with mean zero and standard deviation `Ο`

by projecting a random Matrix onto the tangent space at `vector_at`

.

`Manifolds.uniform_distribution`

β Method`uniform_distribution(M::Grassmann{n,k,β}, p)`

Uniform distribution on given (real-valued) `Grassmann`

`M`

. Specifically, this is the normalized Haar 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.2(iii) in ^{[Chikuse2003]}.

`ManifoldsBase.distance`

β Method`distance(M::Grassmann, p, q)`

Compute the Riemannian distance on `Grassmann`

manifold `M`

$= \mathrm{Gr}(n,k)$.

Let $USV = p^\mathrm{H}q$ denote the SVD decomposition of $p^\mathrm{H}q$, where $\cdot^{\mathrm{H}}$ denotes the complex conjugate transposed or Hermitian. Then the distance is given by

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

where

\[b_{i}=\begin{cases} 0 & \text{if} \; S_i β₯ 1\\ \arccos(S_i) & \, \text{if} \; S_i<1. \end{cases}\]

`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

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

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

`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.,

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

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

`ManifoldsBase.inverse_retract`

β Method`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 $\cdot^{\mathrm{H}}$ denotes the complex conjugate transposed or Hermitian.

`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

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

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

`ManifoldsBase.representation_size`

β Method`representation_size(M::Grassmann{n,k})`

Return the represenation size or matrix dimension of a point on the `Grassmann`

`M`

, i.e. $(n,k)$ for both the real-valued and the complex value case.

`ManifoldsBase.retract`

β Method`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 $\cdot^{\mathrm{H}}$ denotes the complex conjugate transposed or Hermitian.

`ManifoldsBase.retract`

β Method`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).\]

`ManifoldsBase.riemann_tensor`

β Method`riemann_tensor(::Grassmann{n,k,β}, p, X, Y, Z) where {n,k}`

Compute the value of Riemann tensor on the real `Grassmann`

manifold. The formula reads^{[Rentmeesters2011]} $R(X,Y)Z = (XY^\mathrm{T} - YX^\mathrm{T})Z + Z(Y^\mathrm{T}X - X^\mathrm{T}Y)$.

`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.

`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`

.

## 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 β \mathbb F^{n,n}$ projecting onto a $k$-dimensional subspace.

`Manifolds.ProjectorTVector`

β Type`ProjectorTVector <: TVector`

A type to represent tangent vectors to points on a `Grassmann`

manifold that are orthogonal projectors.

`Base.exp`

β Method`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 ^{[BendokatZimmermannAbsil2020]}.

`Manifolds.canonical_project!`

β Method`canonical_project!(M::Grassmann{n,k}, q::ProjectorPoint, p)`

Compute the canonical projection $Ο(p)$ from the `Stiefel`

manifold onto the `Grassmann`

manifold when represented as `ProjectorPoint`

, i.e.

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

`Manifolds.differential_canonical_project!`

β Method`canonical_project!(M::Grassmann{n,k}, q::ProjectorPoint, p)`

Compute the canonical projection $Ο(p)$ from the `Stiefel`

manifold onto the `Grassmann`

manifold when represented as `ProjectorPoint`

, i.e.

\[ DΟ^{\mathrm{SG}}(p)[X] = Xp^{\mathrm{T}} + pX^{\mathrm{T}}\]

`Manifolds.horizontal_lift`

β Method`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)\]

`ManifoldsBase.check_point`

β Method`check_point(::Grassmann{n,k}, 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`

.

`ManifoldsBase.check_size`

β Method`check_size(M::Grassmann{n,k,π½}, p::ProjectorPoint; kwargs...) where {n,k}`

Check that the `ProjectorPoint`

is of correct size, i.e. from $\mathbb F^{nΓn}$

`ManifoldsBase.check_vector`

β Method`check_vector(::Grassmann{n,k,π½}, p::ProjectorPoint, X::ProjectorTVector; kwargs...) where {n,k,π½}`

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.

`ManifoldsBase.get_embedding`

β Method`get_embedding(M::Grassmann{n,k,π½}, p::ProjectorPoint) where {n,k,π½}`

Return the embedding of the `ProjectorPoint`

representation of the `Grassmann`

manifold, i.e. the Euclidean space $\mathbb F^{nΓn}$.

`ManifoldsBase.parallel_transport_direction`

β Method```
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 ^{[BendokatZimmermannAbsil2020]} 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.

`ManifoldsBase.representation_size`

β Method`representation_size(M::Grassmann{n,k}, p::ProjectorPoint)`

Return the represenation size or matrix dimension of a point on the `Grassmann`

`M`

when using `ProjectorPoint`

s, i.e. $(n,n)$.

- BendokatZimmermannAbsil2020
T. Bendokat, R. Zimmermann, and P. -A. Absil:

*A Grassmann Manifold Handbook: Basic Geometry and Computational Aspects*, arXiv preprint 2011.13699, 2020. - Chikuse2003
Y. Chikuse: "Statistics on Special Manifolds", Springer New York, 2003, doi: 10.1007/978-0-387-21540-2.

- Rentmeesters2011
Q. Rentmeesters, βA gradient method for geodesic data fitting on some symmetric Riemannian manifolds,β in 2011 50th IEEE Conference on Decision and Control and European Control Conference, Dec. 2011, pp. 7141β7146. doi: 10.1109/CDC.2011.6161280.