low_rank_toolbox.krylov.utils

Core Krylov Iteration Algorithms.

This submodule implements the fundamental iteration algorithms for constructing Krylov subspaces, including Arnoldi and Lanczos methods and their block and rational variants.

Functions

Arnoldi : Standard Arnoldi iteration rational_Arnoldi : Rational Arnoldi iteration block_Arnoldi : Block Arnoldi iteration block_rational_Arnoldi : Block rational Arnoldi iteration shift_and_invert_Arnoldi : Shift-and-invert Arnoldi block_shift_and_invert_Arnoldi : Block shift-and-invert Arnoldi Lanczos : Lanczos iteration for symmetric matrices block_Lanczos : Block Lanczos iteration

low_rank_toolbox.krylov.utils.Arnoldi(A, x, m)

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:

tuple[ndarray, ndarray]

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.Lanczos(A, x, m)

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:

tuple[ndarray, spmatrix]

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.block_Arnoldi(A, X, m)

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:

tuple[ndarray, ndarray]

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.block_Lanczos(A, X, m)

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:

tuple[ndarray, spmatrix]

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)

low_rank_toolbox.krylov.utils.block_rational_Arnoldi(A, X, poles, inverse_only=False, inverses=None)

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:

tuple[ndarray, ndarray]

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.block_shift_and_invert_Arnoldi(A, X, m, shift=0, invA=None)

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:

tuple[ndarray, ndarray]

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.rational_Arnoldi(A, x, poles, invert_only=False, inverses=None)

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:

• A (ndarray | spmatrix) – Matrix of shape (n,n)

• x (ndarray) – Vector of shape (n,)

• poles (list) – List of poles

• invert_only (bool) – If True, only invert the matrices. Default is False.

• inverses (list) – List of inverses of the matrices (optional). Default is None.

Return type:

tuple[ndarray, ndarray]

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.shift_and_invert_Arnoldi(A, x, m, shift=0)

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:

• A (ndarray | spmatrix) – Matrix of shape (n,n)

• x (ndarray) – Vector of shape (n,)

• m (int) – Size of the Krylov space

• shift (float) – Shift of the matrix. Default is 0 (only invert the matrix)

Return type:

tuple[ndarray, ndarray]

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)

Modules

arnoldi	Arnoldi iteration algorithms for Krylov subspaces.
lanczos	Lanczos iteration algorithms for symmetric matrices.