Logo
ScienceStack
Table of Contents
Fetching ...
Sign In

A subspace-conjugate gradient method for linear matrix equations

Davide Palitta, Martina Iannacito, Valeria Simoncini

TL;DR

This work addresses solving large-scale multiterm linear matrix equations $\mathcal{L}(X)=C$ with symmetric positive definite $\mathcal{L}$ by introducing a preconditioned subspace-conjugate gradient method (ss-cg). Unlike standard matrix-oriented CG that updates along single-vector directions, ss-cg leverages wide subspaces generated by low-rank factor information and enforces $\mathcal{L}$-orthogonality over these subspaces, aided by a randomized range finder to control memory. The method extends to multiterm Sylvester equations and integrates memory-saving truncations, inexact coefficient solves, and preconditioning (including LR-ADI-based approaches) to handle very large problems. Numerical experiments on Lyapunov and Sylvester problems demonstrate strong convergence, reduced memory usage, and competitiveness against existing methods across diverse applications. The approach provides a scalable framework for large-scale matrix equations in control, stochastic PDEs, and related areas.

Abstract

The efficient solution of large-scale multiterm linear matrix equations is a challenging task in numerical linear algebra, and it is a largely open problem. We propose a new iterative scheme for symmetric and positive definite operators, significantly advancing methods such as truncated matrix-oriented Conjugate Gradients (CG). The new algorithm capitalizes on the low-rank matrix format of its iterates by fully exploiting the subspace information of the factors as iterations proceed. The approach implicitly relies on orthogonality conditions imposed over much larger subspaces than in CG, unveiling insightful connections with subspace projection methods. The new method is also equipped with memory-saving strategies. In particular, we show that for a given matrix $\mathbf{Y}$, the action $\mathcal{L}(\mathbf{Y})$ in low rank format may not be evaluated exactly due to memory constraints. This problem is often underestimated, though it will eventually produce Out-of-Memory breakdowns for a sufficiently large number of terms. We propose an ad-hoc randomized range-finding strategy that appears to fully resolve this shortcoming. Experimental results with typical application problems illustrate the potential of our approach over various methods developed in the recent literature.

A subspace-conjugate gradient method for linear matrix equations

TL;DR

This work addresses solving large-scale multiterm linear matrix equations with symmetric positive definite by introducing a preconditioned subspace-conjugate gradient method (ss-cg). Unlike standard matrix-oriented CG that updates along single-vector directions, ss-cg leverages wide subspaces generated by low-rank factor information and enforces -orthogonality over these subspaces, aided by a randomized range finder to control memory. The method extends to multiterm Sylvester equations and integrates memory-saving truncations, inexact coefficient solves, and preconditioning (including LR-ADI-based approaches) to handle very large problems. Numerical experiments on Lyapunov and Sylvester problems demonstrate strong convergence, reduced memory usage, and competitiveness against existing methods across diverse applications. The approach provides a scalable framework for large-scale matrix equations in control, stochastic PDEs, and related areas.

Abstract

The efficient solution of large-scale multiterm linear matrix equations is a challenging task in numerical linear algebra, and it is a largely open problem. We propose a new iterative scheme for symmetric and positive definite operators, significantly advancing methods such as truncated matrix-oriented Conjugate Gradients (CG). The new algorithm capitalizes on the low-rank matrix format of its iterates by fully exploiting the subspace information of the factors as iterations proceed. The approach implicitly relies on orthogonality conditions imposed over much larger subspaces than in CG, unveiling insightful connections with subspace projection methods. The new method is also equipped with memory-saving strategies. In particular, we show that for a given matrix , the action in low rank format may not be evaluated exactly due to memory constraints. This problem is often underestimated, though it will eventually produce Out-of-Memory breakdowns for a sufficiently large number of terms. We propose an ad-hoc randomized range-finding strategy that appears to fully resolve this shortcoming. Experimental results with typical application problems illustrate the potential of our approach over various methods developed in the recent literature.
Paper Structure (16 sections, 7 theorems, 57 equations, 1 figure, 6 tables, 3 algorithms)

This paper contains 16 sections, 7 theorems, 57 equations, 1 figure, 6 tables, 3 algorithms.

Table of Contents

  1. Introduction
  2. Notation
  3. Truncated matrix-oriented CG
  4. The s s-- cg method for the multiterm Lyapunov equation
  5. A first version of the algorithm for the multiterm Lyapunov equation
  6. Discussion on the developed procedure
  7. The iteration for the multiterm Sylvester equation
  8. Advanced implementation devices
  9. Low-rank truncations
  10. Inaccurate coefficients
  11. Preconditioning
  12. The complete algorithm
  13. Numerical results
  14. Numerical experiments for the multiterm Lyapunov equation
  15. Numerical experiments for the multiterm Sylvester equation
  16. ...and 1 more sections

Key Result

Proposition 3.1

\newlabelprop:alpha0 Assume that $\bm A_i$ and $\bm B_i$ are symmetric, and that ${\@fontswitch{}{\mathcal{}} L}$ is positive definite. The minimizer $\mathbf{\alpha} _k\in\mathbb{R}^{s_k\times s_k}$ of eq:min:alpha is the unique solution of or, equivalently, of $P_k^T {\@fontswitch{}{\mathcal{}} L}(P_k \mathbf{\alpha} P_k^T) P_k = P_k^T {\bm{{R}}}_k P_k$.

Figures (1)

  • Figure 1: Example \ref{['ex:diffreac']}. Left: Singular value distribution of the approximate solution matrix $\bm X$ for $\gamma_0(z)=\sin(z\pi)$ and $\gamma_0(z)=\exp(z\pi)$. Right: Convergence history of s s-- cg for $\gamma_0(z)=\exp(z\pi)$, different stopping tolerances tol, and fixed maxrank=20. \newlabelfig:svdX0

Theorems & Definitions (23)

  • Proposition 3.1
  • Proof 1
  • Remark 3.2
  • Remark 3.3
  • Proposition 3.4
  • Proof 2
  • Remark 4.1
  • Proposition 4.2
  • Proof 3
  • Proposition 4.3
  • ...and 13 more