Rank-one Riemannian Subspace Descent for Nonlinear Matrix Equations

Problem Statement
The paper addresses the computational challenge of finding Symmetric Positive Definite (SPD) solutions to nonlinear matrix equations (NMEs) arising in control theory, dynamic programming, and stochastic filtering. These equations, which include Continuous- and Discrete-time Algebraic Riccati Equations (CARE/DARE) and Lyapunov equations, are often formulated as finding $X \in \mathcal{P}_n$ such that $G(X, A, C, \dots) = 0$. This is recast as a nonlinear least-squares problem: $\min_{X \in \mathcal{P}_n} \|G(X, A, C, \dots)\|_F^2$.

The primary bottleneck is scalability. Standard solvers for dense matrix equations typically incur $O(n^3)$ floating-point operations (flops) per iteration. While $O(n^2)$ methods exist, they generally rely on assumptions of sparsity or low-rank solutions, which do not hold for general dense problems. Consequently, for large state dimensions (e.g., $n=10,000$), standard dense solvers become computationally infeasible due to time and memory constraints.

Methodology: Rank-one Riemannian Subspace Descent (R1RSD)
The authors propose a Rank-one Riemannian Subspace Descent (R1RSD) algorithm that operates directly on the manifold of SPD matrices ($\mathcal{P}_n$) equipped with the affine-invariant Riemannian metric.

1. Subspace Selection via Spectral Decomposition: Unlike traditional subspace descent methods that project gradients onto fixed coordinate subspaces, R1RSD decomposes the transformed Riemannian gradient into a sum of rank-one eigen-components. The algorithm selects the dominant eigen-component (corresponding to the largest eigenvalue magnitude) to define the descent direction.
2. Power Method for Direction Identification: To avoid the $O(n^3)$ cost of a full eigen-decomposition, the dominant component is identified using the power method. The authors show that a small number of power iterations (specifically $r = O(\log n)$) is sufficient to approximate the dominant direction effectively.
3. Rank-one Updates and Complexity: The update rule is derived using the matrix exponential of a rank-one matrix, which admits a closed-form rank-one update. By maintaining the Cholesky factor $B_t$ of $X_t$ and its inverse via rank-one updates (using the Woodbury identity and cholupdate), the algorithm avoids dense matrix multiplications.
4. Gradient-Vector Products: For objectives defined as compositions of standard matrix operations (addition, multiplication, inverse, log-determinant), the gradient-vector product $\nabla f(X)v$ is computed in $O(n^2)$ time by traversing the computational graph and applying matrix-vector products, avoiding explicit formation of the full gradient matrix.
5. Step-size Selection: The rank-one structure enables exact line search in many cases with minimal additional cost, though fixed or Armijo step sizes are also viable.

Key Contributions

• Algorithmic Innovation: The paper introduces a method that updates iterates along a single, carefully selected rank-one direction derived from the Riemannian gradient, differing from methods that project onto fixed subspaces.
• Complexity Reduction: The proposed algorithm achieves $O(n^2 \log n)$ arithmetic cost per iteration (dominated by $r$ power iterations of $O(n^2)$ each) and requires only $O(n)$ iterations under standard assumptions. This contrasts with state-of-the-art subspace descent methods (e.g., [43]) that have $O(n^2)$ per-iteration cost but require $O(n^2)$ total iterations, and standard solvers requiring $O(n^3)$ per iteration.
• Theoretical Guarantees: The authors provide convergence analysis for both geodesically smooth (non-convex) and geodesically strongly convex cases. They establish an $O(n)$ iteration bound for the general case and an improved bound for the strongly convex case, with total flop counts of $O(n^3 \log n)$ and $O(n^3 \kappa \log n)$ respectively.
• Practical Implementation: The work includes a publicly available MATLAB implementation that utilizes low-rank updates to maintain efficiency, ensuring the $O(n^2)$ per-iteration complexity is realized in practice.

Numerical Results
The authors benchmark R1RSD against state-of-the-art solvers, including MATLAB's icare, Structure-Preserving Doubling Algorithms (SDA), and the Burer-Monteiro Factorization-based Coordinate (BMFC) subspace descent method.

• Scalability: Experiments on CARE, DARE, and other NMEs with dimensions up to $n=10,000$ demonstrate that standard solvers (SDA, icare) and BMFC fail to return solutions or run out of memory for large $n$. R1RSD successfully solves these instances.
• Efficiency: For $n=100$, R1RSD converges in significantly fewer iterations (1–2 orders of magnitude) compared to BMFC.
• Robustness: The algorithm remains effective even when the optimal step size varies significantly, and the greedy subspace selection (using $\approx 10$ power iterations) provides a strong balance between convergence speed and per-iteration cost.

Significance and Claims
The paper claims that R1RSD provides a practical approach for high-dimensional and dense SPD-constrained matrix equations where traditional dense solvers are prohibitive. By reducing the per-iteration complexity to $O(n^2)$ while maintaining a reasonable iteration count, the method bridges the gap between theoretical scalability and practical applicability for large-scale control problems. The authors note that this approach is particularly relevant for time-varying and adaptive control scenarios (such as State-Dependent Riccati Equation frameworks) where solving Riccati equations at every time step is otherwise infeasible. The work emphasizes that rank-one manifold updates can effectively handle problems that standard solvers cannot, without relying on sparsity or low-rank assumptions on the solution or coefficients.