Riemannian Multinomial Logistics Regression for SPD Neural Networks
Deep neural networks for learning Symmetric Positive Definite (SPD) matrices
are gaining increasing attention in machine learning. Despite the significant
progress, most existing SPD networks use traditional Euclidean classifiers on
an approximated space rather than intrinsic classifiers that accurately capture
the geometry of SPD manifolds. Inspired by Hyperbolic Neural Networks (HNNs),
we propose Riemannian Multinomial Logistics Regression (RMLR) for the
classification layers in SPD networks. We introduce a unified framework for
building Riemannian classifiers under the metrics pulled back from the
Euclidean space, and showcase our framework under the parameterized
Log-Euclidean Metric (LEM) and Log-Cholesky Metric (LCM). Besides, our
framework offers a novel intrinsic explanation for the most popular LogEig
classifier in existing SPD networks. The effectiveness of our method is
demonstrated in three applications: radar recognition, human action
recognition, and electroencephalography (EEG) classification. The code is
available at https://github.com/GitZH-Chen/SPDMLR.git.
Ziheng Chen, Yue Song, Gaowen Liu, Ramana Rao Kompella, Xiaojun Wu, Nicu Sebe
# Riemannian Multiclass Logistics Regression for SPD Neural Networks

###### Abstract

Deep neural networks for learning symmetric positive definite (SPD) matrices are gaining increasing attention in machine learning. Despite the significant progress, most existing SPD networks use traditional Euclidean classifiers on approximated spaces rather than intrinsic classifiers that accurately capture the geometry of SPD manifolds. Inspired by the success of hyperbolic neural networks (HNNs), we propose Riemannian multiclass logistics regression (RMLR) for SPD networks. We introduce a general unified framework for a family of Riemannian metrics on SPD manifolds and showcase the specific \(\mathrm{O}(n)\)-invariant Log-Euclidean Metrics for SPD networks. Moreover, we encompass the most popular classifier in existing SPD networks as a special case of our framework. Extensive experiments on popular SPD learning benchmarks demonstrate the superiority of our classifiers.

## 1 Introduction

Symmetric Positive Definite (SPD) matrices are commonly encountered in a diverse range of scientific fields, such as medical imaging [1; 2], signal processing [3; 4; 5; 6], elasticity [7; 8], question answering [9; 10], and computer vision [11; 12; 13; 14; 15; 14; 16; 17; 18; 19; 20]. Despite their ubiquitous presence, traditional learning algorithms are ineffective in handling the non-Euclidean geometry of SPD matrices. To address this limitation, researchers have proposed different manifold learning techniques based on Riemannian geometry, including Affine-Invariant Metric (AIM), Log-Euclidean Metric (LEM) [21], Bures-Wasserstein Metric (BWM) [22], and Log-Cholesky Metric (LCM) [23]. More recently, Chen _et al._[24] and Thanwerdas _et al._[25] developed two different generalizations of LEM and introduced Adaptive Log-Euclidean Metrics (ALEMs) and \(\mathrm{O}(n)\)-Invariant Log-Euclidean Metrics (OILEMs), respectively.

Inspired by the great success of deep learning [26; 27; 28], several deep networks have been developed on SPD manifolds, exhibiting promising performance in various machine learning applications [11; 1; 13; 6; 2; 15; 14; 17; 18; 29; 30]. However, their classification layers are usually designed in an approximated space or are constrained in a special type of manifold-valued data. For instance, the most commonly used approach is a traditional Euclidean classifier in the tangent space at the identity matrix, realized by stacking matrix logarithm, fully connected (FC) layer, and softmax layer [11; 6; 31; 32; 17; 29; 30]. Other approaches have also been explored, such as the Euclidean classifier in the local coordinate domain [1] and direct classification on SPD matrices [15]. However, all the above three types of classifiers rely on some approximation spaces or tricks that are not intrinsic. More recently, Chakraborty _et al._[2] introduced an invariant layer for manifold-valued data that mimics the invariant FC layer in CNNs which is more intrinsic compared to the previous classifiers, as no approximation spaces or tricks are required. However, it is designed for gridded
manifold-valued data, which is not the main type of data encountered in many other SPD networks. Following the convention of most existing SPD networks, we also only focus on non-gridded cases.

In this paper, we address the issue of ill-intrinsic classifiers by designing Riemannian classifiers based on the geometric reinterpretation of Euclidean classifiers, which has been successful in the design of classifiers on hyperbolic neural networks (HNNs) [33]. One main obstacle is that the formulae of Riemannian classifiers can vary under different Riemannian metrics. Fortunately, we find that on SPD manifolds, several Riemannian metrics can be unified discussed. We first introduce our general Riemannian multiclass logistics regression (RMLR) in this paper, which incorporates different cases under several Riemannian metrics. We then proceed to discuss specific cases under OILEMs, which consist of a family of metrics and can be viewed as parameterized LEM. Our general RMLR allows us to propose unified classifiers for all variants of OILEMs, suiting datasets of different characteristics. To learn the SPD parameters in our RMLR, we discuss in detail three kinds of Riemannian Stochastic Gradient Descent (RSGD) based on OILEMs, AIM, and BWM, respectively. Our framework also allows for an intrinsic explanation for the commonly used Euclidean classifier on SPD manifolds which consists of successive matrix logarithm, FC, and softmax layers. Finally, extensive experiments demonstrate that our proposed Riemannian classifiers exhibit consistent performance gains across widely-used SPD benchmarks. The **contributions** of our work are summarized as follows:

* We are the **first** to introduce a general RMLR framework on SPD manifolds, design specific RMLRs under all OILEMs for SPD networks, and discuss different optimization strategies for training the SPD parameters.
* Our framework encompasses the most popular classifier in SPD networks (_i.e._,stacking matrix logarithm, FC layer, and softmax) as a special case and gives an intrinsic explanation.
* Extensive experiments on widely used SPD learning benchmarks demonstrate the superiority of our proposed classifiers over the previous baselines.

## 2 Geometry of SPD manifolds

