low_rank_toolbox.krylov.solvers.sylvester_solvers

Krylov-based solvers for the Sylvester equation.

Author: Benjamin Carrel, University of Geneva, 2022-2023

Functions

`solve_small_sylvester`(A, B, C)	Solve the Sylvester equation for small systems.
`solve_sparse_low_rank_sylvester`(A, B, C[, ...])	Low-rank solver for the Sylvester equation.
`solve_sylvester`(A, B, C[, tol, max_iter, ...])	Efficient low-rank compatible solver for the Sylvester equation.
`solve_sylvester_large_A_small_B`(A, B, C)	Solve the Sylvester equation when A is large and B is small.

low_rank_toolbox.krylov.solvers.sylvester_solvers.solve_small_sylvester(A, B, C)[source]

Solve the Sylvester equation for small systems.

Find X such that AX + XB = C. Wrapper to scipy.linalg.solve_sylvester. For larger systems, use solve_sparse_low_rank_sylvester or solve_sylvester_large_A_small_B.

Parameters:

A (ndarray | spmatrix) – Matrix of shape (m, m)
B (ndarray | spmatrix) – Matrix of shape (n, n)
C (ndarray) – Matrix of shape (m, n)

Returns:

Dense solution of shape (m, n)

Return type:

ndarray

low_rank_toolbox.krylov.solvers.sylvester_solvers.solve_sparse_low_rank_sylvester(A, B, C, tol=1e-12, max_iter=None, krylov_kwargs=None, is_A_symmetric=None, is_B_symmetric=None)[source]

Low-rank solver for the Sylvester equation.

Find X such that AX + XB = C.

Uses two Krylov spaces (left for A, right for B). The projected problem is solved as a small Sylvester equation: Ak Y + Y Bk = Ck.

Algorithm Note: Each Krylov space independently uses Lanczos if the corresponding matrix is symmetric. If A is symmetric, the left Krylov space uses Lanczos. If B is symmetric, the right Krylov space uses Lanczos. Both can use Lanczos independently.

Special Case: When A = B and both are symmetric, the problem becomes a Lyapunov equation (AX + XA = C). Use solve_lyapunov() instead to exploit the symmetry of X and use only one Krylov space instead of two.

Reference: Simoncini, V., 2016. Computational Methods for Linear Matrix Equations. SIAM Rev. 58, 377-441. https://doi.org/10.1137/130912839

Parameters:

A (spmatrix) – Sparse matrix of shape (m, m)
B (spmatrix) – Sparse matrix of shape (n, n)
C (LowRankMatrix) – Low-rank right-hand side of shape (m, n)
tol (float, optional) – Convergence tolerance, by default 1e-12
max_iter (int, optional) – Maximum iterations, by default m / 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 - ‘invB’ (callable): Custom inverse function for B - ‘poles_A’ (array_like): Poles for rational Krylov space for A - ‘poles_B’ (array_like): Poles for rational Krylov space for B
is_A_symmetric (bool, optional) – If True, A is symmetric and Lanczos is used for left Krylov space. If None (default), symmetry is auto-detected.
is_B_symmetric (bool, optional) – If True, B is symmetric and Lanczos is used for right Krylov space. If None (default), symmetry is auto-detected.

Returns:

Low-rank solution X of shape (m, n)

Return type:

QuasiSVD

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_sylvester
>>> # Create sparse matrices A and B (larger sizes)
>>> 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')
>>> # Create a low-rank right-hand side
>>> 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)
>>> # Solve AX + XB = C
>>> X = solve_sparse_low_rank_sylvester(A, B, C, tol=1e-10)
>>> # Verify the solution
>>> residual = A @ X + X @ B - C
>>> residual.norm() < 1e-8
True
>>> # X is low-rank
>>> X.rank <= 50
True

low_rank_toolbox.krylov.solvers.sylvester_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

low_rank_toolbox.krylov.solvers.sylvester_solvers.solve_sylvester_large_A_small_B(A, B, C)[source]

Solve the Sylvester equation when A is large and B is small.

Find X such that AX + XB = C using eigenvalue decomposition of B.

Reference: Simoncini, V., 2016. Computational Methods for Linear Matrix Equations. SIAM Rev. 58, 377-441. https://doi.org/10.1137/130912839

Parameters:

A (spmatrix) – Sparse matrix of shape (m, m)
B (ndarray) – Dense matrix of shape (n, n)
C (ndarray) – Dense matrix of shape (m, n)

Returns:

Dense solution of the Sylvester equation of shape (m, n)

Return type:

ndarray