loggdet#
- detkit.loggdet(A, X, Xp=None, method='legacy', sym_pos=False, X_orth=False, flops=False, use_scipy=False)#
Compute 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.
- 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'
ormethod='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'
ormethod='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
anduse_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'
andX_orth=True
.method='proj'
andsym_pos=True
.
Warning
When
method='proj'
and A is singular, no error is raised and an incorrect value is returned.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 bysudo 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)