The set of SPD matrices, denoted as \(\mathcal{S}_{++}^{n}\), constitutes a smooth manifold known as the SPD manifold [21; 34]. Endowed with a Riemannian metric, it forms a Riemannian manifold [21; 34]. Several Riemannian metrics have yielded impressive results in machine learning, namely LEM [21], AIM [34], LCM [23], BWM [22], and ALEMs [24]. Recently, the authors of [25] have generalized LEM and AIM into two-parameter \(\mathrm{O}(n)\)-invariant Log-Euclidean Metrics (OILEMs) and \((\alpha,\beta)\)-AIM. This section briefly reviews LEM and OILEMs.

Literally, OILEMs enjoy \(\mathrm{O}(n)\)-invariance, which generalizes LEM by two parameters,

\[g_{S}^{\mathrm{LE}(\alpha,\beta)}(V,V)=\alpha\|\mathrm{mlog}_{\bullet,S}(V)\|_ {\mathrm{F}}+\beta\operatorname{tr}\left(S^{-1}V\right)^{2}, \tag{1}\]

where \(S\in\mathcal{S}_{++}^{n},V\in T_{S}\mathcal{S}_{++}^{n},\alpha>0,\alpha+n\beta>0\), \(\mathrm{mlog}_{\bullet,S}\) denotes the differentail map of matrix logarithm at \(S\), and \(\|\cdot\|_{\mathrm{F}}\) is the standard Frobenius norm. When \(\alpha=1\) and \(\beta=0\), OILEM is reduced back to LEM, which is known for its fast and simple computation. Intuitively, an OILEM is a two-parameter variant of LEM. For this reason, we simply use \((\alpha,\beta)\)-LEM to represent OILEM.

As shown in [25], \((\alpha,\beta)\)-LEM is actually the pullback metric from the standard LEM by

\[f_{p,q}:S\in\mathcal{S}_{++}^{n}\longmapsto\mathrm{mexp}\left(F_{p,q}(\mathrm{ mlog}S)\right)=\det(S)^{\frac{p-q}{n}}S^{q}\in\mathcal{S}_{++}^{n}, \tag{2}\]

where \(\mathrm{mexp}\) is matrix exponential, and \(F_{p,q}(X)=qX+\frac{p-q}{n}\operatorname{tr}(X)I_{n}\) with \(p=\sqrt{\alpha+n\beta}\) and \(q=\sqrt{\alpha}\). With varying \(\alpha,\beta\), OILEMs constitute a family of Riemannian metrics.

We recall an excerpt from Theorem 4.2 of [24], which generally characterize a huge family of Riemannian metrics on SPD manifolds

**Theorem 1** (Pullback Euclidean Metrics (PEMs)).: _Let \(S_{1},S_{2}\in\mathcal{S}_{++}^{n}\), \(\phi:\mathcal{S}_{++}^{n}\rightarrow\mathcal{S}^{n}\) is a diffeomorphism. We define the following operations,_

\[S_{1}\odot_{\phi}S_{2} =\phi^{-1}(\phi(S_{1})+\phi(S_{2})), \tag{3}\] \[g_{S}^{\phi}(V_{1},V_{2}) =\langle\phi_{*,S}(V_{1}),\phi_{*,S}(V_{2})\rangle,\forall S\in \mathcal{S}_{++}^{n},\forall V_{i}\in T_{S}\mathcal{S}_{++}^{n}, \tag{4}\]

_where \(\phi_{*,S}:T_{S}\mathcal{S}_{++}^{n}\to T_{\phi(S)}\mathcal{S}^{n}\) is the differential map of \(\phi\) at \(S\), and \(<\cdot,\cdot>\) is a Euclidean metric. Then, we have the following conclusions: \(\{\mathcal{S}_{++}^{n},\odot_{\phi}\}\) is an abelian Lie group, \(\{\mathcal{S}_{++}^{n},g^{\phi}\}\) is
a Riemannian manifold, and \(g^{\phi}\) is a bi-invariant metric, called Pullback Euclidean Metric (**PEM**). The associated geodesic distance is_

\[d^{\phi}(S_{1},S_{2})=\|\phi(S_{1})-\phi(S_{2})\|, \tag{5}\]

_where \(\|\cdot\|\) is the norm induced by \(\langle\cdot,\cdot\rangle\). The Riemannian operators are as follows_

\[\operatorname{Exp}_{S_{1}}V =\phi^{-1}(\phi(S_{1})+\phi_{*,S_{1}}V), \tag{6}\] \[\operatorname{Log}_{S_{1}}S_{2} =\phi_{*,\phi(S_{1})}^{-1}(\phi(S_{2})-\phi(S_{1})),\] (7) \[\Gamma_{S_{1}\to S_{2}}(V) =\phi_{*,\phi(S_{2})}^{-1}\circ\phi_{*,S_{1}}(V), \tag{8}\]

_where \(V\in T_{S_{1}}\mathcal{S}_{++}^{n}\), \(\operatorname{Exp}\), \(\operatorname{Log}\), and \(\Gamma\) are Riemannian exponential map, logarithmic map, and parallel transportation respectively, and \(\phi_{*}\) & \(\phi_{*}^{-1}\) are the differential maps of \(\phi\) & \(\phi^{-1}\)._

As discussed in [24], both LEM and LCM belongs to PEMs. In the following section, we will further present that all OILEMs actually belong to PEMs. Besides, Theorem 1 also allows us to exempt from various concrete computations. We present the Venn diagram in Figure 1 to illustrate the relationship of some popular Riemannian metrics on SPD manifolds.

_Remark 1_.: The identity element of the Lie group \(\{\mathcal{S}_{++}^{n},\odot_{\phi}\}\) induced by \(\phi\) is not necessarily the identity matrix by definition. However, the current Lie groups on SPD manifolds [21, 23, 35, 24] all incorporate the identity matrix as the identity element. Therefore, by abuse of notation, we use \(I\) to represent the identity matrix or identity element alternatively, according to the context. Besides, the inner product \(\langle\cdot,\cdot\rangle\) in Eq. (4) does not need to be the standard Euclidean inner product. However, without loss of generality, \(\langle\cdot,\cdot\rangle\) is assumed to be the standard one, as \(n\)-dimensional Euclidean spaces are naturally linearly isometric.

