ldl_solve

detkit.ldl_solve(ldu, piv, B, shape=None, lower=True, overwrite=False)

Solve symmetric positive-definite linear system given a sub-matrix of LDL decomposition.

Parameters:

ldunumpy.ndarray

The LDL decomposing of a symmetric matrix with either float64 or float32 data type stored column-major ordering. This matrix can be obtained from detkit.ldl_factor().

pivnumpy.array

Row pivoting indices (zero-based indexing). This array can be obtained from detkit.ldl_factor().

Bnumpy.ndarray

Matrix of the right-hand side. A rectangular matrix with either float64 or float32 data type stored with row-major or column-major ordering. However, if overwrite=True, then B should have column-major ordering only.

shapetuple, default=None

A tuple of size two, determining the shape of an upper-left sub-matrices of ldu and B to be referenced. Namely, if shape is given as the tuple (m, n), the sub-matrices ldu[:m, :m] and B[:m, :n] are used. If None, the full shape of ldu[:, :] and B[:, :] are considered.

lowerbool, default=True

If True, ldu is assumed to contain the lower-triangular LDL factor \(\mathbf{L}\) such that \(\mathbf{A} = \mathbf{L} \mathbf{D} \mathbf{L}^{\intercal}\). If False, ldu is assumed to contain the upper-triangular LDL factor \(\mathbf{U}\) such that \(\mathbf{A} = \mathbf{U} \mathbf{D} \mathbf{U}^{\intercal}\).

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. If overwrite is set to True, the output matrix X is becomes a view for the matrix B.

See also

detkit.ldl_factor

detkit.lu_solve

detkit.cho_solve

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]}.\]

LDL Factor Input Argument:

The input to this function is the LDL factor of \(\mathbf{A}\) ( and not \(\mathbf{A}\) itself). Namely, if lower=True, it is assumed that the input argument ldu is the factor \(\mathbf{L} \mathbf{D} \mathbf{L}^{\intercal}\) where

\[\mathbf{A}_{[:m, :m]} = \mathbf{P}_{[:m, :m]} \mathbf{L}_{[:m, :m]} \mathbf{D}_{[:m, :m]} \mathbf{L}_{[:m, :m]}^{\intercal} \mathbf{P}_{[:m, :m]}^{\intercal}.\]

Similarly, if lower=False, it is assumed that the input argument ldu is the factor \(\mathbf{U} \mathbf{D} \mathbf{U}^{\intercal}\) where

\[\mathbf{A}_{:m, :m} = \mathbf{P}_{[:m, :m]} \mathbf{U}_{[:m, :m]} \mathbf{D}_{[:m, :m]} \mathbf{U}_{[:m, :m]}^{\intercal} \mathbf{P}_{[:m, :m]}^{\intercal}.\]

The shape of ldu 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.

The LDL factor of a matrix can be obtained from detkit.ldl_factor() and its output can be passed to the parameters ldu and piv in detkit.ldl_solve().

Note

When computing LDL decomposition with detkit.ldl_factor() to be used for detkit.ldl_solve(), set the argument return_as_lapack=True in detkit.ldl_factor().

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 slice X[: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:

Scipy does not provide a similar function.

Implementation:

This function is a wrapper around LAPACK’s sgetrf (for 32-bit precision) and dgetrf (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 ldl_factor, ldl_solve, Memory
>>> import numpy

>>> # Create a symmetric matrix with 32-bit precision and column-major
>>> # ordering
>>> A = numpy.random.randn(1000, 1000) + 100 * numpy.eye(1000, 1000)
>>> A = A + A.T
>>> 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 LDL
>>> # decomposition or solving linear system is now creating any new
>>> # memory.
>>> mem = Memory()
>>> mem.set()

>>> # LDL factorization of the upper-left sub-matrix
>>> p, q = 800, 700
>>> ldu, piv = ldl_factor(A, m=p, lower=True, overwrite=True,
...                       return_as_lapack=True)

>>> # Solve a linear system
>>> X = ldl_solve(ldu, piv, 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
>>> numpy.may_share_memory(X, B)
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.