Nested Grassmannians for Dimensionality Reduction with Applications

Yang, Chun-Hao; Vemuri, Baba C.

doi:https://doi.org/10.59275/j.melba.2022-234f

1 Introduction

Riemannian manifolds are often used to model the sample space in which features derived from the raw data encountered in many medical imaging applications live. Common examples include the diffusion tensors (DTs) in diffusion tensor imaging (DTI) (Basser et al., 1994), the ensemble average propagator (EAP) (Callaghan, 1993). Both DTs and EAP are used to capture the diffusional properties of water molecules in the central nervous system by non-invasively imaging the tissue via diffusion weighted magnetic resonance imaging. In DTI, diffusion at a voxel is captured by a DT, which is a $3\times 3$ symmetric positive-definite matrix, whereas EAP is a probability distribution characterizing the local diffusion at a voxel, which can be parametrized as a point on the Hilbert sphere. Another example is the shape space used to represent shapes in shape analysis. There are many ways to represent a shape, and the most simple one is to use landmarks. For the landmark-based representation, the shape space is called Kendall’s shape space (Kendall, 1984). Kendall’s shape space is in general a stratified space (Goresky and MacPherson, 1988; Feragen et al., 2014), but for the special case of planar shapes, the shape space is the complex projective space, which is a complex Grassmann manifold. The examples mentioned above are often high-dimensional: a DTI scan usually contains half a million DTs; the shape of the Corpus Callosum (which is used in our experiments) is represented by a several hundreds of boundary points in $\mathbb{R}^{2}$ . Thus, in these cases, dimension reduction techniques, if applied appropriately, can benefit the subsequent statistical analysis.

For data on Riemannian manifolds, the most widely used dimensionality reduction method is the principal geodesic analysis (PGA) (Fletcher et al., 2004), which generalizes the principal component analysis (PCA) to manifold-valued data. In fact, there are many variants of PGA. Fletcher et al. (2004) proposed to find the geodesic submanifold of a certain dimension that maximizes the projected variance and computationally, it was achieved via a linear approximation, i.e., applying PCA on the tangent space at the intrinsic mean. This is sometimes referred to as the tangent PCA. Note that this approximation requires the data to be clustered around the intrinsic mean, otherwise the tangent space approximation to the manifold leads to inaccuracies. Later on, Sommer et al. (2010) proposed the Exact PGA (EPGA), which does not use any linear approximation. However, EPGA is computationally expensive as it requires two non-linear optimizations steps per iteration (projection to the geodesic submanifold and finding the new geodesic direction such that the loss of information is minimized). Chakraborty et al. (2016) partially solved this problem for manifolds with constant sectional curvature (spheres and hyperbolic spaces) by deriving closed form formulae for the projection. Other variants of PGA include but are not limited to sparse exact PGA (Banerjee et al., 2017), geodesic PCA (Huckemann et al., 2010), and probabilistic PGA (Zhang and Fletcher, 2013). All these methods focus on projecting data to a geodesic submanifold as in PCA where one projects data to a vector subspace. Instead, one can also project data to a submanifold that minimizes the reconstruction error without any further restrictions, e.g. being geodesic. This is the generalization of the principal curve (Hastie and Stuetzle, 1989) to Riemannian manifolds presented in Hauberg (2016).

Another feature of PCA is that it produces a sequence of nested vector subspaces. From this observation, Jung et al. (2012) proposed the principal nested spheres (PNS) by embedding an $n-1$ -sphere into an $n$ -sphere, and the embedding is not necessarily isometric. Hence PNS is more general than PGA in that PNS is not required to be geodesic. Similarly, for the manifold $P_{n}$ of $n\times n$ SPD matrices, Harandi et al. (2018) proposed a geometry-aware dimension reduction by projecting data in $P_{n}$ to $P_{m}$ for some $m\ll n$ . They also applied the nested $P_{n}$ for the supervised dimensionality reduction problem. Damon and Marron (2014) considered a nested sequence of relations which determine a nested sequence of submanifolds that are not necessarily geodesic. They showed various examples, including Euclidean space and the $n$ -sphere, depicting how the nested relations generalized PCA and PNS. However, for an arbitrary Riemannian manifold, it is not clear how to construct a nested submanifold. Another generalization of PGA was proposed by Pennec et al. (2018), called the exponential barycentric subspace (EBS). A $k$ -dimensional EBS is defined as the locus of weighted exponential barycenters of $(k+1)$ affinely independent reference points. The EBSs are naturally nested by removing or adding reference points.

Unlike PGA which can be applied to any Riemannian manifolds, the construction of the nested manifolds relies heavily on the geometry of the specific manifold, and there is no general principle for such a construction. All the examples described above (spheres and $P_{n}$ ) and several others such as the Grassmannian, Stiefel etc. are homogeneous Riemannian manifolds (Helgason, 1979), which intuitively means that all points on the manifold ’look’ the same. In this work, we propose a general framework for constructing a nested sequence of homogeneous Riemannian manifolds, and, via some simple algebraic computation, show that the nested sphere and the nested $P_{n}$ can indeed be derived within this framework. We will then apply this framework to the Grassmann manifolds, called the nested Grassmann manifolds (NG). The Grassmann manifold $\text{Gr}(p,\mathbb{V})$ is the manifold of all $p$ -dimensional subspaces of the vector space $\mathbb{V}$ where $1\leq p\leq\dim\mathbb{V}$ . Usually $\mathbb{V}=\mathbb{R}^{n}$ or $\mathbb{V}=\mathbb{C}^{n}$ . An important example is Kendall’s shape space of 2D shapes. The space of all shapes determined by $k$ landmarks in $\mathbb{R}^{2}$ is denoted by $\Sigma^{k}_{2}$ , and Kendall (1984) showed that it is isomorphic to the complex projective space $\mathbb{C}P^{k-2}\cong\text{Gr}(1,\mathbb{C}^{k-1})$ . In many applications, the number $k$ of landmarks is large, and so is the dimension of $\text{Gr}(1,\mathbb{C}^{k-1})$ . The core of the proposed dimensionality reduction involves projecting data on $\text{Gr}(p,\mathbb{V})$ to $\text{Gr}(p,\widetilde{\mathbb{V}})$ with $\dim\widetilde{\mathbb{V}}\ll\dim\mathbb{V}$ . The main contributions of this work are as follows: (i) We propose a general framework for constructing a nested sequence of homogeneous Riemannian manifolds unifying the recently proposed nested spheres (Jung et al., 2012) and nested SPD manifolds (Harandi et al., 2018). (ii) We present novel dimensionality reduction techniques based on the concept of NG in both supervised and unsupervised settings respectively. (iii) We demonstrate via several simulation studies and real data experiments, that the proposed NG can achieve a higher ratio of expressed variance compared to PGA.

The rest of the paper is organized as follows. In Section 2, we briefly review the definition of homogeneous Riemannian manifolds and present the recipe for the construction of nested homogeneous Riemannian manifolds. In Section 3, we first review the geometry of the Grassmannian. By applying the procedure developed in Section 2, we present the nested Grassmann manifolds representation and discuss some of its properties in details. Then we describe algorithms for our unsupervised and supervised dimensionality reduction techniques, called the Principal Nested Grassmanns (PNG), in Section 4. In Section 5, we present several simulation studies and experimental results showing how the PNG technique performs compared to PGA under different settings. Finally, we draw conclusions in Section 6.

2 Nested Homogeneous Spaces

In this section, we introduce the structure of nested homogeneous Riemannian manifolds. A Riemannian manifold $(M,\tau)$ is homogeneous if the group of isometries $G=I(M)$ admitted by the manifold acts transitively on $M$ (Helgason, 1979), i.e., for $x,y\in M$ , there exists $g\in G$ such that $g(x)=y$ . In this case, $M$ can be identified with $G/H$ where $H$ is an isotropy subgroup of $G$ at some point $p\in M$ , i.e. $H=\{g\in G:g(p)=p\}$ . Examples of homogeneous Riemannian manifolds include but are not limited to, the $n$ -spheres $S^{n-1}=\text{SO}(n)/\text{SO}(n-1)$ , the SPD manifolds $P_{n}=\text{GL}(n)/\text{O}(n)$ , the Stiefel manifolds $\text{St}(m,n)=\text{SO}(n)/\text{SO}(n-m)$ , and the Grassmann manifolds $\text{Gr}(p,n)=\text{SO}(n)/S(\text{O}(p)\times\text{O}(n-p))$ .

In this paper, we focus on the case where $G$ is either a real or a complex matrix Lie group, i.e. $G$ is a subgroup of $\text{GL}(n,\mathbb{R})$ or $\text{GL}(n,\mathbb{C})$ . The main idea behind the construction of nested homogeneous spaces is simple: augmenting the matrix in $G$ in an appropriate way. With an embedding of the isometry group $G$ , the embedding of the homogeneous space $G/H$ follows naturally from the quotient structure.

Let $G$ and $\tilde{G}$ be two connected Lie groups such that $\dim G<\dim\tilde{G}$ and $\tilde{\iota}:G\to\tilde{G}$ be an embedding. For a closed connected subgroup $H$ of $G$ , let $\tilde{H}=\tilde{\iota}(H)$ . Since $\tilde{\iota}$ is an embedding, $\tilde{H}$ is also a closed subgroup of $\tilde{G}$ . Now the canonical embedding of $G/H$ in $\tilde{G}/\tilde{H}$ is defined by $\iota(gH)=\tilde{\iota}(g)\tilde{H}$ for $g\in G$ . It is easy to see that $\iota$ is well-defined. Let $g_{1},g_{2}\in G$ be such that $g_{1}=g_{2}h$ for some $h\in H$ . Then

