diff --git a/curvlinops/diagonal/hutchinson.py b/curvlinops/diagonal/hutchinson.py index 1d82a52d..49c11d45 100644 --- a/curvlinops/diagonal/hutchinson.py +++ b/curvlinops/diagonal/hutchinson.py @@ -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 diff --git a/curvlinops/trace/hutchinson.py b/curvlinops/trace/hutchinson.py index 4765dcbd..a2d1506c 100644 --- a/curvlinops/trace/hutchinson.py +++ b/curvlinops/trace/hutchinson.py @@ -7,7 +7,7 @@ class HutchinsonTraceEstimator: - """Class to perform trace estimation with Hutchinson's method. + r"""Class to perform trace estimation with Hutchinson's method. For details, see @@ -15,6 +15,25 @@ class HutchinsonTraceEstimator: 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 diff --git a/curvlinops/trace/meyer2020hutch.py b/curvlinops/trace/meyer2020hutch.py index f4fe7f8b..e37f32e3 100644 --- a/curvlinops/trace/meyer2020hutch.py +++ b/curvlinops/trace/meyer2020hutch.py @@ -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. @@ -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