Indefinite Stiefel manifold, Riemannian gradient descent, tractable metric, Cayley retraction, trace minimization, eigenvalue problem, matrix least square

1 Introduction↩︎

Optimization of smooth functions of variable \(X\in{\mathbb{R}^{n\times k}}\) with \(k\leq n\) under the condition \[\label{eq:gen95constraint} X^TAX = B,\tag{1}\] in which \(A\in{\mathbb{R}^{n\times n}}, B \in {\mathbb{R}^{k\times k}}\), possibly with some conditions imposed, is popular thanks to its wide range of applications. In many circumstances when the feasible set, \[\mathcal{F}:= \{X\in{\mathbb{R}^{n\times k}}\;:\;X^TAX = B\},\] constitutes a differentiable manifold, various algorithms were developed by extending popular methods for unconstrained optimization in the Euclidean space to the Riemannian manifold once necessary geometric tools such as the tangent space, metric, projection, and retraction are efficiently realized. Such an approach is generally referred to as Riemannian optimization.

The most well known situation is when \(A = I_n\), the \(n\times n\) identity matrix, and \(B=I_k\), where \(\mathcal{F}\) is named as the (orthogonal) Stiefel manifold. Numerous practical issues can be (re)formulated as an optimization problem on the Stiefel manifold, i.e., with the orthogonality constraint. For instance, in dimensionality reduction of data, principal component analysis (PCA) requires solving \[\label{eq:PCAProb} \max_{X^T X = I_k}\mathrm{tr}\left(X^T\bar{S}^T\bar{S}X\right),\tag{2}\] where \(\bar{S}\) is the centered data set, i.e., each of its columns is of the form \(\bar{s}_j = s_j - \mu \in \mathbb{R}^n\) with \(\mu = \frac{1}{m}\sum_{j=1}^ms_j\). Ky-Fan’s theorem is a similar situation in numerical linear algebra which states that given a symmetric matrix \(M\) whose eigenvalues are denoted by \(\lambda_1 \leq \ldots \leq \lambda_n\), we have \[\label{eq:KyFan} \min_{X^TX = I_k}\mathrm{tr}\left(X^TMX\right) = \sum_{j=1}^k\lambda_j,for eachk = 1,\ldots,n.\tag{3}\] Riemannian optimization is among the methods of choice for solving 2 and 3 . The reader is referred to the classic work [@EdelAS98] and excellent monographs [@Udri1994; @AbsiMS08; @Sato2021; @Boumal23] for details. A slightly different context is when \(A\in \mathrm{SPD}(n),\) the set of \(n\times n\) symmetric positive-definite (spd) matrices, and \(B = I_k\) occurring in the study of locality preserving projections where one wants to find a projection matrix that preserves some affinity graph from a given data set; it leads to the optimization problem \[\label{eq:LPP} \min_{X^T(\hat{S}^T\hat{S})X = I_k} \mathrm{tr}\left(X^T\hat{S}^T(I-\hat{W})\hat{S}X\right),\tag{4}\] where \(\hat{S}\) is generally of full rank and \(\hat{W}\) is a modified weighted matrix [@HeN03]. A similar problem arises in canonical correlation analysis whose feasible set is the Cartesian product of two sets as that of 4. Meanwhile, an extension of Ky-Fan’s theorem [@SameW82] for the symmetric matrix pencil \((M,A)\), where \(A\) is spd, also leads to \[\label{eq:KyFan2} \min_{X^{T}AX = I_k}\mathrm{tr}\left(X^{T}MX\right) = \sum_{j=1}^k\lambda_j,for each k = 1,\ldots,n,\tag{5}\] in which \(\lambda_j, j = 1,\ldots,n,\) are the (generalized) eigenvalues of the matrix pencil \((M,A)\). Solutions to the optimization problems with constraints akin to 4 and 5 motivate the investigation of Riemannian optimization methods on the generalized Stiefel manifold, denoted by \(\mathrm{St}_A(k,n)\), [@YgerBGR12; @manopt; @MishS2016; @SatoA19; @ShusA2023; @WangDPY2024]. These results provide us with tools to work directly on \(\mathrm{St}_A(k,n)\) without having to turn back to \(\mathrm{St}(k,n)\) via a Cholesky factorization of \(A\) which should be avoided, especially for large and/or sparse \(A\). More problems on minimization of the trace can be found in [@KokiS11].

Other authors [@GSAS21; @GSAS21a; @BendZ21; @OviH23; @YamaS2023; @GSS24; @JenZ24; @GaoSS2024] concentrated on developing optimization methods on the symplectic Stiefel manifold, which corresponds to the case where both \(A\) and \(B\) are the Poisson matrices, i.e., of the form \(\left[\begin{smallmatrix} 0&I_p\\-I_p&0 \end{smallmatrix}\right]\). This problem arises in model reduction of Hamiltonian systems [@PengM16; @BendZ22; @BuchGH22; @GSS24], computing the symplectic eigenvalues of (semi)spd matrices [@BhatJ15; @SonAGS21; @SonS22].

In [@KovaV95], Ky-Fan’s theorem was further investigated for \(n\times n\) symmetric positive-definite matrix pencils \(M-\lambda A\), i.e., the ones such that the matrix \(M-\lambda_0 A\) is spd for some \(\lambda_0\). In this extension, \(A\) is nonsingular and may be indefinite, i.e., \(A\) has \(p\) positive and \(m\) negative eigenvalues, \(p>0, m > 0, p+m = n\). Namely, let \(J = \mathrm{diag}(I_{k_p},-I_{k_m}),\) the block diagonal matrix, with \(0\leq k_p \leq p, 0\leq k_m \leq m, k_p + k_m = k\), and denote the eigenvalues of \((M,A)\) by \[\lambda^-_m\leq \ldots\leq \lambda^-_1 <\lambda^+_1 \leq \ldots \leq \lambda^+_p,\] where the minus/plus sign means the negative/positive value. Then we have \[\label{eq:KovacStrikoVeselic} \min_{X^TAX = J}\mathrm{tr}\left(X^TMX\right) = \sum_{j=1}^{k_p}\lambda^+_j - \sum_{j=1}^{k_m}\lambda^-_j.\tag{6}\] Once a minimizer \(X_*\) is obtained, one can determine the eigenvalues appearing on the right hand side of 6 and their associated normalized eigenvectors, i.e., an \(n\times k\) matrix \(\hat{X}\) satisfying the constraint \(\hat{X}^TA\hat{X}=J\) and its \(j\)-th column \(\hat{x}_j\) verifies \(M\hat{x}_j = \lambda_jA\hat{x}_j\) and \(\lambda_j = \lambda_j^+\) for \(j=1,\ldots,k_p\) and \(\lambda_{j+k_p} = \lambda_j^-\) for \(j=1,\ldots,k_m\), by solving two standard eigenvalue problems for the matrices of sizes \(k_p\times k_p\) and \(k_m\times k_m\) extracted from \(X_*^TMX_*\), see [@KovaV95] and section 6 for details. The need for computing several smallest positive eigenvalues of symmetric/Hermitian matrices/matrix pencils arises in practice, for instance, in the design of electromagnetic accelerators, see [@Geus2002] and references therein. In the case that \(A=J=\mathrm{diag}(\pm 1,\ldots,\pm 1),\) the signature matrix, the problem of determining the eigenvalues and eigenvectors mentioned above is referred to as the hyperbolic eigenvalue problem [@SlapT2000]. A solution to this problem is a main step for hyperbolic PCA which is arguably more informative than the conventional PCA [@TabaKWM2024] and the computation of hyperbolic singular value decomposition [@OnnSB1989]. Another issue is the linear response eigenvalue problem (LREVP) where one has to compute a few smallest eigenvalues of the matrix \(\left[ \begin{smallmatrix} 0&K\\M&0 \end{smallmatrix}\right]\) with symmetric positive-semidefinite matrices \(K, M\) in which at least one of them is spd. This problem originates from the studies of density functions, random phase approximation, see [@BaiL12] and references therein. Various numerical methods have been proposed, see [@BaiL2017] for a recent survey, among which the ones based on trace minimization [@BaiL13; @BaiL14; @ChenSL2023] are of our interest. It is direct to see that the mentioned LREVP is equivalent to the eigenvalue problem for the matrix pencil \[\label{eq:HGpencil} H - \lambda G := \begin{bmatrix} K&0\\0&M \end{bmatrix} - \lambda \begin{bmatrix} 0&I\\I&0 \end{bmatrix}.\tag{7}\] When both \(K\) and \(M\) are spd, 7 is equivalent to the product eigenvalue problem where one has to find the \(k\) smallest eigenvalues of a product matrix \(KM\) [@KresPS2014; @Pand2019]. Indeed, it can be verified that \(\lambda^2, \lambda\in \mathbb{R},\) is an eigenvalue of \(KM\) if and only if \(\pm\lambda\) are the eigenvalues of the matrix pencil 7 . The statement 6 suggests that the LREVP can be solved via finding a solution \(X\in \mathbb{R}^{2n\times k}\) to \[\label{eq:IndefTraceMin} \min_{X^TGX = J }\mathrm{tr}\left(X^THX\right)\tag{8}\] with \(H\) and \(G\) defined as in 7 and \(J = I_k\).

