pylops.utils.trace_nahutchpp

pylops.utils.trace_nahutchpp(Op, neval=None, sampler='rademacher', c1=0.16666666666666666, c2=0.3333333333333333, backend='numpy')[source]

Trace of linear operator using the NA-Hutch++ method.

Returns an estimate of the trace of a linear operator using the Non-Adaptive variant of 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.
c1 : float, optional

Fraction of neval for sketching matrix \(\mathbf{S}\).

c2 : float, optional

Fraction of neval for sketching matrix \(\mathbf{R}\). Must be larger than c2, ideally by a factor of at least 2.

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 not large enough to accomodate c1 and c2.

NotImplementedError

If the sampler is not one of the available samplers.

Notes

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

  1. Fix constants \(c_1\), \(c_2\), \(c_3\) such that \(c_1 < c_2\) and \(c_1 + c_2 + c_3 = 1\).
  2. Sample sketching matrices \(\mathbf{S} \in \mathbb{R}^{m \times c_1 k}\), \(\mathbf{R} \in \mathbb{R}^{m \times c_2 k}\), and \(\mathbf{G} \in \mathbb{R}^{m \times c_3 k}\) from sub-Gaussian distributions.
  3. Compute \(\mathbf{Z} = \mathbf{Op}\,\mathbf{R}\), \(\mathbf{W} = \mathbf{Op}\,\mathbf{S}\), and \(\mathbf{Y} = (\mathbf{S}^T \mathbf{Z})^+\), where \(+\) denotes the Moore–Penrose inverse.
  4. Return \(\operatorname{tr}(\mathbf{Y} \mathbf{W}^T \mathbf{Z}) + \frac{1}{c_3 k} \left[ \operatorname{tr}(\mathbf{G}^T\,\mathbf{Op}\,\mathbf{G}) - \operatorname{tr}(\mathbf{G}^T\mathbf{Z}\mathbf{Y}\mathbf{W}^T\mathbf{G})\right]\)

The default values for \(c_1\) and \(c_2\) are set to \(1/6\) and \(1/3\), respectively, but [1] suggests \(1/4\) and \(1/2\).

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

[1](1, 2, 3) 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