\iota(g_{1}H)=\tilde{\iota}(g_{1})\tilde{H}=\tilde{\iota}(g_{2}h)\tilde{H}=\tilde{\iota}(g_{2})\tilde{\iota}(h)\tilde{H}=\tilde{\iota}(g_{2})\tilde{H}=\iota(g_{2}H).

Now for the homogeneous Riemannian manifolds $(M=G/H,\tau_{1})$ and $(\tilde{M}=\tilde{G}/\tilde{H},\tau_{2})$ , denote the left- $G$ -invariant, right- $H$ -invariant metric on $G$ (resp. left- $\tilde{G}$ -invariant, right- $\tilde{H}$ -invariant metric on $\tilde{G}$ ) by $\bar{\tau}_{1}$ and $\bar{\tau}_{2}$ , respectively (see Cheeger and Ebin (1975, Prop. 3.16(4))).

Proposition 1

If $\tilde{\iota}:G\to\tilde{G}$ is isometric, then so is $\iota:G/H\to\tilde{G}/\tilde{H}$ .

Proof Denote the Riemannian submersions by $\pi:G\to G/H$ and $\tilde{\pi}:\tilde{G}\to\tilde{G}/\tilde{H}$ . Let $X$ and $Y$ be vector fields on $G/H$ and $\bar{X}$ and $\bar{Y}$ be their horizontal lifts respectively, i.e. $d\pi(\bar{X})=X$ and $d\pi(\bar{Y})=Y$ . By the definition of Riemannian submersions, $d\pi$ is isometric on the horizontal spaces, i.e. $\bar{\tau}_{1}(\bar{X},\bar{Y})=\tau_{1}(d\pi(\bar{X}),d\pi(\bar{Y}))=\tau_{1}(X,Y)$ . Since $\tilde{\iota}$ is isometric, we have $\bar{\tau}_{1}(\bar{X},\bar{Y})=\bar{\tau}_{2}(d\tilde{\iota}(\bar{X}),d\tilde{\iota}(\bar{Y}))$ . By the definition of $\iota$ , we also have $\iota\circ\pi=\tilde{\pi}\circ\tilde{\iota}$ , which implies $d\iota\circ d\pi=d\tilde{\pi}\circ d\tilde{\iota}$ . Hence,

	$\displaystyle\tau_{1}(X,Y)$	$\displaystyle=\bar{\tau}_{1}(\bar{X},\bar{Y})=\bar{\tau}_{2}(d\tilde{\iota}(\bar{X}),d\tilde{\iota}(\bar{Y}))$
		$\displaystyle=\tau_{2}((d\tilde{\pi}\circ d\tilde{\iota})(\bar{X}),(d\tilde{\pi}\circ d\tilde{\iota})(\bar{Y}))$
		$\displaystyle=\tau_{2}((d\iota\circ d\pi)(\bar{X}),(d\iota\circ d\pi)(\bar{Y}))$
		$\displaystyle=\tau_{2}(d\iota(X),d\iota(Y))$

where the third equality follows from the isometry of $d\tilde{\pi}$ .

Proposition 1 simply says that if the isometry group is isometrically embedded, then the associated homogeneous Riemannian manifolds will also be isometrically embedded. Conversely, if we have a Riemannian submersion $\tilde{f}:\tilde{G}\to G$ , it can easily be shown that the induced map $f:\tilde{G}/\tilde{H}\to G/H$ would also be a Riemannian submersion where $H=\tilde{f}(\tilde{H})$ . The construction above can be applied to a sequence of homogeneous spaces $\{M_{m}\}_{m=1}^{\infty}$ , i.e. the embedding $\iota_{m}:M_{m}\to M_{m+1}$ can be induced from the embedding of the isometry groups $\tilde{\iota}_{m}:G_{m}\to G_{m+1}$ where $G_{m}=I(M_{m})$ provided $\dim G_{i}<\dim G_{j}$ for $i<j$ . See Figure 1 for the structure of nested homogeneous spaces.

Figure 1: Commutative diagram of the induced embedding for homogeneous spaces.

3 Nested Grassmann Manifolds

In this section, we will apply the theory of nested homogeneous space from the previous section to the Grassmann manifolds. First, we briefly review the geometry of the Grassmann manifolds in Section 3.1. With the theory in Section 2, we derive the nested Grassmann manifolds in Section 3.2, and the derivation for nested spheres and nested SPD manifolds are carried out in Section 3.3.

3.1 The Riemannian Geometry of Grassmann Manifolds

To simplify the notation, we assume $\mathbb{V}=\mathbb{R}^{n}$ and write $\text{Gr}(p,n)\coloneqq\text{Gr}(p,\mathbb{R}^{n})$ . All the results presented in this section can be easily extended to the case of $\mathbb{V}=\mathbb{C}^{n}$ by replacing transposition with conjugate transposition and orthogonal groups with unitary groups. The Grassmann manifold $\text{Gr}(p,n)$ is the manifold of all $p$ -dimensional subspaces of $\mathbb{R}^{n}$ , and for a subspace $\mathcal{X}\in\text{Gr}(p,n)$ , we write $\mathcal{X}=\operatorname{span}(X)$ where the columns of $X$ form an orthonormal basis for $\mathcal{X}$ . The space of all $n\times p$ matrices $X$ such that $X^{T}X=I_{p}$ is called the Stiefel manifold, denoted by $\text{St}(p,n)$ . Special cases of Stiefel manifolds are the Lie group of all orthogonal matrices, $\text{O}(n)=\text{St}(n,n)$ , and the $n$ -sphere, $S^{n-1}=\text{St}(1,n)$ . The Stiefel manifold with the induced Euclidean metric (i.e. for $U,V\in T_{X}\text{St}(p,n)$ , $\langle U,V\rangle_{X}=\operatorname{tr}(U^{T}V)$ ) is a homogeneous Riemannian manifold, $\text{St}(p,n)=\text{SO}(n)/\text{SO}(n-p)$ . A canonical Riemannian metric on the Grassmann manifold can be inherited from the metric on $\text{St}(p,n)$ since it is invariant to the left multiplication by elements of $O(n)$ (Absil et al., 2004; Edelman et al., 1998). The Grassmann manifold with this metric is also homogeneous, $\text{Gr}(p,n)=\text{SO}(n)/S(\text{O}(p)\times\text{O}(n-p))$ .

With this canonical metric on the Grassmann manifolds, the geodesic can be expressed in closed form. Let $\mathcal{X}=\operatorname{span}(X)\in\text{Gr}(p,n)$ where $X\in\text{St}(p,n)$ and $H$ be an $n\times p$ matrix. Then the geodesic $\gamma(t)$ with $\gamma(0)=\mathcal{X}$ and $\gamma^{\prime}(0)=H$ is given by $\gamma_{\mathcal{X},H}(t)=\operatorname{span}(XV\cos\Sigma t+U\sin\Sigma t)$ where $H=U\Sigma V^{T}$ is the compact singular value decomposition (Edelman et al., 1998, Theorem 2.3). The exponential map at $\mathcal{X}$ is a map from $T_{\mathcal{X}}\text{Gr}(p,n)$ to $\text{Gr}(p,n)$ defined by $\text{Exp}_{\mathcal{X}}H=\gamma_{\mathcal{X},H}(1)=\operatorname{span}(XV\cos\Sigma+U\sin\Sigma)$ . If $X^{T}Y$ is invertible, the geodesic distance between $\mathcal{X}=\operatorname{span}(X)$ and $\mathcal{Y}=\operatorname{span}(Y)$ is given by $d_{g}^{2}(\mathcal{X},\mathcal{Y})=\operatorname{tr}\Theta^{2}=\sum_{i=1}^{p}\theta_{i}^{2}$ where $(I-XX^{T})Y(X^{T}Y)^{-1}=U\Sigma V^{T}$ , $U\in\text{St}(p,n)$ , $V\in\text{O}(p)$ , and $\Theta=\tan^{-1}\Sigma$ . The diagonal entries $\theta_{1},\ldots,\theta_{k}$ of $\Theta$ are known as the principal angles between $\mathcal{X}$ and $\mathcal{Y}$ .

Figure 2: Illustration of the embedding of

\text{Gr}(p,m)

in

\text{Gr}(p,n)

parametrized by

A\in\text{St}(m,n)

and

B\in\mathbb{R}^{n\times p}

such that

A^{T}B=0

.

3.2 Embedding of $\text{Gr}(p,m)$ in $\text{Gr}(p,n)$

Let $\mathcal{X}=\operatorname{span}(X)\in\text{Gr}(p,m)$ , $X\in\text{St}(p,m)$ . The map $\iota:\text{Gr}(p,m)\to\text{Gr}(p,n)$ , for $m<n$ , defined by

\iota(\mathcal{X})=\operatorname{span}\left(\left[\begin{array}[]{c}X\\ 0_{(n-m)\times p}\end{array}\right]\right)