Obviously, the existing Riemannian optimization methods on the generalized Stiefel manifold and the symplectic Stiefel manifold do not help solve the optimization problems mentioned in 6 and 8 . The main reason lies in the difference in the feasible sets. This fact inspires us to consider the optimization problem \[\label{eq:opt95prob} \min_{X\in \mathbb{R}^{n\times k}} f(X)\; s.t. \;X^{T}AX = J,\tag{9}\] where \(f\) is continuously differentiable, \(A\in \mathbb{R}^{n\times n}\) is symmetric, nonsingular and \(J\in\mathrm{Sym}(k)\) satisfying \(J^2 = I_k\). As will be shown in section 2, the feasible set of this problem \[\label{eq:indefStiefel} \mathrm{iSt}_{A,J}(k,n) = \{X\in \mathbb{R}^{n\times k}\;:\;X^{T}AX = J\}\tag{10}\] constitutes a differentiable manifold. We will refer to it as the indefinite Stiefel manifold to emphasize the possible indefiniteness of \(A\) and \(J\)³. If the two matrices \(A,J\) or the size \(n\times k\) of the ambient space are either clear from the context and/or inessential, we also write \(\mathrm{iSt}(k,n)\) or \(\mathrm{iSt}_{A,J}\). This manifold covers several sets considered in the literature. If \(A\) is additionally positive-definite and \(J = I_k\), the constraint set in our problem is a generalized Stiefel manifold which naturally reduces to the Stiefel manifold if \(A=I_n\). For \(A = \mathrm{diag}(-1,1,\ldots,1)\), if \(J = I_k\), we obtain the hyperbolic manifold mentioned in [@XiaoLT2024; @HanJM2024]; if \(J = -1\), the Lorentz model or the hyperboloid model of hyperbolic space, see, e.g., [@NickK2018; @WilsL2018], is a half of \(\mathrm{iSt}_{A,J}\). Moreover, if \(A = J = \mathrm{diag}(\pm 1,\ldots,\pm 1)\), we are having the \(J\)-orthogonal group [@High03; @HeYWX2024].

Despite the similarity, there are essential differences between the indefinite Stiefel manifold and other Stiefel-type manifolds. The generalized Stiefel manifold, including the Stiefel manifold as a special case, is compact while the indefinite Stiefel manifold is generally not so. This can be easily seen by considering the hyperboloid model \(x_1^2 - x_2^2 = 1\) in \(\mathbb{R}^{2\times 1}\). Moreover, \(\mathrm{iSt}_{A,J}\) in this paper is defined generally, both \(J\) and \(A\) are symmetric (with appropriate restrictions), and therefore completely different from the symplectic Stiefel manifold defined in [@GSAS21] which is attached to two fixed (skew-symmetric) Poisson matrices \(J_{2n}\) and \(J_{2k}\).

1.0.0.1 Contribution and outline

To the best of our knowledge, our present work is the first attempt to approach the optimization problem with possibly indefinite quadratic constraint \(X^TAX = J\) in the Riemannian framework. Namely, in section 2 we show that \(\mathrm{iSt}_{A,J}\) is an embedded submanifold of \({\mathbb{R}^{n\times k}}\) and describe its tangent vectors. We then equip the manifold with a family of tractable metrics, construct the associated orthogonal projections onto the tangent and normal spaces, and compute the Riemannian gradient of a general cost function in section 3. A retraction based on the Cayley transformation is studied in section 4 and a gradient descent algorithm is proposed in section 5. The obtained results are then verified by various numerical examples in section 6. Since the results in this work are not necessarily restricted to the indefinite case, they also hold for the generalized Stiefel manifold; with suitable restrictions, some of them reduce to the results derived for generalized Stiefel manifold, see remarks 1, 2, and 3. Moreover, our result enables a Riemannian approach for optimization with \(J\)-orthogonality constraint.

1.0.0.2 Notation

First, \(\dim(\cdot)\) will denote the dimension of a linear space, \(\mathrm{diag}(\cdot,\ldots,\cdot)\) is the (block) diagonal matrix whose diagonal (block) entries are the listed objects. For any matrix \(\Omega\), the notations \(\Omega^T, \mathrm{tr}\,\Omega,\, \|\Omega\|_F,\) and \(\|\Omega\|\) are respectively its transpose, trace, Frobenius norm, and spectral norm, if not specified otherwise. If \(\Omega\) is a square matrix, we denote by \(\mathrm{skew}(\Omega) \mathrel{\vcenter{:}}= \frac{1}{2}(\Omega-\Omega^{T})\) and \(\mathrm{sym}(\Omega) \mathrel{\vcenter{:}}= \frac{1}{2}(\Omega + \Omega^{T})\) its skew-symmetric and symmetric parts, respectively. Meanwhile, \(\mathrm{Skew}(p)\), \(\mathrm{Sym}(p)\), and \(\mathrm{O}(p)\) stand for the set of skew-symmetric, the set of symmetric, and the set of orthogonal matrices of size \(p\times p,\) respectively. The Fréchet derivative of a map \(F:\mathcal{V}_1\longrightarrow\mathcal{V}_2\) at \(X\in\mathcal{V}_1\) between two normed vector spaces \(\mathcal{V}_1\) and \(\mathcal{V}_2\) is defined to be a linear operator \(\mathrm{\small D}F(X):\mathcal{V}_1\longrightarrow\mathcal{V}_2\) such that \(F(X+Z) = F(X)+\mathrm{\small D}F(X)[Z]+o(\| Z\|),\) for any \(Z\in\mathcal{V}_1,\) where \(o(\cdot)\) is the little-\(o\) notation. Finally, for a smooth mapping \(\gamma\) of one variable \(t\), \(\frac{\mathrm{d}}{\mathrm{d}t}(\gamma)\) or \(\dot{\gamma}(t)\) denotes the derivative of \(\gamma\) with respect to \(t.\)

2 The indefinite Stiefel manifolds↩︎

