cho_factor

detkit.cho_factor(A, m=None, lower=True, overwrite=False)

Cholesky decomposition of a sub-matrix.

Parameters:

Anumpy.ndarray

A symmetric positive-definite matrix with either float64 or float32 data type stored with either row-major or column major ordering.

mint, default=None

An integer determining the shape of an upper-left sub-matrix the A[:m, :m] to be considered. If None, the full matrix A[:, :] is used.

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 A will be overwritten as the output, allowing to save memory. In this case, the output matrix lu will point to the same location in the memory as the input matrix A.

Note

When overwrite is set to True, the matrix A should have column-major (Fortran) ordering.

Returns:

lunumpy.ndarray

A 2D matrix of the same shape as the input matrix A (and not the shape of the sub-matrix). The upper-left sub-matrix of lu contains the Cholesky factor of the sub-matrix of A. Namely, the upper or lower triangular parts of the sub-matrix of lu consist of the matrix L (if lower=true) or U (if lower=False) matrix. If overwrite is set to True, the output matrix lu is becomes a view for the matrix A.

See also

detkit.ldl_factor

detkit.lu_factor

detkit.cho_solve

Notes

Cholesky decomposition of sub-matrix:

Let \(\mathbf{A}_{[:p,:p]}\) denote the sub-matrix of the size \(p \times p\) to be the upper-left corner of matrix \(\mathbf{A}\). Its Cholesky decomposition is given by:

\[\mathbf{A}_{[:p, :p]} = \mathbf{L} \mathbf{L}^{\intercal},\]

or

\[\mathbf{A}_{[:p, :p]} = \mathbf{U}^{\intercal} \mathbf{U},\]

where \(\mathbf{L}\) is a lower-triangular matrix of the size \(p \times p\) and \(\mathbf{U}\) is an upper-triangular matrix of the size \(p \times p\).

Interpreting L or U factors from the output cho:

Note that the output matrix cho is given as the same shape as A, and only its upper-left sub-matrix should be used to extract the matrices \(\mathbf{L}\) and \(\mathbf{U}\) as follows.

Let p, p = shape be the given shape of the upper-left sub-matrix, which may be smaller than the shape of A. Then, L = cho[:p, :p] and U = cho[:p, :p].

Comparison with scipy.linalg.cho_factor:

To compute the Cholesky factorization of a sub-matrix of the input matrix using scipy.linalg.cho_factor 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_factor together with the shape 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 spotrf (for 32-bit precision) and dpotrf (for 64-bit precision).

This function is internally used for detkit.memdet() for efficient computation of matrix determinant under memory constraint.

References

• LAPACK

Examples

>>> from detkit import cho_factor, Memory
>>> import numpy

>>> # Create a symmetrc positive-definite matrix with 32-bit precision
>>> # and column-major ordering
>>> n = 1000
>>> A = numpy.random.randn(n, n) + 2*n * numpy.eye(n, n)
>>> A = A.T @ A
>>> A = numpy.asfortranarray(A)
>>> A = A.astype(numpy.float32)

>>> # Get a copy of A (for the purpose of comparison) since we will
>>> # overwrite A
>>> A_copy = numpy.copy(A)

>>> # Track memory allocation to check if Chplesky decomposition is not
>>> # creating any new memory.
>>> mem = Memory()
>>> mem.set()

>>> # Choleksy factorization of the upper-left sub-matrix with a
>>> # smaller shape
>>> m = 800
>>> cho = cho_factor(A, m=m, lower=True, overwrite=True)

>>> # Check peak memory allocation (compared to memory of a sum-matrix)
>>> slice_nbytes = m**2 * A.dtype.itemsize
>>> print(mem.peak() / slice_nbytes)
0.004

>>> # Extract L factor from cho matrix
>>> L = numpy.tril(cho[:m, :m])

>>> # Check if A_copy = L@L.T holds.
>>> atol = numpy.finfo(A.dtype).resolution
>>> print(numpy.allclose(A_copy[:m, :m], L @ L.T, atol=10*atol))
True

>>> # When overwrite is set to True, check if cho is indeed a view of A
>>> print(cho.base == A.base)
True

In the above example, the object mem of class detkit.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.