A Riemannian approach to low-rank algebraic Riccati equations - implementation -

Somehow related to neutron transport theory and a different kind of matrix factorization, here is: A Riemannian approach to low-rank algebraic Riccati equations by Bamdev Mishra, B. Vandereycken

We propose a Riemannian optimization approach for computing low-rank solutions of the algebraic Riccati equation. The scheme alternates between fixed-rank optimization and rank-one updates. The fixed-rank optimization is on the set of fixed-rank symmetric positive definite matrices which is endowed with a particular Riemannian metric (and geometry) that is tuned to the structure of the objective function. We specifically discuss the implementation of a Riemannian trust-region algorithm that is potentially scalable to large-scale problems. The rank-one update is based on a descent direction that ensures a monotonic decrease of the cost function. Preliminary numerical results on standard small-scale benchmarks show that we obtain solutions to the Riccati equation at lower ranks than the standard approaches.

from the code page of Bamdev Mishra

We solve the equation $A^TX + XA + X BB^T X = C^TC$ for X on the low rank manifold of symmetric and positive definite matrices, where A is large and sparse, B is a tall matrix, and C is a fat matrix. For $B = 0$, this boils down to the Lyapunov equation.

The code has been tuned for large scale problems. Computationally, the numerical bottleneck involes solving the linear system $Px = b$, where P is a sparsely-structured matrix.

In the present code, $Px = b$ for smaller to moderate size problems (about 10 000 x 10 000) is solved with the backslash operator of Matlab without any tuning. For larger size matrices, we employ the Matlab's built-in iterative pcg solver without any tuning, which may not be efficient for all problems. Specific problems call for specific solvers! If you are interestd in a particular setup and have a solver, then do look for the function "use_iterative_solver" in the file "Main_files/Riemannian_lowrank_riccati.m" and modify it accordingly.

The code is built on the Manopt platform. Feel free to contact me on any issue.

It will be added to the Advanced Matrix Factorization Jungle page.

Join the CompressiveSensing subreddit or the Google+ Community and post there !
Liked this entry ? subscribe to Nuit Blanche's feed, there's more where that came from. You can also subscribe to Nuit Blanche by Email, explore the Big Picture in Compressive Sensing or the Matrix Factorization Jungle and join the conversations on compressive sensing, advanced matrix factorization and calibration issues on Linkedin.