Without further restriction on \(A\) and \(J\), the set \(\mathrm{iSt}(k,n)\) can be empty. For instance, there is no real number pair \((x_1,x_2)\) such that \(x_1^2 + x_2^2 = -1\). The following statement helps avoid working with such a set.

From now on, it is essential to assume that the condition ?? is always fulfilled. The following result gives the foundation for \(\mathrm{iSt}(k,n).\)

A proof for a special case when both \(A\) and \(J\) are of the form \(\mathrm{diag}(\pm 1,\ldots,\pm 1)\) was given in [@KovaV95]. It can be used to obtain the statement made in Proposition [theo:PropertiesH] via a suitable diffeomorphism but the given direct proof better prepares the reader for later development.

Let us denote by \(\langle\cdot,\cdot\rangle_{\mathbb{R}^{p\times q}}\) the standard Euclidean inner product on \(\mathbb{R}^{p\times q}\). That is, \(\langle U,V\rangle_{\mathbb{R}^{p\times q}} =\mathrm{tr}(U^TV)\) for \(U,V\in\mathbb{R}^{p\times q}\). Given \(X\in\mathrm{iSt}_{A,J},\) we are going to formulate the adjoint operator \(\mathrm{\small D}F^{*}(X)\) of \(\mathrm{\small D}F(X)\) with respect to this inner product, which is the linear operator \[\mathrm{\small D}F^{*}(X) :\mathrm{Sym}(k) \longrightarrow {\mathbb{R}^{n\times k}}\] satisfying the condition \[\label{eq:DF95adjoint} \langle\mathrm{\small D}F(X)[Z],\Omega\rangle_{{\mathbb{R}^{k\times k}}} = \langle Z,\mathrm{\small D}F^{*}(X)[\Omega]\rangle_{{\mathbb{R}^{n\times k}}},\;\forall\;\Omega\in \mathrm{Sym}(k),\; Z\in{\mathbb{R}^{n\times k}},\tag{12}\] and show that \(\mathrm{\small D}F^{*}(X)\) is injective. Indeed, since \[\langle\mathrm{\small D}F(X)[Z],\Omega\rangle_{{\mathbb{R}^{k\times k}}} = \langle X^{T}AZ+Z^{T}AX,\Omega\rangle_{{\mathbb{R}^{k\times k}}} = 2\mathrm{tr}(Z^{T}AX\Omega)\\ = \langle Z,2AX\Omega\rangle_{{\mathbb{R}^{n\times k}}},\] in view of 12 , we deduce that \[\label{eq:DF95adjoint2} \mathrm{\small D}F^{*}(X)[\Omega] = 2AX\Omega, for\Omega\in\mathrm{Sym}(k).\tag{13}\] Combining with the fact that \(X^{T}AX =J\) is invertible, it follows that \(\mathrm{\small D}F^{*}(X)\) is injective.

For \(X\in\mathrm{iSt}_{A,J}\), let us denote by \(X_{\perp}\in\mathbb{R}^{n\times(n-k)}\) a full rank matrix such that \(X^{T}X_{\perp}=0\in\mathbb{R}^{k\times(n-k)}.\) Note that \(X_{\perp}\) is not unique. The following technical lemma, similar to [@GSAS21], will be useful later.

Next, we determine the tangent space at a point on the manifold \(\mathrm{iSt}(k,n).\) This is crucial for optimization with manifold constraint since most iterative methods update the iterate along a vector belonging to this space. Recall that the tangent space at \(X\in \mathrm{iSt}(k,n)\) is the set \[T_{X}\mathrm{iSt}(k,n)= \left\lbrace \dot{\varphi}(0) \; |\; \varphi: I\ni 0 \longrightarrow \mathrm{iSt}(k,n)\text{ is a smooth curve and }\varphi(0)=X \right\rbrace.\] The following proposition characterizes the elements of this tangent space.

3 Metric and geometric objects on \(\mathrm{iSt}_{A,J}\)↩︎

A metric \(\mathrm{g}_X\) on the manifold \(\mathrm{iSt}_{A,J}\) is an inner product defined on each tangent space \(T_X\mathrm{iSt}_{A,J}\) which is required to depend smoothly on \(X\). In this case, the pair \((\mathrm{iSt}_{A,J}, \mathrm{g}_X)\) becomes a Riemannian manifold and the metric is referred to as a Riemannian metric. On a Riemannian manifold, a metric allows quantification such as length, distance, angle and, particularly, the definition of orthogonality which is essential in designing many optimization methods. There are numerous ways to select a metric for a given smooth manifold. This freedom is exploited to get advantages in computation in the framework of preconditioning, see [@MishS2016; @DongGGG2022; @ShusA2023; @GaoPY2024], to name a few. To this end, instead of a single one, in this section we will equip \(\mathrm{iSt}_{A,J}\) with a family of tractable metrics, which was recently introduced in [@NguyD23], used for the generalized Stiefel manifold [@ShusA2023] and the symplectic Stiefel manifold [@GaoSS2024], and present the resulting geometric structure.

Let \(\mathbf{M}_{X}\in\mathrm{SPD}(n)\) depend smoothly on \(X\in\mathrm{iSt}(k,n)\). We define an inner product \(\mathrm{g}_{\mathbf{M}_{X}}\) on \(T_{X}\mathrm{iSt}(k,n)\) by \[\label{eq:def95gMx} \mathrm{g}_{\mathbf{M}_{X}}(Z_1,Z_2)=\mathrm{tr}(Z_1^{T}\mathbf{M}_{X}Z_2), \text{ for all } Z_1, Z_2\in T_X\mathrm{iSt}(k,n).\tag{17}\] It can be seen that \(\mathrm{g}_{\mathbf{M}_{X}}\) is indeed a metric not only on \(\mathrm{iSt}(k,n)\) but also on \({\mathbb{R}^{n\times k}}\) by letting \(Z_1, Z_2\) vary in \({\mathbb{R}^{n\times k}}\). As a result, we can define the normal space \(T_{X}\mathrm{iSt}(k,n)^{\perp}_{\mathrm{g}_{\mathbf{M}_{X}}}\) with respect to the inner product \(\mathrm{g}_{\mathbf{M}_{X}}\) of \({\mathbb{R}^{n\times k}}\) as \[T_{X}\mathrm{iSt}(k,n)^{\perp}_{\mathrm{g}_{\mathbf{M}_{X}}} \mathrel{\vcenter{:}}= \left\{ N\in{\mathbb{R}^{n\times k}}: \mathrm{g}_{\mathbf{M}_{X}}(N,Z)=0,\text{ for all } Z\in T_{X}\mathrm{iSt}(k,n)\right\}.\] The following characterization of \(T_{X}\mathrm{iSt}(k,n)^{\perp}_{\mathrm{g}_{\mathbf{M}_{X}}}\) helps formulate the orthogonal projection onto \(T_{X}\mathrm{iSt}(k,n)\).

