low_rank_toolbox.krylov.solvers.lyapunov_solvers
Krylov-based solvers for the Lyapunov equation.
Author: Benjamin Carrel, University of Geneva, 2022-2023
Functions
|
Efficient low-rank compatible solver for the Lyapunov equation. |
|
Solve the Lyapunov equation AX + XA = C for small matrices. |
Low-rank solver for the general Lyapunov equation: |
|
Low-rank solver for the symmetric Lyapunov equation: |
- low_rank_toolbox.krylov.solvers.lyapunov_solvers.solve_lyapunov(A, C, tol=1e-12, max_iter=None, is_symmetric=None, krylov_kwargs=None)[source]
Efficient low-rank compatible solver for the Lyapunov equation.
Find X such that AX + XA^H = C.
This function is a wrapper that selects the appropriate solver based on the types of A and C.
For lower computational cost, if A and C are symmetric, set is_symmetric=True. This enables the use of the Lanczos algorithm (3-term recurrence) instead of Arnoldi (full orthogonalization), approximately halving the computational cost and memory usage per iteration.
- Parameters:
A (ndarray | spmatrix) – The matrix A of shape (n, n)
C (ndarray | LowRankMatrix) – The matrix C of shape (n, n)
tol (float, optional) – Convergence tolerance for Krylov solver, by default 1e-12
max_iter (int, optional) – Maximum iterations for Krylov solver, by default None
is_symmetric (bool, optional) – If True, use symmetric solver (only for sparse A and low-rank C). Exploits symmetry using Lanczos algorithm for ~2x speedup. If None, symmetry is auto-detected (with efficiency warning).
krylov_kwargs (dict, optional) – Krylov space configuration (see solve_sparse_low_rank_symmetric_lyapunov or solve_sparse_low_rank_non_symmetric_lyapunov for details)
- Returns:
The solution X, either dense or low-rank and symmetric if and only if C is symmetric.
- Return type:
ndarray | LowRankMatrix
Examples
>>> import numpy as np >>> from scipy.sparse import csr_matrix >>> from low_rank_toolbox import LowRankMatrix >>> from low_rank_toolbox.krylov import solve_lyapunov >>> # Small dense case >>> A_small = np.array([[4, 1], [1, 3]]) >>> C_small = np.array([[1, 0], [0, 1]]) >>> X_small = solve_lyapunov(A_small, C_small) >>> # Verify: AX + XA = C >>> np.allclose(A_small @ X_small + X_small @ A_small, C_small) True >>> # Large sparse case with low-rank RHS >>> from scipy.sparse import diags >>> n = 100 >>> A = diags([1, 4, 1], [-1, 0, 1], shape=(n, n), format='csr') >>> U = np.zeros((n, 1)) >>> U[0, 0] = 1.0 >>> C = LowRankMatrix(U, U.T) >>> X = solve_lyapunov(A, C, is_symmetric=True, tol=1e-10) >>> # Solution is low-rank and symmetric >>> type(X).__name__ 'SVD' >>> X.rank <= n # Much lower rank than n True
- low_rank_toolbox.krylov.solvers.lyapunov_solvers.solve_small_lyapunov(A, C)[source]
Solve the Lyapunov equation AX + XA = C for small matrices.
- low_rank_toolbox.krylov.solvers.lyapunov_solvers.solve_sparse_low_rank_non_symmetric_lyapunov(A, C, tol=1e-12, max_iter=None, krylov_kwargs=None)[source]
Low-rank solver for the general Lyapunov equation:
Find X such that A X + X A^H = C.
Uses two Krylov spaces (left for A, right for A^H). The projected problem is solved as a Sylvester equation: Ak Y + Y Bk^H = Ck.
- Parameters:
A (spmatrix) – Sparse matrix of shape (n, n)
C (LowRankMatrix) – Low-rank right-hand side of shape (n, n)
tol (float, optional) – Convergence tolerance, by default 1e-12
max_iter (int, optional) – Maximum iterations, by default n / rank(C)
krylov_kwargs (dict, optional) – Krylov space configuration with keys: - ‘extended’ (bool): Use extended Krylov space (default True) - ‘invA’ (callable): Custom inverse function for A - ‘invAH’ (callable): Custom inverse function for A^H - ‘poles_A’ (array_like): Poles for rational Krylov space for A - ‘poles_AH’ (array_like): Poles for rational Krylov space for A^H
- Returns:
Low-rank solution X (not necessarily symmetric)
- Return type:
- low_rank_toolbox.krylov.solvers.lyapunov_solvers.solve_sparse_low_rank_symmetric_lyapunov(A, C, tol=1e-12, max_iter=None, krylov_kwargs=None)[source]
Low-rank solver for the symmetric Lyapunov equation:
Find X such that A X + X A = C.
NOTE: Matrices A and C must be symmetric, which is exploited to halve the computational cost. Uses a single Krylov space since the solution X is also symmetric.
Performance: Uses the Lanczos algorithm (3-term recurrence) instead of Arnoldi (full orthogonalization), reducing computational cost by approximately 2x and storage by O(m²) → O(m) for the recurrence coefficients.
- Parameters:
A (spmatrix) – Symmetric sparse matrix of shape (n, n)
C (LowRankMatrix) – Symmetric low-rank matrix of shape (n, n)
tol (float, optional) – Convergence tolerance, by default 1e-12
max_iter (int, optional) – Maximum iterations, by default n / rank(C)
krylov_kwargs (dict, optional) –
Krylov space configuration with keys: - ‘extended’ (bool): Use extended Krylov space (default True)
Extended Krylov uses both A and A^(-1), providing faster convergence. Uses Lanczos algorithm for symmetric A (efficient 3-term recurrence).
’invA’ (callable): Custom inverse function for A
’poles’ (array_like): Poles for rational Krylov space WARNING: Rational Krylov currently uses Arnoldi even for symmetric A (Lanczos not yet implemented for rational case).
- Returns:
Symmetric low-rank solution X
- Return type:
Examples
>>> import numpy as np >>> from scipy.sparse import diags >>> from low_rank_toolbox import LowRankMatrix >>> from low_rank_toolbox.krylov import solve_sparse_low_rank_symmetric_lyapunov >>> # Create a symmetric sparse matrix (100x100 tridiagonal) >>> n = 100 >>> A = diags([1, 4, 1], [-1, 0, 1], shape=(n, n), format='csr') >>> # Create a symmetric low-rank right-hand side >>> U = np.zeros((n, 1)) >>> U[0, 0] = 1.0 >>> C = LowRankMatrix(U, U.T) # Rank-1 symmetric matrix >>> # Solve AX + XA = C >>> X = solve_sparse_low_rank_symmetric_lyapunov(A, C, tol=1e-10) >>> # Verify the solution >>> residual = A @ X + X @ A - C >>> residual.norm() < 1e-8 True >>> # X is low-rank >>> X.rank <= n True