imate.traceinv(method=’cholesky’)#

imate.traceinv(A, gram=False, p=1, return_info=False, method='cholesky', B=None, invert_cholesky=True, cholmod=None)

Trace of inverse of non-singular matrix using Cholesky method.

See also

This page describes only the cholesky method. For other methods, see imate.traceinv().

Given the matrices \(\mathbf{A}\) and \(\mathbf{B}\) and the integer exponent \(p\), the following is computed:

\[\mathrm{trace} \left(\mathbf{A}^{-p} \mathbf{B} \right).\]

If B is None, it is assumed that \(\mathbf{B}\) is the identity matrix.

If gram is True, then \(\mathbf{A}\) in the above is replaced by the Gramian matrix \(\mathbf{A}^{\intercal} \mathbf{A}\), and the following is instead computed:

\[\mathrm{trace} \left((\mathbf{A}^{\intercal}\mathbf{A})^{-p} \mathbf{B} \right).\]
Parameters:
Anumpy.ndarray, scipy.sparse

A positive-definite sparse or dense matrix. If gram is True, the matrix can be non-square.

Warning

This function does not pre-check whether the input matrix is positive-definite.

Note

In the Cholesky method, the matrix cannot be a type of Matrix or imate.AffineMatrixFunction classes.

grambool, default=False

If True, the trace of the Gramian matrix, \((\mathbf{A}^{\intercal}\mathbf{A})^{-p}\), is computed. The Gramian matrix itself is not directly computed. If False, the trace of \(\mathbf{A}^{-p}\) is computed.

pfloat, default=1.0

The integer exponent \(p\) in \(\mathbf{A}^p\).

return_infobool, default=False

If True, this function also returns a dictionary containing information about the inner computation, such as process time, algorithm settings, etc.

Bnumpy.ndarray, scipy.sparse

A positive-definite sparse or dense matrix. B should be the same size and type of A. If B is None, is it assumed that the matrix B is identity.

Warning

This function does not pre-check whether the input matrix is positive-definite.

invert_choleskybool, default=True

If True, the inverse of Cholesky decomposition is computed. This approach is fast but it is only suitable for small matrices. If set to False, it uses the inverse of Cholesky matrix indirectly (see Notes). This approach is suitable for larger matrices but slower.

cholmodbool, default=None

If set to True, it uses the Cholmod library from scikit-sparse package to compute the Cholesky decomposition. If set to False, it uses scipy.sparse.cholesky method. If set to None, first, it tries to use Cholmod library, but if Cholmod is not available, then it uses scipy.sparse.cholesky method.

Returns:
traceinvfloat or numpy.array

Trace of inverse of A.

infodict

(Only if return_info is True) A dictionary of information with the following keys.

  • matrix:

    • data_type: str, {float32, float64, float128}. Type of the matrix data.

    • gram: bool, whether the matrix A or its Gramian is considered.

    • exponent: float, the exponent p in \(\mathbf{A}^{-p}\).

    • size: (int, int), the size of matrix A.

    • sparse: bool, whether the matrix A is sparse or dense.

    • nnz: int, if A is sparse, the number of non-zero elements of A.

    • density: float, if A is sparse, the density of A, which is the nnz divided by size squared.

    • num_inquiries: int, for the cholesky method, this is always 1.

  • device:

    • num_cpu_threads: int, number of CPU threads used in shared memory parallel processing.

    • num_gpu_devices: int, for the cholesky method, this is always 0.

    • num_gpu_multiprocessors: int, for the cholesky method, this is always 0.

    • num_gpu_threads_per_multiprocessor: int, for the cholesky method, this is always 0.

  • time:

    • tot_wall_time: float, total elapsed time of computation.

    • alg_wall_time: float, elapsed time of computation during only the algorithm execution.

    • cpu_proc_time: float, CPU processing time of computation.

  • solver:

    • cholmod_used: bool, whether the Cholmod from SparseSuite library was used.

    • version: str, version of imate.

    • method: ‘cholesky’

