ldl_factor#
- detkit.ldl_factor(A, m=None, lower=True, overwrite=False, return_as_lapack=False)#
LDL decomposition of a sub-matrix.
- Parameters:
- Anumpy.ndarray
A symmetric matrix with either
float64
orfloat32
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 matrixA[:, :]
is used.- lowerbool, default=True
If True, ldu is assumed to contain the lower-triangular Cholesky factor \(\mathbf{L}\) such that \(\mathbf{A} = \mathbf{L} \mathbf{D} \mathbf{L}^{\intercal}\). If False, ldu is assumed to contain the upper-triangular Cholesky factor \(\mathbf{U}\) such that \(\mathbf{A} = \mathbf{U} \mathbf{D} \mathbf{U}^{\intercal}\).
- 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 ldu 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.- return_as_lapackboolean, default=False
If True, the outputs are encoded as given in LAPACK’s
dsytrf
function. This format is not intuitively interpretable, but can be directly used to pass as input todetkit.ldl_solve()
. If False, the outputs with LAPACK’s encoded format are post-processed to become the matriceslu
andd
and arrayperm
, wherelu
stores the matrix L (iflower
is True) or the matrix `U (iflower
is False), andperm
is the row permutations.Note
When using
return_as_lapack=False
, the operation is not memory-efficient as new memory will be allocated.
- Returns:
- If
return_as_lapack
is False: - lunumpy.ndarray
2D array 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 lower-triangular matrix L (iflower=True
) or the upper-triangular matrix U (iflower=False
).- dnumpy.array
2D array of the size of
m
containing the diagonal blocks of the size 1 by 1 or 2 by 2.- permnumpy.array
The row-permutation index array that brings
lu[:m, :m]
into triangular form. The size of this array ism
.
- If
return_as_lapack
is True: - ldunumpy.ndarray
2D matrix of the same shape as the input matrix A (and not the shape of the sub-matrix), containing the LDL decomposition. This matrix stores both L (or U) and D, however, the storage format of this matrix is not intuitively interpretable.
- pivnumpy.array
Row pivoting indices (zero-based indexing). The size of this array is
m
.
- If
Notes
LU decomposition of sub-matrix:
Let \(\mathbf{A}_{[:p,:q]}\) denote the sub-matrix of the size \(p \times q\) to be the upper-left corner of matrix \(\mathbf{A}\). Its LU decomposition is given by:
\[\mathbf{A}_{[:p, :q]} = \mathbf{P} \mathbf{L} \mathbf{U},\]where \(\mathbf{P}\) is a permutation matrix of the size \(p \times p\), \(\mathbf{L}\) is a lower-triangular matrix of the size \(p \times k\) with diagonals \(1\), where \(k = \min(p, q)\), and \(\mathbf{U}\) is an upper-triangular matrix of the size \(k \times q\).
Interpreting L and U factors from the output ldu:
Note that the output matrix
ldu
is given as the same shape asA
, and only its upper-left sub-matrix should be used to extract the matrices \(\mathbf{L}\) and \(\mathbf{U}\) as follows.Let
p, q = shape
be the given shape of the upper-left sub-matrix, which may be smaller than the shape ofA
. Also, letk = min(p, q)
. Then,L = ldu[:p, :k]
(excluding the diagonals) andU = ldu[:k, :q]
.The diagonals of
L
are not stored in the above matrix, and you should assume the diagonals ofL
are1
(see example below to adjust the diagonals ofL
to one manually).Interpreting permutation array:
The output
perm
is a 1D array of the size of the number of rows of the input matrix A (and not the number of rows of the sub-matrix). The arrayA[perm, :]
is equivalent of \(\mathbf{P}^{\intercal} \mathbf{A}\). As such, the relation relation \(\mathbf{P}^{\intercal} \mathbf{A} = \mathbf{L} \mathbf{U}\) can be achieved byA[perm[:p], :] = L @ U
.Comparison with scipy.linalg.ldl:
To compute the LDL factorization of a sub-matrix of the input matrix using
scipy.linalg.ldl
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.ldl_factor
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
ssytrf
(for 32-bit precision) anddsytrf
(for 64-bit precision).This function is internally used for
detkit.memdet()
for efficient computation of matrix determinant under memory constraint.References
Examples
Using
return_as_lapack=True
:>>> from detkit import ldl_factor, Memory >>> import numpy >>> # Create a symmetric matrix with 32-bit precision and column-major >>> # ordering >>> A = numpy.random.randn(1000, 1000) >>> A = A + A.T + 100 * numpy.eye(1000, 1000) >>> A = numpy.asfortranarray(A) >>> A = A.astype(numpy.float32) >>> # Track memory allocation to check if LU decomposition is not >>> # creating any new memory. >>> mem = Memory() >>> mem.set() >>> # LU factorization of the upper-left sub-matrix of smaller shape >>> m = 800 >>> ldu, piv = ldl_factor(A, m=800, overwrite=True, ... return_as_lapack=True) >>> # Check peak memory allocation (compared to memory of a sum-matrix) >>> slice_nbytes = m * m * A.dtype.itemsize >>> print(mem.peak() / slice_nbytes) 0.001 >>> # When overwrite is set to True, check if lu is indeed a view of A >>> print(ldu.base == A.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.Repeating the above example, but using
return_as_lapack=True
:>>> from detkit import ldl_factor >>> import numpy >>> # Create a symmetric matrix with 32-bit precision and column-major >>> # ordering >>> A = numpy.random.randn(1000, 1000) >>> A = A + A.T + 100 * numpy.eye(1000, 1000) >>> 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) >>> # LU factorization of the upper-left sub-matrix of smaller shape >>> m = 800 >>> lu, d, perm = ldl_factor(A, m=m, overwrite=False, ... return_as_lapack=False) >>> # Check if A = LDLt holds. >>> atol = numpy.finfo(A.dtype).resolution >>> print(numpy.allclose(lu[:m, :m] @ d[:m, :m] @ lu[:m, :m].T, ... A[:m, :m], atol=10*atol)) True