## 3 Riemannian multiclass logistic regression under general PEMs

In this section, we first reformulate the Euclidean multiclass logistic regression (MLR). Then we proceed to deal with the general RMLR under arbitrary PEM on SPD manifolds.

### Reformulation of Euclidean MLR

In HNNs [33], hyperbolic MLR is designed based on the reformulation of Euclidean MLR from the perspective of distances to margin hyperplanes. Lebanon _et al._[36] first showcase this framework on spherical geometry. We now briefly review the reformulation of Euclidean MLR and then present our RMLR on SPD manifolds.

Given \(C\) classes, each margin hyperplane in \(\mathbb{R}^{n}\) can be represented by softmax probabilities:

\[\forall k\in\{1,\dots,C\},\quad p(y=k\mid x)\propto\exp\left((\langle a_{k}, x\rangle-b_{k})\right),\quad\text{ where }b_{k}\in\mathbb{R},x,a_{k}\in\mathbb{R}^{n}. \tag{9}\]

Every hyperplane in \(\mathbb{R}^{n}\) can be parameterized by a normal vector \(a\) and a scalar shift \(b\):

\[H_{a,b}=\{x\in\mathbb{R}^{n}:\langle a,x\rangle-b=0\},\quad\text{ where }a\in\mathbb{R}^{n}\backslash\{\mathbf{0}\},\text{ and }b\in\mathbb{R}. \tag{10}\]

As in [36, Sec. 5] and [33, Sec. 3.1], by a proper \(p\in\mathbb{R}^{n}\), we have \(\langle a,x\rangle-b=\langle a,x-p\rangle\). Hyperplane \(H_{a,b}\) can be reformulated as

\[H_{a,p}=\{x\in\mathbb{R}^{n}:\langle a,x-p\rangle=0\},\quad\text{ where }a\in\mathbb{R}^{n}\backslash\{\mathbf{0}\},\text{ and }p\in\mathbb{R}^{n}. \tag{11}\]

In view of \(\langle a,x-p\rangle=\operatorname{sign}(\langle a,x-p\rangle)\|a\|d(x,H_{a,p})\), Eq. (9) can be rewritten as :

\[p(y=k\mid x)\propto\exp(\operatorname{sign}(\langle a_{k},x-p_{k}\rangle)\|a_ {k}\|d(x,H_{a_{k},p_{k}})),p_{k},x\in\mathbb{R}^{n},\text{ and }a_{k}\in\mathbb{R}^{n}\backslash\{\mathbf{0}\}. \tag{12}\]

In geometry, \(\operatorname{Log}_{p}x\) is the natural generalization of the directional vector \(\vec{p}=x-p\) starting at \(p\) and ending at \(x\). The inner product can also be replaced by the Riemannian metric at \(p\). More detail can be found in [34, Table 1]. Therefore, the hyperplane in Eq. (11) and MLR in Eq. (12) can be readily generalized into SPD manifolds \(\{\mathcal{S}_{++}^{n},g\}\).

Figure 1: Conceptual illustration of some popular Riemannian metrics on SPD manifolds.
**Definition 3.1** (SPD hyperplanes).: Given \(P\in\mathcal{S}^{n}_{++},A\in T_{P}\mathcal{S}^{n}_{++}\backslash\{\mathbf{0}\}\), we define the SPD hyperplane as

\[\tilde{H}_{A,P}=\{S\in\mathcal{S}^{n}_{++}:g_{P}(\operatorname{Log}_{P}S,A)= \langle\operatorname{Log}_{P}S,A\rangle_{P}=0\} \tag{13}\]

**Definition 3.2** (SPD multiclass logistics regression).: SPD multiclass logistics regression is defined as

\[p(y=k\mid S)\propto\exp(\operatorname{sign}(\langle A_{k},\operatorname{Log} _{P_{k}}(S)\rangle_{P_{k}})\|A_{k}\|_{P_{k}}d(S,\tilde{H}_{A_{k},P_{k}})), \tag{14}\]

where \(P_{k}\in\mathcal{S}^{n}_{++}\), \(A_{k}\in T_{P_{k}}\mathcal{S}^{n}_{++}\backslash\{\mathbf{0}\}\), \(\langle\cdot,\cdot\rangle_{P_{k}}=g_{P_{k}}\), and \(\|\|_{P_{k}}\) is the norm on \(T_{P_{k}}\mathcal{S}^{n}_{++}\) induced by \(g\) at \(P_{k}\), and \(\tilde{H}_{A_{k},P_{k}}\) is a margin hyperplane in \(\mathcal{S}^{n}_{++}\) as defined in Eq. (13). \(d(S,\tilde{H}_{A_{k},P_{k}})\) denotes the distance between \(S\) and SPD hyperplane \(\tilde{H}_{A_{k},P_{k}}\), which is formulated as:

\[d(S,\tilde{H}_{A_{k},P_{k}}))=\inf_{Q\in H_{A_{k},P_{k}}}d(S,Q), \tag{15}\]

where \(d(S,Q)\) is the geodesic distance induced by \(g\).

_Remark 2_.: Note that since \(g\) could be any kind of the existing Riemannian metrics on SPD manifolds, the specific formulae of Eq. (13) and Eq. (14) would vary with \(g\). Also, simple computation shows that SPD hyperplane is actually a submanifold of \(\mathcal{S}^{n}_{++}\) but we still follow the nomenclature of [33, 36]. Lastly, Definition 3.1 and Definition 3.2 can also be literally applied to other matrix manifolds, as matrix manifolds are usually geodesically complete.

### General SPD MLR under PEMs

As stated in Theorem 1, PEMs are not a single Riemannian metric. Instead, they denote a family of Riemannian metrics pulled back from Euclidean space. In this subsection, we follow the notation in Theorem 1 and will conclude that RMLR under any PEM can be uniformly expressed.

