low_rank_toolbox.krylov.solvers

Matrix Equation Solvers.

This submodule provides iterative solvers for large-scale matrix equations using Krylov subspace methods. These solvers are efficient for large sparse matrices where direct methods would be computationally prohibitive.

Functions

solve_sylvester : Solve the Sylvester equation AX + XB = C solve_lyapunov : Solve the Lyapunov equation AX + XA^T = C

Author: Benjamin Carrel, University of Geneva

low_rank_toolbox.krylov.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.solve_sylvester(A, B, C, tol=1e-12, max_iter=None, krylov_kwargs=None, is_A_symmetric=None, is_B_symmetric=None)[source]

Efficient low-rank compatible solver for the Sylvester equation.

Find X such that AX + XB = C.

This function is a wrapper that selects the appropriate solver based on the types of A, B, and C.

Parameters:

A (ndarray | spmatrix) – Matrix of shape (m, m)
B (ndarray | spmatrix) – Matrix of shape (n, n)
C (ndarray | LowRankMatrix) – Right-hand side of shape (m, 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
krylov_kwargs (dict, optional) – Krylov space configuration for sparse low-rank solver
is_A_symmetric (bool, optional) – If True, A is symmetric. Auto-detected if None.
is_B_symmetric (bool, optional) – If True, B is symmetric. Auto-detected if None.

Returns:

Solution X of shape (m, n), either dense or low-rank depending on inputs

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_sylvester
>>> # Small dense case
>>> A_small = np.array([[4, 1], [1, 3]])
>>> B_small = np.array([[2, 1], [1, 1]])
>>> C_small = np.array([[1, 0], [0, 1]])
>>> X_small = solve_sylvester(A_small, B_small, C_small)
>>> # Verify: AX + XB = C
>>> np.allclose(A_small @ X_small + X_small @ B_small, C_small)
True
>>> # Large sparse case with low-rank RHS
>>> from scipy.sparse import diags
>>> n, m = 100, 80
>>> A = diags([1, 4, 1], [-1, 0, 1], shape=(n, n), format='csr')
>>> B = diags([1, 2, 1], [-1, 0, 1], shape=(m, m), format='csr')
>>> U = np.zeros((n, 1))
>>> U[0, 0] = 1.0
>>> V = np.zeros((m, 1))
>>> V[0, 0] = 1.0
>>> C = LowRankMatrix(U, V.T)
>>> X = solve_sylvester(A, B, C, tol=1e-10)
>>> # Solution is low-rank
>>> type(X).__name__
'SVD'
>>> X.rank <= 50
True

Modules

`lyapunov_solvers`	Krylov-based solvers for the Lyapunov equation.
`sylvester_solvers`	Krylov-based solvers for the Sylvester equation.