pylops.utils.trace_hutchinson

pylops.utils.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:
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 exceed shape[1]. Defaults to 100 or neval.

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 or cupy). Note that this must be consistent with how the operator has been created.

Returns:
trace : self.dtype

Operator 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 to shape[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.