is an embedding and it is easy to check that this embedding is isometric (Ye and Lim, 2016, Eq. (8)). However, for the dimensionality reduction problem, the above embedding is insufficient as it is not flexible enough to encompass other possible embeddings. To design flexible embeddings, we apply the framework proposed in Section 2. We consider $M_{m}=\text{Gr}(p,m)$ for which the isometry groups are $G_{m}=\text{SO}(m)$ and $H_{m}=\text{S}(\text{O}(p)\times\text{O}(m-p))$ .

In this paper, we consider the embedding $\tilde{\iota}_{m}:\text{SO}(m)\to\text{SO}(m+1)$ given by,

\displaystyle\tilde{\iota}_{m}(O)=\text{GS}\left(R\left[\begin{array}[]{cc}O&a\\ b^{T}&c\end{array}\right]\right)

(3)

where $O\in\text{SO}(m)$ , $R\in\text{SO}(m+1)$ , $a,b\in\mathbb{R}^{m}$ , $c\in\mathbb{R}$ , $c\neq b^{T}O^{-1}a$ , and $\text{GS}(\cdot)$ is the Gram-Schmidt process. Since the Riemannian submersion $\pi:\text{SO}(m)\to\text{Gr}(p,m)$ is defined by $\pi(O)=\operatorname{span}(O_{p})$ where $O\in\text{SO}(m)$ and $O_{p}$ is the $m\times p$ matrix containing the first $p$ columns of $O$ , the induced embedding $\iota_{m}:\text{Gr}(p,m)\to\text{Gr}(p,m+1)$ is given by,

\iota_{m}(\mathcal{X})=\operatorname{span}\left(R\left[\begin{array}[]{c}X\\ b^{T}\end{array}\right]\right)=\operatorname{span}(\tilde{R}X+vb^{T}),

where $b\in\mathbb{R}^{p}$ , $R\in\text{SO}(m+1)$ , $\tilde{R}$ contains the first $m$ columns of $R$ (which means $\tilde{R}\in\text{St}(m,m+1)$ ), $v$ is the last column of $R$ , and $\mathcal{X}=\operatorname{span}(X)\in\text{Gr}(p,m)$ . It is easy to see that for $R=I$ and $b=0$ , this gives the natural embedding described in Ye and Lim (2016) and at the beginning of this section.

Proposition 2

If $b=0$ , then $\iota_{m}$ is an isometric embedding.

Proof With Proposition 1, it suffices to show that $\tilde{\iota}_{m}$ is isometric when $b=0$ . Note that as $\iota_{m}$ is independent of $a$ and $c$ in the definition of $\tilde{\iota}_{m}$ , we can assume $a=0$ and $c=1$ without loss of generality. If $b=0$ , $\tilde{\iota}_{m}$ simplifies to

\tilde{\iota}_{m}(O)=R\left[\begin{array}[]{cc}O&0\\ 0&1\end{array}\right]

where $R\in\text{SO}(m+1)$ . The Riemannian distance on $\text{SO}(n)$ given the induced Euclidean metric is $d_{\text{SO}}(O_{1},O_{2})=\frac{1}{\sqrt{2}}\|\log O_{1}^{T}O_{2}\|_{F}$ . Then for $O_{1},O_{2}\in\text{SO}(m)$ ,

d_{\text{SO}}(\tilde{\iota}_{m}(O_{1}),\tilde{\iota}_{m}(O_{2}))=\frac{1}{\sqrt{2}}\left\|\log\left(\left[\begin{array}[]{cc}O_{1}^{T}O_{2}&0\\ 0&1\end{array}\right]\right)\right\|_{F}=d_{\text{SO}}(O_{1},O_{2}).

Therefore $\tilde{\iota}_{m}$ is an isometric embedding, and so is $\iota_{m}$ by Proposition 1.

With the embedding $\iota_{m}$ , we can construct the corresponding projection $\pi_{m}:\text{Gr}(p,m+1)\to\text{Gr}(p,m)$ using the following proposition.

Proposition 3

The projection $\pi_{m}:\text{Gr}(p,m+1)\to\text{Gr}(p,m)$ corresponding to $\iota_{m}(\mathcal{X})=\operatorname{span}(\tilde{R}X+vb^{T})$ is given by $\pi_{m}(\mathcal{X})=\operatorname{span}(\tilde{R}^{T}X)$ .

Proof First, let $\mathcal{Y}=\operatorname{span}(Y)\in\text{Gr}(p,m)$ and $\mathcal{X}=\operatorname{span}(X)\in\text{Gr}(p,m+1)$ be such that $\mathcal{X}=\operatorname{span}(\tilde{R}Y+vb^{T})$ . Then $XL=\tilde{R}Y+vb^{T}$ for some $L\in\text{GL}(p)$ . Therefore, $Y=\tilde{R}^{T}(XL-vb^{T})=\tilde{R}^{T}XL$ and $\mathcal{Y}=\operatorname{span}(Y)=\operatorname{span}(\tilde{R}^{T}XL)=\operatorname{span}(\tilde{R}^{T}X)$ . Hence, the projection is given by $\pi_{m}(\mathcal{X})=\operatorname{span}(\tilde{R}^{T}X)$ . This completes the proof.

Note that for $\mathcal{X}=\operatorname{span}(X)\in\text{Gr}(p,m+1)$ , $\iota_{m}(\pi_{m}(\mathcal{X}))=\operatorname{span}(RR^{T}X+vb^{T})=\operatorname{span}((I-vv^{T})X+vb^{T})$ where $v\in\mathbb{R}^{m+1}$ and $\|v\|=1$ . The nested relation can be extended inductively and we refer to this construction as the nested Grassmann structure:

\text{Gr}(p,m)\stackrel{{\scriptstyle\iota_{m}}}{{\hookrightarrow}}\text{Gr}(p,m+1)\stackrel{{\scriptstyle\iota_{m+1}}}{{\hookrightarrow}}\ldots\stackrel{{\scriptstyle\iota_{n-2}}}{{\hookrightarrow}}\text{Gr}(p,n-1)\stackrel{{\scriptstyle\iota_{n-1}}}{{\hookrightarrow}}\text{Gr}(p,n).

Thus the embedding from $\text{Gr}(p,m)$ into $\text{Gr}(p,n)$ can be constructed inductively by $\iota\coloneqq\iota_{n-1}\circ\ldots\circ\iota_{m-1}\circ\iota_{m}$ and similarly for the corresponding projection. The explicit forms of the embedding and the projection are given in the following proposition.

Proposition 4

The embedding of $\text{Gr}(p,m)$ into $\text{Gr}(p,n)$ for $m<n$ is given by $\iota_{A,B}(\mathcal{X})=\operatorname{span}(AX+B)$ where $A\in\text{St}(m,n)$ and $B\in\mathbb{R}^{n\times p}$ such that $A^{T}B=0$ . The corresponding projection from $\text{Gr}(p,n)$ to $\text{Gr}(p,m)$ is given by $\pi_{A}=\operatorname{span}(A^{T}X)$ .

Proof By the definition, $\iota\coloneqq\iota_{n-1}\circ\ldots\circ\iota_{m-1}\circ\iota_{m}$ and thus the embedding $\iota:\text{Gr}(p,m)\to\text{Gr}(p,n)$ can be simplified as

\iota_{A,B}(\mathcal{X})=\operatorname{span}\left(\Bigg{(}\prod_{i=m}^{n-1}R_{i}\Bigg{)}X+\sum_{i=m}^{n-1}\Bigg{(}\prod_{j=i+1}^{n-1}R_{j}\Bigg{)}v_{i}b_{i}^{T}\right)=\operatorname{span}(AX+B)

where $R_{i}\in\text{St}(i,i+1)$ , $v_{i}$ is such that $[R_{i}\;v_{i}]\in\text{O}(i+1)$ , $b_{i}\in\mathbb{R}^{p}$ , $A=R_{n-1}R_{n-2}\cdots R_{m}\in\text{St}(m,n)$ , and $B=\sum_{i=m}^{n-1}\big{(}\prod_{j=i+1}^{n-1}R_{j}\big{)}v_{i}b_{i}^{T}$ is an $n\times p$ matrix. It is easy to see that $A^{T}B=0$ . Similar to Proposition 3, the projection $\pi_{A}:\text{Gr}(p,n)\to\text{Gr}(p,m)$ is then given by $\pi_{A}(\mathcal{X})=\operatorname{span}(A^{T}X)$ . This completes the proof.

From Proposition 2, if $B=0$ then $\iota_{A}$ is an isometric embedding. Hence, our nested Grassmann structure is more flexible than PGA as it allows one to project the data onto a non-geodesic submanifold. An illustration is shown in Figure 2. The results discussed in this subsection can be generalized to any homogeneous space in principle. For a given homogeneous space, once the embedding of the groups of isometries (e.g., Eq. (3)) is determined, the induced embedding and the corresponding projection can be derived akin to the case of Grassmann manifolds.

3.3 Connections to Other Nested Structures

The nested homogeneous spaces proposed in this work (see Figure 1) actually provides a unified framework within which, the nested spheres (Jung et al., 2012) and the nested SPD manifolds (Harandi et al., 2018) are special cases.

The $n$ -Sphere Example: Since the $n$ -sphere can be identified with a homogeneous space $S^{n-1}\cong\text{O}(n)/\text{O}(n-1)$ , with the embedding (3), the induced embedding of $S^{n-1}$ into $S^{n}$ is

