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
orcupy
). 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
.- 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.
- Compute reduced QR decomposition of \(\mathbf{Op}\,\mathbf{S}\), retaining only \(\mathbf{Q}\).
- 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 - neval :