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 thanc2
, ideally by a factor of at least 2.- 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
not large enough to accomodatec1
andc2
.- 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
.- Fix constants \(c_1\), \(c_2\), \(c_3\) such that \(c_1 < c_2\) and \(c_1 + c_2 + c_3 = 1\).
- 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.
- 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.
- 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 - neval :