\iota(\boldsymbol{x})=\text{GS}\left(R\left[\begin{array}[]{c}x\\ b\end{array}\right]\right)=\frac{1}{\sqrt{1+b^{2}}}R\left[\begin{array}[]{c}x\\ b\end{array}\right]=R\left[\begin{array}[]{c}\sin(r)x\\ \cos(r)\end{array}\right]

where $x\in S^{n-1}$ , $b\in\mathbb{R}$ , and $r=\cos^{-1}\big{(}\frac{b}{\sqrt{1+b^{2}}}\big{)}$ . This is precisely the nested sphere proposed in Jung et al. (2012, Eq. (2)).

The SPD Manifold Example: For the $m$ -dimensional SPD manifold denoted by $P_{m}$ , $G_{m}=\text{GL}(m)$ and $H_{m}=\text{O}(m)$ . Consider the embedding $\tilde{\iota}_{m}:\text{GL}(m)\to\text{GL}(m+1)$ given by

\tilde{A}=\tilde{\iota}_{m}(A)=R\left[\begin{array}[]{cc}A&0\\ 0&1\end{array}\right]R^{T},

where $A\in\text{GL}(m)$ , $R\in\text{O}(m+1)$ and the corresponding projection $\tilde{\pi}_{m}:\text{GL}(m+1)\to\text{GL}(m)$ is

\tilde{\pi}_{m}(\tilde{A})=W^{T}\tilde{A}W

where $W$ contains the first $m$ columns of $R=[W\;v]\in\text{O}(m+1)$ (i.e., $W\in\text{St}(m,m+1)$ and $W^{T}v=0$ ). The submersion $\psi\circ f:\text{GL}(m)\to P_{m}$ is given by $\psi\circ f(A)=A^{T}A$ . Hence the induced embedding $\iota_{m}:P_{m}\to P_{m+1}$ and projection $\pi_{m}:P_{m+1}\to P_{m}$ are

\iota_{m}(X)=WXW^{T}+vv^{T}\text{ and }\pi_{m}(X)=W^{T}XW

which is the projection map used in Harandi et al. (2018, Eq. (13)). However, Harandi et al. (2018) did not perform any embedding or construct a nested family of SPD manifolds. Recently, it came to our attention that Jaquier and Rozo (2020) derived a similar nested family of SPD manifolds based on the projection maps in Harandi et al. (2018) described above.

4 Dimensionality Reduction with Nested Grassmanns

In this section, we discuss how to apply the nested Grassmann structure to the problem of dimension reduction. In Section 4.1 and 4.2, we describe the unsupervised and supervised dimension reduction based on the nested Grassmann manifolds. In Section 4.3, we will discuss the choice of distance metrics required by the dimensionality reduction algorithm and present some technical details regarding the implementation. Analysis of principal nested Grassmann (PNG) will be introduced and discussed in Section 4.4 and Section 4.5.

4.1 Unsupervised Dimensionality Reduction

We can now apply the nested Grassmann structure to the problem of unsupervised dimensionality reduction. Suppose that we are given the points, $\mathcal{X}_{1},\ldots,\mathcal{X}_{N}\in\text{Gr}(p,n)$ . We would like to have lower dimensional representations in $\text{Gr}(p,m)$ for $\mathcal{X}_{1},\ldots,\mathcal{X}_{N}$ with $m\ll n$ . The desired projection map $\pi_{A}$ that we seek is obtained by the minimizing the reconstruction error, i.e.

L_{u}(A,B)=\frac{1}{N}\sum_{i=1}^{N}d^{2}(\mathcal{X}_{i},\hat{\mathcal{X}}_{i})=\frac{1}{N}\sum_{i=1}^{N}d^{2}(\operatorname{span}(X_{i}),\operatorname{span}(AA^{T}X_{i}+B))

where $d(\cdot,\cdot)$ is a distance metric on $\text{Gr}(p,n)$ . It is clear that $L_{u}$ has a $\text{O}(m)$ -symmetry in the first argument, i.e. $L_{u}(AO,B)=L_{u}(A,B)$ for $O\in\text{O}(m)$ . Hence, the optimization is performed over the space $\text{St}(m,n)/\text{O}(m)\cong\text{Gr}(m,n)$ when optimizing with respect to this particular loss function. Now we can apply the Riemannian gradient descent algorithm (Edelman et al., 1998) to obtain $A$ and $B$ by optimizing $L_{u}(A,B)$ over $\operatorname{span}(A)\in\text{Gr}(m,n)$ and $B\in\mathbb{R}^{n\times p}$ such that $A^{T}B=0$ . Note that the restriction $A^{T}B=0$ simply means that the columns of $B$ are in the null space of $A^{T}$ , denoted $N(A^{T})$ . Hence in practice this restriction can be handled as follows. For arbitrary $\tilde{B}\in\mathbb{R}^{n\times p}$ , project $\tilde{B}$ on to $N(A^{T})$ , i.e. $B=P_{N(A^{T})}\tilde{B}$ where $P_{N(A^{T})}=I-AA^{T}$ is the projection from $\mathbb{R}^{n}$ to $N(A^{T})$ . Thus, the loss function can be written as

L_{u}(A,B)=\frac{1}{N}\sum_{i=1}^{N}d^{2}(\operatorname{span}(X_{i}),\operatorname{span}(AA^{T}X_{i}+(I-AA^{T})B))

and it is optimized over $\text{Gr}(m,n)\times\mathbb{R}^{n\times p}$ . Note that since we represent a subspace by its orthonormal basis, when $m>n/2$ , we can use the isomorphism $\text{Gr}(m,n)\cong\text{Gr}(n-m,n)$ to further reduce the computational burden. This will be particularly useful when $m=n-1$ as in Section 4.4. Under this isomorphism $\text{Gr}(m,n)\cong\text{Gr}(n-m,n)$ , the corresponding subspace of $\operatorname{span}(A)\in\text{Gr}(m,n)$ is $\operatorname{span}(A_{\perp})\in\text{Gr}(n-m,n)$ where $A_{\perp}$ is an $n\times(n-m)$ matrix such that $[A\;A_{\perp}]$ is an orthogonal matrix. Hence the loss function $L_{u}$ can alternatively be expressed as

L_{u}(A,B)=\frac{1}{N}\sum_{i=1}^{N}d^{2}(\operatorname{span}(X_{i}),\operatorname{span}((I-A_{\perp}A_{\perp}^{T})X_{i}+A_{\perp}A_{\perp}^{T}B)).

Remark 1

The reduction of the space of all possible projections from $\text{St}(m,n)$ to $\text{Gr}(m,n)$ is a consequence of the choice of the loss function $L_{u}$ . With a different loss function, one might have a different space of possible projections.

4.2 Supervised Dimensionality Reduction

If in addition to $\mathcal{X}_{1},\ldots,\mathcal{X}_{N}\in\text{Gr}(p,n)$ , we are given the associated labels $y_{1},\ldots,y_{N}\in\{1,\ldots,k\}$ , then we would like to use this extra information to sharpen the result of dimensionality reduction. Specifically, we expect that after reducing the dimension, points from the same class preserve their proximity while points from different classes are well separated. We use an affinity function $a:\text{Gr}(p,n)\times\text{Gr}(p,n)\to\mathbb{R}$ to encode the structure of the data as suggested by Harandi et al. (2018, Sec 3.1, Eq. (14)-(16)). For completeness, we repeat the definition of the affinity function here. The affinity function is defined as $a(\mathcal{X}_{i},\mathcal{X}_{j})=g_{w}(\mathcal{X}_{i},\mathcal{X}_{j})-g_{b}(\mathcal{X}_{i},\mathcal{X}_{j})$ where

	$\displaystyle g_{w}(\mathcal{X}_{i},\mathcal{X}_{j})$	$\displaystyle=\begin{cases}1&\text{if }\mathcal{X}_{i}\in N_{w}(\mathcal{X}_{j})\text{ or }\mathcal{X}_{j}\in N_{w}(\mathcal{X}_{i})\\ 0&\text{\text{Otherwise}}\end{cases}$
	$\displaystyle g_{b}(\mathcal{X}_{i},\mathcal{X}_{j})$	$\displaystyle=\begin{cases}1&\text{if }\mathcal{X}_{i}\in N_{b}(\mathcal{X}_{j})\text{ or }\mathcal{X}_{j}\in N_{b}(\mathcal{X}_{i})\\ 0&\text{\text{Otherwise}}\end{cases}.$

The set $N_{w}(\mathcal{X}_{i})$ consists of $\nu_{w}$ nearest neighbors of $\mathcal{X}_{i}$ that have the same labels as $y_{i}$ , and the set $N_{b}(\mathcal{X}_{i})$ consists of $\nu_{b}$ nearest neighbors of $\mathcal{X}_{i}$ that have different labels from $y_{i}$ . The nearest neighbors can be computed using the geodesic distance.

The desired projection map $\pi_{A}$ that we seek is obtained by the minimizing the following loss function

L_{s}(A)=\frac{1}{N^{2}}\sum_{i,j=1}^{N}a(\mathcal{X}_{i},\mathcal{X}_{j})d^{2}(\operatorname{span}(A^{T}X_{i}),\operatorname{span}(A^{T}X_{j}))

