imate.trexp(method=’eigenvalue’)#

imate.trexp(A, gram=False, p=1.0, coeff=1.0, return_info=False, method='eigenvalue', eigenvalues=None, assume_matrix='gen', non_zero_eig_fraction=0.9)

Trace of the exponential of matrix using eigenvalue method.

See also

This page describes only the eigenvalue method. For other methods, see imate.trexp().

Given the matrix \(\mathbf{A}\) and the real exponent \(p\), the following is computed:

\[\mathrm{trace} \; \mathrm{exp} \left(c \mathbf{A}^p \right),\]

where \(c\) is a given real coefficient.

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} \; \mathrm{exp} \left(c (\mathbf{A}^{\intercal}\mathbf{A})^p \right).\]

Parameters:

Anumpy.ndarray, scipy.sparse: A sparse or dense matrix. If gram is True, the matrix can be non-square.

Note

In the eigenvalue method, the matrix cannot be a type of Matrix or imate.AffineMatrixFunction classes.
grambool, default=False: If True, the log-determinant of the Gramian matrix, \((\mathbf{A}^{\intercal}\mathbf{A})^p\), is computed. The Gramian matrix itself is not directly computed. If False, the trace-exponential of \(\mathbf{A}^p\) is computed.
pfloat, default=1.0: The real exponent \(p\) in \(\mathbf{A}^p\).
coefffloat, default=1.0: The coefficient \(c\) in \(\mathrm{trace} \mathrm{exp} \left(c \mathbf{A}^p \right)\).
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.
eigenvaluesarray_like [float], default=None: The array of all eigenvalues of A, if available. The size of this array must be the same as the size of A. If None, the eigenvalues of A are computed.
assume_matrixstr {‘gen’, ‘sym’}, default: ‘gen’: Type of matrix. gen assumes generic matrix, while sym assumes A is symmetric.
non_zero_eig_fractionfloat, default=0.9: A fraction (between 0 and 1) of eigenvalues assumed to be non-zero. For large matrices, it is not possible to compute all eigenvalues, and only the largest eigenvalues can be computed and the rest are assumed to be negligible. By setting this parameter, a fraction of non-negligible eigenvalues is determined.

Returns:

trexpfloat or numpy.array

Trace of exponential 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\).
- assume_matrix: str, {gen, sym}, determines whether matrix is generic or symmetric.
- 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 eigenvalue 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 eigenvalue method, this is always 0.
- num_gpu_multiprocessors: int, for the eigenvalue method, this is always 0.
- num_gpu_threads_per_multiprocessor: int, for eigenvalue 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, the CPU processing time of computation.
solver:
- version: str, version of imate.
- method: ‘eigenvalue’

See also

imate.logdet
imate.trace
imate.traceinv
imate.schatten

Notes

Computational Complexity:

The eigenvalue method uses spectral decomposition. The computational complexity of this method is \(\mathcal{O}(n^3)\) where \(n\) is the matrix size. This method is only suitable for small matrices (\(n < 2^{12}\)). The solution is exact and can be used as a benchmark to test randomized methods of computing trace-exponential.

Warning

It is not recommended to use this method for sparse matrices, as not all eigenvalues of sparse matrices can be computed.

Examples

Dense matrix:

Compute the trace-exponential of a sample sparse Toeplitz matrix created by imate.toeplitz() function:

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

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

>>> # Convert the sparse matrix to a dense matrix
>>> A = A.toarray()

>>> # Compute trace-exponential with Cholesky method (default method)
>>> trexp(A, method='eigenvalue', assume_matrix='sym')
138.62943611198907

Precomputed Eigenvalues:

Alternatively, compute the eigenvalues of A in advance, and pass it to this function:

>>> # Compute eigenvalues of symmetric matrix A.
>>> from scipy.linalg import eigh
>>> eigenvalues = eigh(A, eigvals_only=True)

>>> # Pass the eigenvalues to trace-exponential function
>>> trexp(A, method='eigenvalue', eigenvalues=eigenvalues)
138.62943611198907

Pre-computing eigenvalues can be useful if imate.trexp() function should be called repeatedly for the same matrix A but other parameters may change, such as p.

Print Information:

Print information about the inner computation:

>>> ld, info = trexp(A, method='eigenvalue', return_info=True)
>>> print(ld)
138.6294361119891

>>> # Print dictionary neatly using pprint
>>> from pprint import pprint
>>> pprint(info)
{
    'matrix': {
        'assume_matrix': 'gen',
        'data_type': b'float64',
        'density': 1.0,
        'exponent': 1.0,
        'gram': False,
        'nnz': 10000,
        'num_inquiries': 1,
        'size': (100, 100),
        'sparse': False
    },
    'solver': {
        'method': 'eigenvalue',
        'version': '0.15.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.007327683997573331,
        'cpu_proc_time': 0.014451992999999774,
        'tot_wall_time': 0.007327683997573331
    }
}

Sparse matrix:

Using a large matrix and ignoring 10% of its eigenvalues:

>>> # Generate a symmetric sparse matrix
>>> A = toeplitz(2, 1, size=2000, gram=True)

>>> # Assume only 80% of eigenvalues of A are non-zero
>>> trexp(A, method='eigenvalue', assume_matrix='sym',
...       non_zero_eig_fraction=0.9)
2451.2640192906174