cho_solve#
- detkit.cho_solve(cho, B, shape=None, lower=True, overwrite=False)#
Solve symmetric positive-definite linear system given a sub-matrix output of Cholesky decomposition.
- Parameters:
- chonumpy.ndarray
The Cholesky decomposing of a symmetric positive-definite matrix with either
float64
orfloat32
data type stored column-major ordering. This matrix can be obtained fromdetkit.cho_factor()
.- Bnumpy.ndarray
Matrix of the right-had side. A rectangular matrix with either
float64
orfloat32
data type stored with row-major or column-major ordering.- shapetuple, default=None
A tuple of size two, determining the shape of an upper-left sub-matrices of cho and B to be referenced. Namely, if
shape
is given as the tuple(m, n)
, the sub-matricescho[:m, :m]
andB[:m, :n]
are used. If None, the full shape ofcho[:, :]
andB[:, :]
are considered.- lowerbool, default=True
If True, cho is assumed to be a lower-triangular Cholesky factor \(\mathbf{L}\) such that \(\mathbf{A} = \mathbf{L} \mathbf{L}^{\intercal}\). If False, cho is assumed to be an upper-triangular Cholesky factor \(\mathbf{U}\) such that \(\mathbf{A} = \mathbf{U}^{\intercal} \mathbf{U}\).
- overwritebool, default=False
If True, the input matrix B will be overwritten as the output, allowing to save memory. In this case, the output matrix X will point to the same location in the memory as the input matrix B.
Note
When
overwrite
is set to True, the matrix B should have column-major (Fortran) ordering.
- Returns:
- Xnumpy.ndarray
A 2D matrix of the same shape as the input matrix B (and not the shape of the sub-matrix). The upper-left sub-matrix of B contains the solution to the linear system of equations corresponding to the sub-matrices determined by the
shape
argument. Ifoverwrite
is set to True, the output matrix X is becomes a view for the matrix B.
See also
Notes
Linear system of equations for sub-matrix:
Let \(\mathbf{B}_{[:m,:n]}\) denote the sub-matrix of the size \(m \times n\) to be the upper-left corner of matrix \(\mathbf{B}\). Given matrices \(\mathbf{A}\) and \(\mathbf{B}\), this function solves
\[\mathbf{A}_{[:m, :m]} \mathbf{X}_{[:m, :n]} = \mathbf{B}_{[:m, :n]}.\]Cholesky Factor Input Argument:
The input to this function is the Cholesky factor of \(\mathbf{A}\) ( and not \(\mathbf{A}\) itself). Namely, if
lower=True
, it is assumed that the input argumentcho
is the matrix \(\mathbf{L}\) where\[\mathbf{A}_{[:m, :m]} = \mathbf{L}_{[:m, :m]} \mathbf{L}_{[:m, :m]}^{\intercal}.\]Similarly, if
lower=False
, it is assumed that the input argumentcho
is the matrix \(\mathbf{U}\) where\[\mathbf{A}_{:m, :m} = \mathbf{U}_{[:m, :m]}^{\intercal} \mathbf{U}_{[:m, :m]}.\]The shape of
cho
can be larger than the shape of the \(m \times m\) sub-matrix, however, its upper-left \(m \times m\) sub-matrix will be referenced during the computation.Interpreting the output matrix:
The shape of the output variable X is the same as the shape of B, even if a smaller sub-matrix is considered. Regardless, only the corresponding upper-left sub-matrix of X has meaningful data. Namely, if
shape=(m, n)
, the sliceX[:m, :n]
should be considered as the referenced. As such, the following relation should hold:A[:m, :m] @ X[:m, :n] = B[:m, :n]
Comparison with scipy.linalg.cho_solve:
To solve a linear system for a sub-matrix of the input matrices using
scipy.linalg.cho_solve
function, you should pass a slice of the matrix to the function. This approach is not memory-efficient since the sliced array allocates new memory.In contrast, using
detkit.cho_solve
together with theshape
argument, no memory slice is created during the inner computation, rather, the data from the original input matrix is accessed efficiently.Implementation:
This function is a wrapper around LAPACK’s
spotrs
(for 32-bit precision) anddpotrs
(for 64-bit precision).This function is internally used for
detkit.memdet()
for efficient computation of matrix determinant under memory constraint.References
Examples
>>> from detkit import cho_factor, cho_solve, Memory >>> import numpy >>> # Create a symmetric positive-definite matrix with 32-bit precision >>> # and column-major ordering >>> A = numpy.random.randn(1000, 1000) >>> A = A.T @ A + 1000 * numpy.eye(1000, 1000) >>> A = numpy.asfortranarray(A) >>> A = A.astype(numpy.float32) >>> # Matrix of the right-hand side >>> B = numpy.random.randn(1000, 2000) >>> B = numpy.asfortranarray(B) >>> B = B.astype(numpy.float32) >>> # Get a copy of A and B (for the purpose of comparison) since we >>> # will overwrite A and B >>> A_copy = numpy.copy(A) >>> B_copy = numpy.copy(B) >>> # Track memory allocation to check if either of the Cholesky >>> # decomposition or solving linear system is now creating any new >>> # memory. >>> mem = Memory() >>> mem.set() >>> # Cholesky factorization of the upper-left sub-matrix >>> p, q = 800, 700 >>> cho = cho_factor(A, m=p, lower=True, overwrite=True) >>> # Solve a linear system >>> X = cho_solve(cho, B, shape=(p, q), lower=True, overwrite=True) >>> # Check peak memory allocation (compared to memory of a sum-matrix) >>> slice_nbytes = p * q * B.dtype.itemsize >>> print(mem.peak() / slice_nbytes) 0.001 >>> # Check if A_copy @ X = B_copy holds. >>> atol = numpy.finfo(A.dtype).resolution >>> print(numpy.allclose(A_copy[:p, :p] @ X[:p, :q], B_copy[:p, :q], ... atol=10*atol)) True >>> # When overwrite is set to True, check if X is indeed a view of B >>> print(X.base == B.base) True
In the above example, the object
mem
of classdetkit.Memory
tracks memory allocation. The peak of allocated memory during the matrix multiplication is three orders of magnitude smaller than the size of one of the matrices slices, confirming that no new array slice was created during the operation.