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