PyLops
pylops.utils.estimators.trace_hutchpp¶
- pylops.utils.estimators.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:
- Op
pylops.LinearOperator Linear operator to compute trace on.
- 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 (
numpyorcupy). Note that this must be consistent with how the operator has been created.
- Op
- Returns:
- trace
float Operator trace.
- trace
- Raises:
- ValueError
If
nevalis smaller than 3.- NotImplementedError
If the
sampleris 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.