solve_triangular

detkit.solve_triangular(A, B, shape=None, trans=False, lower=False, unit_diagonal=False, overwrite=False)

Solve triangular linear system given a sub-matrix.

Parameters:

Anumpy.ndarray

Matrix of coefficients. A triangular matrix with either float64 or float32 data type stored with column-major ordering.

Bnumpy.ndarray

Matrix of the right-had side. A rectangular matrix with either float64 or float32 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 A and B to be referenced. Namely, if shape is given as the tuple (m, n), the sub-matrices A[:m, :m] and B[:m, :n] are used. If None, the full shape of A[:, :] and B[:, :] are considered.

transbool, default=False

If False, the system \(\mathbf{A} \mathbf{X} = \mathbf{B}\) is solved. If True, the system \(\mathbf{A}^{\intercal} \mathbf{X} = \mathbf{B}\) is solved.

lowerbool, default=False

If True, A is assumed to be lower-triangular. If False, A is assumed to be upper-triangular.

unit_diagonalbool, default=False

If True, the diagonals of A are assumed to be 1, even though a different value of diagonals are stored on the memory.

overwritebool, default=False

If True, the output X is overwritten to B, hence, X and B would share the same memory. If False, a new memory will be allocated for the output X.

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.lu_factor

detkit.ldl_solve

detkit.cho_solve

detkit.cho_factor

detkit.lu_factor

detkit.ldl_factor

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

if trans is False, or

\[\mathbf{A}_{[:m, :m]}^{\intercal} \mathbf{X}_{[:m, :n]} = \mathbf{B}_{[:m, :n]},\]

if trans is True.

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 relation

A[:m, :m] @ X[:m, :n] = B[:m, :n]

or

A[:m, :m].T @ X[:m, :n] = B[:m, :n]

(when trans is True) should hold.

Comparison with scipy.linalg.solve_triangular:

To solve a linear system for a sub-matrix of the input matrices using scipy.linalg.solve_triangular 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.solve_triangular 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 strtrs (for 32-bit precision) and dtrtrs (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 solve_triangular, Memory
>>> import numpy

>>> # Create a lower-triangular matrix with 32-bit precision and
>>> # column-major ordering
>>> A = numpy.random.randn(1000, 900) + 100 * numpy.eye(1000, 900)
>>> A = numpy.tril(A)
>>> A = numpy.asfortranarray(A)
>>> A = A.astype(numpy.float32)

>>> # Create the matrix of right-hand side
>>> B = numpy.random.randn(900, 800)
>>> B = numpy.asfortranarray(B)
>>> B = B.astype(numpy.float32)

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

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

>>> # Solve the system A X = B for the sub-matrix A[:m, :m] and
>>> # B[:m, :n]
>>> m, n = (800, 700)
>>> X = solve_triangular(A, B, shape=(m, n), lower=True,
...                      overwrite=True)

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

>>> # Check if A @ X = B_copy holds.
>>> atol = numpy.finfo(A.dtype).resolution
>>> print(numpy.allclose(A[:m, :m] @ X[:m, :n], B_copy[:m, :n],
...                      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 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.