pylops.utils.estimators.trace_nahutchpp

pylops.utils.estimators.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
Oppylops.LinearOperator

Linear operator to compute trace on.

nevalint, optional

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

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

c1float, optional

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

c2float, optional

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

backendstr, optional

Backend used to densify matrix (numpy or cupy). Note that this must be consistent with how the operator has been created.

Returns
tracefloat

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