# 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. trace : self.dtype Operator trace. 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