Next, for each \(X\in \mathrm{iSt}_{A,J}\), considering \(\mathrm{g}_{\mathbf{M}_{X}}\) as an inner product on the ambient space \({\mathbb{R}^{n\times k}}\), it holds that \[{\mathbb{R}^{n\times k}}= T_X\mathrm{iSt}(k,n)\; {\bigoplus}\;T_X\mathrm{iSt}(k,n)_{\mathrm{g}_{\mathbf{M}_{X}}}^{\perp}.\] So, there exists a unique pair \((Y_1,Y_2)\in T_X \mathrm{iSt}(k,n)\times T_X\mathrm{iSt}(k,n)_{\mathrm{g}_{\mathbf{M}_{X}}}^{\perp}\) for each \(Y\) in \({\mathbb{R}^{n\times k}}\) such that \(Y = Y_1 + Y_2.\) This enables us to define two orthogonal projections \[\begin{align} \mathcal{P}_{X}^{\mathbf{M}_{X}}: {\mathbb{R}^{n\times k}}&\longrightarrow T_{X}\mathrm{iSt}(k,n)\\ Y& \longmapsto \mathcal{P}_{X}^{\mathbf{M}_{X}}(Y) = Y_1 \end{align}\] and \[\begin{align} \mathcal{P}_{X}^{\perp_{\mathbf{M}_{X}}}: {\mathbb{R}^{n\times k}}&\longrightarrow T_{X}\mathrm{iSt}(k,n)_{\mathbf{M}_{X}}^{\perp}\\ Y&\longmapsto \mathcal{P}_{X}^{\perp_{\mathbf{M}_{X}}}(Y) = Y_2. \end{align}\] The following proposition provides formulae for the two projections above.

Based on the presented derivation of the orthogonal projection onto the tangent space, we examine the Riemannian gradient. Let \(f:\mathrm{iSt}(k,n)\longrightarrow\mathbb{R}\) be a continuously differentiable function. The Riemannian gradient of \(f\) with respect to the tractable metric \(\mathrm{g}_{\mathbf{M}_{X}}\) is the unique vector \(\mathrm{grad}_{\mathbf{M}_{X}}f(X)\) of \(T_X\mathrm{iSt}(k,n)\) satisfying \[\mathrm{g}_{\mathbf{M}_{X}}\left( \mathrm{grad}_{\mathbf{M}_{X}}f(X),Z\right) = \mathrm{\small D}\bar{f}(X)[Z]= \langle\nabla \bar{f}(X),Z\rangle_{{\mathbb{R}^{n\times k}}},\;\forall\;Z\in T_X\mathrm{iSt}(k,n),\] where \(\bar{f}\) is any smooth extension of \(f\) on a neighborhood of \(X\) in \({\mathbb{R}^{n\times k}},\) and \(\nabla \bar{f}(X)\) is the Euclidean gradient of \(\bar{f}\) at \(X.\) This definition is also used for a function defined on the ambient space \({\mathbb{R}^{n\times k}}\) when one wants to refer to its gradient associated with a metric other than the conventional Euclidean one. Next, we will obtain an explicit formula for the Riemannian gradient.

In many practical cases, the matrix \(\mathbf{M}_{X}\) representing the tractable metric is chosen to be constant. We will refer to this case as the weighted Euclidean metric. When \(\mathbf{M}_{X}\) is further set to be the identity matrix \(I_n\), \(\mathrm{g}_{\mathbf{M}_{X}}\) reduces to the popular Euclidean metric.

A specific choice for the metric that can help avoid having to solve the Lyapunov equation and leads to a different geometric structure and optimization tools will be addressed in an ongoing work.

4 Cayley retraction↩︎

In this section, we present a tool which allows us to update a solution to the problem 9 along a given direction in the tangent space \(T_{X}\mathrm{iSt}(k,n).\) Namely, we define a retraction [@AdleDMMS02; @AbsiMS08] on \(\mathrm{iSt}(k,n)\) to be a smooth mapping \(\mathcal{R}\) from the tangent bundle \(T\mathrm{iSt}(k,n)\mathrel{\vcenter{:}}= \bigcup_{X\in\mathrm{iSt}(k,n)}T_X\mathrm{iSt}(k,n)\) to \(\mathrm{iSt}(k,n)\) satisfying the following conditions for all \(X\in\mathrm{iSt}(k,n)\): \[\begin{align}\; \;\; &\mathcal{R}_X(0) = X, \text{ where } 0 \text{ is the origin of } T_X\mathrm{iSt}(k,n),\tag{22}\\\; \; \;& \left.\frac{\mathrm{d}}{\mathrm{d}t}\mathcal{R}_X\left(tZ\right)\right|_{t=0}= Z \text{ for all } Z\in T_X\mathrm{iSt}(k,n),\tag{23} \end{align}\] where \(\mathcal{R}_X\) is the restriction of \(\mathcal{R}\) to \(T_X\mathrm{iSt}(k,n).\)

In an embedded submanifold, retraction can be formulated based on matrix decomposition [@AbsiMS08] or distance projection [@Boumal23]. In the case of indefinite Stiefel manifold, however, deriving such a retraction does not seem straightforward due to the generality of \(A\) and \(J\). Here, we will construct a retraction based on the Cayley transform for quadratic Lie group [@HairLW06], which will be referred to as the Cayley retraction. This type of retraction has been established for the generalized and the symplectic Stiefel manifolds, see [@WenY2013; @GSAS21; @ShusA2023; @WangDPY2024] for details. To this end, for \(X\in\mathrm{iSt}(k,n)\) and \(Z\in T_{X}\mathrm{iSt}(k,n),\) let us define the matrix \[\label{eq:2identicalSxz} S_{X,Z} = XJZ^{T}AXJX^{T} - XJZ^{T} +ZJX^{T}.\tag{24}\] Taking 14 into account, it can be rewritten as \[\begin{align} \notag S_{X,Z} &=\dfrac{1}{2}XJZ^{T}AXJX^{T}- \dfrac{1}{2}XJX^{T}AZJX^{T}- XJZ^{T} +ZJX^{T}\\ &= \left( I_n - \dfrac{1}{2}XJX^{T}A\right)ZJX^{T} - XJZ^{T}\left( I_n - \dfrac{1}{2}AXJX^{T}\right). \label{eq:2identicalSxz952} \end{align}\tag{25}\] By definition, \(S_{X,Z}\) is skew-symmetric.

It is seen that \(\mathcal{R}_X^{\mathrm{cay}}(Z)\) is defined whenever \(\left( I_n - \frac{1}{2}S_{X,Z}A\right)\) is invertible. When applied to an iterative optimization scheme, the update is performed in the direction of the tangent vector \(Z\in T_{X}\mathrm{iSt}(k,n)\) along the curve \[\label{eq:Cayley95curve} Y^{\mathrm{cay}}(t,Z) \mathrel{\vcenter{:}}= \left(I_n - \dfrac{t}{2}S_{X,Z}A\right)^{-1}\left( I_n + \dfrac{t}{2}S_{X,Z}A\right) X =: \mathrm{\mathrm{cay}}\left( \dfrac{t}{2}S_{X,Z}A\right)X,\tag{26}\] for \(t\geq 0\), on \(\mathrm{iSt}(k,n).\) The following proposition verifies Definition 1 for the Cayley retraction.

The well-definedness of the so-called Cayley curve \(Y^{\mathrm{cay}}(t,Z)\) in 26 depends on the fact whether \(2/t\) with \(t>0\) is an eigenvalue of \(S_{X,Z}A\). When \(A = I_n, J = I_k\), i.e., the Stiefel manifold case, since \(S_{X,Z}A\) is skew-symmetric, it admits only pure imaginary eigenvalues. Therefore, \(2/t\) cannot be an eigenvalue of \(S_{X,Z}A\) for any \(t > 0\). While [@GSAS21] provides a condition for the well-definedness of the Cayley curve for the symplectic Stiefel manifold, we cannot find such a discussion for the generalized Stiefel manifold in the literature. Therefore, we present it here as a special case of the indefinite Stiefel manifold for completeness.