where $d$ is a distance metric on $\text{Gr}(p,m)$ . Note that if the distance metric $d$ has $\text{O}(m)$ -symmetry, e.g. the geodesic distance, so does $L_{s}$ . In this case the optimization can be done on $\text{St}(m,n)/\text{O}(m)\cong\text{Gr}(m,n)$ . Otherwise it is on $\text{St}(m,n)$ . This supervised dimensionality reduction is termed as, supervised nested Grassmann (sNG).

4.3 Choice of the distance function

The loss functions $L_{u}$ and $L_{s}$ depend on the choice of the distance $d:\text{Gr}(p,n)\times\text{Gr}(p,n)\to\mathbb{R}_{\geq 0}$ . Besides the geodesic distance, there are many widely used distances on the Grassmann manifold, see, for example, Edelman et al. (1998, p. 337) and Ye and Lim (2016, Table 2). In this work, we use two different distance metrics: (1) the geodesic distance $d_{g}$ and (2) the projection distance, which is also called the chordal distance in Ye and Lim (2016) and the projection $F$ -norm in Edelman et al. (1998). The geodesic distance was defined in Section 3.1 and the projection distance is defined as follows. For $\mathcal{X},\mathcal{Y}\in\text{Gr}(p,n)$ , denote the projection matrices onto $\mathcal{X}$ and $\mathcal{Y}$ by $P_{\mathcal{X}}$ and $P_{\mathcal{Y}}$ respectively. Then, the distance between $\mathcal{X}$ and $\mathcal{Y}$ is given by $d_{p}(\mathcal{X},\mathcal{Y})=\|P_{\mathcal{X}}-P_{\mathcal{Y}}\|_{F}/\sqrt{2}=\big{(}\sum_{i=1}^{p}\sin^{2}\theta_{i}\big{)}^{1/2}$ where $\theta_{1},\ldots,\theta_{p}$ are the principal angles of $\mathcal{X}$ and $\mathcal{Y}$ . If $\mathcal{X}=\operatorname{span}(X)$ , then $P_{\mathcal{X}}=X(X^{T}X)^{-1}X^{T}$ . It is also easy to see the the projection distance has $\text{O}(n)$ -symmetry. We choose the projection distance mainly for its computational efficiency as it involves only matrix multiplication which has a time complexity $O(n^{2})$ whereas the geodesic distance requires an SVD which has a time complexity of $O(n^{3})$ .

4.4 Analysis of Principal Nested Grassmannians

To determine the dimension of the nested submanifold that fits the data well enough, we can choose $p<m_{1}<\ldots<m_{k}<n$ and estimate the projection onto these nested Grassmann manifolds. The ratio of expressed variance for each projection is the ratio of the variance of the projected data and the variance of the original data. With these ratios, we can choose the desired dimension according to the pre-specified percentage of expressed variance as one would do for choosing the number of principal components in PCA.

Alternatively, one can have a full analysis of principal nested Grassmanns (PNG) as follows. Starting from $\text{Gr}(p,n)$ , one can reduce the dimension down to $\text{Gr}(p,p+1)$ . Using the diffeomorphism between $\text{Gr}(p,n)$ and $\text{Gr}(p,n-p)$ , we have $\text{Gr}(p,p+1)\cong\text{Gr}(1,p+1)$ , and hence we can continue reducing the dimension down to $\text{Gr}(1,2)$ . The resulting sequence will be

\text{Gr}(p,n)\to\text{Gr}(p,n-1)\to\cdots\to\text{Gr}(p,p+1)=\text{Gr}(1,p+1)\to\text{Gr}(1,p)\to\cdots\to\text{Gr}(1,2).

Furthermore, we can reduce the points on $\text{Gr}(1,2)$ , which is a 1-dimensional manifold, to a 0-dimensional manifold, which is simply a point, by computing the Fréchet mean (FM). We call this FM the nested Grassmannian mean (NGM) of $\mathcal{X}_{1},\ldots,\mathcal{X}_{N}\in\text{Gr}(p,n)$ . The NGM is unique since $\text{Gr}(1,2)\cong\mathbb{R}P^{1}$ can be identified as the half circle in $\mathbb{R}^{2}$ and the FM is unique in this case. Note that in general, the NGM will not be the same as the FM of $\mathcal{X}_{1},\ldots,\mathcal{X}_{N}$ since the embedding/projection need not be isometric. The supervised PNG (sPNG) can be obtained similarly by replacing each projection with it supervised counterpart.

4.5 Principal Scores

Whenever we apply a projection $\pi_{m}:\text{Gr}(p,m+1)\to\text{Gr}(p,m)$ to the data, we might lose some information contained in the data. More specifically, since we project data on a $p(m+1-p)$ -dimensional manifold to a $p(m-p)$ -dimensional manifold, we need to describe this $p(m+1-p)-p(m-p)=p$ dimensional information loss during the projection. In PCA, this is done by computing the scores of each principal component, which are the transformed coordinates of each sample in the eigenspace of the covariance matrix. We can generalize the notion of principal scores to the nested Grassmanns as follows: For each $\mathcal{X}\in\text{Gr}(p,m+1)$ , denote by $M_{\pi_{m}(\mathcal{X})}$ , the fiber of $\pi_{m}(\mathcal{X})$ , i.e. $M_{\pi_{m}(\mathcal{X})}=\pi_{m}^{-1}(\pi_{m}(\mathcal{X}))=\{\mathcal{Y}\in\text{Gr}(p,m+1):\pi_{m}(\mathcal{Y})=\pi_{m}(\mathcal{X})\}$ , which is a $p$ -dimensional submanifold of $\text{Gr}(p,m+1)$ . An illustration of this fiber is given in Figure 3(a). Let $\widetilde{\mathcal{X}}=\iota_{m}(\pi_{m}(\mathcal{X}))$ and let the unit tangent vector $V\in T_{\widetilde{\mathcal{X}}}M_{\pi_{m}(\mathcal{X})}$ be the geodesic direction from $\widetilde{\mathcal{X}}$ to $\mathcal{X}$ . Given a suitable basis on $T_{\widetilde{\mathcal{X}}}M_{\pi_{m}(\mathcal{X})}$ , $V$ can be realized as a $p$ -dimensional vector, and this will be the score vector of $\mathcal{X}$ associated with the projection $\pi_{m}$ .

Refer to caption — (a) The fiber $M_{\pi_{m}(\mathcal{X})}$ in $\text{Gr}(p,m+1)$ .

By the definition of $M_{\pi_{m}(\mathcal{X})}$ , we have the following decomposition of the tangent space of $\text{Gr}(p,m+1)$ at $\widetilde{\mathcal{X}}$ into the horizontal space and the vertical space induced by $\pi_{m}$ ,

T_{\widetilde{\mathcal{X}}}\text{Gr}(p,m+1)=T_{\widetilde{\mathcal{X}}}M_{\pi_{m}(\mathcal{X})}\oplus(d\iota_{m})_{\pi_{m}(\mathcal{X})}(T_{\pi_{m}(\mathcal{X})}\text{Gr}(p,m)).

An illustration of this decomposition is given in Figure 3(b). A tangent vector of $M_{\pi_{m}(\mathcal{X})}$ at $\widetilde{\mathcal{X}}$ is of the form $\Delta=A_{\perp}b^{T}$ where $A_{\perp}$ is any $(m+1)$ -dim vector such that $[A\;A_{\perp}]$ is orthogonal and $b\in\mathbb{R}^{p}$ . It is easy to check that $\pi_{m}(\operatorname{span}(AA^{T}X+A_{\perp}b^{T}))=\pi_{m}(\operatorname{span}(X))=\operatorname{span}(A^{T}X)$ . Hence a natural coordinate for the tangent vector $\Delta=A_{\perp}b^{T}$ is $b\in\mathbb{R}^{p}$ , and the geodesic direction from $\widetilde{\mathcal{X}}$ to $\mathcal{X}$ would be $V=X^{T}A_{\perp}$ . It is easy to see that $\|V\|_{F}=1$ since $X$ has orthonormal columns. To reflect the distance between $\widetilde{\mathcal{X}}$ and $\mathcal{X}$ , i.e. the reconstruction error, we define $d(\widetilde{\mathcal{X}},\mathcal{X})V$ as the score vector for $\mathcal{X}$ associated with $\pi_{m}$ . In the case of $\text{Gr}(1,2)\to\text{NGM}$ , we use the sign distance to the NGM as the score. For complex nested Grassmanns however, the principal score associated with each projection is a $p$ -dimensional complex vector. For the sake of visualization, we transform this $p$ -dimensional complex vector to a $2p$ -dimensional real vector. The procedure for computing the PNG and the principal scores is summarized in Algorithm 1.

Remark 2

Note that this definition of principal score is not intrinsic as it depends on the choice of basis. Indeed, it is impossible to choose a $p$ -dimensional vector for the projection $\pi_{m}$ in an intrinsic way, since the only property of a map that is independent of bases is the rank of the map. A reasonable choice of basis is made by viewing the Grassmann $\text{Gr}(p,m)$ as a quotient manifold of $\text{St}(p,m)$ , which is a submanifold in $\mathbb{R}^{m\times p}$ . This is how we define the principal score for nested Grassmanns.

Input :

\mathcal{X}_{1},\ldots\mathcal{X}_{N}\in\text{Gr}(p,n)

Output : An

N\times p(n-p)

score matrix.

Let