Before embarking on technical details, let us first clarify why we choose PEMs as our starting metrics. Several Riemannian metrics, including LEM, LCM, and ALEMs, all end up as PEMs. We will further show that all OILEMs are PEMs. Besides, when formulating Eq. (12) into manifolds, the core concerns lie in the optimization problem of calculating the distance \(d(x,H_{a_{k},b_{k}})\). The calculation of this distance under PEMs enjoys theoretical convenience, while other metrics like AIM would be complicated to obtain the distances to hyperplanes.

Now, we start by calculating Eq. (15), the distance to an SPD hyperplane.

**Lemma 2**.: _The distance of \(S\in\mathcal{S}^{n}_{++}\) to the SPD hyperplane \(\tilde{H}_{A_{k},P_{k}}\) is reduced to the distance of \(\phi(S)\) to the Euclidean hyperplane \(H_{\phi_{*,P_{k}}(A_{k}),\phi(P_{k})}\) in the Euclidean space of \(T_{P_{k}}\mathcal{S}^{n}_{++}\):_

\[d(S,\tilde{H}_{A_{k},P_{k}}))=d(\phi(S),H_{\phi_{*,P_{k}}(A_{k}),\phi(P_{k})}) =\frac{|\langle\phi(S)-\phi(P_{k}),\phi_{*,P_{k}}(A_{k})\rangle|}{\|A_{k}\|_{P _{k}}}, \tag{16}\]

_where \(|\cdot|\) is the absolute value._

Putting Eq. (16) into Eq. (14), we obtain our SPD classifiers under any PEMs:

\[p(y=k\mid S)\propto\exp(\langle A_{k},\operatorname{Log}_{P_{k}}(S)\rangle_{P_ {k}})=\exp(\langle\phi(S)-\phi(P_{k}),\phi_{*,P_{k}}(A_{k})\rangle). \tag{17}\]

One might have observed that \(A_{k}\in T_{P_{k}}\mathcal{S}^{n}_{++}\) in Eq. (17) are a non-Euclidean parameter, as \(P_{k}\) would vary during training. Fortunately, we have different tricks to avoid this issue. The first solution is the parallel transportation from a fixed tangent space, writing \(A_{k}=\Gamma_{Q\to P_{k}}(\tilde{A}_{k})\) with \(\tilde{A}_{k}\in T_{Q}\mathcal{S}^{n}_{++}\) as a Euclidean parameter. This is the solution adopted by HNNs [33], where the tangent point is the identity element. Alternatively, one can also rely on the differential of a Lie group translation, which is widely used in differential manifolds [37, SS 20]. More interestingly, under PEMs, the above two solutions are equivalent and anchor points can be arbitrarily chosen (See Appendix D for technical details). Therefore, without loss of generality, we generate \(A_{k}\) from tangent space at the identity by parallel transportation, _i.e._,\(A_{k}=\Gamma_{I\to P_{k}}(\tilde{A}_{k})\) with \(\tilde{A}_{k}\in T_{I}\mathcal{S}^{n}_{++}\cong\mathcal{S}^{n}\). Together with Eq. (8), Eq. (17) can be further simplified.
**Theorem 3** (SPD MLR under PEMs).: _Under any PEMs, SPD multiclass logistics regression and SPD hyperplane is_

\[p(y =k\mid S)\propto\exp(\langle\phi(S)-\phi(P_{k}),\phi_{*,I}(\tilde{A} _{k})\rangle), \tag{18}\] \[\tilde{H}_{\tilde{A}_{k},P_{k}} =\{S\in\mathcal{S}^{n}_{++}:\langle\phi(S)-\phi(P_{k}),\phi_{*,I} (\tilde{A}_{k})\rangle=0\}, \tag{19}\]

One can observe that Eq. (18) and Eq. (19) are very similar to a Euclidean MLR. However, since \(\phi\) is normally non-linear and \(P_{k}\) is an SPD parameter, Eq. (18) can not hastily be identified with a Euclidean MLR. However, under some special circumstances, SPD MLR can be reduced to the familiar Euclidean MLR. To show this result, we first need to present the Riemannian Stochastic Gradient Descent (RSGD) under PEMs. General RSGD [38] is formulated as

\[W_{t+1}=\exp_{W_{t}}(-\gamma_{t}\Pi_{W_{t}}(\nabla_{W}f|_{W_{t}})) \tag{20}\]

where \(\Pi_{W_{t}}\) denotes the projection map mapping Euclidean gradient \(\nabla_{W}f|_{W_{t}}\) to Riemannian gradient. We have already obtained the formula for the Riemannian exponential map as shown in Eq. (7). We proceed to formulate \(\Pi\).

**Lemma 4**.: _For a smooth function \(f:\mathcal{S}^{n}_{++}\rightarrow\mathbb{R}\) on \(\mathcal{S}^{n}_{++}\) endowed with any kind of PEMs, the projection map \(\Pi_{P}:\mathcal{S}^{n}\to T_{P}\mathcal{S}^{n}_{++}\) at \(P\in\mathcal{S}^{n}_{++}\) is_

\[\Pi_{P}(\nabla_{P}f)=\phi_{*,P}^{-1}(\phi_{*,P}^{-*})(\nabla_{P}f), \tag{21}\]

_where \(\phi_{*,P}^{-*}\) is the adjoint operator of \(\phi_{*,P}^{-1}\), i.e.,\((V_{1},\phi_{*,P}^{-1}V_{2})_{P}=\langle\phi_{*,P}^{-*}V_{1},V_{2}\rangle_{P}\), for all \(V_{i}\in T_{P}\mathcal{S}^{n}_{++}\)._

Together with the above lemma, we can describe the special case we mentioned.

