imate.trace#
- imate.trace(A, gram=False, p=1.0, return_info=False, method='exact', **options)#
Trace of matrix or linear operator.
Given the matrix or the linear operator \(\mathbf{A}\) and the real non-negative exponent \(p \geq 0\), the following is computed:
\[\mathrm{trace} \left(\mathbf{A}^p \right).\]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 \right).\]If \(\mathbf{A} = \mathbf{A}(t)\) is a linear operator of the class
imate.AffineMatrixFunction
with the parameter \(t\), then for an input tuple \(t = (t_1, \dots, t_q)\), an array output of the size \(q\) is returned, namely:\[\mathrm{trace} \left((\mathbf{A}(t_i))^p \right), \quad i=1, \dots, q.\]- Parameters:
- Anumpy.ndarray, scipy.sparse,
imate.Matrix
, orimate.AffineMatrixFunction
A non-singular sparse or dense matrix or linear operator. The linear operators
imate.Matrix
andimate.AffineMatrixFunction
can be used only ifmethod=slq
. See details in slq method. Ifmethod=slq
andgram=False
, the input matrix A should be symmetric. Ifgram=True
, the matrix can be non-square.- 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 exponent \(p\) in \(\mathbf{A}^p\).
If
method=exact
, \(p\) should be a non-negative integer.If
method=eigenvalue
, \(p\) can be any real number.If
method=slq
, \(p\) should be non-negative real number.
Note
If \(p < 0\), use
imate.traceinv()
function.- 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. See the documentation for each method for details.
- method{‘exact’, ‘eigenvalue’, ‘slq’}, default=’exact’
The method of computing trace. See documentation for each method:
- options**kwargs
Extra arguments that are specific to each method. See the documentation for each method for details.
- Anumpy.ndarray, scipy.sparse,
- Returns:
- tracefloat or numpy.array
Trace of matrix. If
method=slq
and if A is of typeimate.AffineMatrixFunction
with an array ofparameters
, then the output is an array.- infodict
(Only if
return_info
is True) A dictionary of information with at least the following keys. Further keys specific to each method can be found in the documentation of each method.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, 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, The size of inquiries of each parameter of the linear operator A. If A is a matrix, this is always 1. For more details see slq method.
device
:num_cpu_threads
: int, number of CPU threads used in shared memory parallel processing.num_gpu_devices
: int, number of GPU devices used in the multi-GPU (GPU farm) computation.num_gpu_multiprocessors
: int, number of GPU multi-processors.num_gpu_threads_per_multiprocessor
: int, number of GPU threads on each GPU multi-processor.
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
:version
: str, version of imate.method
: str, method of computation.
- Raises:
- ImportError
If the package has not been compiled with GPU support, but
gpu
is True. Either setgpu
to False to use the existing installed package. Alternatively, export the environment variableUSE_CUDA=1
and recompile the source code of the package.
See also
Notes
Method of Computation:
See documentation for each method below.
exact: Computes trace directly from its diagonal entries. This is used when \(p\) is integer. If \(p=1\), this method is preferred.
eigenvalue: uses spectral decomposition. Suitable for small matrices (\(n < 2^{12}\)). The solution is exact.
slq: uses stochastic Lanczos quadrature (SLQ), which is a randomized algorithm. Can be used on very large matrices (\(n > 2^{12}\)). The solution is an approximation.
Note
If \(p=1\) and
gram
is False, always use exact method. If \(p\) is non-integer, you may use eigenvalue or slq method, though, for large matrices, the slq method is preferred.Input Matrix:
The input A can be either of:
A matrix, such as numpy.ndarray, or scipy.sparse.
A linear operator representing a matrix using
imate.Matrix
( only ifmethod=slq
).A linear operator representing a one-parameter family of an affine matrix function \(t \mapsto \mathbf{A} + t\mathbf{B}\), using
imate.AffineMatrixFunction
(only ifmethod=slq
).
Output:
The output is a scalar. However, if A is the linear operator of the type
imate.AffineMatrixFunction
representing the matrix function \(\mathbf{A}(t) = \mathbf{A} + t \mathbf{B}\), then if the parameter \(t\) is given as the tuple \(t = (t_1, \dots, t_q)\), then the output of this function is an array of size \(q\) corresponding to the trace of each \(\mathbf{A}(t_i)\).Note
When A represents \(\mathbf{A}(t) = \mathbf{A} + t \mathbf{I}\), where \(\mathbf{I}\) is the identity matrix, and \(t\) is given by a tuple \(t = (t_1, \dots, t_q)\), by setting
method=slq
, the computational cost of an array output of size q is the same as computing for a single \(t_i\). Namely, the trace of only \(\mathbf{A}(t_1)\) is computed, and the trace of the rest of \(q=2, \dots, q\) are obtained from the result of \(t_1\) immediately.Examples
Sparse matrix:
Compute the trace of a sample sparse Toeplitz matrix created by
imate.toeplitz()
function.>>> # Import packages >>> from imate import toeplitz, trace >>> # Generate a sample matrix (a toeplitz matrix) >>> A = toeplitz(2, 1, size=100) >>> # Compute trace with the exact method (default method) >>> trace(A) 200.0
Compute the trace of \((\mathbf{A}^{\intercal} \mathbf{A})^3\):
>>> # Compute trace of the Gramian of A to the power of 3. >>> trace(A, p=3, gram=True) 24307.0
Output information:
Print information about the inner computation:
>>> tr, info = trace(A, return_info=True) >>> print(tr) 200.0 >>> # Print dictionary neatly using pprint >>> from pprint import pprint >>> pprint(info) { 'matrix': { 'data_type': b'float64', 'density': 0.0199, 'exponent': 1.0, 'gram': False, 'nnz': 199, 'num_inquiries': 1, 'size': 100, 'sparse': True }, 'solver': { 'method': 'exact', 'version': '0.14.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.00013329205103218555, 'cpu_proc_time': 0.00017459900000016404, 'tot_wall_time': 0.00013329205103218555 } }
Large matrix:
Compute the trace of a very large sparse matrix using SLQ method. This method does not compute the trace exactly, rather, the result is an approximation using Monte-Carlo sampling. The following example uses at least 100 samples.
>>> # Generate a matrix of size one million >>> A = toeplitz(2, 1, size=1000000) >>> # Approximate trace using stochastic Lanczos quadrature >>> # with at least 100 Monte-Carlo sampling >>> tr, info = trace(A, method='slq', min_num_samples=100, ... max_num_samples=200, return_info=True) >>> print(tr) 4999741.080000001 >>> # Find the time it took to compute the above >>> print(info['time']) { 'tot_wall_time': 16.221865047933534, 'alg_wall_time': 16.20779037475586, 'cpu_proc_time': 116.213995219 }
Compare the result of the above approximation with the exact solution of the trace using the analytic relation for Toeplitz matrix. See
imate.sample_matrices.toeplitz_trace()
for details.>>> from imate.sample_matrices import toeplitz_trace >>> toeplitz_trace(2, 1, size=1000000) 4999999
It can be seen that the error of approximation is \(0.0018 \%\). This accuracy is remarkable considering that the computation on such a large matrix took only a 16 seconds. Computing the trace of such a large matrix using any of the exact methods (such as
exact
oreigenvalue
) is infeasible.Matrix operator:
The following example uses an object of
imate.Matrix
. Note that this can be only applied tomethod=slq
. See further details in slq method.>>> # Import matrix operator >>> from imate import toeplitz, trace, Matrix >>> # Generate a sample matrix (a toeplitz matrix) >>> A = toeplitz(2, 1, size=100, gram=True) >>> # Create a matrix operator object from matrix A >>> Aop = Matrix(A) >>> # Compute trace of Aop >>> trace(Aop, method='slq') 495.0
Affine matrix operator:
The following example uses an object of
imate.AffineMatrixFunction
to create the linear operator:\[t \mapsto \mathbf{A} + t \mathbf{I}\]Note that this can be only applied to
method=slq
. See further details in slq method.>>> # Import affine matrix function >>> from imate import toeplitz, trace, AffineMatrixFunction >>> # Generate a sample matrix (a toeplitz matrix) >>> A = toeplitz(2, 1, size=100, gram=True) >>> # Create a matrix operator object from matrix A >>> Aop = AffineMatrixFunction(A) >>> # A list of parameters t to pass to Aop >>> t = [-1.0, 0.0, 1.0] >>> # Compute trace of Aop with non-integer power for all parameters t >>> trace(Aop, method='slq', parameters=t) array([398.04, 498.04, 598.04])