\mathcal{X}=(\mathcal{X}_{1},\ldots,\mathcal{X}_{N})

.

/*

\text{Gr}(p,n)\to\text{Gr}(p,n-1)\to\cdots\to\text{Gr}(p,p+1)

*/

for $m=n-1,\ldots,p+1$ do

\hat{v},\hat{b}=\operatorname*{arg\,min}_{v,b}\sum^{N}_{i=1}d^{2}(\mathcal{X}_{i},\operatorname{span}((I-vv^{T})X_{i}+vb^{T}))

;

/*

v\in\mathbb{R}^{m+1}

,

\|v\|=1

,

b\in\mathbb{R}^{p}

*/

for $i=1,\ldots,N$ do

Let

\widehat{\mathcal{X}}_{i}=\operatorname{span}((I-\hat{v}\hat{v}^{T})X_{i}+\hat{v}\hat{b}^{T})

where

\mathcal{X}_{i}=\operatorname{span}(X_{i})

.

V_{i,m}=d_{g}(\mathcal{X}_{i},\widehat{\mathcal{X}}_{i})X_{i}^{T}v

;

// the score vector for

\mathcal{X}_{i}

\mathcal{X}_{i}\leftarrow\operatorname{span}(R^{T}X_{i})

;

//

R\in\text{St}(m,m+1)

such that

R^{T}\hat{v}=0

.

if $p>1$ then

for $i=1,\ldots,N$ do

\mathcal{X}_{i}\leftarrow\mathcal{X}_{i}^{\perp}

;

//

\text{Gr}(p,p+1)\cong\text{Gr}(1,p+1)

/*

\text{Gr}(1,p+1)\to\text{Gr}(1,p)\to\cdots\to\text{Gr}(1,2)

*/

for $m=p,\ldots,2$ do

\hat{v},\hat{b}=\operatorname*{arg\,min}_{v,b}\sum^{N}_{i=1}d^{2}(\mathcal{X}_{i},\operatorname{span}((I-vv^{T})X_{i}+vb^{T}))

;

/*

v\in\mathbb{R}^{m+1}

,

\|v\|=1

,

b\in\mathbb{R}^{p}

*/

for $i=1,\ldots,N$ do

Let

\widehat{\mathcal{X}}_{i}=\operatorname{span}((I-\hat{v}\hat{v}^{T})X_{i}+\hat{v}\hat{b}^{T})

where

\mathcal{X}_{i}=\operatorname{span}(X_{i})

.

V_{i,m}=d_{g}(\mathcal{X}_{i},\widehat{\mathcal{X}}_{i})X_{i}^{T}v

;

// the score vector for

\mathcal{X}_{i}

\mathcal{X}_{i}\leftarrow\operatorname{span}(R^{T}X_{i})

;

//

R\in\text{St}(m,m+1)

such that

R^{T}v=0

.

/*

\text{Gr}(1,2)\to\text{NGM}

*/

/* Note that

\widehat{\mathcal{X}}_{i}\in\text{Gr}(1,2)

*/

\text{NGM}\leftarrow\text{FM}(\widehat{\mathcal{X}}_{1},\ldots,\widehat{\mathcal{X}}_{N})

;

Let

U_{i}\in\mathbb{R}^{2}

be the geodesic direction from the NGM to

\mathcal{X}_{i}

such that

\|U_{i}\|=d_{g}(\text{NGM},\mathcal{X}_{i})

.;

for $i=1,\ldots,N$ do

V_{i,1}\in\mathbb{R}

is such that

U_{i}=V_{i,1}\times\frac{U_{1}}{\|U_{1}\|}

;

Return the score matrix

S=[V_{i,m}]

. (Note that

V_{i,m}\in\mathbb{R}

for

m=1,\ldots,p

and

V_{i,m}\in\mathbb{R}^{p}

for

m=p+1,\ldots,n-1

.)

Algorithm 1 Principal Nested Grassmanns

5 Experiments

In this section, we will demonstrate the performance of the proposed dimensionality reduction technique, i.e. PNG and sPNG, via experiments on synthetic and real data. The implementation¹¹1Our code is available at https://github.com/cvgmi/NestedGrassmann. is based on the python library pymanopt (Townsend et al., 2016) and we use the steepest descent algorithm for the optimization (with default parameters in pymanopt). The optimization was performed on a desktop with 3.6GHz Intel i7 processors and took about 30 seconds to converge.

5.1 Synthetic Data

In this subsection, we compare the performance of the projection and the geodesic distances respectively. The questions we will answer are the following. (1) From Section 4.3, we see that using projection distance is more efficient than using the geodesic distance. But how do they perform compared to each other under varying dimension $n$ and variance level $\sigma^{2}$ ? (2) Is our method of dimensionality reduction ”better” than PGA? Under what conditions does our method outperform PGA?

5.1.1 Projection and Geodesic Distance Comparisons

The procedure we used to generate random points on $\text{Gr}(p,n)$ for the synthetic data experiments is as follows: First, we generate $N$ points from a uniform distribution on $\text{St}(p,m)$ (Chikuse, 2003, Ch. 2.5), generate $A$ from the uniform distribution on $\text{St}(m,n)$ , and generate $B$ as an $n\times p$ matrix with i.i.d entries from $N(0,0.1)$ . Then we compute $\widetilde{\mathcal{X}}_{i}=\operatorname{span}(AX_{i}+(I-AA^{T})B)\in\text{Gr}(p,n)$ . Finally, we compute $\mathcal{X}_{i}=\text{Exp}_{\widetilde{\mathcal{X}}_{i}}(\sigma U_{i})$ , where $U_{i}=\widetilde{U}_{i}/\|\widetilde{U}_{i}\|$ and $\widetilde{U}_{i}\in T_{\widetilde{\mathcal{X}}_{i}}\text{Gr}(p,n)$ , to include some perturbation.

This experiment involves comparing the performance of the NG representation in terms of the explained variance, under different levels of data variance. In this experiment, we set $N=50$ , $n=10$ , $m=3$ , and $p=1$ and $\sigma$ is ranging from 0.5 to 1. The results are averaged over $100$ repetitions and are shown in Figure 4. From these results, we can see that the explained variance for the projection distance and the geodesic distance are indistinguishable but using projection distance leads to much faster convergence than when using the geodesic distance. The reason is that when two points on the Grassmann manifold are close, the geodesic distance can be well-approximated by the projection distance. When the algorithm converges, the original point $\mathcal{X}_{i}$ and the reconstructed point $\hat{\mathcal{X}}_{i}$ should be close and the geodesic distance can thus be well-approximated by the projection distance. Therefore, for all the experiments in the next section, we use the projection distance for the sake of efficiency.

5.1.2 Comparison of PNG and PGA

Now we compare the proposed PNG to PGA. In order to have a fair comparison between PNG and PGA, we define the principal components of PNG as the principal components of the scores obtained as in Section 4.5. Similar to the previous experiment, we set $N=50$ , $n=10$ , $m=5$ , $p=2$ , and $\sigma=0.01,0.05,0.1,0.5$ and apply the same procedure to generate synthetic data. The results are averaged over $100$ repetitions and are shown in Figure 5.

From Figure 5, we can see that our method outperforms PGA by virtue of the fact that it is able to capture a larger amount of variance contained in the data. We can see that when the variance is small, the improvement of PNG over PGA is less significant, whereas, our method is significantly better for the large data variance case (e.g. comparing $\sigma=0.5$ and $\sigma=0.01$ ). Note that when the variance in the data is small, i.e. the data are tightly clustered around the FM, and PGA captures the essence of the data well. However, the requirement in PGA on the geodesic submanifold to pass through the anchor point, namely the FM, is not meaningful for data with large variance as explained through the following simple example. Consider, a few data points spread out on the equator of a sphere. The FM in this case is likely to be the north pole of the sphere if we restrict ourselves to the upper hemisphere. Thus, the geodesic submanifold computed by PGA will pass through this FM. However, what is more meaningful is a submanifold corresponding to the equator, which is what a nested spheres representation (Jung et al., 2012) in this case yields. In similar vein, for data with large variance on a Grassmann manifold, our NG representation will yield a more meaningful representation than PGA.

5.2 Application to Planar Shape Analysis

We now apply our method to planar (2-dimensional) shape analysis. A planar shape $\sigma$ can be represented as an ordered set of $k>2$ points in $\mathbb{R}^{2}$ , called a $k$ -ad or a configuration. Here we assume that these $k$ points are not all identical. Denote the configuration by $X$ which is a $k\times 2$ matrix. Let $H$ be the $(k-1)\times k$ Helmert submatrix (Dryden and Mardia, 2016, Ch. 2.5). Then $Z=HX/\|HX\|_{F}$ is called the pre-shape of $X$ from which the information about translation and scaling is removed. The space of all pre-shapes is called the pre-shape space, denoted by $\mathcal{S}^{k}_{2}$ . By definition, the pre-shape space is a $(2k-3)$ -dimensional sphere. The shape is obtained by removing the rotation from the pre-shape, and thus the shape space is $\Sigma^{k}_{2}=\mathcal{S}^{k}_{2}/\text{O}(2)$ . It was shown by Kendall (1984) that $\Sigma^{k}_{2}$ is a smooth manifold and, when equipped with the quotient metric, is isometric to the complex projective space $\mathbb{C}P^{k-2}$ equipped with the Fubini-Study metric (up to a scale factor) which is a special case of the complex Grassmannians, i.e. $\mathbb{C}P^{k-2}\cong\text{Gr}(1,\mathbb{C}^{k-1})$ . Hence, we can apply the proposed PNG to planar shapes. For planar shapes, we also compare with the recently proposed principal nested shape spaces (PNSS) (Dryden et al., 2019), which is an application of PNS on the pre-shape space. We will now demonstrate how the PNG performs compared to PGA and PNSS using some simple examples of planar shapes and the OASIS dataset.

