logpdet#
- detkit.logpdet(A, X, Xp=None, method='legacy', sym_pos=False, X_orth=False, flops=False, use_scipy=False)#
Compute 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.
- 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:
- 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
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='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 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 to1
, this function is parallelized. The compile-time default is0
. To see a list of compile-time definitions, seedetkit.get_config()
function.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 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)