**Theorem 5**.: _Supposing the differential map \(\phi_{*,I}:T_{I}\mathcal{S}^{n}_{++}\to T_{0}\mathcal{S}^{n}\) is the identity map, and \(P_{k}\) in Eq. (18) is optimized by PEM-based RSGD, then Eq. (18) can be reduced to a Euclidean MLR in the codomain of \(\phi\) updated by Euclidean SGD._

## 4 Riemannian multiclass logistic regression under OILEMs

Although Theorem 3 can deal with Riemannian classifiers under all PEMs, we focus on OILEMs to establish our SPD MLR in this paper.

### OILEMs-based SPD MLR

LEM enjoys fast and simple computation on SPD manifolds [21] and has shown many successful applications in dealing with SPD data [39; 17; 29; 10; 18]. The nascent OILEMs [25] are natural generalizations of LEM. Furthermore, as shown in the following lemma, all the OILEMs are PEMs.

**Lemma 6**.: \((\alpha,\beta)\)_-LEM is a pullback metric by \(\phi_{p,q}=F_{p,q}\circ\mathrm{mlog}\) from the standard Euclidean space \(\mathcal{S}^{n}\)._

Therefore, the previous discussion of PEMs-based MLR can be directly applied to OILEMs-based ones as well. For the above two reasons, we choose OILEMs as the Riemannian metrics to build our SPD MLR.

Calculating the differential of \(\phi_{p,q}\) at \(I\) directly gives the final formulation of SPD MLR.

**Proposition 7** (Final formulation of SPD MLR under \((\alpha,\beta)\)-LEM).: _On the SPD manifold with \((\alpha,\beta)\)-LEM, the SPD MLR and SPD hyperplane can be written as_

\[p(y =k\mid S)\propto\exp(\langle\phi_{p,q}(S)-\phi_{p,q}(P_{k}),F_{p, q}(\tilde{A_{k}})\rangle), \tag{22}\] \[\tilde{H}_{\tilde{A}_{k},P_{k}} =\{S\in\mathcal{S}^{n}_{++}:\langle\phi_{p,q}(S)-\phi_{p,q}(P_{k} ),F_{p,q}(\tilde{A_{k}})=0\}, \tag{23}\]

_where \(\tilde{A}_{k}\in\mathcal{S}^{n}\) is the normal matrix, \(P_{k}\in\mathcal{S}^{n}_{++}\) is the shift matrix, and \(\phi_{p,q}=F_{p,q}\circ\mathrm{mlog}\)._

Figure 2: Conceptual illustration of SPD hyperplanes on \(\mathcal{S}^{2}_{++}\). The black dots are symmetric positive semi-definite matrices, denoting the boundary of \(\mathcal{S}^{2}_{++}\). The blue, red, and yellow dots denote three SPD hyperplanes.
Figure 2 illustrates three OILEMs-based hyperplanes on SPD manifolds. As a submanifold of \(\mathbb{R}^{3}\), \(\mathcal{S}^{2}_{++}\) can be visualize in \(\mathbb{R}^{3}\), by \(\forall P=\left(\begin{array}{cc}x&y\\ y&z\end{array}\right)\in\mathcal{S}^{2}\) is positive definite iff \(x,z>0\wedge xz>y^{2}\).

_Remark 3_.: When \(p=1,q=1\) (\(\alpha=1,\beta=0\)), OILEM is exactly the familiar LEM. Eq. (22) and Eq. (23) are also reduced to SPD MLR under LEM.

### Understanding LogEig classifier in existing SPD networks

Many of the existing SPD neural networks [11; 6; 31; 32; 17; 29; 30] rely on a Euclidean MLR in the codomain of matrix logarithm. _i.e._,matrix logarithm followed by an FC layer and a softmax layer. For simplicity, we call this classifier as LogEig MLR. The existing explanation of LogEig MLR is approximating manifolds by tangent space. Here, we claim that it is a special case of our SPD MLR.

Although Eq. (22) is very similar to the LogEig MLR, as we stated before, due to the nonlinearity of \(\mathrm{mlog}\) and non-Euclideanness of SPD parameter \(P_{k}\), SPD MLR cannot be hastily viewed equivalent to LogEig MLR. That being said, under special circumstances, as a direct corollary of Theorem 5, Eq. (22) is equivalent to a LogEig MLR.

**Corollary 8**.: _Endowing SPD manifolds with LEM, optimizing SPD parameter \(P_{k}\) in Eq. (22) by LEM-based RSGD and Euclidean parameter \(A_{k}\) by Euclidean SGD, the LEM-based SPD MLR is equivalent to a LogEig MLR with parameters in FC layer optimized by Euclidean SGD._

The widely used LogEig MLR can be thus geometrically explained as a special case of our approach.

### Learning SPD parameters

Similar to Corollary 8, optimizing \(P_{k}\) in SPD MLR by OILEM-based RSGD would be equivalent to a Euclidean MLR in the codomain of \(\phi_{p,q}\). However, both theoretical analysis [40] and empirical experiments [41; 5] have demonstrated the benefits of AIM-based optimization. Therefore, we choose AIM to optimize the SPD parameter \(P_{k}\) in our SPD MLR. Besides, we also adopt the RSGD based on the less explored Bures-Wasserstein Metric (BWM) [22], which has shown promising performance for ill-conditioned matrices [42].

The required operators for RSGD under these two metrics are well studied in the existing literature [34; 40; 22]. We summarize them in Table 1, where \(\mathcal{L}_{P}(V)\) is known as Lyapunov operator, _i.e._,\(\mathcal{L}_{P}(V)P+P\mathcal{L}_{P}(V)=V\), and \((A)_{\mathrm{sym}}=\frac{A+A^{\top}}{2}\). For fast computation of BWM-based RSGD, we further adopt the Newton-Schulz method to calculate the Lyapunov operator. More detail is exposed in Appendix E.

As for the gradient computation and backpropagation of SPD MLR, please refer to Appendix C.

### Final SPD MLR algorithm

