logpdet#

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

Computes the logpdet of a matrix.

The logpdet function is defined by

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

where

\[\mathrm{pdet}(\mathbf{A}, \mathbf{X}) := \det(\mathbf{X}^{\intercal} \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.

The value of logpdet is independent of whether \(\mathbf{X}\) is orthogonal or not.

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

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 logpdet directly by the equation given in the above.

  • ‘proj’: Computes logpdet 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:
logpdetfloat

logpdet 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='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.

loggdet

Log of determinant terms used in Gaussian process regression.

Notes

When the method is legacy, the logpdet 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:

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

  • The value of logpdet function is independent of whether X is orthogonal or not. However, when X is orthogonal, the option X_orth=True can enhance the performance of both proj and legacy method.

The inner computation of the function logpdet 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.

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

If the compile-time variable USE_OPENMP is set to 1, this function is parallelized. The compile-time default is 0. To see a list of compile-time definitions, see detkit.get_config() function.

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 logpdet, 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 logpdet of matrix
>>> logpdet(A, X)
(-680.9183141420649, -1)

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

>>> # Compute logpdet when X is orthogonal
>>> orthogonalize(X)
>>> logpdet(B, X, sym_pos=True, X_orth=True)
(-2059.6208046883685, 1)

>>> # Compute logpdet of a singular matrix
>>> A[:, 0] = 0
>>> logpdet(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
>>> logpdet(B, X, sym_pos=True, X_orth=True, flops=True)
(-2059.6208046883685, 1, 8520537034)