pylops.utils.estimators.trace_hutchinson¶
- pylops.utils.estimators.trace_hutchinson(Op, neval=None, batch_size=None, sampler='rademacher', backend='numpy')[source]¶
Trace of linear operator using the Hutchinson method.
Returns an estimate of the trace of a linear operator using the Hutchinson method [1].
- Parameters
- Op
pylops.LinearOperator
Linear operator to compute trace on.
- neval
int
, optional Maximum number of matrix-vector products compute. Defaults to 10% of
shape[1]
.- batch_size
int
, optional Vectorize computations by sampling sketching matrices instead of vectors. Set this value to as high as permitted by memory, but there is no guarantee of speedup. Coerced to never exceed
neval
. When using “unitvector” as sampler, is coerced to not exceedshape[1]
. Defaults to 100 orneval
.- sampler
str
, optional Sample sketching matrices from the following distributions:
“gaussian”: Mean zero, unit variance Gaussian.
“rayleigh”: Sample from mean zero, unit variance Gaussian and normalize the columns.
“rademacher”: Random sign.
“unitvector”: Samples from the unit vectors \(\mathrm{e}_i\) without replacement.
- backend
str
, optional Backend used to densify matrix (
numpy
orcupy
). Note that this must be consistent with how the operator has been created.
- Op
- Returns
- trace
float
Operator trace.
- trace
- Raises
- ValueError
If
neval
is smaller than 3.- NotImplementedError
If the
sampler
is not one of the available samplers.
Notes
Let \(m\) =
shape[1]
and \(k\) =neval
. This algorithm estimates the trace via\[\frac{1}{k}\sum\limits_{i=1}^k \mathbf{z}_i^T\,\mathbf{Op}\,\mathbf{z}_i\]where vectors \(\mathbf{z}_i\) are sampled according to the sampling function. See [2] for a description of the variance and \(\epsilon\)-approximation of different samplers.
Prefer the Rademacher sampler if the goal is to minimize variance, but the Gaussian for a better probability of approximating the correct value. Use the Unit Vector approach if you are sampling a large number of
neval
(compared toshape[1]
), especially if the operator is highly-structured.- 1
Hutchinson, M. F. (1990). A stochastic estimator of the trace of the influence matrix for laplacian smoothing splines. Communications in Statistics - Simulation and Computation, 19(2), 433–450.
- 2
Avron, H., and Toledo, S. (2011). Randomized algorithms for estimating the trace of an implicit symmetric positive semi-definite matrix. Journal of the ACM, 58(2), 1–34.