Now we write the full algorithm of our SPD MLR. Recalling Eq. (22), for each class \(k\in\{1,\cdots,C\}\), we have a normal matrix \(\tilde{A}_{k}\in\mathcal{S}^{n}\) and a biasing matrix \(P_{k}\in\mathcal{S}^{n}_{++}\). Computationally, Eq. (22) first applies \(\phi_{p,q}\) to each \(P_{k}\) and input SPD feature \(S_{i}\), and apply \(F_{p,q}\) to each \(\tilde{A}_{k}\). Then, the associated inner products are calculated and softmax is applied. The inner products can be further efficiently carried on by matrix product. Therefore, we can concatenate all the Euclidean data \(F_{p,q}(\tilde{A}_{k})\) into \(A\) and apply traditional linear operation with dimensionality reduction matrix \(A\) and biasing vector \(b=(b_{i},\cdots,b_{C})\) with \(b_{k}=\langle F_{p,q}(\tilde{A}_{k}),\phi_{p,q}(P_{k})\rangle)\). In practice, our MLR can be applied as a classifier to any SPD networks. We present the above process in Algorithm 1.

\begin{table}
\begin{tabular}{c c c} \hline \hline Operators & AIM & BWM \\ \hline Exponential map \(\mathrm{Exp}_{P}(V)\) & \(P^{\frac{1}{2}}\mathrm{mexp}(P^{-\frac{1}{2}}SP^{-\frac{1}{2}})P^{\frac{1}{2}}\) & \(P+V+\mathcal{L}_{P}(V)P\mathcal{L}_{P}(V)\) \\ Riemannian gradient \(\Pi_{P}(\nabla_{P}f)\) & \(P(\nabla_{P}f)_{\mathrm{sym}}P\) & \(\Pi_{P}(\nabla_{P}f)=4(\nabla_{P}fP)_{\mathrm{sym}}\) \\ \hline \hline \end{tabular}
\end{table}
Table 1: Riemannian optimization operators for AIM and BWM.
```
Hyper-parameters: \(\alpha,\beta,\mathrm{s.t.}p=\sqrt{\alpha+n\beta}>0,q=\sqrt{\alpha}>0\); Parameters: \(n\times n\) SPD parameters \(\{P_{k}\}_{k\leq C}\), \(n\times n\) symmetric matrices \(\{\tilde{A}_{k}\}_{k\leq C}\); Inputs: a batch of \(n\times n\) SPD matrices \(\{S_{i}\}_{i<N}\); Step 1: mapping SPD features and parameters: \[\forall k\leq C,i\leq N,\bar{P}_{k}\leftarrow\phi_{p,q}(P_{k}),\bar{A}_{k} \gets F_{p,q}(\tilde{A}_{k}),\bar{S}_{i}\leftarrow\phi_{p,q}(S_{i});\] Step 2: concatenation: \[A\leftarrow\mathrm{concat}(\mathrm{vec}(\bar{A}_{1}),\cdots,\mathrm{vec}(\bar{A }_{C})),S\leftarrow\mathrm{concat}(\mathrm{vec}(\bar{S}_{1}),\cdots,\mathrm{ vec}(\bar{S}_{N}));\] Step 3: calculating multinomial probabilities: \[v\leftarrow\mathrm{softmax}(SA^{\top}-(b,\cdots,b)^{\top}),\] where bias vector \(b\in\mathbb{R}^{C}\) and \(b_{i}=\langle\bar{P}_{i},\bar{A}_{i}\rangle\); Output: \(v\in\mathbb{R}^{N\times C}\); Updating \(P_{k}\) by AIM/BWM-based RSGD, \(\tilde{A}_{k}\) by Euclidean optimization.
```

**Algorithm 1**Training \((\alpha,\beta)\)-LEM-based SPD multiclass logistics regression on SPD networks

### Further explanations on hyper-parameters

As we have discussed, SPD MLR and the corresponding hyperplane would vary with different \((\alpha,\beta)\), respecting various OILEMs. Here, we make further explanations about the effect of \(\alpha,\beta\) or equivalently \((p,q)\).

Recalling \(F_{p,q}(X)=qX+\frac{p-q}{n}\operatorname{tr}(X)I_{n}\) in Lemma 6, the second term can be viewed as \((p-q)\frac{\operatorname{tr}(X)I_{n}}{n}\). Intuitively, \(p-q\) is related to the importance of the average of trace. To see this effect more clearly, we present the general expression of OILEMs as follows, which is not covered in [25].

**Proposition 9** (\((\alpha,\beta)\)-LEM).: _The \((\alpha,\beta)\)-LEM is formulated as_

\[\langle V_{1},V_{2}\rangle_{S}=q^{2}\langle\tilde{V}_{1},\tilde{V}_{2}\rangle_ {\mathbb{F}}+(p-q)^{2}\frac{\operatorname{tr}(\tilde{V}_{1})\operatorname{ tr}(\tilde{V}_{2})}{n}+q(p-q)\frac{2(\operatorname{tr}(\tilde{V}_{1}) \operatorname{tr}(\tilde{V}_{2}))}{n}, \tag{24}\]

_where \(S\) is an SPD matrix, \(V_{i}\in T_{S}\mathcal{S}_{++}^{n}\), \(\tilde{V}_{i}=\mathrm{mlog}_{*,S}V_{i},p=\sqrt{\alpha+n\beta}\), and \(q=\sqrt{\alpha}\)._

Now we can see clearly the effect of \(p,q\) by Eq. (24). Compared with the vanilla LEM (\(p=q=1\)), the second and third terms in Eq. (24) indicates that _the absolute value of \(p-q\) control the magnitude of trace in the Riemannian metric._ The third term indicates that _the sign of \(p-q\)_ (\(+,-,0\)) _determines whether the trace is amplified, suppressed, or neutralized._ Correspondingly, there are three basic ratios between \(q\) and \((p-q)\), _i.e._,\(1:1\), \(1:-0.99\), and \(1:0\). Note that since \(1:-1\) cannot ensure \(\mathrm{O}(n)\)-invariance, we lower the ratio and use \(1:-0.99\) instead. Although there are more fine-grained ratios between \(q\) and \((p-q)\), we will focus on the above three most basic ones.