When \(A\) is indefinite, the following counterexample tells us that the Cayley curve \(Y^{\mathrm{cay}}(t,Z)\) is generally not globally defined. Indeed, consider the case when \(A = \left[\begin{smallmatrix} -1&0\\0& 1 \end{smallmatrix}\right], J = -1\). For any \(X = \left[\begin{smallmatrix} \pm\sqrt{\xi^2 + 1}\\\xi \end{smallmatrix}\right]\) with \(\xi\in \mathbb{R}\) and \(Z = \left[\begin{smallmatrix} \xi\\ \pm\sqrt{\xi^2 + 1} \end{smallmatrix}\right] \in T_X\mathrm{iSt}_{A,J}\). Calculation shows that \(S_{X,Z}A = \left[\begin{smallmatrix} 0&1\\1& 0 \end{smallmatrix}\right]\) which admits the eigenvalues \(\pm 1\).

Although the Cayley curve is not defined globally, the following feature is enough for establishing the global convergence of the gradient descent algorithm in section 5.

The main computational burden for the Cayley retraction \(\mathcal{R}^{\mathrm{cay}}\) on \(\mathrm{iSt}(k,n)\) is contributed by the term \(\left( I_n - \frac{1}{2}S_{X,Z}A\right)^{-1}\) which means solving a linear system of equations whose coefficient matrix is of size \(n\times n.\) This task principally requires \(O(n^3)\) flops. In some cases, \(k\) is considerably smaller than \(n\) which makes that cost unreasonably high. It is therefore expected to develop a more efficient alternative. In the same framework as in [@OviH23; @BendZ21] for the symplectic Stiefel manifold and [@WangDPY2024] for the generalized Stiefel manifold, we obtain the following result.

To close this section, we discuss the possibility that the Cayley retraction ?? is of second order. It can be interpreted from [@Boumal23] that a retraction \(\mathcal{R}\) is of second order if, for any \(X\) and \(Z\in T_X\mathrm{iSt}_{A,J}\), the second derivative of its curve \(\frac{\mathrm{d}^2}{dt^2}(\mathcal{R}_X(tZ))\) belongs to the normal space \(T_X\mathrm{iSt}_{A,J}^{\perp}\). Second-order retraction helps compute the Riemannian Hessian of the cost function at any point without invoking the Riemannian connection [@Boumal23]. Differentiating the identity 27 yields \[\dfrac{\mathrm{d}^2}{\mathrm{d}t^2}\mathcal{R}_X^{\mathrm{cay}}(tZ) = 2 \left(\left(I_{n}-\dfrac{t}{2}S_{X,Z}A\right)^{-1}\dfrac{1}{2}S_{X,Z}A\right)^2\left(\mathcal{R}_X^{\mathrm{cay}}(tZ) + X\right).\] Particularly at \(t=0\), from 24 and the fact \(S_{X,Z}AX=Z\), we have \[\begin{align} \left.\dfrac{\mathrm{d}^2}{\mathrm{d}t^2}\mathcal{R}_X^{\mathrm{cay}}(tZ)\right|_{t=0} &= S_{X,Z}AS_{X,Z}AX =S_{X,Z}AZ\\ &=\left(XJZ^{T}AXJX^{T} - XJZ^{T} +ZJX^{T}\right)AZ. \end{align}\]

With the example below, we point out that \(\left.\frac{\mathrm{d}^2}{\mathrm{d}t^2}\mathcal{R}_X^{\mathrm{cay}}(tZ)\right|_{t=0}\) is generally not a normal vector at \(X\) which in turn implies that the Cayley retraction is not of second order. Let us consider a point \(X\in \mathrm{iSt}_{A,J}(2,3)\) with \[A = \begin{bmatrix}\begin{array}{ccc} \dfrac{-7}{3} & \dfrac{-2}{3} & \dfrac{4}{3}\\[0.8em] \dfrac{-2}{3} & \dfrac{-23}{15} & \dfrac{-14}{15}\\[0.8em] \dfrac{4}{3} & \dfrac{-14}{15} & \dfrac{-2}{15} \end{array} \end{bmatrix},\; J=\begin{bmatrix} 0&-1\\ -1&0 \end{bmatrix},\; X = \begin{bmatrix}\begin{array}{cc} \dfrac{\sqrt{5}-\sqrt{3}}{6} & \dfrac{\sqrt{5}+\sqrt{3}}{6}\\[0.8em] \dfrac{\sqrt{5}+5\sqrt{3}}{30} & \dfrac{\sqrt{5}-5\sqrt{3}}{30}\\[0.8em] -\dfrac{\sqrt{5}+5\sqrt{3}}{15} & \dfrac{-\sqrt{5}+5\sqrt{3}}{15} \end{array} \end{bmatrix}.\] It follows from 15 that, for \(W= \left[\begin{smallmatrix} 1&0\\ 0&-1 \end{smallmatrix}\right],\,X_{\perp}=\begin{bmatrix} 0&2&1 \end{bmatrix}^T, \text{ and } K=\begin{bmatrix} 1&0 \end{bmatrix},\) \[Z \mathrel{\vcenter{:}}= XW + A^{-1}X_{\perp}K = \begin{bmatrix}\begin{array}{cc} \dfrac{\sqrt{5}-\sqrt{3}}{6} & -\dfrac{\sqrt{5}+\sqrt{3}}{6}\\[1em] \dfrac{\sqrt{5}+5\sqrt{3}-30}{30} & \dfrac{-\sqrt{5}+5\sqrt{3}}{30}\\[1em] -\dfrac{2\sqrt{5}+10\sqrt{3}+15}{30}& \dfrac{\sqrt{5}-5\sqrt{3}}{15} \end{array} \end{bmatrix}\in T_X\mathrm{iSt}_{A,J}.\] Therefore, \[S_{X,Z} AZ= \begin{bmatrix}\begin{array}{cc} \dfrac{-3\sqrt{5}-7\sqrt{3}}{12} & \dfrac{\sqrt{5}+\sqrt{3}}{6}\\[1em] \dfrac{-3\sqrt{5}+35\sqrt{3}-60}{60} & \dfrac{\sqrt{5}-5\sqrt{3}}{30}\\[1em] \dfrac{3\sqrt{5}-35\sqrt{3}-15}{30}& \dfrac{-\sqrt{5}+5\sqrt{3}}{15} \end{array} \end{bmatrix}.\] Since \(JX^TS_{X,Z} AZ= \left[\begin{smallmatrix} 2&-2/3\\ -3/2&1/3 \end{smallmatrix}\right]\) is not symmetric, by contradiction it results in the fact that \(S_{X,Z}AZ\neq AX\hat{W},\) for all \(\hat{W}\in\mathrm{Sym}(2).\) That is, in view of Proposition [prop:normalspace95gM], \(S_{X,Z}AZ\) does not belong to the normal space \(T_X\mathrm{iSt}_{A,J}^{\perp}\) with respect to the Euclidean metric.

5 A gradient descent method↩︎

As the necessary tools are available, we are now ready to present a Riemannian gradient descent method for solving the problem 9 . Starting with an \(X_0\in \mathrm{iSt}_{A,J}\), this method generates a sequence \(\{X_j\}\) using a line search along the direction \(-\mathrm{grad}_{\mathbf{M}_{X}}f(X_j)\) and updates the solution as \[X_{j+1} = \mathcal{R}_{X_j}^{\mathrm{cay}}(-\tau_j \mathrm{grad}_{\mathbf{M}_{X}}f(X_j)),\] where \(\tau_j >0\) is a carefully chosen step size. The choice for step size undoubtedly has a great impact on the method’s performance. Here, we adopt a combination of a nonmonotone line search technique [@ZhanH2004] with the alternating Barzilai-Borwein [@BarB88] which has been adapted to Riemannian optimization in [@IannP18] and simplified in [@HuLWY2020] for embedded submanifolds, see Algorithm 1 for details and [@GSAS21] for more choices of step size. One can see from [eq:nonmonotone95cond] that if \(\alpha\) is chosen to be zero, this scheme reduces to the popular monotonically decreasing line search with Amijo’s rule.

