Rank-one Riemannian Subspace Descent for Nonlinear Matrix Equations

TL;DR

This work tackles the challenge of solving large-scale symmetric positive definite nonlinear matrix equations by recasting the problem as residual minimization on the SPD manifold. It introduces Rank-one Riemannian Subspace Descent (R1RSD), which greedily updates along a dominant eigen-component of a transformed gradient, identified via a small number of power iterations, and exploits rank-one matrix updates for efficiency. The method achieves ${\mathcal{O}}(n^2)$ per-iteration cost with ${\mathcal{O}}(n)$ iterations under standard assumptions, and shows strong scalability in CARE, DARE, and other NMEs up to $n=10{,}000$, outperforming or matching state-of-the-art solvers under practical budgets. The combination of low per-iteration cost, favorable convergence in practice, and public MATLAB code positions R1RSD as a viable tool for high-dimensional SPD-constrained matrix equations in control and dynamic filtering contexts.

Abstract

We propose a rank-one Riemannian subspace descent algorithm for computing symmetric positive definite (SPD) solutions to nonlinear matrix equations arising in control theory, dynamic programming, and stochastic filtering. For solution matrices of size $n\times n$, standard approaches for dense matrix equations typically incur $\mathcal{O}(n^3)$ cost per-iteration, while the efficient $\mathcal{O}(n^2)$ methods either rely on sparsity or low-rank solutions, or have iteration counts that scale poorly. The proposed method entails updating along the dominant eigen-component of a transformed Riemannian gradient, identified using at most $\mathcal{O}(\log(n))$ power iterations. The update structure also enables exact step-size selection in many cases at minimal additional cost. For objectives defined as compositions of standard matrix operations, each iteration can be implemented using only matrix--vector products, yielding $\mathcal{O}(n^2)$ arithmetic cost. We prove an $\mathcal{O}(n)$ iteration bound under standard smoothness assumptions, with improved bounds under geodesic strong convexity. Numerical experiments on large-scale CARE, DARE, and other nonlinear matrix equations show that the proposed algorithm solves instances (up to $n=10{,}000$ in our tests) for which the compared solvers, including MATLAB's \texttt{icare}, structure-preserving doubling, and subspace-descent baselines fail to return a solution. These results demonstrate that rank-one manifold updates provide a practical approach for high-dimensional and dense SPD-constrained matrix equations. MATLAB code implementation is publicly available on GitHub : \href{https://github.com/yogeshd-iitk/nonlinear_matrix_equation_R1RSD}{\textcolor{blue}{https://github.com/yogeshd-iitk/nonlinear\_matrix \_equation\_R1RSD}}

Figures (5)

• Figure 1: Per-Iteration Complexity of R1RSD, BMFC, and SDA for DARE. The number of power iterations performed at each iteration is denoted by $r$.
• Figure 2: BMFC coordinate selection: cyclic vs sampling with replacement for problems of size $n = 100$.
• Figure 3: Performance of R1RSD for different number of power iterations ($r$)
• Figure 4: R1RSD step size $\beta$ selected via line search at each iteration for two instances of NME
• Figure 5: R1RSD step size selection: line search vs. fixed

Theorems & Definitions (7)

• Definition 1: Directional derivative
• Definition 2: Riemannian gradient
• Lemma 1
• Theorem 1
• proof
• Theorem 2
• proof