## 5 Experiments

To fairly validate the effectiveness of our methods, we adopt the two most classic SPD networks, _i.e._,SPDNet [11] and SPDNetBN [6], as backbones and validate the performance on Radar [6], HDM05 [43], and AFEW [44] datasets. A brief review of their basic layers is presented in Appendix B.3. We implement our SPD MLR based on the official code by SPDNetBN1. We apply our methods to these two networks, by substituting their LogEig MLR with our SPD MLR. We focus on three basic cases, where \(q:(p-q)\) is \(1:0,1:1\), and \(1:-0.99\), respectively. We call SPDNet-RMLR-(1,0) an SPDNet with our Riemannian classifier with \(q=1,p-q=0\). So does SPDNet-RMLR-(1,1) and SPDNet-RMLR-(1,-0.99). We initialize our normal matrices \(\tilde{A}_{k}\) by Kaiming uniform strategies [45], and \(P_{k}\) as identity matrix. Note that, with this initialization, the
initial status of our classifier under LEM (\(q=1,p-q=0\)) is exactly equivalent to LogEig MLR in the PyTorch implementation of SPDNetBN. We denote \(\{d_{0},d_{1},\cdots,d_{L}\}\) as the dimensions of each transformation layer in the SPDNet backbone. The batch size is 30, and the optimizer is SGD. We will first validate SPDNet-RMLR optimized by AIM RSGD. In Section 5.2, we will further implement our classifier into SPDNetBN, and implement our classifier by BWM-based RSGD, respectively.

**Impact of \(p\) and \(p-q\).** For the hyper-parameters \(p\) and \(p-q\), as we stated before, the variants of \(p,q\) respect different kinds of OILEMs. The best setting depends on the characteristics of the datasets. Our general observations are that when the dimension is relatively low (\(8\times 8\) on Radar), extra attention to trace might be beneficial, while under relatively high dimension (\(30\times 30\) on HDM05 and \(50\times 50\) on AFEW), \(1:0\) tends to be already saturated. This matches our intuition that with the growing dimensions, the proportion of diagonal elements gets smaller. The trace is therefore relatively less important when calculating the Riemannian metric in Eq. (24). Below are our detailed results.

### Experiments on SPDNet

**Drone recognition.** We adopt the Radar dataset2[6], composed of \(3,000\) synthetic radar signals, each of which is split into windows of length 20, resulting in 3000 \(20\times 20\) SPD matrices equally distributed in 3 classes. In line with [6], we assign 50%, 25%, and 25% of data for training, validation, and testing. We test our classifiers under two suggested architectures [6] on this dataset, {20, 16, 8} and {20, 16, 14, 12, 10, 8}. Different learning rates are also tested. The 10-fold results are presented in Table 2. Our classifier with ratios (1,0) and (1,1) can bring consistent performance gain to SPDNet under different settings. In particular, SPDNet-RMLR-(1,1) achieves much better performance and less variance, compared with the vanilla SPDNet. In contrast, suppressing the trace part, _i.e._,RMLR-\((1,-1)\), is harmful to this task. The accuracy curve versus epochs is also presented in Figure 3. The performance gain is consistent throughout the training.

