Rank-one Riemannian Subspace Descent for Nonlinear Matrix Equations

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}}

💡 Research Summary

The paper introduces a novel algorithm called Rank‑One Riemannian Subspace Descent (R1RSD) for solving large‑scale symmetric positive definite (SPD) nonlinear matrix equations (NMEs) that arise in control theory, dynamic programming, and stochastic filtering. The authors start by reformulating a generic NME G(X)=0 as a nonlinear least‑squares problem ‖G(X)‖²_F, where the variable X belongs to the SPD manifold Pₙ. This casts the problem as a Riemannian optimization task on a curved space, allowing the use of geometric tools that respect the SPD constraint.

A key observation is that the Riemannian gradient grad_R f(X) can be expressed, after a congruence transformation C_X(·)=B⁻¹(·)B⁻ᵀ with B = L(X) (the Cholesky factor of X), as a symmetric matrix that admits an eigen‑decomposition. The authors write
C_X(grad_R f(X)) = Σ_{i=1}^n λ_i x_i x_iᵀ,
which translates back to the original tangent space as
grad_R f(X) = Σ_{i=1}^n λ_i H_i, H_i = B x_i x_iᵀ Bᵀ.
Each H_i is a rank‑one tangent direction and the set {H_i} is orthonormal with respect to the affine‑invariant metric on Pₙ. This decomposition is the foundation of the algorithm.

Instead of moving along the full gradient (which would require O(n³) operations due to dense matrix products), R1RSD selects a single rank‑one component H_k that promises the largest first‑order decrease. To identify the dominant eigenpair (λ₁, x₁), the algorithm employs the power method, requiring only matrix‑vector multiplications. By running O(log n) power iterations, the dominant eigenvector is approximated with sufficient accuracy, and the cost per iteration becomes O(n² log n). Because each update is rank‑one, the Cholesky factor B and its inverse can be updated in O(n) time using classic rank‑one update formulas, preserving the SPD property without recomputing a full factorization.

The step size α can often be obtained in closed form. For many NMEs of interest (e.g., algebraic Riccati equations), the objective function along the chosen direction reduces to a low‑degree polynomial or rational function, allowing exact line search at negligible extra cost. When a closed‑form step size is unavailable, a simple backtracking line search adds only a few extra matrix‑vector products.

Theoretical analysis proceeds under two standard assumptions. First, assuming Lipschitz‑continuous Riemannian gradient (smoothness), the authors prove that the algorithm reaches an ε‑optimal solution in O(n) iterations, yielding an overall complexity of O(n³ log n) in the worst case, but because each iteration costs O(n² log n), the total cost is O(n³ log n) which is dominated by the quadratic per‑iteration cost for large n. Second, when the objective is geodesically strongly convex, a linear convergence rate is established, reducing the iteration count to O(log(1/ε)). In both regimes, the per‑iteration cost remains O(n² log n), a dramatic improvement over traditional dense solvers that require O(n³) per iteration.

Empirical evaluation focuses on continuous‑time algebraic Riccati equations (CARE), discrete‑time algebraic Riccati equations (DARE), and other representative NMEs. Experiments scale the problem dimension from 1 000 up to 10 000. The authors compare R1RSD against MATLAB’s built‑in icare, structure‑preserving doubling algorithms, and recent subspace‑descent methods such as RRSD‑multi. Results show that for dimensions beyond roughly 5 000, the traditional solvers either run out of memory or exceed reasonable time limits, while R1RSD consistently converges within a few hundred iterations. Moreover, the algorithm attains high accuracy (residual norms below 10⁻⁸) and exhibits robust behavior even when the underlying NME is ill‑conditioned or slightly infeasible; the least‑squares formulation still returns a meaningful SPD approximation.

The paper also discusses practical implementation details. The authors provide MATLAB code and a public GitHub repository, including scripts to reproduce all reported experiments. They highlight that the rank‑one update framework naturally extends to time‑varying Riccati equations, enabling warm‑starts and rapid tracking of slowly changing solutions. In the context of the State‑Dependent Riccati Equation (SDRE) framework for nonlinear control, R1RSD could reduce the per‑time‑step computational burden by an order of magnitude, making real‑time embedded control of high‑dimensional systems feasible.

In summary, R1RSD offers a scalable, geometry‑aware algorithm for dense SPD‑constrained matrix equations. By exploiting the spectral structure of the transformed Riemannian gradient and restricting updates to a single dominant rank‑one direction, the method achieves O(n² log n) per‑iteration arithmetic, O(n) (or O(log n) under strong convexity) iteration complexity, and provable convergence guarantees. The experimental evidence confirms that the approach outperforms state‑of‑the‑art dense solvers on problems of size up to 10 000, opening the door to large‑scale applications in control, signal processing, and statistical estimation where traditional cubic‑cost methods are prohibitive.