Given Lemma 3, according to [@GSAS21], the global convergence of Algorithm 1 is guaranteed. Namely, regardless of the initial guess \(X_0\), every accumulation point \(X_*\) of the sequence \(\{X_j\}\) generated by Algorithm 1 is a critical point of the cost function \(f\) in 9 , i.e., \(\mathrm{grad}_{\mathbf{M}_{X}}f(X_*)=0\). The material presented in this section is theoretically known. Nevertheless, to the best of our knowledge, it is the first algorithmic development for Riemannian optimization on the indefinite Stiefel manifold. Despite of simplicity, the algorithm will prove its efficiency via numerical tests in the next section. To finish this section, we give a rough estimate for the computational complexity of Algorithm 1. To this end, we provide in Table 1 the complexity of its main steps. In the computation of the Riemannian gradient, if we choose the Euclidean metric, the term \(\frac{4}{3}n^3\) will disappear, but it does not necessarily result in a shorter execution time since the total runtime also depends on the number of iterations niter and the number of the backtracking line searches nls in each of them. If \(A\) is sparse, e.g., diagonal, the \(4n^2k\) becomes \(O(n)\). Meanwhile, the retraction formula 26 always costs \(O(n^3)\) and is obviously more expensive than 31 and 28 . From this table, one can give a rough statement on the cost of Algorithm 1. For instance, using the retraction 28 , the total cost is \[\begin{align} \begin{aligned} \label{eq:roughcost} &\texttt{niter}\left[\left(\frac{4}{3}n^3 + 4n^2k +8nk^2 -2nk + 18k^3\right) + 6nk\right.\\ &\left. \qquad+ (8nk^2 -nk + 4k^3 - 4k^2)+ \texttt{nls}\left(nk^2 + 4nk +\frac{2}{3}k^3 + 4k^2 + \texttt{fval}\right)\right], \end{aligned} \end{align}\tag{32}\] where fval is the cost to evaluate the cost function.

6 Numerical examples↩︎

In this section, we run several examples to verify the theoretical results obtained previously and to better understand the overall numerical behavior of the method. First, the parameters in Algorithm 1 need to be set up. Most of them are fixed for all examples. Namely, for the backtracking search, we set \(\beta = 1\mathrm{e}{-4}\), \(\delta = 5\mathrm{e}{-1}\), \(\gamma_0 = 1\mathrm{e}{-3}\), \(\gamma_{\min} = 1\mathrm{e}{-15}\), \(\gamma_{\max} =1\mathrm{e}{+5}\). Moreover, \(\alpha = 0.85\) is used in the nonmonotone condition. The stopping tolerance \(\texttt{rstop}\), which indicates a convergence when \(\|\mathrm{grad}_{\mathbf{M}_{X}}f(X_j)\|_{X_j}\leq\texttt{rstop}\|\mathrm{grad}_{\mathbf{M}_{X}}f(X_0)\|_{X_0}\), is usually set to be either \(1\mathrm{e}{-9}\) or \(1\mathrm{e}{-6}\) upon example.

