low_rank_toolbox.krylov.utils.lanczos
Lanczos iteration algorithms for symmetric matrices.
Author: Benjamin Carrel, University of Geneva, 2022-2023
Functions
|
Lanczos algorithm. Computes orthogonal basis of a Krylov space: K_m(A,x) = span{x, A x, A^2 x, ..., A^(m-1) x} where $A$ is a symmetric matrix, and $x$ is a vector. If $A$ is non-symmetric, use the Arnoldi algorithm instead. |
|
Block Lanczos algorithm. |
- low_rank_toolbox.krylov.utils.lanczos.Lanczos(A, x, m)[source]
Lanczos algorithm. Computes orthogonal basis of a Krylov space:
K_m(A,x) = span{x, A x, A^2 x, …, A^(m-1) x}
where $A$ is a symmetric matrix, and $x$ is a vector. If $A$ is non-symmetric, use the Arnoldi algorithm instead.
Inspired from Martin J. Gander’s lecture.
- Parameters:
A (ndarray | spmatrix) – Matrix of shape (n,n)
x (ndarray) – Vector of shape (n,)
m (int) – Size of the Krylov space
- Return type:
- Returns:
Q (ndarray) – Matrix of shape (n,m) containing the basis of the Krylov space.
T (spmatrix) – Tridiagonal matrix of shape (m,m). It is also the projection of A on the Krylov space.
Examples
>>> import numpy as np >>> from scipy.sparse import csr_matrix >>> from low_rank_toolbox.krylov import Lanczos >>> # Create a symmetric matrix >>> A = csr_matrix([[4, 1, 0], [1, 3, 1], [0, 1, 2]]) >>> x = np.array([1.0, 0.0, 0.0]) >>> Q, T = Lanczos(A, x, m=3) >>> # Verify orthogonality >>> np.allclose(Q.T @ Q, np.eye(3)) True >>> # Verify Lanczos relation: A @ Q = Q @ T (up to numerical error) >>> np.allclose(A @ Q, Q @ T.toarray()) True >>> # T is tridiagonal (more efficient than Arnoldi's Hessenberg) >>> T.toarray() array([[...]])
- low_rank_toolbox.krylov.utils.lanczos.block_Lanczos(A, X, m)[source]
Block Lanczos algorithm. Initialize a Krylov Space where X is a matrix
- Parameters:
A (ndarray) – Matrix of shape (n,n)
X (ndarray) – Matrix of shape (n,r)
m (int) – Size of the Krylov space
- Return type:
- Returns:
Q (ndarray) – Matrix of shape (n,m*r) containing the basis of the Krylov space
T (spmatrix) – Tridiagonal matrix of shape (m*r,m*r). It is also the projection of A on the Krylov space.
Examples
>>> import numpy as np >>> from scipy.sparse import csr_matrix >>> from low_rank_toolbox.krylov import block_Lanczos >>> # Create a symmetric matrix >>> A = csr_matrix([[4, 1, 0, 0], [1, 3, 1, 0], [0, 1, 2, 1], [0, 0, 1, 1]]) >>> X = np.array([[1.0, 0.0], [0.0, 1.0], [0.0, 0.0], [0.0, 0.0]]) >>> # Build block Krylov space with m=2 blocks >>> Q, T = block_Lanczos(A, X, m=2) >>> Q.shape (4, 4) >>> # Verify orthogonality >>> np.allclose(Q.T @ Q, np.eye(4)) True >>> # T is block tridiagonal >>> T.shape (4, 4)