loggdet#

detkit.loggdet(A, X, Xp=None, method='legacy', sym_pos=False, X_orth=False, flops=False, use_scipy=False)#

Computes the loggdet of a matrix.

The loggdet function is defined by

\[\mathrm{loggdet}(\mathbf{A}, \mathbf{X}) := \log_{e} |\mathrm{gdet}(\mathbf{A}, \mathbf{X})|,\]

where

\[\mathrm{gdet}(\mathbf{A}, \mathbf{X}) := \det(\mathbf{A}) \det(\mathbf{X}^{\intercal} \mathbf{A}^{-1} \mathbf{X}),\]

and

  • \(\mathbf{A}\) is an \(n \times n\) invertible matrix.

  • \(\mathbf{X}\) is an \(n \times m\) full column-rank matrix.

Parameters:
A(n, n) array_like

Invertible matrix. The matrix type can be float32, float64, or float128. If a matrix of the type int32 or int64 is given, the type is cast to float64.

X(n, m) array_like

Rectangular matrix with full column-rank.

Xp(n, n-m) array_like, default: None

Rectangular matrix with full column-rank. Xp is the orthonormal complement of X. If None, this matrix will be generated. Only used if method is comp.

method{‘legacy’, ‘proj’, ‘comp’}, default=’legacy’

Method of computing, and can be one of legacy or proj.

  • ‘legacy’: Computes loggdet directly by the equation given in the above.

  • ‘proj’: Computes loggdet using Bott-Duffin inverse (See [1]).

  • ‘comp’: Computes loggdet using compression matrix (See [1]).

sym_posbool, default=False

If True, the matrix A is assumed to be symmetric and positive-definite (SPD). The computation can be twice as fast as when the matrix is not SPD. This function does not verify whether A is symmetric or positive-definite. This option is applicable when method='legacy' or method='comp'.

X_orthbool, default=False

If True, the matrix X is assumed to have orthogonal columns. The computation in this case is faster. This function will not verify whether X is orthogonal. This option is only applicable when method='proj' or method='comp'.

flopsbool, default=False

If True, the count of the retired hardware instructions is returned using perf tool. This functionality is only available on Linux operating system and only newer models of Intel CPUs. This option is only relevant when use_scipy=False.

use_scipybool, default=False

If True, it uses scipy functions which are the wrappers around Fortran routines in BLAS and LAPACK. If False, it uses a C++ library developed in this package.

Returns:
loggdetfloat

loggdet of A.

signint

Sign of gdet function and can be +1 for positive or -1 for negative determinant.

flopsint

(returned if flops=True and use_scipy=False)

Count of the retired hardware instructions of the processor during the runtime of this function.

Raises:
RuntimeError

Error raised when:

  • sym_pos=True and matrix A is not symmetric positive-definite.

  • method='legacy' and matrix A is degenerate.

  • flops=True and either Perf tool is not installed (Linux only), or the user permission for the Perf tool is not set, or the performance counter is not supported on the user’s CPU.

ValueError

Error raised when

  • method='legacy' and X_orth=True.

  • method='proj' and sym_pos=True.

Warning

When method='proj' and A is singular, no error is raised and an incorrect value is returned.

See also

logdet

Computes the logdet of a matrix.

logpdet

Log of pseudo-determinant of the precision matrix in Gaussian process regression.

Notes

When the method is legacy, the loggdet function is computed using the equation given in the above. However, when the method is set to ‘proj’ or ‘comp’, an alternative formulation is used. Note that:

  • legacy method is independent of whether X is orthogonal or not, thus, cannot take advantage of an orthogonal X.

  • proj method is independent of whether A is symmetric positive-definite or not, thus, cannot take advantage of an SPD matrix A.

The inner computation of the function loggdet applies the following algorithms:

  • When sym_pos=False, the logdet function is computed using the PLU decomposition of A.

  • When sym_pos=True, the logdet function is computed using the Cholesky decomposition of A.

This function is not parallelized and does not accept sparse matrices.

The loggdet function is used in the likelihood function of the Gaussian process regression (see [2]).

Counting Flops:

When flops is set to True, make sure perf tool is installed (Linux only). On Ubuntu, install perf by

sudo apt-get install linux-tools-common linux-tools-generic \
        linux-tools-`uname -r`

sudo sh -c 'echo -1 >/proc/sys/kernel/perf_event_paranoid'

Test if the perf tool works by

perf stat -e instructions:u dd if=/dev/zero of=/dev/null count=1000

Alternatively, you may also test perf tool with detkit by

>>> import detkit
>>> detkit..get_instructions_per_task()

If the installed perf tool is configured properly, the output of either of the above commands should not be empty.

References

[1] (1,2)

Ameli, S., Shadden, S. C. (2022) A Singular Woodbury and Pseudo-Determinant Matrix Identities and Application to Gaussian Process Regression.

[2]

Ameli, S., Shadden, S. C. (2022) Noise Estimation in Gaussian Process Regression.

Examples

>>> import numpy
>>> from detkit import loggdet, orthogonalize

>>> # Generate a random matrix
>>> n, m = 1000, 500
>>> rng = numpy.random.RandomState(0)
>>> A = rng.rand(n, n)
>>> X = rng.rand(n, m)

>>> # Compute loggdet of matrix
>>> loggdet(A, X)
(1710.9576831500378, -1)

>>> # Compute loggdet of a symmetric and positive-definite matrix
>>> B = A.T @ A
>>> loggdet(B, X, sym_pos=True)
(3421.9153663693114, 1)

>>> # Compute loggdet when X is orthogonal
>>> orthogonalize(X)
>>> loggdet(B, X, sym_pos=True, X_orth=True)
(3421.9153663693114, 1)

>>> # Compute loggdet of a singular matrix
>>> A[:, 0] = 0
>>> loggdet(A, X)
RuntimeError: LUP factorization failed since matrix "A" is degenerate.

The count of the hardware instructions is only supported on the Linux operating system and on recent Intel processors.

>>> # Count the processor instructions
>>> loggdet(B, X, sym_pos=True, X_orth=True, flops=True)
(3421.9153663693114, 1)