For every test, if the Euclidean Hessian of the extended cost function \(\nabla^2 \bar{f}(X)\) is positive-definite, we will take it to define the tractable metric as an adoption of the suggestion in [@DongGGG2022; @GaoPY2024]. The hint from [@ShusA2023] employing \(A\) is not applied here because \(A\) in the context of this paper is generally not spd. We will examine three schemes for the Cayley retraction to see the difference in efficiency and feasibility. We will also compare our proposed method with the augmented Lagrangian method (ALM) [@NoceW06]. This method has recently been reported to deliver good results when applied to the optimization problem under the symplecticity constraint [@OviedoLara2025]. Specifically applied to our present problem, at step \(j\) of the main iteration, the method approximately solves the subproblem with the stopping criterion \(\|X_j^TAX_j-J\|_F/(j + 10)\) using the spectral gradient method as in [@OviedoLara2025]⁴. After this step, we update the multiplier and penalty parameter as \(\Lambda_{j+1} = \Lambda_{j}+\mu_j(X_{j+1}^TAX_{j+1}-J)\) and \[\mu_{j+1} = \left\{ \begin{array}{ll} &\mu_j, if\|X_{j+1}^TAX_{j+1}-J\|_F\leq \tau\|X_j^TAX_j-J\|_F or\mu_j\geq 1\mathrm{e}+14,\\ &\rho\mu_{j}, \end{array}\right.\] respectively, where \(\tau\in (0,1), \rho > 1\). Usually, we initialize the parameters as: Lagrange multiplier \(\Lambda_0 = 0, \mu_0 = 1\), and \(X_0\) is the normalized matrix of 1’s and set \(\tau = 0.25, \rho = 10\) as recommended in [@SunYuan2006]. We stop the main iteration when the Euclidean norm of the associated Lagrangian function meets the relative convergence condition, specified by \(\texttt{rstop}_{ALM}\), in the same manner as the proposed Riemannian gradient descent (RGD) method and the feasibility violation does not exceed the tolerance feastol. Additionally, we also point out that the step size strategy described in Algorithm 1 is helpful by comparing the result with that using a fixed step size. The numerical experiments are performed on a laptop with Intel(R) Core(TM) i7-4500U (at 1.8GHz, 5MB Cache, 8GB RAM) running MATLAB R2023a under Windows 10 Home. The main code is made available at https://sites.google.com/view/ntson/code.

6.1 Trace minimization for positive-definite matrix pencil↩︎

We consider in this section the minimization problem 9 where \(f(X) = \mathrm{tr}\left(X^TMX\right)\), \(M\in \mathrm{SPD}(n)\), \(A\in \mathrm{Sym}(n)\) is nonsingular with \(\mathrm{i}_+(A)=p, \mathrm{i}_-(A) = m\), and \(p+m = n\), \(J = \mathrm{diag}(I_{k_p},-I_{k_m})\) and \(k_p+k_m = k, k_p\leq p, k_m \leq m\). In view of [@KovaV95], recalling 6 , a global solution \(X_*\) can be used to compute \(k_p\) positive and \(k_m\) negative eigenvalues having the smallest magnitude of the matrix pencil \((M,A)\). Indeed, it follows that \(MX_* = AX_*\Lambda\) where \(\Lambda = \mathrm{diag}(\Lambda_+,\Lambda_-)\) is determined as in the theorem mentioned. This equality leads to \(X_*^TMX_* = \mathrm{diag}(\Lambda_+,-\Lambda_-)\). Finally, the sought eigenvalues can be determined as that of the \(k_p\times k_p\) matrix \(\Lambda_+\) and \(k_m\times k_m\) matrix \(\Lambda_-\). In each test, we present the value of the objective function (obj.), the norm in the chosen metric of the Riemannian gradient (grad.), the relative error in computing the mentioned eigenvalues (eig. rel. err.) \(\|MV - AVD\|_F/\|AVD\|_F\) in which \(V\) is a matrix of (generalized) eigenvectors associated with the computed eigenvalues listed in the diagonal of the diagonal matrix \(D\), the feasibility (feas.) \(\|X_*^TAX_*-J\|_F\), the number of iterations (\(\#\)iter.), the number of function evaluations (\(\#\)eval.), and the CPU time (CPU) needed to obtain \(X_*\).

In the first context, we choose the Lehmer matrix of size \(n = 200\) as the matrix \(M\), and set \(A = \mathrm{diag}(1,\ldots,p,-m,\ldots,-1)\) with \(p=150\) and \(m=50\). We consider both Euclidean and weighted Euclidean metrics for two triplets \((k,k_p,k_m) = (5,3,2)\) and \((20,15,5)\) which results in four instances. With rstop\(=1\mathrm{e}-9\), the obtained numerical results, noted as RGD, are given in Table 2. It is obvious to see that using an appropriately chosen metric can significantly speed up the iterative scheme. Surprisingly, it also helps reduce the error of eigenvalue and eigenvector computations. In the forthcoming tests, we only use an individually chosen weighted Euclidean metric. We also run the ALM for this model with \(\texttt{rstop}_{ALM} = 1\mathrm{e}-7\) and \(\texttt{feastol} = 1\mathrm{e}-10\). This specific convergence condition and the maximal numbers of the main and inner iterations, 10 and 5000, respectively, are set so that the accuracy of the obtained result can be comparable to that of the RGD method. Note that the norm of the Riemannian gradient of the cost function, which is indeed quite the same as the norm of the gradient of the Lagrangian function, is displayed to make the comparison appropriate.

We continue with this model to check the scalability of the proposed algorithm. To this end, we let \(n=200, 400,\ldots,6400\), \(p = 3n/4, m = n/4\), while we fix \(k=10, k_p=k_m=5\). In Table 3, we report some measures versus the problem size. One can observe that while the gradient norm and the feasibility violation are quite stable, the number of iterations as well as the function evaluations increase predictably. Specifically, the CPU time indicates that the computational complexity of Algorithm 1 is approximately \(O(n^3)\), which basically matches the estimate given in 32 .

Next, we numerically examine the constant step size strategy for the proposed algorithm using the same model problem with \(n = 5000\) and \(\texttt{rstop} = 1\mathrm{e}-6\). For unconstrained optimization of quadratic functions, it was confirmed in [@Bertsekas1999] that a convergence to a solution is guaranteed for the gradient descent method if the fixed step size is smaller than \(2/\lambda_{\max}\), where \(\lambda_{\max}\) is the largest eigenvalue of the Hessian. Our situation is a bit more involved, but this hint can provide a range for the tests. Since the recommended upper bound in our cases is approximately \(7\mathrm{e}-4\), we run the RGD method with fixed step sizes, namely, \(5\mathrm{e}-6, 1\mathrm{e}-5, 5\mathrm{e}-5, 1\mathrm{e}-4, 5\mathrm{e}-4, 1\mathrm{e}-3, 5\mathrm{e}-3\) and compare the results with our Algorithm 1. As expected, only Algorithm 1 delivers a convergence after 242 iterations. We report the running history in Figure [fig:fixstep]. One can see that, except for a slight loss of feasibility and about 1.5 times slower per iteration, the non-monotone line search strategy in Algorithm 1 is still much more efficient than a fixed step size in reaching a solution.

of the step size strategy used in Algorithm :non-monotone32gradient with different fixed step sizes (specified on the top) on the model problem described in Table :scalability with \(n=5000\) and \(\textttrstop = 1\mathrm{e}-6\). :fixstep

Now, we proceed with different positive-definite matrices in the MATLAB gallery with a similar setting. Namely, we consider the cost function 9 determined by different matrices \(M\) with \((n,p,m) = (2000,1000,1000)\) and \((k,k_p,k_m)=(10,5,5)\), \(A = \mathrm{diag}(1,\ldots,p,-1,\ldots,-m)\) and \(\texttt{rstop}=1\mathrm{e}-9\). For ALM, we set \(\texttt{rstop}_{ALM} = 1\mathrm{e}-7\), \(\texttt{feastol} = 1\mathrm{e}-8\), with maximal 20 and 1000 iterations for the main and inner loops, respectively. More details and the obtained results are described in Table 4. One can easily see that the Cayley retraction formulae 31 and 28 are faster in yielding the solution, while the conventional one 26 is always the best in preserving the feasibility. Compared to RGD, the ALM method delivers less acurate result, especially the gradient norm and the relative error in eigenvalue computation, although it consumes more time. Particularly, it does not yield a solution in the case of the \(\texttt{tridiag}\) matrix.

Now, we move our attention to the LREVP in terms of the minimization problem 8 . Here, we use data of larger sizes, stemming from practice. In the first test, following [@BaiL14], we make use of the University of Florida sparse matrix collection [@DavisHu2011]. Namely, the matrix bcsstk21, which is of size \(3600\times 3600\), and the leading principal submatrix of the same size of the matrix sts4098 are assigned to \(K\) and \(M\), respectively. The block structure makes the size of the problem 7200. We aim at computing the smallest four positive eigenvalues, and therefore the matrix \(J\) is chosen as \(I_4\). For the RGD method, parameter setting is as before, while an initial guess \(X_0\) can be simply chosen as \(\left[\begin{smallmatrix} \tilde{V}\\\tilde{V} \end{smallmatrix}\right]\), where \(\tilde{V}\) is any \(n\times k\) random matrix with orthonormal columns divided by \(\sqrt{2}\). In our test, we start with an \(n\times k\) random matrix, preset by rng(1,‘philox’), and then orthogonalize the obtained matrix. With \(\texttt{rstop} = 1\mathrm{e}-9\), the convergence is reached after 45 iterations. The condition \(\texttt{rstop} = 1\mathrm{e}-11\) is fulfilled at the 108-th iteration. Meanwhile, after trying various sets of paramters and initial guesses, we observe that ALM helps reduce the norm of the associated Lagrangian function very fast at the beginning and usually stagnates when it reaches \(1\mathrm{e}+1\), which is too far from yielding a reliable approximation for the eigenvalues. Therefore, in Figure [fig:Florida] and Table 5, we report the results obtained using the proposed RGD only. One can see that, with the tolerance \(\texttt{rstop}=1\mathrm{e}-11\), our method can deliver a comparably accurate solution as in [@BaiL14]⁵

history of the LREVP with data specified in Table :Florida. :Florida

Another model for the LREVP is considered next where the matrices \(K\) and \(M\) are derived from electronic structure calculations, which were also used in [@KresPS2014; @Pand2019; @BaiL13]. Here, for the problem of size \(n=11320\), the stopping tolerance rstop is set as \(1\mathrm{e}-6\) while other parameters are the same as before. Results are given in Table 6. Note however that our use of this example is merely as a model with practical data; the eigenvalue computation scheme proposed in [@KresPS2014] is apparently more efficient. Indeed, a test using the retraction 26 with rstop \(=1\mathrm{e}-8\) delivers four approximate eigenvalues that can match about the first ten digits of the reference eigenvalues \[\begin{align} &\lambda^+_1 \approx 0.541812517132466,\quad \lambda^+_2 \approx 0.541812517132473,\\ &\lambda^+_3 \approx 0.541812517132498,\quad \lambda^+_4 \approx 0.615143209274579, \end{align}\] which is a bit worse than the approximation reported in [@KresPS2014].

6.2 Matrix least square problems↩︎

We consider minimization of the cost function of the form \[\label{eq:fitting} f(X) = \|GX-B\|^2_F.\tag{33}\] This form is an approach for various problems. First, with suitable specification, it is the well known Procrustes problem [@Scho1966; @High1988; @AndeE1997; @BojaL1999; @BhatC2019]. In this subsection, to intimate the conventional situation [@Scho1966], we examine the problem on the \(J\)-orthogonal group. That is, given \(G, B\in {\mathbb{R}^{l\times n}},\, l\geq n\), find a matrix \(X\in \mathrm{iSt}_{J,J}\) with \(J = \mathrm{diag}(\pm 1,\ldots,\pm 1) \in {\mathbb{R}^{n\times n}}\) that minimizes \(f\) defined in 33 . It is worth noting that in a recent work [@HeYWX2024], a block coordinate descent method was developed. To set up, let us consider \(J = (I_p,-I_m)\) with \(p+m=n\). Then \(G\) is randomly generated with rng default in MATLAB. A global solution is prescribed as \(V = \mathrm{diag}(V_1,V_2),\) where \(V_1\in \mathrm{O}(p), V_2\in \mathrm{O}(m)\) are initialized randomly with rng(1,‘twister’) and then orthogonalized by orth. This results in the setting that \(B = GV\). The identity matrix \(I_n\) is chosen as the initial guess and \(\texttt{rstop} = 1\mathrm{e}-6\) for the RGD method; naturally, only the Cayley retraction 26 is used. For the ALM method, we set \(\texttt{rstop}_{ALM} = 1\mathrm{e}-9\) and \(\texttt{feastol} = 1\mathrm{e}-12\) while the maximal number of the main and the subproblem iterations are 20 and 100, respectively. Table 7 displays the results for different sizes of problem with \(p = 3n/4\) in which diff. represents the difference between the prescribed and the numerically obtained solutions. Apparently, the iterate generated by the RGD method converges to the prescribed global minimizer only in the cases \(l=n=1000\) and \((l,n)=(1500,1000)\); this is normal since the problem can have more than one global solution and/or the obtained solution might not be a global one. Overall, the generated iterate converges as expected. Surprisingly, the ALM method offers better results in the case of \((l,n)=(1000,500)\).

As the second application, we consider a matrix equation \[\label{eq:mateq} GX = B, G\in \mathrm{SPD}(n),\tag{34}\] on the indefinite Stiefel manifold \(\mathrm{iSt}_{A,I_k}\) with nonsingular \(A\in \mathrm{Sym}(n)\). If there is no information about the existence of a solution, a natural idea is to minimize the associated cost function of the form 33 . One can see that if \(B \in \mathrm{iSt}_{G^{-1}AG^{-1},I_k}\) then 34 admits a unique solution \(G^{-1}B\) which is also the unique global solution of the associated optimization problem. When \(G = I_n\), it is the nearest matrix problem considered in [@GSAS21]; some minimal distance and averaging problems can also be converted into this form, see, e.g., [@FioriT2009; @Fior11].

In the first scenario, most of the numerical data are fixed as follows. Let \((n,p,m) = (4000,3000,1000)\) and \(k=10\). First, let MATLAB generate an \(n\times n\) random matrix with default setting and then this matrix is orthogonalized to yield \(V\) whose columns are denoted by \(v_j,j=1,\ldots,n\). Next, we set \(A = V\mathrm{diag}(1,\ldots,p,-m,\ldots,-1)V^T\). The right hand side is set as \(B = G[v_1,\ldots,v_k]\mathrm{diag}(1/\sqrt{1},\ldots,1/\sqrt{k})\). A set of \(k\) orthogonal eigenvectors corresponding to \(k\) positive eigenvalues of \(A\) divided by the square root of the associated eigenvalue is used as an initial guess \(X_0.\) Table 8 reports the result of running different matrices \(G\) with the initial guess to be the scaled orthogonal eigenvectors of \(A\) associated with \(k\) eigenvalues of the largest magnitudes obtained by eigs in MATLAB. The same quantities as in the previous examples are displayed.

We try various sets of parameter values and initial guesses with the ALM method, but the obtained result is unsatisfactory and therefore not reported here.

In order to practically verify the global convergence of Algorithm 1, we reuse the above model setting with random data. Namely, \(V\) is randomly generated without any prescribed generator; \(X_0\) is the set of \(k\) randomly picked columns among the first \(p\) columns of \(V\) using randsample\((p,k)\) in MATLAB. Then, for each model matrix, we run ten tests. The average of quantities obtained from these tests are given in Table 9. Apparently, the iterate generated by the proposed algorithm always converges to the global minimizer of the cost function which is also the solution to the matrix equation.

From the examples presented in both two subsections 6.1 and 6.2, it can be concluded that (i) the algorithm with all of its components proposed for minimization on the indefinite Stiefel manifold works as expected: it converges to a solution regardless of the initial guess. In most cases, it apparently delivers better results than that obtained by the classical augmented Lagrangian method for the equality constrained problem. Moreover, (ii) an appropriately chosen metric, enabled by the consideration of a tractable metric, always results in a faster solution: fewer iterations and CPU time. And finally, (iii) the three formulae for the Cayley retraction all work properly: the conventional one 26 is suitable for large \(k\), e.g., the indefinite orthogonal group, while the two others 31 and 28 are advisable for small \(k\) as they help replace inversion of an \(n\times n\) matrix by that of a \(k\)-sized matrix. Using these two schemes can sometimes result in a slight loss of feasibility. This is most probably due to the occasional instability of Sherman-Morrison-Woodbury formula as reported in, e.g., [@FineS2001; @HaoS2021]; it can be restored using a suitable Gram-Schmidt algorithm, upon the individual case.

7 Conclusion↩︎

We addressed the minimization problem subject to a quadratic constraint of the form \(X^TAX = J\) with symmetric, nonsingular matrix \(A\) and the matrix \(J\) satisfying \(J^2 = I\), in the Riemannian framework. Our result provided a unified framework for Riemannian optimization on the orthogonal and generalized Stiefel manifolds, \(J\)-orthogonal group, and more than that.

The numerical examples showed that, on the one hand, the proposed Riemannian gradient descent scheme is a highly capable method, especially for the case that an appropriate metric is used, but on the other hand, it may suffer from slow convergence. More robust algorithms such as conjugate gradient or an efficient way to exploit the second-order information of the cost function beyond the tractable metric might be considered in future work.

Acknowledgement↩︎

The authors would like to thank Dr. Marija Miloloža Pandur for providing the data used in the last example of subsection 6.1. Dinh Van Tiep’s research was supported by Thai Nguyen University of Technology.

---

1. Faculty of Fundamental and Applied Sciences, Thai Nguyen University of Technology, 24131 Thai Nguyen, Vietnam (tiep.dv@tnut.edu.vn).↩︎

2. Department of Mathematics and Informatics, Thai Nguyen University of Sciences, 24118 Thai Nguyen, Vietnam (ntson@tnus.edu.vn).↩︎

3. There is a mixed use of terms for “Stiefel-type" manifolds. The set \(\mathrm{iSt}_{A,I_k}\), which is a special case of \(\mathrm{iSt}_{A,J}\) in 10 , was named the indefinite Stiefel manifold in [@Kobayashi1992], pseudo-Stiefel manifold in [@Chikuse2003] and generalized Stiefel manifold in [@Li1993; @SenadoMendoza2022]. At the same time, as mentioned previously, the generalized Stiefel manifold means \(\mathrm{iSt}_{A,I_k}\) with \(A\) is spd in the optimization community and in this work. The set \(\{X\in {\mathbb{R}^{n\times k}}:\;X^TX = C\}\), where \(C\) is spd, was called the Stiefel C-manifold in [@Downs1972].↩︎

4. In fact, we adapt the MATLAB code, provided by the authors at https://www.mathworks.com/matlabcentral/fileexchange/181008, to our problem.↩︎

5. Note that [@BaiL14] considered the generalized LREVP including the LREVP as a special case.↩︎