See also

imate.trace
imate.logdet
imate.schatten

Notes

Algorithm:

The trace of inverse is computed from the Cholesky decompositions \(\mathbf{A}^{p} = \mathbf{L}_{\vert p \vert} \mathbf{L}_{\vert p \vert}^{\intercal}\) and \(\mathbf{B} = \mathbf{L}_{\mathbf{B}} \mathbf{L}_{\mathbf{B}}^{\intercal}\) as follows:

\[\mathrm{trace} \left( \mathbf{A}^{-p} \mathbf{B} \right) = \Vert \mathbf{L}_{\vert p \vert}^{-1} \mathbf{L}_{\mathbf{B}} \Vert_F^2.\]

where \(\Vert \cdot \Vert_F\) is the Frobenius norm. If inverst_cholesky is True, the inverse \(\mathbf{L}_{\vert p \vert}^{-1}\) is computed directly. Inverting this matrix directly is only feasible for small matrices.

For larger matrices, set invert_cholesky to False. This approach is slower than setting invert_cholesky to True, however, it can process larger matrices. In this case, \(\Vert \mathbf{L}_{\vert p \vert}^{-1} \Vert_F^2\) is computed indirectly by

\[\Vert \mathbf{L}_{\vert p \vert}^{-1} \Vert_F^2 = \sum_{i=1}^n \Vert \boldsymbol{x}_i \Vert^2,\]

where \(\boldsymbol{x}_i\) is solved by the lower-triangular system

\[\mathbf{L}_{\vert p \vert} \boldsymbol{x}_i = \boldsymbol{b}_i,\]

where \(\boldsymbol{b}_i\) is the \(i\)-th column of \(\mathbf{L}_{\mathbf{B}}\). If B is None, then \(\mathbf{B}\) is assumed to be the identity and hence, \(\boldsymbol{b}_i = (0, \dots, 0, 1, 0, \dots, 0)\) is a vector of zeross, except its \(i\)-th element is 1.

Computational Complexity:

The computational complexity of this method is \(\mathcal{O}(\frac{1}{3}n^3)\) for dense matrices and \(\mathcal{O}(\rho n^2)\) for sparse matrices where \(n\) is the matrix size and \(\rho\) is the sparse matrix density.

Implementation:

This function is essentially a wrapper for the Cholesky function of the scipy and scikit-sparse packages and is primarily used for testing and comparison (benchmarking) against the randomized methods that are implemented in this package. If cholmod is set to True, this function uses the Suite Sparse package to compute the Cholesky decomposition.

Examples

Compute the trace of inverse of a sparse positive-definite Toeplitz matrix:

>>> # Import packages
>>> from imate import toeplitz, traceinv

>>> # Generate a sample symmetric and positive-definite matrix
>>> A = toeplitz(2, 1, size=100, gram=True)

>>> # Compute with Cholesky method (default method)
>>> traceinv(A, method='cholesky')
33.22222222222223

Print information about the inner computation:

>>> ti, info = traceinv(A, method='cholesky', return_info=True)
>>> print(ti)
33.22222222222223

>>> # Print dictionary neatly using pprint
>>> from pprint import pprint
>>> pprint(info)
{
    'matrix': {
        'data_type': b'float64',
        'density': 0.0298,
        'exponent': 1,
        'gram': False,
        'nnz': 298,
        'num_inquiries': 1,
        'size': (100, 100),
        'sparse': True
    },
    'solver': {
        'cholmod_used': True,
        'invert_cholesky': True,
        'method': 'cholesky',
        'version': '0.16.0'
    },
    'device': {
        'num_cpu_threads': 8,
        'num_gpu_devices': 0,
        'num_gpu_multiprocessors': 0,
        'num_gpu_threads_per_multiprocessor': 0
    },
    'time': {
        'alg_wall_time': 0.031367766903713346,
        'cpu_proc_time': 0.03275534099998367,
        'tot_wall_time': 0.031367766903713346
    }
}