low_rank_toolbox.krylov.utils.arnoldi
Arnoldi iteration algorithms for Krylov subspaces.
Author: Benjamin Carrel, University of Geneva, 2022-2023
Functions
|
Arnoldi algorithm. Computes orthogonal basis of a Krylov space: PK_m(A,x) = span{x, A x, A^2 x, ..., A^(m-1) x} where A is a non-symmetric matrix, and x is a vector. If A is symmetric, the Lanczos algorithm will be more efficient. |
|
Block Arnoldi algorithm. Compute orthogonal basis of a Krylov space: PK_m(A,X) = span{X, A X, A^2 X, ..., A^(m-1) X} where A is a (non-symmetric) matrix, and X is a (tall) matrix. |
|
Block Arnoldi algorithm with rational Krylov space. Compute orthogonal basis of a Krylov space: RK_m(A,X) = q_m(A) PK_m(A,X) q_m(A) = (A - p_1 I)^(-1) ... (A - p_m I)^(-1) where A is a matrix, and X is a (tall) matrix. p_i are the poles. |
|
Block Arnoldi algorithm with shift and invert. Compute orthogonal basis of a Krylov space: SK_m(A,X) = span{X, (A - sI)^(-1) X, ..., (A - sI)^(-m+1) X} where A is a matrix, and X is a (tall) matrix. s is the shift. |
|
Arnoldi algorithm with rational Krylov space. Computes orthogonal basis of a Krylov space: RK_m(A,x) = q_m(A) PK_m(A,x) q_m(A) = (A - p_1 I)^(-1) ... (A - p_m I)^(-1) where A is a matrix, and x is a vector. p_i are the poles. |
|
Arnoldi algorithm with shift and invert. Computes orthogonal basis of a Krylov space: SK_m(A,x) = span{x, (A - sI)^(-1) x, ..., (A - sI)^(-m+1) x} where A is a matrix, and x is a vector. s is the shift. |
- low_rank_toolbox.krylov.utils.arnoldi.Arnoldi(A, x, m)[source]
Arnoldi algorithm. Computes orthogonal basis of a Krylov space:
PK_m(A,x) = span{x, A x, A^2 x, …, A^(m-1) x}
where A is a non-symmetric matrix, and x is a vector. If A is symmetric, the Lanczos algorithm will be more efficient.
- 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) – Orthogonal Matrix of shape (n,m) containing the basis of the Krylov space
H (ndarray) – Hessenberg Matrix of shape (m,m) containing 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 Arnoldi >>> # Create a non-symmetric matrix >>> A = csr_matrix([[1, 2, 0], [0, 3, 1], [1, 0, 2]]) >>> x = np.array([1.0, 0.0, 0.0]) >>> Q, H = Arnoldi(A, x, m=3) >>> # Verify orthogonality >>> np.allclose(Q.T @ Q, np.eye(3)) True >>> # Verify Arnoldi relation: A @ Q = Q @ H (up to numerical error) >>> np.allclose(A @ Q, Q @ H[:3, :]) True
- low_rank_toolbox.krylov.utils.arnoldi.block_Arnoldi(A, X, m)[source]
Block Arnoldi algorithm. Compute orthogonal basis of a Krylov space:
PK_m(A,X) = span{X, A X, A^2 X, …, A^(m-1) X}
where A is a (non-symmetric) matrix, and X is a (tall) matrix.
See Y. Saad, Iterative methods for sparse linear systems, SIAM, 2003.
- Parameters:
A (ndarray) – Matrix of shape (n,n)
X (ndarray) – Matrix of shape (n,r), r > 1
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
H (ndarray) – Matrix of shape (m*r,m*r) containing the Hessenberg matrix
Examples
>>> import numpy as np >>> from scipy.sparse import csr_matrix >>> from low_rank_toolbox.krylov import block_Arnoldi >>> # Create a matrix and block starting vector >>> A = csr_matrix([[1, 2, 0, 0], [0, 3, 1, 0], [1, 0, 2, 1], [0, 1, 0, 3]]) >>> 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, H = block_Arnoldi(A, X, m=2) >>> Q.shape (4, 4) >>> # Verify orthogonality >>> np.allclose(Q.T @ Q, np.eye(4)) True
- low_rank_toolbox.krylov.utils.arnoldi.block_rational_Arnoldi(A, X, poles, inverse_only=False, inverses=None)[source]
Block Arnoldi algorithm with rational Krylov space. Compute orthogonal basis of a Krylov space:
RK_m(A,X) = q_m(A) PK_m(A,X) q_m(A) = (A - p_1 I)^(-1) … (A - p_m I)^(-1)
where A is a matrix, and X is a (tall) matrix. p_i are the poles.
- Parameters:
A (ndarray | spmatrix) – Matrix of shape (n,n)
X (ndarray) – Matrix of shape (n,r), r > 1
poles (list) – List of poles
inverse_only (bool) – If True, only invert the matrices. Default is False.
inverses (list) – List of functions that compute the matrix-vector product with the inverse of (A - pole*I). Optional, faster if provided.
- Return type:
- Returns:
Q (ndarray) – Matrix of shape (n,m*r) containing the basis of the Krylov space
H (ndarray) – Matrix of shape (m*r,m*r) containing the Hessenberg matrix
- low_rank_toolbox.krylov.utils.arnoldi.block_shift_and_invert_Arnoldi(A, X, m, shift=0, invA=None)[source]
Block Arnoldi algorithm with shift and invert. Compute orthogonal basis of a Krylov space:
SK_m(A,X) = span{X, (A - sI)^(-1) X, …, (A - sI)^(-m+1) X}
where A is a matrix, and X is a (tall) matrix. s is the shift.
- Parameters:
A (ndarray) – Matrix of shape (n,n)
X (ndarray) – Matrix of shape (n,r), r > 1
m (int) – Size of the Krylov space
shift (float) – Shift of the matrix. Default is 0 (only invert the matrix)
invA (callable) – Function that computes the matrix-vector product with the inverse of A - shift*I. Optional, faster if provided.
- Return type:
- Returns:
Q (ndarray) – Matrix of shape (n,m*r) containing the basis of the Krylov space
H (ndarray) – Matrix of shape (m*r,m*r) containing the Hessenberg matrix
- low_rank_toolbox.krylov.utils.arnoldi.rational_Arnoldi(A, x, poles, invert_only=False, inverses=None)[source]
Arnoldi algorithm with rational Krylov space. Computes orthogonal basis of a Krylov space:
RK_m(A,x) = q_m(A) PK_m(A,x) q_m(A) = (A - p_1 I)^(-1) … (A - p_m I)^(-1)
where A is a matrix, and x is a vector. p_i are the poles.
- Parameters:
- Return type:
- Returns:
Q (ndarray) – Orthogonal Matrix of shape (n, m) containing the basis of the Krylov space
H (ndarray) – Hessenberg Matrix of shape (m, m) containing 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 rational_Arnoldi >>> # Create a sparse matrix >>> A = csr_matrix([[4, 1, 0], [1, 3, 1], [0, 1, 2]]) >>> x = np.array([1.0, 0.0, 0.0]) >>> # Use poles to emphasize specific spectral regions >>> poles = [1.0, 2.0] >>> Q, H = rational_Arnoldi(A, x, poles) >>> Q.shape (3, 3) >>> # Verify orthogonality >>> np.allclose(Q.T @ Q, np.eye(3)) True
- low_rank_toolbox.krylov.utils.arnoldi.shift_and_invert_Arnoldi(A, x, m, shift=0)[source]
Arnoldi algorithm with shift and invert. Computes orthogonal basis of a Krylov space:
SK_m(A,x) = span{x, (A - sI)^(-1) x, …, (A - sI)^(-m+1) x}
where A is a matrix, and x is a vector. s is the shift.
- Parameters:
- Return type:
- Returns:
Q (ndarray) – Orthogonal Matrix of shape (n, m) containing the basis of the Krylov space
H (ndarray) – Hessenberg Matrix of shape (m, m) containing 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 shift_and_invert_Arnoldi >>> # Create a sparse matrix >>> A = csr_matrix([[4, 1, 0], [1, 3, 1], [0, 1, 2]]) >>> x = np.array([1.0, 0.0, 0.0]) >>> # Use shift-and-invert to find eigenvectors near shift=3 >>> Q, H = shift_and_invert_Arnoldi(A, x, m=2, shift=3.0) >>> # The Krylov space emphasizes eigenvalues close to 3 >>> Q.shape (3, 2)