pylops.utils.trace_hutchpp

pylops.utils.trace_hutchpp(Op, neval=None, sampler='rademacher', backend='numpy')[source]

Trace of linear operator using the Hutch++ method.

Returns an estimate of the trace of a linear operator using the Hutch++ method [1].

Parameters:
neval : int, optional

Maximum number of matrix-vector products compute. Defaults to 10% of shape[1].

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.
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

This function follows Algorithm 1 of [1]. Let \(m\) = shape[1] and \(k\) = neval.

  1. Sample sketching matrices \(\mathbf{S} \in \mathbb{R}^{m \times \lfloor k/3\rfloor}\), and \(\mathbf{G} \in \mathbb{R}^{m \times \lfloor k/3\rfloor}\), from sub-Gaussian distributions.
  2. Compute reduced QR decomposition of \(\mathbf{Op}\,\mathbf{S}\), retaining only \(\mathbf{Q}\).
  3. Return \(\operatorname{tr}(\mathbf{Q}^T\,\mathbf{Op}\,\mathbf{Q}) + \frac{1}{\lfloor k/3\rfloor}\operatorname{tr}\left(\mathbf{G}^T(\mathbf{I} - \mathbf{Q}\mathbf{Q}^T)\,\mathbf{Op}\,(\mathbf{I} - \mathbf{Q}\mathbf{Q}^T)\mathbf{G}\right)\)

Use the Rademacher sampler unless you know what you are doing.

[1](1, 2) Meyer, R. A., Musco, C., Musco, C., & Woodruff, D. P. (2021). Hutch++: Optimal Stochastic Trace Estimation. In Symposium on Simplicity in Algorithms (SOSA) (pp. 142–155). Philadelphia, PA: Society for Industrial and Applied Mathematics. link