Examples of Planar Shapes

We take three datasets: digit3, gorf, and gorm, from the R package shapes (Dryden, 2021). The digit3 dataset consists of 30 shapes of the digit 3, each of which is represented by 13 points in $\mathbb{R}^{2}$ ; the gorf dataset consists of 30 shapes of female gorilla skull, each of which is represented by 8 points in $\mathbb{R}^{2}$ ; the gorm dataset consists of 29 shapes of male gorilla skull, each of which is represented by 8 points in $\mathbb{R}^{2}$ . Example shapes from these three datasets are shown in Figure 6. The cumulative ratios of variance explained by the first 5 principal components²²2Here the principal components in PNG and PGA are complex whereas the principal components in PNSS are real. Hence, we choose 5 principal components in PNG and PGA and 10 principal components in PNSS, so that the reduced (real) dimension is 10 in all three cases. of PNG, PGA, and PNSS are shown in Figure 7. It can be seen from Figure 7 that the proposed PNG achieves higher explained variance than PGA and PNSS respectively in most cases.

OASIS Corpus Callosum Data Experiment

The OASIS database (Marcus et al., 2007) is a publicly available database that contains T1-MR brain scans of subjects of age ranging from 18 to 96. In particular, it includes subjects that are clinically diagnosed with mild to moderate Alzheimer’s disease. We further classify them into three groups: young (aged between 10 and 40), middle-aged (aged between 40 and 70), and old (aged above 70). For demonstration, we randomly choose 4 brain scans within each decade, totalling 36 brain scans. From each scan, the Corpus Callosum (CC) region is segmented and 250 points are taken on the boundary of the CC region. See Figure 8 for samples of the segmented corpus callosi. In this case, the shape space is $\Sigma^{248}_{2}\cong\mathbb{C}P^{248}\cong\mathsf{Gr}(1,\mathbb{C}^{249})$ . Results of application of the three methods to this data are shown in Figure 9.

Since the data are divided into three groups (young, middle-aged, and old), we can apply the sPNG described in Section 4.2 to reduce the dimension. The purpose of this experiment is not to demonstrate state-of-the-art classification accuracy for this dataset. Instead, our goal here is to demonstrate that the proposed nested Grassmann representation in a supervised setting is much more discriminative than the competition, namely the supervised PGA. Hence, we choose a simple classifier such as the support vector machine (SVM) Vapnik (1995) to highlight the aforementioned discriminative power of the nested Grassmann over PGA.

For comparison, the PGA can be easily extended to supervised PGA (sPGA) by first diffeomorphically mapping all the data to the tangent space anchored at the FM and then performing supervised PCA Bair et al. (2006); Barshan et al. (2011) on the tangent space. However, generalizing PNSS to the supervised case is nontrivial and is beyond the scope of this paper. Therefore, we limit our comparison to the unsupervised PNSS. In this demonstration, we apply an SVM on the scores obtained from different dimension reduction algorithms, and we choose only the first three principal scores to show that even with the 3-dimensional representation of the original shapes, we can still achieve good classification results. The results are shown in Table 1. These results are in accordance with our expectation since in sPNG, we seek a projection that minimizes the within-group variance while maximizing the between-group variance. However, as we observed earlier, the constraint of requiring the geodesic submanifold to pass through the FM is not well suited for this dataset which has a large variance across the data. This accounts for why the sPNG exhibits far superior performance compared to sPGA in accuracy.

Method	Accuracy
sPNG	83.33%
PNG	75%
sPGA	66.67%
PGA	63.89%
PNSS	80.56%

Table 1: Classification accuracies for sPGA and sPNG respectively.

6 Conclusion

In this work, we proposed a novel nested geometric structure for homogeneous spaces and used this structure to achieve dimensionality reduction for data residing in Grassmann manifolds. We also discuss how this nested geometric structure served as a natural generalization of other existing nested geometric structures in literature namely, spheres and the manifold of SPD matrices. Specifically, we showed that a lower dimensional Grassmann manifold can be embedded into a higher dimensional Grassmann manifold and via this embedding we constructed a sequence of nested Grassmann manifolds. Compared to the PGA, which is designed for general Riemannian manifolds, the proposed method can capture a higher percentage of data variance after reducing the dimensionality. This is primarily because our method, unlike the PGA, does not require the submanifold to be a geodesic submanifold and to pass through the Fréchet mean of the data. Succinctly, the nested Grassmann structure allows us to fit the data to a larger class of submanifolds than PGA. We also proposed a supervised dimensionality reduction technique which simultaneously differentiates data classes while reducing dimensionality. Efficacy of our method was demonstrated on the OASIS Corpus Callosi data for dimensionality reduction and classification. We showed that our method outperforms the widely used PGA and the recently proposed PNSS by a large margin.

Acknowledgments

This research was in part funded by the NSF grant IIS-1724174, the NIH NINDS and NIA via R01NS121099 to Vemuri and the MOST grant 110-2118-M-002-005-MY3 to Yang.

Ethical Standards

The work follows appropriate ethical standards in conducting research and writing the manuscript, following all applicable laws and regulations regarding treatment of animals or human subjects.

Conflicts of Interest

We declare we don’t have conflicts of interest.

References

Absil et al. (2004) P-A Absil, Robert Mahony, and Rodolphe Sepulchre. Riemannian geometry of Grassmann manifolds with a view on algorithmic computation. Acta Applicandae Mathematica, 80(2):199–220, 2004.
Bair et al. (2006) Eric Bair, Trevor Hastie, Debashis Paul, and Robert Tibshirani. Prediction by supervised principal components. Journal of the American Statistical Association, 101(473):119–137, 2006.
Banerjee et al. (2017) Monami Banerjee, Rudrasis Chakraborty, and Baba C Vemuri. Sparse exact PGA on Riemannian manifolds. In Proceedings of the IEEE International Conference on Computer Vision, pages 5010–5018, 2017.
Barshan et al. (2011) Elnaz Barshan, Ali Ghodsi, Zohreh Azimifar, and Mansoor Zolghadri Jahromi. Supervised principal component analysis: Visualization, classification and regression on subspaces and submanifolds. Pattern Recognition, 44(7):1357–1371, 2011.
Basser et al. (1994) Peter J Basser, James Mattiello, and Denis LeBihan. MR diffusion tensor spectroscopy and imaging. Biophysical journal, 66(1):259–267, 1994.
Callaghan (1993) Paul T Callaghan. Principles of Nuclear Magnetic Resonance Microscopy. Oxford University Press on Demand, 1993.
Chakraborty et al. (2016) Rudrasis Chakraborty, Dohyung Seo, and Baba C Vemuri. An efficient exact-PGA algorithm for constant curvature manifolds. In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 3976–3984, 2016.
Cheeger and Ebin (1975) Jeff Cheeger and David Gregory Ebin. Comparison theorems in Riemannian geometry, volume 9. North-holland Amsterdam, 1975.
Chikuse (2003) Yasuko Chikuse. Statistics on Special Manifolds, volume 174. Springer Science & Business Media, 2003.
Damon and Marron (2014) James Damon and JS Marron. Backwards principal component analysis and principal nested relations. Journal of Mathematical Imaging and Vision, 50(1-2):107–114, 2014.
Dryden (2021) I. L. Dryden. shapes package. R Foundation for Statistical Computing, Vienna, Austria, 2021. URL http://www.R-project.org. Contributed package, Version 1.2.6.
Dryden and Mardia (2016) Ian L Dryden and Kanti V Mardia. Statistical shape analysis: with applications in R, volume 995. John Wiley & Sons, 2016.
Dryden et al. (2019) Ian L Dryden, Kwang-Rae Kim, Charles A Laughton, Huiling Le, et al. Principal nested shape space analysis of molecular dynamics data. Annals of Applied Statistics, 13(4):2213–2234, 2019.
Edelman et al. (1998) Alan Edelman, Tomás A Arias, and Steven T Smith. The geometry of algorithms with orthogonality constraints. SIAM Journal on Matrix Analysis and Applications, 20(2):303–353, 1998.
Feragen et al. (2014) Aasa Feragen, Mads Nielsen, Eva Bjørn Vedel Jensen, Andrew du Plessis, and François Lauze. Geometry and statistics: Manifolds and stratified spaces. Journal of Mathematical Imaging and Vision, 50(1):1–4, 2014.
Fletcher et al. (2004) P Thomas Fletcher, Conglin Lu, Stephen M Pizer, and Sarang Joshi. Principal geodesic analysis for the study of nonlinear statistics of shape. IEEE Transactions on Medical Imaging, 23(8):995–1005, 2004.
Goresky and MacPherson (1988) Mark Goresky and Robert MacPherson. Stratified Morse Theory. Springer, 1988.
Harandi et al. (2018) Mehrtash Harandi, Mathieu Salzmann, and Richard Hartley. Dimensionality reduction on SPD manifolds: The emergence of geometry-aware methods. IEEE Transactions on Pattern Analysis and Machine Intelligence, 40(1):48–62, 2018.
Hastie and Stuetzle (1989) Trevor Hastie and Werner Stuetzle. Principal curves. Journal of the American Statistical Association, 84(406):502–516, 1989.
Hauberg (2016) Søren Hauberg. Principal curves on Riemannian manifolds. IEEE Transactions on Pattern Analysis and Machine Intelligence, 38(9):1915–1921, 2016.
Helgason (1979) Sigurdur Helgason. Differential geometry, Lie groups, and symmetric spaces. Academic press, 1979.
Huckemann et al. (2010) Stephan Huckemann, Thomas Hotz, and Axel Munk. Intrinsic shape analysis: Geodesic PCA for Riemannian manifolds modulo isometric Lie group actions. Statistica Sinica, pages 1–58, 2010.
Jaquier and Rozo (2020) Noémie Jaquier and Leonel Rozo. High-dimensional bayesian optimization via nested riemannian manifolds. Advances in Neural Information Processing Systems, 33:20939–20951, 2020.
Jung et al. (2012) Sungkyu Jung, Ian L Dryden, and JS Marron. Analysis of principal nested spheres. Biometrika, 99(3):551–568, 2012.
Kendall (1984) David G Kendall. Shape manifolds, Procrustean metrics, and complex projective spaces. Bulletin of the London Mathematical Society, 16(2):81–121, 1984.
Marcus et al. (2007) Daniel S Marcus, Tracy H Wang, Jamie Parker, John G Csernansky, John C Morris, and Randy L Buckner. Open Access Series of Imaging Studies (OASIS): cross-sectional MRI data in young, middle aged, nondemented, and demented older adults. Journal of Cognitive Neuroscience, 19(9):1498–1507, 2007.
Pennec et al. (2018) Xavier Pennec et al. Barycentric subspace analysis on manifolds. The Annals of Statistics, 46(6A):2711–2746, 2018.
Sommer et al. (2010) Stefan Sommer, François Lauze, Søren Hauberg, and Mads Nielsen. Manifold valued statistics, exact principal geodesic analysis and the effect of linear approximations. In European Conference on Computer Vision, pages 43–56. Springer, 2010.
Townsend et al. (2016) James Townsend, Niklas Koep, and Sebastian Weichwald. Pymanopt: A python toolbox for optimization on manifolds using automatic differentiation. Journal of Machine Learning Research, 17(137):1–5, 2016. URL http://jmlr.org/papers/v17/16-177.html.
Vapnik (1995) Vladimir Vapnik. The Nature of Statistical Learning Theory. Springer Science & Business Media, 1995.
Ye and Lim (2016) Ke Ye and Lek-Heng Lim. Schubert varieties and distances between subspaces of different dimensions. SIAM Journal on Matrix Analysis and Applications, 37(3):1176–1197, 2016.
Zhang and Fletcher (2013) Miaomiao Zhang and Tom Fletcher. Probabilistic principal geodesic analysis. In Advances in Neural Information Processing Systems, pages 1178–1186, 2013.

