Skip to content

Commit

Permalink
[DOC] Short summary for each trace/diagonal estimation method
Browse files Browse the repository at this point in the history
  • Loading branch information
f-dangel committed Oct 17, 2023
1 parent 4efa2b3 commit e32f649
Show file tree
Hide file tree
Showing 3 changed files with 61 additions and 3 deletions.
21 changes: 20 additions & 1 deletion curvlinops/diagonal/hutchinson.py
Original file line number Diff line number Diff line change
Expand Up @@ -8,13 +8,32 @@


class HutchinsonDiagonalEstimator:
"""Class to perform diagonal estimation with Hutchinson's method.
r"""Class to perform diagonal estimation with Hutchinson's method.
For details, see
- Martens, J., Sutskever, I., & Swersky, K. (2012). Estimating the hessian by
back-propagating curvature. International Conference on Machine Learning (ICML).
Let :math:`\mathbf{A}` be a square linear operator. We can approximate its diagonal
:math:`\mathrm{diag}(\mathbf{A})` by drawing a random vector :math:`\mathbf{v}`
which satisfies :math:`\mathbb{E}[\mathbf{v} \mathbf{v}^\top] = \mathbf{I}` and
sample from the estimator
.. math::
\mathbf{a}
:= \mathbf{v} \odot \mathbf{A} \mathbf{v}
\approx \mathrm{diag}(\mathbf{A})\,.
This estimator is unbiased,
.. math::
\mathbb{E}[a_i]
= \sum_j \mathbb{E}[v_i A_{i,j} v_j]
= \sum_j A_{i,j} \mathbb{E}[v_i v_j]
= \sum_j A_{i,j} \delta_{i, j}
= A_{i,i}\,.
Example:
>>> from numpy import diag, mean, round
>>> from numpy.random import rand, seed
Expand Down
21 changes: 20 additions & 1 deletion curvlinops/trace/hutchinson.py
Original file line number Diff line number Diff line change
Expand Up @@ -7,14 +7,33 @@


class HutchinsonTraceEstimator:
"""Class to perform trace estimation with Hutchinson's method.
r"""Class to perform trace estimation with Hutchinson's method.
For details, see
- Hutchinson, M. (1989). A stochastic estimator of the trace of the influence
matrix for laplacian smoothing splines. Communication in Statistics---Simulation
and Computation.
Let :math:`\mathbf{A}` be a square linear operator. We can approximate its trace
:math:`\mathrm{Tr}(\mathbf{A})` by drawing a random vector :math:`\mathbf{v}`
which satisfies :math:`\mathbb{E}[\mathbf{v} \mathbf{v}^\top] = \mathbf{I}` and
sample from the estimator
.. math::
a
:= \mathbf{v}^\top \mathbf{A} \mathbf{v}
\approx \mathrm{Tr}(\mathbf{A})\,.
This estimator is unbiased,
.. math::
\mathbb{E}[a]
= \mathrm{Tr}(\mathbb{E}[\mathbf{v}^\top\mathbf{A} \mathbf{v}])
= \mathrm{Tr}(\mathbf{A} \mathbb{E}[\mathbf{v} \mathbf{v}^\top])
= \mathrm{Tr}(\mathbf{A} \mathbf{I})
= \mathrm{Tr}(\mathbf{A})\,.
Example:
>>> from numpy import trace, mean, round
>>> from numpy.random import rand, seed
Expand Down
22 changes: 21 additions & 1 deletion curvlinops/trace/meyer2020hutch.py
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,7 @@


class HutchPPTraceEstimator:
"""Class to perform trace estimation with the Huch++ method.
r"""Class to perform trace estimation with the Huch++ method.
In contrast to vanilla Hutchinson, Hutch++ has lower variance, but requires more
memory.
Expand All @@ -20,6 +20,26 @@ class HutchPPTraceEstimator:
- Meyer, R. A., Musco, C., Musco, C., & Woodruff, D. P. (2020). Hutch++:
optimal stochastic trace estimation.
Let :math:`\mathbf{A}` be a square linear operator whose trace we want to
approximate. First, we compute an orthonormal basis :math:`\mathbf{Q}` of a
sub-space spanned by :math:`\mathbf{A} \mathbf{S}` where :math:`\mathbf{S}` is a
tall random matrix with i.i.d. elements. Then, we compute the trace in the sub-space
and apply Hutchinson's estimator in the remaining space spanned by
:math:`\mathbf{I} - \mathbf{Q} \mathbf{Q}^\top`: We can draw a random vector
:math:`\mathbf{v}` which satisfies
:math:`\mathbb{E}[\mathbf{v} \mathbf{v}^\top] = \mathbf{I}` and sample from the
estimator
.. math::
a
:= \mathrm{Tr}(\mathbf{Q}^\top \mathbf{A} \mathbf{Q})
+ \mathbf{v}^\top (\mathbf{I} - \mathbf{Q} \mathbf{Q}^\top)^\top
\mathbf{A} (\mathbf{I} - \mathbf{Q} \mathbf{Q}^\top) \mathbf{v}
\approx \mathrm{Tr}(\mathbf{A})\,.
This estimator is unbiased, :math:`\mathbb{E}[a] = \mathrm{Tr}(\mathbf{A})`, as the
first term is constant and the second part is Hutchinson's estimator in a sub-space.
Example:
>>> from numpy import trace, mean, round
>>> from numpy.random import rand, seed
Expand Down

0 comments on commit e32f649

Please sign in to comment.