Footnote 2: [https://www.dropbox.com/s/dfnlx2bnyh3kjwy/data.zip?dl=0](https://www.dropbox.com/s/dfnlx2bnyh3kjwy/data.zip?dl=0)

**Action recognition.** HDM05 dataset [43] contains 2,273 skeleton-based motion capture sequences executed by various actors. Each frame consists of 3D coordinates of 31 joints of the subjects and each sequence can be therefore modeled by a \(93\times 93\) covariance matrix. For a fair comparison, we adopt the pre-processed covariance features3 released by SPDNet [11], where datasets are augmented

\begin{table}
\begin{tabular}{c|c c|c c} \hline \hline Learning Rate & \multicolumn{2}{c|}{\(1e^{-2}\)} & \multicolumn{2}{c}{\(2.5\)e\({}^{-2}\)} \\ \hline Architecture & {20, 16, 8} & {20, 16, 14, 12, 10, 8} & {20, 16, 8} & {20, 16, 14, 12, 10, 8} \\ \hline SPDNet & 89.89\(\pm\)1.21 & 90.88\(\pm\)0.61 & 91.93\(\pm\)0.84 & 92.32\(\pm\)0.50 \\ SPDNet-RMLR-(1,0) & 92.74\(\pm\)1.01 & 93.25\(\pm\)0.98 & 94.19\(\pm\)0.91 & 93.33\(\pm\)0.79 \\ SPDNet-RMLR-(1,-0.99) & 85.15\(\pm\)1.34 & 85.02\(\pm\)0.76 & 85.65\(\pm\)1.34 & 85.84\(\pm\)0.78 \\ SPDNet-RMLR-(1,1) & **94.57\(\pm\)1.08** (\(\pm\)4.68) & **95.08\(\pm\)0.48** (\(\uparrow\)4.22) & **94.97\(\pm\)0.70** (\(\uparrow\)3.04) & **95.33\(\pm\)0.61** (\(\uparrow\)3.01) \\ \hline \hline \end{tabular}
\end{table}
Table 2: Results of SPDNet with and without Riemannian MLR on the Radar dataset.

Figure 3: Test accuracy on Radar (left) and HDM05 (right) datasets. The architectures are {20, 16, 8} and {93, 30}. The learning rates are \(1e^{-2}\) and \(5e^{-2}\).
to 18,294 for training, and 1,197 for testing. Following [11; 6], we adopt the architecture of {93, 30} and the learning rate of \(5e^{-2}\). The 10-fold results are presented in Table 3. For this task, RMLR-\((1,0)\) is the best classifier. Compared with the vanilla LogEig MLR, RMLR-\((1,0)\) shows less variance as well. In fact, RMLR-\((1,0)\) respects the geometry induced by LEM, indicating that for this task, no extra attention should be paid to trace when computing Riemannian metric tensor in Eq. (24). We also present the accuracy curve in Figure 3.

**Emotion recognition.** AFEW dataset [44] covers 7 kinds of emotion, consisting of 2,118 video clips with 1,747 for training and 371 for testing. Following [11; 5], each video is modeled by a \(400\times 400\) covariance matrix. Following [11], we also validate our classifier under different network architectures. The learning rate is set to be \(5e^{-2}\). The 10-fold results are presented in Table 4. Note that since SPDNet shows relatively large fluctuation on this dataset, we also present the best result among 10 folds. Similar to HDM05, SPDNet-RMLR-\((1,0)\) shows the most promising results. The best result of SPDNet-RMLR-\((1,0)\) reaches 37.2, while the best result of vanilla SPDNet is 35.8.

### Ablation studies

**Experiments on SDPNetBN.** We further apply our classifier to SPDNetBN. We adopt the best settings of \(q\) and \(p-q\) obtained in our above experiments, _i.e._,(1,1) for Radar, (1,0) for HDM05 and AFEW. The learning rate is \(1e^{-2},1e^{-2},\) and \(5e^{-2}\), respectively. The results are presented in Table 5. Our classifier brings better performance and less variance on Radar and HDM05 datasets. Although the average performance on the AFEW dataset is slightly worse than SPDNetBN, the best performance is achieved by our approach. Nevertheless, we admit that the general performance gain of our SPDNetBN-MLR against SPDNetBN might not be as obvious as SPDNet-MLR against SPDNet. We think the underlying reason comes from the used metrics: the Riemannian batch normalization (RBN) in SPDNetBN is based on AIM, while our classifiers are based on OILEMs. This inconsistency could undermine the effectiveness of our methods.

**BWM-based RSGD.** We focus on the SPDNet backbone. We call SPDNet-RMLR-BWM a network where RMLR is updated by BWM-based RSGD. So does SPDNet-RMLR-AIM, which is the model we have tested. We carry 10-fold experiments on all three datasets, with learning rates \(1e^{-2},5e^{-2}\), and \(5e^{-2}\), respectively. We adopt the optimal value of \((p,p-q)\) obtained before. We set the maximal iteration in the Newton-Schulz method as 8 while other settings remain the same. The average performance and training time (s/epoch) are presented in Table 6. We observe that although BWM shows obviously better performance on the HDM05 dataset against AIM, on the other two datasets, BWM is generally inferior to AIM, similar to the conclusion reached in [42]. However, BWM could be more efficient than AIM, especially under multiple SPD parameters. The main reason is the fast computation of the Lyapunov equation.

\begin{table}
\begin{tabular}{c|c c|c c|c c} \hline \hline Dataset & \multicolumn{3}{c|}{Radar} & \multicolumn{2}{c|}{HDM05} & \multicolumn{2}{c}{AFEW} \\ \hline Architecture & \multicolumn{2}{c|}{[20, 16, 8]} & \multicolumn{2}{c|}{[20, 16, 14, 12, 10, 8]} & \multicolumn{2}{c|}{[93, 30]} & \multicolumn{2}{c}{[400, 200, 100, 50]} \\ \hline Measure & Mean\$STD & Time & Mean\$STD & Time & Mean\$STD & Time & Mean\$STD & Max & Time \\ \hline SPDNet & \$89.89\$1.21 & 1.21 & 0.98\$0.61 & 2.59 & 62.80\$1.01 & 21.98 & 34.16\$1.04 & 35.8 & 74.18 \\ SPDNet-RMLR-AIM & **94.57\$1.08** & 1.39 & **95.08\$0.48** & 2.56 & 63.93\$0.49\$2 & 16.27 & **34.66\$1.45** & **37.2** & 79.15 \\ SPDNet-RMLR-BWM & 93.59\$0.99 & 1.32 & 94.03\$0.48 & 2.51 & **65.26\$0.25** & 44.89 & 34.35\$0.98 & 35.99 & 76.13 \\ \hline \hline \end{tabular}
\end{table}
Table 6: Results of SPDNet-RMLR under BWM RSGD against AIM RSGD.

\begin{table}
\begin{tabular}{c|c c|c c|c c} \hline \hline Architecture & \multicolumn{2}{c|}{[400, 50]} & \multicolumn{2}{c|}{[400, 200, 50]} & \multicolumn{2}{c}{[400, 200, 100, 50]} \\ Measure & Mean\$STD & Max & Mean\$STD & Max & Mean\$STD & Max \\ \hline SPDNet & 31.86\$0.95 & 33.05 & 33.13\$0.55 & 34.08 & 34.16\$1.04 & 35.8 \\ SPDNet-RMLR-(1,0) & **33.14\$1.48\$1.7\$1\$2\$3\$3\$6\$6\$7\$1\$2\$3\$3\$6\$1\$2\$1\$1\$2\$1\$1\$2\$1\$1\$2\$1\$2\$1\$2\$1\$2\$1\$2\$1\$2\$1\$2\$1\$2\$1\$2\$1\$2\$1\$2\$1\$2\$1\$2\$1\$2$1\$2\$1\$2\$1\$2\$1\$2\$1\$2\$1\$2\$1\$2\$1\$2\$1\$2\$1\$2\$1\$2\$1\$2\$1\$2\$1\$2\$1\$2\$1\$2\$1\$2\$2\$1\$2\$1\$2\$1\$2\$2\$1\$2\$2\$1\$2\$2\$1\$2\$2\$2\$1\$2\$2\$2\$2\$1\$2$2\$2\$2$2\$2$2\$2$
Conclusion

In this paper, we are the first to develop Riemannian classifiers for SPD neural networks. We provide a general classifier for all kinds of metrics pulled back from Euclidean space. Furthermore, we showcase our framework under a family of \(\mathrm{O}(n)\)-invariant Riemannian metrics, called OILEMs. To the best of our knowledge, our work is the first one to maneuver with these metrics in machine learning. The consistent superior performance in extensive experiments also supports our claims. As a future avenue, our framework can also be readily applied to other kinds of PEMs.