[bib.bib1] Absil et al. (2004) P-A Absil, Robert Mahony, and Rodolphe Sepulchre. Riemannian geometry of Grassmann manifolds with a view on algorithmic computation. Acta Applicandae Mathematica, 80(2):199–220, 2004.

[bib.bib2] Bair et al. (2006) Eric Bair, Trevor Hastie, Debashis Paul, and Robert Tibshirani. Prediction by supervised principal components. Journal of the American Statistical Association, 101(473):119–137, 2006.

[bib.bib3] Banerjee et al. (2017) Monami Banerjee, Rudrasis Chakraborty, and Baba C Vemuri. Sparse exact PGA on Riemannian manifolds. In Proceedings of the IEEE International Conference on Computer Vision, pages 5010–5018, 2017.

[bib.bib4] Barshan et al. (2011) Elnaz Barshan, Ali Ghodsi, Zohreh Azimifar, and Mansoor Zolghadri Jahromi. Supervised principal component analysis: Visualization, classification and regression on subspaces and submanifolds. Pattern Recognition, 44(7):1357–1371, 2011.

[bib.bib5] Basser et al. (1994) Peter J Basser, James Mattiello, and Denis LeBihan. MR diffusion tensor spectroscopy and imaging. Biophysical journal, 66(1):259–267, 1994.

[bib.bib6] Callaghan (1993) Paul T Callaghan. Principles of Nuclear Magnetic Resonance Microscopy. Oxford University Press on Demand, 1993.

[bib.bib7] Chakraborty et al. (2016) Rudrasis Chakraborty, Dohyung Seo, and Baba C Vemuri. An efficient exact-PGA algorithm for constant curvature manifolds. In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 3976–3984, 2016.

[bib.bib8] Cheeger and Ebin (1975) Jeff Cheeger and David Gregory Ebin. Comparison theorems in Riemannian geometry, volume 9. North-holland Amsterdam, 1975.

[bib.bib9] Chikuse (2003) Yasuko Chikuse. Statistics on Special Manifolds, volume 174. Springer Science & Business Media, 2003.

[bib.bib10] Damon and Marron (2014) James Damon and JS Marron. Backwards principal component analysis and principal nested relations. Journal of Mathematical Imaging and Vision, 50(1-2):107–114, 2014.

[bib.bib11] Dryden (2021) I. L. Dryden. shapes package. R Foundation for Statistical Computing, Vienna, Austria, 2021. URL http://www.R-project.org. Contributed package, Version 1.2.6.

[bib.bib12] Dryden and Mardia (2016) Ian L Dryden and Kanti V Mardia. Statistical shape analysis: with applications in R, volume 995. John Wiley & Sons, 2016.

[bib.bib13] Dryden et al. (2019) Ian L Dryden, Kwang-Rae Kim, Charles A Laughton, Huiling Le, et al. Principal nested shape space analysis of molecular dynamics data. Annals of Applied Statistics, 13(4):2213–2234, 2019.

[bib.bib14] Edelman et al. (1998) Alan Edelman, Tomás A Arias, and Steven T Smith. The geometry of algorithms with orthogonality constraints. SIAM Journal on Matrix Analysis and Applications, 20(2):303–353, 1998.

[bib.bib15] Feragen et al. (2014) Aasa Feragen, Mads Nielsen, Eva Bjørn Vedel Jensen, Andrew du Plessis, and François Lauze. Geometry and statistics: Manifolds and stratified spaces. Journal of Mathematical Imaging and Vision, 50(1):1–4, 2014.

[bib.bib16] Fletcher et al. (2004) P Thomas Fletcher, Conglin Lu, Stephen M Pizer, and Sarang Joshi. Principal geodesic analysis for the study of nonlinear statistics of shape. IEEE Transactions on Medical Imaging, 23(8):995–1005, 2004.

[bib.bib17] Goresky and MacPherson (1988) Mark Goresky and Robert MacPherson. Stratified Morse Theory. Springer, 1988.

[bib.bib18] Harandi et al. (2018) Mehrtash Harandi, Mathieu Salzmann, and Richard Hartley. Dimensionality reduction on SPD manifolds: The emergence of geometry-aware methods. IEEE Transactions on Pattern Analysis and Machine Intelligence, 40(1):48–62, 2018.

[bib.bib19] Hastie and Stuetzle (1989) Trevor Hastie and Werner Stuetzle. Principal curves. Journal of the American Statistical Association, 84(406):502–516, 1989.

[bib.bib20] Hauberg (2016) Søren Hauberg. Principal curves on Riemannian manifolds. IEEE Transactions on Pattern Analysis and Machine Intelligence, 38(9):1915–1921, 2016.

[bib.bib21] Helgason (1979) Sigurdur Helgason. Differential geometry, Lie groups, and symmetric spaces. Academic press, 1979.

[bib.bib22] Huckemann et al. (2010) Stephan Huckemann, Thomas Hotz, and Axel Munk. Intrinsic shape analysis: Geodesic PCA for Riemannian manifolds modulo isometric Lie group actions. Statistica Sinica, pages 1–58, 2010.

[bib.bib23] Jaquier and Rozo (2020) Noémie Jaquier and Leonel Rozo. High-dimensional bayesian optimization via nested riemannian manifolds. Advances in Neural Information Processing Systems, 33:20939–20951, 2020.

[bib.bib24] Jung et al. (2012) Sungkyu Jung, Ian L Dryden, and JS Marron. Analysis of principal nested spheres. Biometrika, 99(3):551–568, 2012.

[bib.bib25] Kendall (1984) David G Kendall. Shape manifolds, Procrustean metrics, and complex projective spaces. Bulletin of the London Mathematical Society, 16(2):81–121, 1984.

[bib.bib26] Marcus et al. (2007) Daniel S Marcus, Tracy H Wang, Jamie Parker, John G Csernansky, John C Morris, and Randy L Buckner. Open Access Series of Imaging Studies (OASIS): cross-sectional MRI data in young, middle aged, nondemented, and demented older adults. Journal of Cognitive Neuroscience, 19(9):1498–1507, 2007.

[bib.bib27] Pennec et al. (2018) Xavier Pennec et al. Barycentric subspace analysis on manifolds. The Annals of Statistics, 46(6A):2711–2746, 2018.