# Monte Carlo Methods for Estimating the Diagonal of a Real Symmetric Matrix

Eric Hallman Department of Mathematics, North Carolina State University,
Raleigh, NC, 27695<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>Ilse
C.F. Ipsen11footnotemark: 1 Arvind K. Saibaba11footnotemark: 1

###### Abstract

For real symmetric matrices that are accessible only through matrix vector
products, we present Monte Carlo estimators for computing the diagonal
elements. Our probabilistic bounds for normwise absolute and relative errors
apply to Monte Carlo estimators based on random Rademacher, sparse Rademacher,
normalized and unnormalized Gaussian vectors, and to vectors with bounded
fourth moments. The novel use of matrix concentration inequalities in our
proofs represents a systematic model for future analyses. Our bounds mostly do
not depend on the matrix dimension, target different error measures than
existing work, and imply that the accuracy of the estimators increases with
the diagonal dominance of the matrix. An application to derivative-based
global sensitivity metrics corroborates this, as do numerical experiments on
synthetic test matrices. We recommend against the use in practice of sparse
Rademacher vectors, which are the basis for many randomized sketching and
sampling algorithms, because they tend to deliver barely a digit of accuracy
even under large sampling amounts.

###### keywords:

Concentration inequalities, Monte Carlo Methods, Relative error, Rademacher
random vectors, Gaussian random vectors

15A15, 65C05, 65F50, 60G50, 68W20

##  1 Introduction

We compute the diagonal elements of symmetric matrices
$\boldsymbol{A}\in\mathbb{R}^{n\times n}$ with Monte Carlo estimators of the
form

$\widehat{\boldsymbol{A}}=\frac{1}{N}\sum_{k=1}^{N}{\boldsymbol{A}\boldsymbol{z}_{k}\boldsymbol{z}_{k}^{\top}}$

where $\boldsymbol{z}_{k}$ are independent random vectors. This approach is
crucial when the elements of $\boldsymbol{A}$ are available only implicitly,
via matrix vector products.

Estimating the diagonal elements of a matrix is important in many areas of
science and engineering: In electronic structure calculations, one computes
the diagonal elements of a projector onto the smallest eigenvectors of a
Hamiltonian matrix [4]. In statistics, leverage scores for column subset
selection can be computed from the diagonals of the projector onto the column
space. In Bayesian inverse problems, the diagonal elements of the posterior
covariance are computed with matrix-free estimators. Diagonal, or Jacobi
preconditioners can accelerate the convergence of iterative linear solvers
[21]. More recently, diagonal estimators have been used to accelerate second
order optimization techniques for machine learning [22]. In network science,
subgraph centrality measures and ranks the importance of the network nodes
based on the diagonal of a scaled exponential of the adjacency matrix. In
sensitivity analysis, Monte Carlo diagonal estimators efficiently compute the
derivative-based global sensitivity metrics [6, 10].

Diagonal estimation is related to trace estimation. Once the diagonal elements
are known, the trace can be computed from their sum. Therefore, estimators for
the diagonal of a matrix can be easily adapted to trace estimators. Monte
Carlo methods were first proposed by Hutchinson [8], and subsequently improved
and expanded to different distributions [2, 7, 16]. Applications of trace
estimators, reviewed in [19], include estimating density of states, log
determinants, and Schatten $p$-norms.

#### Literature review

To our knowledge, Monte Carlo diagonal estimators were first proposed by
Bekas, Kokiopoulou, and Saad [4], and a sufficient condition was given for a
Monte Carlo estimator to be unbiased. However, this paper identified that
large offdiagonal entries can result in large relative errors and developed
probing methods to mitigate the effects of the offdiagonal entries. This idea
is further explored in the following works [11, 9].

We are aware of a recent paper [3] as the only other work to analyze the
number of samples required for a relative $(\epsilon,\delta)$ estimator. In
contrast to [3], our proofs are the first to exploit matrix concentration
inequalities to impose a systematic structure that can serve as a model for
future analyses and allow us to analyze the normwise errors in a different
norm. We analyze more general distributions such as random vectors with
bounded fourth moments and sparse Rademacher vectors with a user-specified
sparsity parameter, and—in contrast to [3]—focus on un-normalized estimators.
Most of our bounds do not show an explicit dependence on the matrix dimension
which is desirable for large-scale problems.

###  1.1 Contributions and overview

After introducing notation, relevant concentration inequalities, and the setup
for our analysis (section 2.1), we derive normwise error bounds for Monte
Carlo estimators based on independent Rademacher vectors (section 3), random
vectors with bounded fourth moments and Gaussian vectors (section 4);
componentwise bounds for Rademacher and Gaussian vectors (section 5); and
apply Monte Carlo estimators to derivative-based global sensitivity metrics
(section 6). Numerical experiments (section 7) illustrate the accuracy of the
Monte Carlo estimators and the bounds. The novel and noteworthy features of
contributions are:

  1. 1.

Most of our bounds do not depend on the matrix dimension $n$, and hold for all
symmetric matrices, whether positive definite or not.

  2. 2.

We extend the concept of relative $(\epsilon,\delta)$ estimators to diagonal
estimation to determine the minimal number of samples $N$ for a user-specified
choice of relative error $\epsilon$ and failure probability $\delta$
(Definitions 2.6, 2.7).

  3. 3.

Our normwise bounds suggest that for Rademacher vectors, the Monte Carlo
estimators are more accurate for matrices that are more strongly diagonally
dominant (Theorem 3.4). In particular, the least number of samples required
for the Monte Carlo estimators to achieve a user-specified relative error
decreases with increasing diagonal dominance of $\boldsymbol{A}$ in the
relative sense (Corollaries 3.6, 4.3).

  4. 4.

For Rademacher vectors parameterized in terms of sparsity levels (Definition
3.8), we show that the Monte Carlo estimators lose accuracy with increasing
sparsity (Theorem 3.10, Corollary 3.13). Numerical experiments (section 7)
confirm that, even for large sampling amounts, the estimators barely achieve a
single digit of accuracy. Therefore we recommend against their use in
practice.

  5. 5.

Our componentwise bounds suggest that the accuracy for computing a diagonal
element $a_{ii}$ depends only on the diagonal dominance of column/row $i$ of
$\boldsymbol{A}$ (Corollaries 5.7, 5.11).

  6. 6.

In the context of derivative-based global sensitivity metrics, we design and
analyze Monte Carlo estimators based on random vectors from a problem-specific
probability distribution (Theorem 6.1, Corollary 6.3).

##  2 Background

After reviewing notation (section 2.1) and relevant concentration inequalities
(section 2.2), we present the setup for our analysis (section 2.3).

###  2.1 Notation

The Schur product (or Hadamard, or elementwise product) of
$\boldsymbol{A},\boldsymbol{B}\in\mathbb{R}^{m\times n}$ is denoted by
$\boldsymbol{C}=\boldsymbol{A}\circ\boldsymbol{B}\in\mathbb{R}^{m\times n}$
and has elements

$c_{ij}=a_{ij}b_{ij}\qquad 1\leq 1\leq m,\quad 1\leq j\leq n.$

For $\boldsymbol{A},\boldsymbol{B},\boldsymbol{C}\in{\mathbb{R}}^{m\times n}$,
the Schur product is commutative and distributive,

$\displaystyle\boldsymbol{A}\circ\boldsymbol{B}=\boldsymbol{B}\circ\boldsymbol{A},\qquad\boldsymbol{A}\circ(\boldsymbol{B}+\boldsymbol{C})=\boldsymbol{A}\circ\boldsymbol{B}+\boldsymbol{A}\circ\boldsymbol{C}.$

Following MATLAB convention, we define
$\operatorname{\mathrm{diag}}(\boldsymbol{A})=\begin{pmatrix}a_{11}&\cdots&a_{nn}\end{pmatrix}^{\top}\in\mathbb{R}^{n}$
as the column vector of diagonal elements of
$\boldsymbol{A}\in{\mathbb{R}}^{n\times n}$. The operator
$\operatorname{\mathrm{diag}}$ is overloaded, and
$\operatorname{\mathrm{diag}}(\boldsymbol{x})\in\mathbb{R}^{n\times n}$
represents a diagonal matrix whose diagonal elements are the elements of the
vector $\boldsymbol{x}\in\mathbb{R}^{n}$. In particular,

(1)
$\displaystyle\mathcal{D}\,(\boldsymbol{A})\equiv\operatorname{\mathrm{diag}}(\operatorname{\mathrm{diag}}(\boldsymbol{A}))=\boldsymbol{I}\circ\boldsymbol{A}\
\in\mathbb{R}^{n\times n}$

represents the diagonal matrix whose diagonal elements are the diagonal
elements of $\boldsymbol{A}$. In other words,
$\boldsymbol{I}\circ\boldsymbol{A}$ zeros out the offdiagonal elements of
$\boldsymbol{A}$.

If the first factor in a Schur product is a square matrix
$\boldsymbol{M}\in\mathbb{R}^{n\times n}$, and the second factor an outer
product involving $\boldsymbol{x},\boldsymbol{y}\in\mathbb{R}^{n}$, then

(2)
$\displaystyle\boldsymbol{M}\circ(\boldsymbol{xy}^{\top})=\operatorname{\mathrm{diag}}(\boldsymbol{x})\,\boldsymbol{M}\,\operatorname{\mathrm{diag}}(\boldsymbol{y}).$

For symmetric matrices $\boldsymbol{A},\boldsymbol{B}\in\mathbb{R}^{n\times
n}$, the partial order $\boldsymbol{A}\preceq\boldsymbol{B}$, or equivalently
$\boldsymbol{B}\succeq\boldsymbol{A}$, says that
$\boldsymbol{B}-\boldsymbol{A}$ is positive semidefinite. If $\boldsymbol{A}$
and $\boldsymbol{B}$ are positive semidefinite, then
$\boldsymbol{A}\preceq\boldsymbol{B}$ implies
$\boldsymbol{A}^{1/2}\preceq\boldsymbol{B}^{1/2}$.

The intrinsic dimension of a nonzero symmetric positive semidefinite matrix
$\boldsymbol{A}\in\mathbb{R}^{n\times n}$ is

$\operatorname{\mathrm{intdim}}(\boldsymbol{A})\equiv\frac{\operatorname{\mathrm{trace}}(\boldsymbol{A})}{\|\boldsymbol{A}\|_{2}},\qquad\text{with}\quad
1\leq\operatorname{\mathrm{intdim}}(\boldsymbol{A})\leq\operatorname{\mathrm{rank}}(\boldsymbol{A})\leq
n.$

If, additionally, $\boldsymbol{A}$ is a diagonal matrix, then

$\operatorname{\mathrm{intdim}}(\boldsymbol{A})=\frac{\sum_{i=1}^{n}a_{ii}}{\max_{1\leq
i\leq n}a_{ii}}.$

The columns of
$\boldsymbol{A}=\begin{bmatrix}\boldsymbol{a}_{1}&\cdots&\boldsymbol{a}_{n}\end{bmatrix}\in\mathbb{R}^{m\times
n}$ are $\boldsymbol{a}_{j}\in\mathbb{R}^{m}$, $1\leq j\leq n$, and the
columns of the identity
$\boldsymbol{I}=\begin{bmatrix}\boldsymbol{e}_{1}&\cdots&\boldsymbol{e}_{n}\end{bmatrix}\in\mathbb{R}^{n\times
n}$ are $\boldsymbol{e}_{j}\in\mathbb{R}^{n}$. The transpose of
$\boldsymbol{A}$ is $\boldsymbol{A}^{\top}$.

###  2.2 Concentration inequalities

We rely on two scalar and two matrix concentration inequalities.

Markov’s inequality [15, Section 3.1] bounds the the probability that a random
variable exceeds a constant.

######  Theorem 2.1 (Markov’s inequality).

If $Z$ is a non-negative random variable, then for $t>0$

Problem 1

Hoeffding’s inequality for general bounded random variables [20, Theorem
2.2.6] bounds the probability that a sum of scalar random variables exceeds
its mean.

######  Theorem 2.2 (Scalar Hoeffding inequality).

Let $Z_{1},\ldots,Z_{N}$ be independent random variables, bounded by
$m_{k}\leq Z_{k}\leq M_{k}$, $1\leq k\leq N$, with sum
$Z\equiv\sum_{k=1}^{N}{Z_{k}}$. Then for $t>0$

$\operatorname{\mathbb{P}}\left[|Z-\operatorname{\mathbb{E}}[Z]|\geq
t\right]\leq
2\exp\left(\frac{-2t^{2}}{\sum_{k=1}^{N}{(M_{k}-m_{k})^{2}}}\right).$

Next are two bounds for sums of independent symmetric matrix-valued random
variables. The first is a matrix Bernstein concentration inequality [18,
Theorems 7.3.1 and 7.7.1] for sums of independent, symmetric, bounded, zero-
mean random matrices.

######  Theorem 2.3 (Matrix Bernstein inequality).

Let $\boldsymbol{S}_{1},\dots,\boldsymbol{S}_{N}\in\mathbb{R}^{n\times n}$ be
independent symmetric random matrices with

$\mathbb{E}[\boldsymbol{S}_{k}]=\boldsymbol{0},\qquad\|\boldsymbol{S}_{k}\|_{2}\leq
L\qquad 1\leq k\leq N.$

Let the sum $\boldsymbol{S}\equiv\sum_{k=1}^{N}\boldsymbol{S}_{k}$ have a
matrix-valued variance that is majorized by
$\boldsymbol{V}\in\mathbb{R}^{n\times n}$,

$\boldsymbol{V}\succeq\operatorname{\mathrm{Var}}(\boldsymbol{S})=\mathbb{E}[\boldsymbol{S}^{2}]=\sum_{k=1}^{N}\mathbb{E}[\boldsymbol{S}_{k}^{2}].$

Abbreviate $\nu\equiv\|\boldsymbol{V}\|_{2}$ and
$d\equiv\operatorname{\mathrm{intdim}}(\boldsymbol{V})$. Then for $t>0$

(3) $\operatorname{\mathbb{P}}\left[\|\boldsymbol{S}\|_{2}\geq t\right]\leq
8d\,\exp\left(\frac{-t^{2}}{2(\nu+Lt/3)}\right).$

######  Proof 2.4.

In [18, Theorems 7.3.1 and 7.7.1] it is shown that (3) holds, provided
$t\geq\sqrt{\nu}+\frac{L}{3}$. We show that (3) always holds and the lower
bound on $t$ is not necessary. To see this, note that

$\frac{-t^{2}}{2(\nu+Lt/3)}$

decreases monotonically as $t$ increases. Therefore we can bound it from below
as long as $t<\sqrt{\nu}+L/3$, by

$\frac{-t^{2}}{2(\nu+Lt/3)}>\frac{-(\sqrt{\nu}+L/3)^{2}}{2\left(\nu+(L/3)(\sqrt{\nu}+L/3)\right)}\geq-\frac{2}{3}.$

The second inequality comes from setting $x=\frac{L}{3\sqrt{\nu}}$ and noting
that $f(x)=-\frac{(1+x)^{2}}{2(1+x+x^{2})}$ has a minimum at $\hat{x}=1$ where
$f(\hat{x})=-2/3$, and $\sqrt{\nu}=L/3$. Substituting the lower bound into
Theorem 2.3 gives

$8d\,\exp\left(\frac{-t^{2}}{2(\nu+Lt/3)}\right)\geq 8\exp(-2/3)>4.$

But now (3) holds trivially since

$\operatorname{\mathbb{P}}\left[\|\boldsymbol{S}\|_{2}\geq t\right]\leq
1<4<8d\,\exp\left(\frac{-t^{2}}{2(\nu+Lt/3)}\right).$

The second matrix concentration inequality [5, Theorem 3.2] bounds the mean of
the squared norm of the sum of symmetric random matrices.

######  Theorem 2.5.

Let $\boldsymbol{S}_{1},\dots,\boldsymbol{S}_{N}\in\mathbb{R}^{n\times n}$
with $n\geq 3$ be independent symmetric random matrices with zero mean. Then

$\mathbb{E}\left[\left\|\sum_{k=1}^{N}\boldsymbol{S}_{k}\right\|_{2}^{2}\right]^{1/2}\leq\sqrt{2e\ln{n}}\>\left\|\left(\sum_{k=1}^{N}\mathbb{E}[\boldsymbol{S}_{k}^{2}]\right)^{1/2}\right\|_{2}+4e\ln{n}\>\left(\mathbb{E}\left[\max_{1\leq
k\leq N}{\|\boldsymbol{S}_{k}\|_{2}^{2}}\right]\right)^{1/2}.$

###  2.3 Setup for the analysis

Our Monte Carlo estimators compute the diagonal elements of a symmetric matrix
$\boldsymbol{A}\in\mathbb{R}^{n\times n}$ by means of matrix vector products
with $\boldsymbol{A}$. It samples $N$ independent random vectors
$\boldsymbol{w}_{k}\in\mathbb{R}^{n}$, and approximates the vector of diagonal
elements $\operatorname{\mathrm{diag}}(\boldsymbol{A})\in\mathbb{R}^{n}$ by
the mean

$\operatorname{\mathrm{diag}}(\widehat{\boldsymbol{A})}=\frac{1}{N}\sum_{k=1}^{N}{\left((\boldsymbol{Aw}_{k})\circ\boldsymbol{w}_{k}\right)}\in\mathbb{R}^{n}\qquad\text{where}\quad\widehat{\boldsymbol{A}}\equiv\frac{1}{N}\sum_{k=1}^{N}{\boldsymbol{Aw}_{k}\boldsymbol{w}_{k}^{\top}}\in\mathbb{R}^{n\times
n}.$

To see this, apply (1) and (2) to

$\displaystyle\mathcal{D}\,(\widehat{\boldsymbol{A}})=\boldsymbol{I}\circ\widehat{\boldsymbol{A}}$
$\displaystyle=\frac{1}{N}\sum_{k=1}^{N}{\operatorname{\mathrm{diag}}(\boldsymbol{Aw}_{k})\boldsymbol{I}\operatorname{\mathrm{diag}}(\boldsymbol{w}_{k})}=\frac{1}{N}\sum_{k=1}^{N}{\operatorname{\mathrm{diag}}((\boldsymbol{Aw}_{k})\circ\boldsymbol{w}_{k})}$
$\displaystyle=\operatorname{\mathrm{diag}}\left(\frac{1}{N}\sum_{k=1}^{N}{((\boldsymbol{Aw}_{k})\circ\boldsymbol{w}_{k})}\right)=\operatorname{\mathrm{diag}}\left(\operatorname{\mathrm{diag}}(\widehat{\boldsymbol{A}})\right)\in\mathbb{R}^{n\times
n}.$

Alternately, the diagonal elements of the estimators can be expressed as

$\displaystyle\widehat{\boldsymbol{A}}_{ii}$
$\displaystyle=\frac{1}{N}\sum_{k=1}^{N}{(\boldsymbol{A}\boldsymbol{w}_{k}\boldsymbol{w}_{k}^{\top})_{ii}}=\frac{1}{N}\sum_{k=1}^{N}{(\boldsymbol{A}\boldsymbol{w}_{k})_{i}\,(\boldsymbol{w}_{k}^{\top})_{i}}$
$\displaystyle=\frac{1}{N}\sum_{k=1}^{N}\left((\boldsymbol{Aw}_{k})\circ\boldsymbol{w}_{k}\right)_{i},\qquad
1\leq i\leq n.$

We measure the cost of a diagonal estimator by the number $N$ of samples. To
assess the accuracy, we introduce a relative error in the form of normwise and
componentwise $(\epsilon,\delta)$ estimators, which extend the notion of
$(\epsilon,\delta)$ trace estimator from [2, 16].

######  Definition 2.6 (Normwise $(\epsilon,\delta)$ diagonal estimator).

Let $\boldsymbol{A}\in{\mathbb{R}}^{n\times n}$ be symmetric. Given user-
specified parameters $0<\epsilon,\delta<1$, we say that
$\mathcal{D}\,(\widehat{\boldsymbol{A}})$ is a normwise $(\epsilon,\delta)$
estimator for the diagonal elements of $\boldsymbol{A}$, if

$\|\mathcal{D}\,(\boldsymbol{A})-\mathcal{D}\,(\widehat{\boldsymbol{A}})\|_{2}\leq\epsilon\,\|\mathcal{D}\,(\boldsymbol{A})\|_{2}$

holds with probability at least $1-\delta$.

In other words, for a user-specified failure probability $\delta>0$ and
tolerance $\epsilon>0$, the normwise relative error of the diagonal estimator
is, with probability at most $1-\delta$, at most $\epsilon$. The two-norm in
Definition 2.6 can be replaced by any matrix-$p$ norm, because the $p$-norm of
a diagonal matrix $\boldsymbol{D}\in\mathbb{R}^{n\times n}$ is
$\|\boldsymbol{D}\|_{p}=\max_{1\leq i\leq n}{|d_{ii}|}$ for $p\geq 1$. Next we
define a componentwise $(\epsilon,\delta)$ estimator.

######  Definition 2.7 (Componentwise $(\epsilon,\delta)$ diagonal estimator).

Let $\boldsymbol{A}\in{\mathbb{R}}^{n\times n}$ be symmetric. Given user-
specified parameters $0<\epsilon,\delta<1$ and diagonal element $a_{ii}\neq 0$
of $\boldsymbol{A}$, we say that
$\widehat{a}_{ii}=(\widehat{\boldsymbol{A}})_{ii}$ is a componentwise
$(\epsilon,\delta)$ estimator for $a_{ii}$, if

$|\widehat{a}_{ii}-a_{ii}|\leq\epsilon\,|a_{ii}|$

holds with probability at least $1-\delta$.

##  3 Normwise bounds for Rademacher random vectors

We present normwise bounds for Monte Carlo estimators based on standard
(section 3.1) and on sparse Rademacher vectors (section 3.2).

###  3.1 Standard Rademacher vectors

After defining Rademacher vectors (De-finition 3.1) and discussing their
properties (Remarks 3.2 and 3.3), we present a normwise absolute error bound
(Theorem 3.4), and a bound on the minimal sampling amount that makes the
Rademacher Monte Carlo estimator a normwise $(\epsilon,\delta)$ diagonal
estimator (Corollary 3.6).

######  Definition 3.1.

A Rademacher random variable takes on the values $\pm 1$ with equal
probability 1/2. A Rademacher vector is a random vector whose elements are
independent Rademacher random variables.

Standard Rademacher vectors have the advantage of cheap matrix vector products
and immediately recovering diagonal matrices.

######  Remark 3.2.

The elements $w_{j}$ of a Rademacher vector $\boldsymbol{w}$ have the
following properties:

  1. 1.

Zero mean: $\mathbb{E}[w_{j}]=0$

  2. 2.

Constant square: $w_{j}^{2}=1$

  3. 3.

Independence: $\mathbb{E}[w_{j}w_{i}]=0$ for $i\neq j$.

######  Remark 3.3.

Standard Rademacher vectors recover a diagonal matrix with a single sample,
$N=1$.

To see this, let
$\boldsymbol{A}=\mathcal{D}\,(\boldsymbol{A})\in\mathbb{R}^{n\times n}$ be
diagonal, and $\boldsymbol{w}\in\mathbb{R}^{n}$ a Rademacher vector. Remark
3.2 implies that $\boldsymbol{Aw}\boldsymbol{w}^{\top}\in\mathbb{R}^{n\times
n}$ has diagonal elements $a_{ii}w_{i}^{2}=a_{ii}$, $1\leq i\leq n$.

As a consequence, we can focus the analysis of standard Rademacher-based
estimators on non-diagonal matrices. The results below are special cases of
those for sparse Rademacher vectors in section 3.2.

######  Theorem 3.4.

Let $\boldsymbol{A}\in\mathbb{R}^{n\times n}$ be non-diagonal symmetric, and

$K_{1}\equiv\|\mathcal{D}\,(\boldsymbol{A}^{2})-\mathcal{D}\,(\boldsymbol{A})^{2}\|_{2},\qquad
K_{2}\equiv\|\boldsymbol{A}-\mathcal{D}\,(\boldsymbol{A})\|_{\infty},\qquad
d\equiv(\|\boldsymbol{A}\|_{F}^{2}-\|\mathcal{D}\,(\boldsymbol{A})\|_{F}^{2})/K_{1}.$

If
$\widehat{\boldsymbol{A}}\equiv\frac{1}{N}\sum_{k=1}^{N}\boldsymbol{Aw}_{k}\boldsymbol{w}_{k}^{\top}$
is a Monte Carlo estimator with independent Rademacher vectors
$\boldsymbol{w}_{k}\in\mathbb{R}^{n}$, $1\leq k\leq N$, then the probability
that the absolute error exceeds $t>0$ is at most

$\operatorname{\mathbb{P}}\left[\|\mathcal{D}\,(\boldsymbol{A})-\mathcal{D}\,(\widehat{\boldsymbol{A}})\|_{2}\geq
t\right]\leq 8d\,\exp\left(\frac{-Nt^{2}}{2(K_{1}+tK_{2}/3)}\right).$

######  Proof 3.5.

This is the special case $s=1$ of Theorem 3.10.

The constants $K_{1}$ and $K_{2}$ represent the absolute deviation of
$\boldsymbol{A}$ from diagonality, and more specifically the degree of
diagonal dominance of $\boldsymbol{A}$ in the absolute sense. Theorem 3.4
implies that the Rademacher estimator has a small absolute error when applied
to strongly diagonally dominant matrices. In other words, the normwise
absolute error in the Rademacher estimator decreases with increasing diagonal
dominance of $\boldsymbol{A}$ in the absolute sense.

We determine the least sampling amount required for the Monte Carlo estimator
with Rademacher vectors to be a normwise $(\epsilon,\delta)$ diagonal
estimator.

######  Corollary 3.6.

Let $\boldsymbol{A}\in\mathbb{R}^{n\times n}$ be non-diagonal symmetric. Let

$\displaystyle K_{1}$
$\displaystyle\equiv\|\mathcal{D}\,(\boldsymbol{A}^{2})-\mathcal{D}\,(\boldsymbol{A})^{2}\|_{2},$
$\displaystyle\Delta_{1}$
$\displaystyle\equiv\frac{K_{1}}{\|\mathcal{D}\,(\boldsymbol{A})\|_{2}^{2}},\qquad\Delta_{2}\equiv\frac{\|\boldsymbol{A}-\mathcal{D}\,(\boldsymbol{A})\|_{\infty}}{\|\mathcal{D}\,(\boldsymbol{A})\|_{\infty}},\qquad
d\equiv\frac{\|\boldsymbol{A}\|_{F}^{2}-\|\mathcal{D}\,(\boldsymbol{A})\|_{F}^{2}}{K_{1}},$

and let
$\widehat{\boldsymbol{A}}\equiv\frac{1}{N}\sum_{k=1}^{N}\boldsymbol{Aw}_{k}\boldsymbol{w}_{k}^{\top}$
a Monte Carlo estimator with independent Rademacher vectors
$\boldsymbol{w}_{k}\in\mathbb{R}^{n}$, $1\leq k\leq N$. Pick $\epsilon>0$. For
any $0<\delta<1$, if the sampling amount is at least

(4)
$N\geq\frac{\Delta_{2}}{3\epsilon^{2}}\left(2\epsilon+6\,\frac{\Delta_{1}}{\Delta_{2}}\right)\ln{(8d/\delta)},$

then
$\|\mathcal{D}\,(\boldsymbol{A})-\mathcal{D}\,(\widehat{\boldsymbol{A}})\|_{2}\leq\epsilon\,\|\mathcal{D}\,(\boldsymbol{A})\|_{2}$
holds with probability at least $1-\delta$.

######  Proof 3.7.

This is the special case $s=1$ of Corollary 3.13.

The constants $\Delta_{1}$ and $\Delta_{2}$ in Corollary 3.6 represent the
respective relative counterparts of $K_{1}$ and $K_{2}$ in Theorem 3.4: they
represent the relative deviation of $\boldsymbol{A}$ from diagonality, and
more specifically the degree of diagonal dominance of $\boldsymbol{A}$ in the
relative sense. Corollary 3.6 implies that if $\boldsymbol{A}$ is strongly
diagonally dominant in the relative sense, then a small sampling amount
suffices to make the Rademacher estimator a normwise $(\epsilon,\delta)$
diagonal estimator. As with many randomized sampling algorithms, the lower
bound for $N$ is proportional to $1/\epsilon^{2}$.

###  3.2 Sparse Rademacher vectors

For Rademacher vectors that are parameterized in terms of sparsity (Definition
3.8), we derive a normwise absolute error bound (Theorem 3.10), followed by
the minimal sampling amount that makes the sparse Rademacher estimator a
normwise $(\epsilon,\delta)$ diagonal estimator (Corollary 3.13).

The random vectors in [1] have elements that assume values from the discrete
distribution $\\{-\sqrt{3}\,,0\,,\sqrt{3}\\}$ with respective probability
$\left\\{\frac{1}{6},\,\frac{2}{3},\,\frac{1}{6}\right\\}$. This concept was
extended in [14, (2)] to Rademacher vectors that are parameterized in terms of
a sparsity parameter $s$.

######  Definition 3.8.

A sparse Rademacher random variable with parameter $s\geq 1$ takes the values
$\\{-\sqrt{s},\,0,\,\sqrt{s}\\}$ with probability
$\\{\frac{1}{2s},\,1-\frac{1}{s},\,\frac{1}{2s}\\}$ respectively.

A Sparse Rademacher vector is a random vector whose elements are independent
sparse Rademacher random variables.

The properties of sparse Rademacher vectors are almost the same as those of
the original Rademacher vectors in Remark 3.2.

######  Remark 3.9.

The elements of a sparse Rademacher vector $\boldsymbol{w}$ with parameter
$s\geq 1$ have the following properties

  1. 1.

Zero mean: $\operatorname{\mathbb{E}}[w_{j}]=0$

  2. 2.

Unit variance $\operatorname{\mathbb{E}}[w_{j}^{2}]=1$

  3. 3.

Independence: For $i\neq j$, and integer $\ell\geq 1$

$\mathbb{E}[w_{i}^{2}w_{j}^{2}]=1,\qquad\mathbb{E}[w_{i}^{\ell}w_{j}]=\mathbb{E}[w_{i}w_{j}^{\ell}]=0.$

The case $s=1$ corresponds to the original Rademacher vectors (Definition
3.1), while $s=3$ corresponds to the choice in [1].

Below is the extension of Theorem 3.4 to sparse Rademacher vectors with
integer parameters $s$.

######  Theorem 3.10.

Let $\boldsymbol{A}\in\mathbb{R}^{n\times n}$ be non-diagonal symmetric, and
and let

$\displaystyle K_{1}(s)$
$\displaystyle\equiv\|\mathcal{D}\,(\boldsymbol{A}^{2})+(s-2)\mathcal{D}\,(\boldsymbol{A})^{2}\|_{2},\qquad
K_{2}(s)\equiv\|s\boldsymbol{A}-\mathcal{D}\,(\boldsymbol{A})\|_{\infty},$
$\displaystyle d(s)$
$\displaystyle\equiv\frac{1}{K_{1}(s)}\|\boldsymbol{A}\|_{F}^{2}+(s-2)\|\mathcal{D}\,(\boldsymbol{A})\|_{F}^{2}.$

Let
$\widehat{\boldsymbol{A}}\equiv\frac{1}{N}\sum_{j=1}^{N}{\boldsymbol{Aw}_{k}\boldsymbol{w}_{k}^{\top}}$
be a Monte Carlo estimator with independent sparse Rademacher random vectors
$\boldsymbol{w}_{k}\in\mathbb{R}^{n}$, with integer parameter $s\geq 1$. Then
the probability that the absolute error exceeds $t>0$ is at most

$\operatorname{\mathbb{P}}\left[\|\mathcal{D}\,(\boldsymbol{A})-\mathcal{D}\,(\widehat{\boldsymbol{A}})\|_{2}\geq
t\right]\leq 8d(s)\,\exp\left(\frac{-Nt^{2}}{2(K_{1}(s)+tK_{2}(s)/3)}\right).$

######  Proof 3.11.

Define the random diagonal matrices

$\displaystyle\boldsymbol{S}_{k}\equiv\frac{1}{N}\left(\boldsymbol{I}\circ(\boldsymbol{A}\boldsymbol{w}_{k}\boldsymbol{w}_{k}^{\top})-\boldsymbol{I}\circ\boldsymbol{A}\right),\qquad
1\leq k\leq N,$

and their sum

(5)
$\displaystyle\boldsymbol{Z}\equiv\sum_{k=1}^{N}\boldsymbol{S}_{k}=\boldsymbol{I}\circ\widehat{\boldsymbol{A}}-\boldsymbol{I}\circ\boldsymbol{A}=\mathcal{D}\,(\boldsymbol{A})-\mathcal{D}\,(\widehat{\boldsymbol{A}}).$

Before applying Corollary 2.3, we need to verify the assumptions for the
Bernstein inequality.

#### 1\. Expectation

With (2), the first summand in $\boldsymbol{S}_{k}$ equals

$\displaystyle\boldsymbol{I}\circ(\boldsymbol{A}\boldsymbol{w}_{k}\boldsymbol{w}_{k}^{\top})=\operatorname{\mathrm{diag}}(\boldsymbol{Aw}_{k})\boldsymbol{I}\operatorname{\mathrm{diag}}(\boldsymbol{w}_{k})=\operatorname{\mathrm{diag}}(\boldsymbol{Aw}_{k})\operatorname{\mathrm{diag}}(\boldsymbol{w}_{k}),$

and with (1) this gives

(6)
$\displaystyle\boldsymbol{S}_{k}=\frac{1}{N}\left(\operatorname{\mathrm{diag}}(\boldsymbol{Aw}_{k})\operatorname{\mathrm{diag}}(\boldsymbol{w}_{k})-\mathcal{D}\,(\boldsymbol{A})\right),\qquad
1\leq k\leq N.$

The linearity of expectation implies

$\mathbb{E}[\boldsymbol{S}_{k}]=\frac{1}{N}\left(\underbrace{\mathbb{E}[\operatorname{\mathrm{diag}}(\boldsymbol{Aw}_{k})\operatorname{\mathrm{diag}}(\boldsymbol{w}_{k})]}_{\boldsymbol{T}}-\mathcal{D}\,(\boldsymbol{A})\right),\qquad
1\leq k\leq N.$

Remark 3.9 implies that element $(i,j)$ of
$\boldsymbol{T}\equiv\mathbb{E}[\operatorname{\mathrm{diag}}(\boldsymbol{Aw}_{k})\operatorname{\mathrm{diag}}(\boldsymbol{w}_{k})]$
equals

(7) $\displaystyle
t_{ij}=\mathbb{E}[(\boldsymbol{Aw}_{k})_{i}\,\boldsymbol{I}_{ij}\,(\boldsymbol{w}_{k})_{j}]=\boldsymbol{I}_{ij}a_{ij}\qquad
1\leq i\leq n,\ 1\leq j\leq n.$

Hence $\boldsymbol{T}=\boldsymbol{I}\circ\boldsymbol{A}$, and
$\mathbb{E}[\boldsymbol{S}_{k}]=\frac{1}{N}(\boldsymbol{T}-\boldsymbol{I}\circ\boldsymbol{A})=\boldsymbol{0}$,
$1\leq k\leq N$. Therefore, the random matrices have zero mean,
$\mathbb{E}[\boldsymbol{S}_{k}]=\boldsymbol{0}$; and so has their sum,
$\mathbb{E}[\boldsymbol{Z}]=\boldsymbol{0}$.

#### 2\. Boundedness

From (6) follows that the diagonal matrices $\boldsymbol{S}_{k}$ have diagonal
elements

$\displaystyle(\boldsymbol{S}_{k})_{ii}$
$\displaystyle=\frac{1}{N}\left((\boldsymbol{A}\boldsymbol{w}_{k})_{i}(\boldsymbol{w}_{k})_{i}-a_{ii}\right)=\frac{1}{N}\left((\boldsymbol{w}_{k})_{i}\,\sum_{j=1}^{n}{a_{ij}(\boldsymbol{w}_{k})_{j}}-a_{ii}\right)$
$\displaystyle=\frac{1}{N}\left(a_{ii}((\boldsymbol{w}_{k})_{i}^{2}-1)+(\boldsymbol{w}_{k})_{i}\,\sum_{j\neq
i}{a_{ij}(\boldsymbol{w}_{k})_{j}}\right),\qquad 1\leq i\leq n,\quad 1\leq
k\leq N.$

Since $s\geq 1$ is an integer,
$|((\boldsymbol{w}_{k})_{ii}^{2}-1)a_{ii}|\leq(s-1)|a_{ii}|$ holds for any
value of $(\boldsymbol{w}_{k})_{i}$, and

$\displaystyle|(\boldsymbol{S}_{k})_{ii}|\leq\frac{1}{N}\left((s-1)|a_{ii}|+s\sum_{j\neq
i}{|a_{ij}|}\right),\qquad 1\leq i\leq n,\quad 1\leq k\leq N.$

Since the matrix infinity norm is absolute, with
$\|\boldsymbol{S}_{k}\|_{\infty}=\|\,|\boldsymbol{S}_{k}|\,\|_{\infty}$ we can
bound the two-norms of the diagonal matrices $\boldsymbol{S}_{k}$ by

$\displaystyle\|\boldsymbol{S}_{k}\|_{2}$ $\displaystyle\leq\max_{1\leq i\leq
n}{\frac{1}{N}\left((s-1)|a_{ii}|+s\sum_{j\neq
i}{|a_{ij}|}\right)}=\max_{1\leq i\leq
n}{\frac{1}{N}\left(-|a_{ii}|+s\sum_{j=1}^{n}{|a_{ij}|}\right)}$
$\displaystyle=\frac{1}{N}\|s\boldsymbol{A}-\mathcal{D}\,(\boldsymbol{A})\|_{\infty}=\frac{K_{2}(s)}{N},\qquad
1\leq k\leq N.$

Set $L(s)\equiv K_{2}(s)/N$ where $K_{2}(s)>0$ since $\boldsymbol{A}$ is not
diagonal.

#### 3\. Variance

In (6) abbreviate
$\boldsymbol{D}_{k}\equiv\operatorname{\mathrm{diag}}(\boldsymbol{Aw}_{k})$,
$\boldsymbol{W}_{k}\equiv\operatorname{\mathrm{diag}}(\boldsymbol{w}_{k})$,
and $\boldsymbol{D}\equiv\mathcal{D}\,(\boldsymbol{A})$, and write the
summands in
$\operatorname{\mathrm{Var}}[\boldsymbol{Z}]=\sum_{k=1}^{N}{\mathbb{E}[\boldsymbol{S}_{k}^{2}]}$
as

$\displaystyle\boldsymbol{S}_{k}^{2}=\frac{1}{N^{2}}\left(\boldsymbol{D}_{k}\boldsymbol{W}_{k}-\boldsymbol{D}\right)^{2}=\frac{1}{N^{2}}\left((\boldsymbol{D}_{k}\boldsymbol{W}_{k})^{2}-2\boldsymbol{D}\boldsymbol{D}_{k}\boldsymbol{W}_{k}+\boldsymbol{D}^{2}\right),\qquad
1\leq k\leq N,$

taking advantage of the fact that diagonal matrices commute. The expectation
is

(8)
$\displaystyle\mathbb{E}\left[\boldsymbol{S}_{k}^{2}\right]=\frac{1}{N^{2}}\left(\mathbb{E}\left[(\boldsymbol{D}_{k}\boldsymbol{W}_{k})^{2}\right]-2\boldsymbol{D}\,\mathbb{E}\left[\boldsymbol{D}_{k}\boldsymbol{W}_{k}\right]+\boldsymbol{D}^{2}\right),\qquad
1\leq k\leq N.$

To determine $\mathbb{E}[\boldsymbol{D}_{k}\boldsymbol{W}_{k}]$, look at the
individual diagonal elements and abbreviate
$\boldsymbol{u}=\boldsymbol{w}_{k}$,

(9)
$\displaystyle(\boldsymbol{D}_{k}\boldsymbol{W}_{k})_{ii}={u}_{i}\,\boldsymbol{e}_{i}^{\top}\boldsymbol{Au}={u}_{i}\,\sum_{j=1}^{n}{a_{ij}{u}_{j}}=a_{ii}{u}_{i}^{2}+\sum_{j\neq
i}{a_{ij}{u}_{i}{u}_{j}},\qquad 1\leq i\leq n.$

Remark 3.9 implies for the expectation

$\displaystyle\mathbb{E}\left[(\boldsymbol{D}_{k}\boldsymbol{W}_{k})_{ii}\right]=a_{ii}\,\mathbb{E}\left[{u}_{i}^{2}\right]+\sum_{j\neq
i}{a_{ij}\,\mathbb{E}\left[{u}_{i}{u}_{j}\right]}=a_{ii},\qquad 1\leq i\leq
n.$

Thus
$\mathbb{E}\left[\boldsymbol{D}_{k}\boldsymbol{W}_{k}\right]=\boldsymbol{D}$.
Insert this into (8)

(10)
$\displaystyle\mathbb{E}\left[\boldsymbol{S}_{k}^{2}\right]=\frac{1}{N^{2}}\left(\mathbb{E}\left[(\boldsymbol{D}_{k}\boldsymbol{W}_{k})^{2}\right]-\boldsymbol{D}^{2}\right),\qquad
1\leq k\leq N.$

To determine
$\mathbb{E}\left[(\boldsymbol{D}_{k}\boldsymbol{W}_{k})^{2}\right]$, start
from the individual diagonal elements in (9)

$\displaystyle(\boldsymbol{D}_{k}\boldsymbol{W}_{k})_{ii}^{2}$
$\displaystyle=\left(a_{ii}{u}_{i}^{2}+{u}_{i}\,\sum_{j\neq
i}{a_{ij}{u}_{j}}\right)^{2}$
$\displaystyle=\underbrace{a_{ii}^{2}{u}_{i}^{4}}_{\alpha}+\underbrace{2a_{ii}{u}_{i}^{3}\,\sum_{j\neq
i}{a_{ij}{u}_{j}}}_{\beta}+\underbrace{{u}_{i}^{2}\left(\sum_{j\neq
i}{a_{ij}{u}_{j}}\right)^{2}}_{\gamma}\qquad 1\leq i\leq n.$

In the expectation

(11)
$\displaystyle\mathbb{E}\left[(\boldsymbol{D}_{k}\boldsymbol{W}_{k})_{ii}^{2}\right]=\mathbb{E}[\alpha]+\mathbb{E}[\beta]+\mathbb{E}[\gamma]$

we inspect each summand in turn.

Definition 3.8 implies for the fourth moment

(12)
$\displaystyle\mathbb{E}[\alpha]=a_{ii}^{2}\mathbb{E}\left[{u}_{i}^{4}\right]=a_{ii}^{2}\left(\frac{1}{2s}\left(-\sqrt{s}\right)^{4}+\left(1-\frac{1}{s}\right)0+\frac{1}{2s}\left(\sqrt{s}\right)^{4}\right)=a_{ii}^{2}s.$

Since ${u}_{i}$ is independent from ${u}_{j}$ for $j\neq i$ Remark 3.9 implies

(13) $\displaystyle\mathbb{E}[\beta]=2a_{ii}\sum_{j\neq
i}{a_{ij}\,\underbrace{\mathbb{E}\left[{u}_{i}^{3}{u}_{j}\right]}_{=0}}=0.$

Collect the independent elements of $\boldsymbol{u}$,

$\displaystyle\gamma$ $\displaystyle={u}_{i}^{2}\left(\sum_{j\neq
i}{a_{ij}{u}_{j}}\right)\left(\sum_{\ell\neq
i}{a_{i\ell}{u}_{\ell}}\right)={u}_{i}^{2}\left(\sum_{j\neq
i}{a_{ij}{u}_{j}\left(a_{ij}{u}_{j}+\sum_{\ell\neq
i,j}{a_{i\ell}{u}_{\ell}}\right)}\right)$ $\displaystyle=\sum_{j\neq
i}{a_{ij}^{2}{u}_{j}^{2}{u}_{i}^{2}}+\sum_{j\neq i}{\sum_{\ell\neq
i,j}{a_{ij}a_{i\ell}{u}_{j}{u}_{\ell}{u}_{i}^{2}}}$

and apply Remark 3.9,

(14) $\displaystyle\mathbb{E}[\gamma]$ $\displaystyle=\sum_{j\neq
i}{a_{ij}^{2}\,\underbrace{\mathbb{E}\left[{u}_{j}^{2}{u}_{i}^{2}\right]}_{=1}}+\sum_{j\neq
i}{\sum_{\ell\neq
i,j}{a_{ij}a_{i\ell}\,\underbrace{\mathbb{E}\left[{u}_{j}{u}_{\ell}{u}_{i}^{2}\right]}_{=0}}}=\sum_{j\neq
i}{a_{ij}^{2}}$

Substitute (12), (13) and (14) into (11). Then Remark 3.9 and the symmetry of
$\boldsymbol{A}$ imply

$\displaystyle\mathbb{E}\left[(\boldsymbol{D}_{k}\boldsymbol{W}_{k})_{ii}^{2}\right]=sa_{ii}^{2}+\sum_{j\neq
i}{a_{ij}^{2}}=\sum_{j=1}^{n}{a_{ij}^{2}}+(s-1)a_{ii}^{2}=\|\boldsymbol{a}_{i}\|_{2}^{2}+(s-1)|a_{ii}|^{2}.$

With
$\mathcal{D}\,(\boldsymbol{A}^{2})=\operatorname{\mathrm{diag}}\begin{pmatrix}(\|\boldsymbol{a}_{1}\|_{2}^{2}&\cdots&\|\boldsymbol{a}_{n}\|_{2}^{2}\end{pmatrix}$
this gives for the whole diagonal matrix

$\displaystyle\mathbb{E}[(\boldsymbol{D}_{k}\boldsymbol{W}_{k})^{2}]=\mathcal{D}\,(\boldsymbol{A}^{2})+(s-1)\boldsymbol{D}^{2},\qquad
1\leq k\leq N.$

At last, substitute this into the variance (10)

$\displaystyle\mathbb{E}\left[\boldsymbol{S}_{k}^{2}\right]$
$\displaystyle=\frac{1}{N^{2}}\left(\mathcal{D}\,(\boldsymbol{A}^{2})+(s-1)\boldsymbol{D}^{2}-\boldsymbol{D}^{2}\right)$
$\displaystyle=\frac{1}{N^{2}}\left(\mathcal{D}\,(\boldsymbol{A}^{2})+(s-2)\boldsymbol{D}^{2}\right),\qquad
1\leq k\leq N.$

Sum up the individual variances,

$\displaystyle\boldsymbol{V}(s)\equiv\operatorname{\mathrm{Var}}[\boldsymbol{Z}]$
$\displaystyle=\sum_{k=1}^{N}{\mathbb{E}[\boldsymbol{S}_{k}^{2}]}=\frac{1}{N^{2}}\,\sum_{k=1}^{N}{\left(\mathcal{D}\,(\boldsymbol{A}^{2})+(s-2)\boldsymbol{D}^{2}\right)}$
$\displaystyle=\frac{1}{N}\left(\mathcal{D}\,(\boldsymbol{A}^{2})+(s-2)\boldsymbol{D}^{2}\right).$

Since $\boldsymbol{V}(s)$ is diagonal, its norm is

$\nu(s)\equiv\|\boldsymbol{V}(s)\|_{2}=\frac{1}{N}\|\mathcal{D}\,(\boldsymbol{A}^{2})+(s-2)\boldsymbol{D}^{2}\|_{2}=\frac{K_{1}(s)}{N},$

where $K_{1}(s)>0$ since $\boldsymbol{A}$ is not diagonal. The intrinsic
dimension of $\boldsymbol{V}(s)$ is

$\displaystyle d(s)$
$\displaystyle\equiv\operatorname{\mathrm{intdim}}(\boldsymbol{V})=\frac{\operatorname{\mathrm{trace}}(\mathcal{D}\,(\boldsymbol{A}^{2})+(s-2)\boldsymbol{D}^{2})}{N\nu(s)}=\frac{1}{K_{1}(s)}\sum_{j=1}^{n}{\left(\|\boldsymbol{a}_{j}\|_{2}^{2}+(s-2)|a_{jj}|^{2}\right)}$
$\displaystyle=\frac{1}{K_{1}(s)}\|\boldsymbol{A}\|_{F}^{2}+(s-2)\|\mathcal{D}\,(\boldsymbol{A})\|_{F}^{2}.$

#### 4\. Apply Theorem 2.3

Substituting $L(s)=K_{2}(s)/N$ and $\nu(s)=K_{1}(s)/N$ into Theorem 2.3 and
remembering that the sum $\boldsymbol{Z}$ in (5) has zero mean gives

$\operatorname{\mathbb{P}}\left[\|\boldsymbol{Z}\|_{2}\geq
t\right]=\operatorname{\mathbb{P}}\left[\|\mathcal{D}\,(\boldsymbol{A})-\mathcal{D}\,(\widehat{\boldsymbol{A}})\|_{2}\geq
t\right]\leq 8d(s)\,\exp\left(\frac{-Nt^{2}}{2(K_{1}(s)+tK_{2}(s)/3)}\right).$

In the special case $s=1$ of standard Rademacher vectors, Theorem 3.10 reduces
to Theorem 3.4. However, sparse Rademacher Monte Carlo estimators with $s>1$
do, in general, not recover a diagonal matrix with a single sample, $N=1$.

As $s$ increases, so do the constants $K_{1}(s)$ and $K_{2}(s)$, and the upper
bound on
$\operatorname{\mathbb{P}}\left[\|\mathcal{D}\,(\boldsymbol{A})-\mathcal{D}\,(\widehat{\boldsymbol{A}})\|_{2}\geq
t\right]$. In other words, the sparser the vectors $\boldsymbol{w}_{k}$, the
less accurate the Monte Carlo estimate
$\mathcal{D}\,(\widehat{\boldsymbol{A}})$.

######  Remark 3.12 (Non-integer sparsity levels).

The restriction to integers $s$ in Theorem 3.10 is relevant only for $1<s<2$.
More generally, Theorem 3.10 holds for $s=1$, and any real number $s\geq 2$.

The extension below of Corollary 3.6 presents the minimal sampling amount that
makes the sparse Rademacher Monte Carlo estimator a normwise
$(\epsilon,\delta)$ diagonal estimator.

######  Corollary 3.13.

Let $\boldsymbol{A}\in\mathbb{R}^{n\times n}$ be non-diagonal symmetric, and
let

$\displaystyle K_{1}(s)$
$\displaystyle\equiv\|\mathcal{D}\,(\boldsymbol{A}^{2})+(s-2)\mathcal{D}\,(\boldsymbol{A})^{2}\|_{2},\qquad\Delta_{1}(s)\equiv\frac{K_{1}(s)}{\|\mathcal{D}\,(\boldsymbol{A})\|_{2}}$
$\displaystyle\Delta_{2}(s)$
$\displaystyle\equiv\frac{\|s\boldsymbol{A}-\mathcal{D}\,(\boldsymbol{A})\|_{\infty}}{\|\mathcal{D}\,(\boldsymbol{A})\|_{\infty}},\qquad
d\equiv\frac{\|\boldsymbol{A}\|_{F}^{2}+(s-2)\|\mathcal{D}\,(\boldsymbol{A})\|_{F}^{2}}{K_{1}(s)}.$

Let
$\widehat{\boldsymbol{A}}\equiv\frac{1}{N}\sum_{j=1}^{N}{\boldsymbol{Aw}_{k}\boldsymbol{w}_{k}^{\top}}$
be a Monte Carlo estimator with independent sparse Rademacher random vectors
$\boldsymbol{w}_{k}\in\mathbb{R}^{n}$, with integer parameter $s\geq 1$. Pick
$\epsilon>0$. For any $0<\delta<1$, if the sampling amount is at least

$N\geq\frac{\Delta_{2}(s)}{3\epsilon^{2}}\left(2\epsilon+6\,\frac{\Delta_{1}(s)}{\Delta_{2}(s)}\right)\ln{(8d(s)/\delta)},$

then
$\|\mathcal{D}\,(\boldsymbol{A})-\mathcal{D}\,(\widehat{\boldsymbol{A}})\|_{2}\leq\epsilon\,\|\mathcal{D}\,(\boldsymbol{A})\|_{2}$
holds with probability at least $1-\delta$.

######  Proof 3.14.

Denote the bound for the failure probability in Theorem 3.10 by

$\delta\equiv
8d(s)\,\exp\left(\frac{-Nt^{2}}{2(K_{1}(s)+tK_{2}(s)/3)}\right),$

and solve it for $t$,

$\displaystyle
t=\frac{K_{2}(s)}{3N}\ln{(8d(s)/\delta)}+\sqrt{\frac{K_{2}(s)^{2}}{9N^{2}}\ln^{2}(8d(s)/\delta)+\frac{2K_{1}(s)}{N}\ln{(8d(s)/\delta)}}.$

We can restate Theorem 3.4 in terms of the failure probability: With
probability at most $\delta$, the normwise absolute error exceeds
$\|\mathcal{D}\,(\boldsymbol{A})-\mathcal{D}\,(\widehat{\boldsymbol{A}})\|_{2}\geq
t$, where

$\displaystyle
t=\frac{K_{2}(s)}{3N}\ln{(8d(s)/\delta)}+\sqrt{\frac{K_{2}(s)^{2}}{9N^{2}}\ln^{2}(8d(s)/\delta)+\frac{2K_{1}(s)}{N}\ln{(8d(s)/\delta)}}.$

And in terms of the success probability: With probability at least $1-\delta$,
the normwise absolute error is bounded above by
$\|\mathcal{D}\,(\boldsymbol{A})-\mathcal{D}\,(\widehat{\boldsymbol{A}})\|_{2}\leq
t$. Converting this absolute error into a relative one requires
$t\leq\epsilon\|\mathcal{D}\,(\boldsymbol{A})\|_{2}$, in other words,

$\frac{t}{\|\mathcal{D}\,(\boldsymbol{A})\|_{2}}=\gamma+\sqrt{\gamma^{2}+\beta\gamma}\leq\epsilon,$

where

$\displaystyle\gamma\equiv\frac{K_{2}(s)\ln{(8d(s)/\delta)}}{3N\|\mathcal{D}\,(\boldsymbol{A})\|_{2}},\qquad\beta\equiv\frac{6K_{1}(s)}{K_{2}(s)\|\mathcal{D}\,(\boldsymbol{A})\|_{2}}$

Solving for $\gamma$ gives $\gamma\leq\frac{\epsilon^{2}}{2\epsilon+\beta}$,
which implies

$\displaystyle N$
$\displaystyle\geq\frac{K_{2}(s)}{3\epsilon^{2}\,\|\mathcal{D}\,(\boldsymbol{A})\|_{2}}\left(2\epsilon+\frac{6K_{1}(s)}{K_{2}(s)\|\mathcal{D}\,(\boldsymbol{A})\|_{2}}\right)\ln{(8d(s)/\delta)}$
$\displaystyle=\frac{1}{3\epsilon^{2}}\,\frac{K_{2}(s)}{\|\mathcal{D}\,(\boldsymbol{A})\|_{2}}\left(2\epsilon+6\,\frac{K_{1}(s)}{\|\mathcal{D}\,(\boldsymbol{A})\|_{2}^{2}}\frac{\|\mathcal{D}\,(\boldsymbol{A})\|_{2}}{K_{2}(s)}\right)\ln{(8d(s)/\delta)}$
$\displaystyle=\frac{\Delta_{2}}{3\epsilon^{2}}\left(2\epsilon+6\,\frac{\Delta_{1}(s)}{\Delta_{2}(s)}\right)\ln{(8d(s)/\delta)},$

where we exploited the norms of diagonal matrices,

$\displaystyle\Delta_{1}(s)$
$\displaystyle=\frac{K_{1}(s)}{\|\mathcal{D}\,(\boldsymbol{A})\|_{2}^{2}},\qquad\Delta_{2}=\frac{\|\boldsymbol{A}-\mathcal{D}\,(\boldsymbol{A})\|_{\infty}}{\|\mathcal{D}\,(\boldsymbol{A})\|_{\infty}}=\frac{K_{2}(s)}{\|\mathcal{D}\,(\boldsymbol{A})\|_{2}}.$

Corollary 3.13 suggests that increasing the sparsity parameter $s$ could on
the one hand lower the computational cost per sample, but on the other hand
increase the sampling amount for the same accuracy.

##  4 Gaussian vectors

We present normwise bounds for random vectors with bounded fourth moment
(section 4.1), and standard Gaussian vectors (section 4.2).

###  4.1 Random vectors with bounded fourth moment

We bound the expectation of the squared absolute error (Theorem 4.1) for Monte
Carlo estimators based on random vectors $\boldsymbol{w}_{k}$ (assumed to have
independent entries with zero mean and variance $1$) with bounded fourth
moment, $1\leq k\leq N$,

(15) $\mathbb{E}\left[\max_{1\leq k\leq
N}\|\boldsymbol{w}_{k}\|_{\infty}^{4}\right]<+\infty.$

These include standard Rademacher (section 3.1) and sparse Rademacher vectors
(section 3.2), as well as standard Gaussian vectors (section 4.2).

######  Theorem 4.1.

Let $\boldsymbol{A}\in\mathbb{R}^{n\times n}$ with $n\geq 3$ be symmetric, and

$\widehat{\boldsymbol{A}}\equiv\frac{1}{N}\sum_{j=1}^{N}{\boldsymbol{Aw}_{k}\boldsymbol{w}_{k}^{\top}}$

be a Monte Carlo estimator with independent random vectors
$\boldsymbol{w}_{k}\in\mathbb{R}^{n}$, $1\leq k\leq N$, that have independent
elements with zero mean and variance 1. If the vectors $\boldsymbol{w}_{k}$
also have a bounded fourth moment (15), then

$\mathbb{E}\left[\|\mathcal{D}\,(\widehat{\boldsymbol{A}})-\mathcal{D}\,(\boldsymbol{A})\|_{2}^{2}\right]^{1/2}\leq\|\boldsymbol{A}\|_{\infty}\,\left(\sqrt{\frac{8e\ln
n}{N}}+\frac{8e\ln n}{N}\right)\left(\mathbb{E}\left[\max_{1\leq k\leq
N}\|\boldsymbol{w}_{k}\|_{\infty}^{4}\right]\right)^{1/2}.$

######  Proof 4.2.

We make use of matrix concentration inequalities but follow the spirit of the
analysis in [5].

#### 1\. Symmetrization

Write the normwise error by exploiting diagonal Schur products (1)

$\displaystyle\mathcal{D}\,(\widehat{\boldsymbol{A}})-\mathcal{D}\,(\boldsymbol{A})=\boldsymbol{I}\circ\widehat{\boldsymbol{A}}-\boldsymbol{I}\circ\boldsymbol{A}=\frac{1}{N}\sum_{k=1}^{N}{\left(\boldsymbol{I}\circ(\boldsymbol{A}\boldsymbol{w}_{k}\boldsymbol{w}_{k}^{\top})-\boldsymbol{I}\circ\boldsymbol{A}\right)},$

and take expectations of the squared norms

$\displaystyle\mathbb{E}\left[\left\|\mathcal{D}\,(\widehat{\boldsymbol{A}})-\mathcal{D}\,(\boldsymbol{A})\right\|_{2}^{2}\right]=\frac{1}{N^{2}}\,\mathbb{E}\left[\left\|\sum_{k=1}^{N}{\left(\boldsymbol{I}\circ(\boldsymbol{A}\boldsymbol{w}_{k}\boldsymbol{w}_{k}^{\top})-\boldsymbol{I}\circ\boldsymbol{A}\right)}\right\|_{2}^{2}\right].$

From the assumption that $\boldsymbol{w}_{k}$ has independent elements with
zero mean and and variance 1 follows
$\mathbb{E}\left[\boldsymbol{I}\circ(\boldsymbol{A}\boldsymbol{w}_{k}\boldsymbol{w}_{k}^{\top})\right]=\boldsymbol{I}\circ\boldsymbol{A}$.
Hence the matrix random variables

$\boldsymbol{X}_{k}\equiv\boldsymbol{I}\circ(\boldsymbol{A}\boldsymbol{w}_{k}\boldsymbol{w}_{k}^{\top})-\boldsymbol{I}\circ\boldsymbol{A},\qquad
1\leq k\leq N$

have zero mean. We use symmetrization [20, Lemma 6.4.2] to create symmetric
random variables $\varepsilon_{k}\,\boldsymbol{X}_{k}$, where
$\varepsilon_{k}$ are independent symmetric Bernoulli random variables, that
is, they are Rademacher variables as in Definition 3.1. The Rademacher
variables $\varepsilon_{k}$ are independent of each other and also independent
of the random vectors $\boldsymbol{w}_{k}$. Remark 3.2 implies
$\operatorname{\mathbb{E}}[\varepsilon_{k}]=0$, hence

(16)
$\mathbb{E}\left[\sum_{k=1}^{N}{\varepsilon_{k}\,\boldsymbol{I}\circ{\boldsymbol{A}}}\right]=\boldsymbol{0}.$

Then [20, Lemma 6.4.2], (16) and (2) imply

(17)
$\displaystyle\mathbb{E}\left[\left\|\sum_{k=1}^{N}\boldsymbol{X}_{k}\right\|_{2}^{2}\right]\leq
2\,\mathbb{E}\left[\left\|\sum_{k=1}^{N}{\varepsilon_{k}\,\boldsymbol{I}\circ(\boldsymbol{A}\boldsymbol{w}_{k}\boldsymbol{w}_{k}^{\top})}\right\|_{2}^{2}\right]=2\,\mathbb{E}\left[\left\|\sum_{k=1}^{N}{\boldsymbol{Y}_{k}}\right\|_{2}^{2}\right],$

where
$\boldsymbol{Y}_{k}\equiv\varepsilon_{k}\,\operatorname{\mathrm{diag}}(\boldsymbol{A}\boldsymbol{w}_{k})\operatorname{\mathrm{diag}}(\boldsymbol{w}_{k})$
are symmetric random matrices, and the second and third expectations range
over all random vectors $\boldsymbol{w}_{k}$ and all Rademacher variables
$\varepsilon_{k}$.

#### 2\. Concentration inequality

Applying Cauchy-Schwartz inequality and Theorem 2.5 to the sum
$\boldsymbol{Z}\equiv\sum_{k=1}^{N}\boldsymbol{Y}_{k}$ gives

(18)
$\displaystyle\mathbb{E}[\|\boldsymbol{Z}\|_{2}^{2}]^{1/2}\leq\sqrt{2e\ln{n}}\>\left\|\left(\sum_{i=1}^{N}{\mathbb{E}[\boldsymbol{Y}_{k}^{2}]}\right)^{1/2}\right\|_{2}+4e\ln{n}\>\mathbb{E}\left[\max_{1\leq
k\leq N}{\|\boldsymbol{Y}_{k}\|_{2}^{2}}\right]^{1/2}.$

We bound the expectations that represent the matrix variance and the maximal
two-norm separately.

#### 3\. Variance

As in item 3 of the proof of Theorem 3.10 abbreviate

$\displaystyle\boldsymbol{D}_{k}\equiv\operatorname{\mathrm{diag}}(\boldsymbol{Aw}_{k}),\qquad\boldsymbol{W}_{k}\equiv\operatorname{\mathrm{diag}}(\boldsymbol{w}_{k}),\qquad\boldsymbol{O}\equiv\operatorname{\mathrm{diag}}\left(\begin{bmatrix}\|\boldsymbol{a}_{1}\|_{1}&\cdots&\|\boldsymbol{a}_{n}\|_{1}\end{bmatrix}\right).$

With this notation,
$\mathbb{E}[\boldsymbol{Y}_{k}^{2}]=\mathbb{E}[\boldsymbol{D}_{k}^{2}\boldsymbol{W}_{k}^{2}]$.
Consider the diagonal term

$(\boldsymbol{D}_{k}^{2}\boldsymbol{W}_{k}^{2})_{ii}=\left(\sum_{j=1}^{n}a_{ij}(\boldsymbol{w}_{k})_{j}\right)^{2}(\boldsymbol{w}_{k})_{i}^{2}\leq\|\boldsymbol{a}_{i}\|_{1}^{2}\|\boldsymbol{w}_{k}\|_{\infty}^{4}\qquad
1\leq i\leq n.$

Therefore,
$\boldsymbol{D}_{k}^{2}\boldsymbol{W}_{k}^{2}\preceq\|\boldsymbol{w}_{k}\|_{\infty}^{4}\boldsymbol{O}^{2}$.
The symmetry of $\boldsymbol{A}$ implies
$\boldsymbol{O}\preceq\|\boldsymbol{A}\|_{\infty}\boldsymbol{I}$. Combining
the two inequalities gives

(19)
$\displaystyle\mathbb{E}[\boldsymbol{Y}_{k}^{2}]=\mathbb{E}[\boldsymbol{D}_{k}^{2}\boldsymbol{W}_{k}^{2}]\preceq\|\boldsymbol{A}\|_{\infty}^{2}\>\mathbb{E}\left[\max_{1\leq
k\leq N}\|\boldsymbol{w}_{k}\|_{\infty}^{4}\right]\boldsymbol{I}.$

Now take square roots,

$\displaystyle\left(\sum_{k=1}^{N}{\mathbb{E}[\boldsymbol{Y}_{k}^{2}]}\right)^{1/2}\preceq\sqrt{N}\>\|\boldsymbol{A}\|_{\infty}\>\mathbb{E}\left[\max_{1\leq
k\leq N}{\|\boldsymbol{w}_{k}\|_{\infty}^{4}}\right]^{1/2}\boldsymbol{I},$

and bound the norm,

(20)
$\displaystyle\left\|\left(\sum_{k=1}^{N}{\mathbb{E}[\boldsymbol{Y}_{k}^{2}]}\right)^{1/2}\right\|_{2}\leq\sqrt{N}\>\|\boldsymbol{A}\|_{\infty}\>\mathbb{E}\left[\max_{1\leq
k\leq N}{\|\boldsymbol{w}_{k}\|_{\infty}^{4}}\right]^{1/2}.$

#### 4: Maximal two-norm

In analogy to (19), we derive

$\displaystyle\mathbb{E}\left[\max_{1\leq k\leq
N}{\|\boldsymbol{Y}_{k}\|_{2}^{2}}\right]\leq\|\boldsymbol{A}\|_{\infty}^{2}\>\mathbb{E}\left[\max_{1\leq
k\leq N}{\|\boldsymbol{w}_{k}\|_{\infty}^{4}}\right],$

and its square root

(21) $\displaystyle\mathbb{E}\left[\max_{1\leq k\leq
N}{\|\boldsymbol{Y}_{k}\|_{2}^{2}}\right]^{1/2}\leq\|\boldsymbol{A}\|_{\infty}\>\mathbb{E}\left[\max_{1\leq
k\leq N}{\|\boldsymbol{w}_{k}\|_{\infty}^{4}}\right]^{1/2}.$

#### 5\. Putting everything together

Substitute the variance bound (20) and the norm bound (21) into the
expectation (18) for the sum

$\mathbb{E}[\|\boldsymbol{Z}\|_{2}]^{1/2}\leq\left(\sqrt{2e\ln{n}}\,N^{1/2}+4e\ln{n}\right)\|\boldsymbol{A}\|_{\infty}\>\mathbb{E}\left[\max_{1\leq
k\leq N}{\|\boldsymbol{w}_{k}\|_{\infty}^{4}}\right]^{1/2},$

substitute this, in turn, into the expectation (17) for the absolute error,

$\displaystyle\mathbb{E}\left[\|\mathcal{D}\,(\widehat{\boldsymbol{A}})-\mathcal{D}\,(\boldsymbol{A})\|_{2}^{2}\right]^{1/2}\leq\frac{2}{N}\>\mathbb{E}\left[\|\boldsymbol{Z}\|_{2}\right]^{1/2}$
$\displaystyle\qquad\qquad\leq\frac{2}{N}\>\left(\sqrt{2e\ln{n}}\,N^{1/2}+4e\ln{n}\right)\|\boldsymbol{A}\|_{\infty}\>\mathbb{E}\left[\max_{1\leq
k\leq N}{\|\boldsymbol{w}_{k}\|_{\infty}^{4}}\right]^{1/2},$

and simplify.

###  4.2 Gaussian vectors

We determine the minimal sampling amount for Gaussian Monte Carlo estimators
to be normwise $(\epsilon,\delta)$ diagonal estimators.

######  Corollary 4.3.

Let $\boldsymbol{A}\in\mathbb{R}^{n\times n}$ with $n\geq 3$ be symmetric, and
let

$\widehat{\boldsymbol{A}}\equiv\frac{1}{N}\sum_{j=1}^{N}{\boldsymbol{Aw}_{k}\boldsymbol{w}_{k}^{\top}}$

be a Monte Carlo estimator with independent Gaussian random vectors
$\boldsymbol{w}_{k}\sim\mathcal{N}(\boldsymbol{0},\boldsymbol{I})$ in
$\mathbb{R}^{n}$, $1\leq k\leq N$. Pick $\epsilon>0$. For any $0<\delta<1$, if
the sampling amount $N$ satisfies ${8e\ln n}\leq N\leq n$, and is at least

(22) $\displaystyle
N\geq\frac{128\,(e\,\ln{n})^{3}}{\epsilon^{2}\,\delta}\>\left(\frac{\|\boldsymbol{A}\|_{\infty}}{\|\mathcal{D}\,(\boldsymbol{A})\|_{\infty}}\right)^{2},$

then
$\|\mathcal{D}\,(\widehat{\boldsymbol{A}})-\mathcal{D}\,(\boldsymbol{A})\|_{2}\leq\epsilon\|\mathcal{D}\,(\boldsymbol{A})\|_{2}$
holds with probability at least $1-\delta$.

######  Proof 4.4.

For Gaussian random vectors
$\boldsymbol{w}_{k}\sim\mathcal{N}(\boldsymbol{0},\boldsymbol{I})$, [5, (3.7)]
implies

$\left(\mathbb{E}\left[\max_{1\leq k\leq
N}{\|\boldsymbol{w}_{k}\|_{\infty}^{4}}\right]\right)^{1/2}\leq
e\ln{(nN)}\max_{1\leq i,j\leq n}{|\boldsymbol{I}_{ij}|}=e\ln{(nN)}.$

Substituting this into Theorem 4.1 gives

$\left(\mathbb{E}\left[\|\mathcal{D}\,(\widehat{\boldsymbol{A}})-\mathcal{D}\,(\boldsymbol{A})\|_{2}^{2}\right]\right)^{1/2}\leq\left(\sqrt{\frac{8e\ln{n}}{N}}+\frac{8e\ln{n}}{N}\right)\>e\ln{(nN})\>\|\boldsymbol{A}\|_{\infty}.$

Square both sides and apply Markov’s inequality (Theorem 2.1) to the random
variable
$Z\equiv\|\mathcal{D}\,(\widehat{\boldsymbol{A}})-\mathcal{D}\,(\boldsymbol{A})\|_{2}$
using $t\equiv\epsilon\|\mathcal{D}\,(\boldsymbol{A})\|_{2}$,

(23)
$\displaystyle\qquad\qquad\leq\left(\sqrt{\frac{8e\ln{n}}{N}}+\frac{8e\ln{n}}{N}\right)^{2}\>(e\ln{(nN)})^{2}\>\frac{\|\boldsymbol{A}\|_{\infty}^{2}}{\epsilon^{2}\,\|\mathcal{D}\,(\boldsymbol{A})\|_{\infty}^{2}}.$

Substituting the assumption ${8e\ln{n}}\leq N\leq n$ into the relevant part of
the above bound gives

$\displaystyle\left(\sqrt{\frac{8e\ln{n}}{N}}+\frac{8e\ln{n}}{N}\right)^{2}\>(e\ln{(nN)})^{2}\leq\left(2\sqrt{\frac{8e\ln{n}}{N}}\right)^{2}(2e\ln{n})^{2}=\frac{128\,(e\,\ln{n})^{3}}{N}.$

Substitute this, in turn, into (23), set the failure probability equal to

$\displaystyle\delta\equiv\frac{128\,(e\,\ln{n})^{3}}{\epsilon^{2}\,N}\>\left(\frac{\|\boldsymbol{A}\|_{\infty}}{\|\mathcal{D}\,(\boldsymbol{A})\|_{\infty}}\right)^{2}$

and solve for the sampling amount $N$.

##  5 Componentwise bounds

We present componentwise bounds for Monte Carlo estimators based on standard
Rademacher vectors (section 5.1), as well as on standard (section 5.2) and
normalized Gaussian vectors (section 5.3).

Our attempt at deriving alternative normwise bounds by applying a union bound
over the componentwise bounds for all diagonal elements did not produce
results that were substantially tighter than our previous normwise bounds.

###  5.1 Standard Rademacher vectors

We present a componentwise worst case absolute error bound (Corollary 5.1),
and a bound on the minimal sampling amount that makes the Rademacher Monte
Carlo estimator a componentwise ($\epsilon,\delta)$ diagonal estimator
(Theorem 5.5).

Theorem 3.4 is identical to the following worst case componentwise bound.

######  Corollary 5.1.

Let $\boldsymbol{A}\in\mathbb{R}^{n\times n}$ be non-diagonal symmetric, and
let

$K_{1}\equiv\|\mathcal{D}\,(\boldsymbol{A}^{2})-\mathcal{D}\,(\boldsymbol{A})^{2}\|_{2},\quad
K_{2}\equiv\|\boldsymbol{A}-\mathcal{D}\,(\boldsymbol{A})\|_{\infty},\quad
d\equiv\frac{1}{K_{1}}\|\boldsymbol{A}-\mathcal{D}\,(\boldsymbol{A})\|_{F}^{2}.$

If
$\widehat{\boldsymbol{A}}\equiv\frac{1}{N}\sum_{k=1}^{N}\boldsymbol{Aw}_{k}\boldsymbol{w}_{k}^{\top}$
is a Monte Carlo estimator with independent Rademacher vectors
$\boldsymbol{w}_{k}\in\mathbb{R}^{n}$, $1\leq k\leq N$, then the probability
that the absolute error exceeds $t>0$ is at most

$\operatorname{\mathbb{P}}\left[\max_{1\leq i\leq
n}{|a_{ii}-\widehat{a}_{ii}|}\geq t\right]\leq
8d\,\exp\left(\frac{-Nt^{2}}{2(K_{1}+tK_{2}/3)}\right).$

######  Proof 5.2.

In Theorem 3.4, the $p$-norm of the diagonal matrix
$\mathcal{D}\,(\boldsymbol{A})-\mathcal{D}\,(\widehat{\boldsymbol{A}})$ is a
largest magnitude diagonal element.

In contrast to Corollary 5.1, the next bound depends on the particular
diagonal element. We first require a lemma on the independence of products of
Rademacher variables.

######  Lemma 5.3.

Let $Z_{1},W_{1},W_{2},\ldots,W_{n}$ be independent Rademacher variables. Then
the products $X_{1}\equiv ZW_{1},\ldots,X_{n}\equiv ZW_{n}$ are also
independent Rademacher variables.

######  Proof 5.4.

For any $x_{1},\ldots,x_{n}\in\\{-1,+1\\}$, using the law of total
probability, the joint probability mass function satisfies

$\displaystyle\mathbb{P}\left[\cap_{i=1}^{n}\,\\{X_{i}=x_{i}\\}\right]$
$\displaystyle=\sum_{z\in\\{-1,+1\\}}\mathbb{P}\left[\cap_{i=1}^{n}\,\\{X_{i}=x_{i}\\}|Z=z\right]\mathbb{P}[Z=z]$
$\displaystyle=\frac{1}{2}\mathbb{P}\left[\cap_{i=1}^{n}\,\\{W_{i}=-x_{i}\\}\right]+\frac{1}{2}\mathbb{P}\left[\cap_{i=1}^{n}\,\\{W_{i}=x_{i}\\}\right]$
$\displaystyle=2^{-n}=\prod_{i=1}^{n}\mathbb{P}\left[X_{i}=x_{i}\right].$

Since the joint PMF factorizes, the variables $X_{1},\ldots,X_{n}$ are
independent.

######  Theorem 5.5.

Let $\boldsymbol{A}\in\mathbb{R}^{n\times n}$ be non-diagonal symmetric, and
let

$\widehat{\boldsymbol{A}}\equiv\frac{1}{N}\sum_{k=1}^{N}\boldsymbol{Aw}_{k}\boldsymbol{w}_{k}^{\top}$

be a Monte Carlo estimator with independent Rademacher vectors
$\boldsymbol{w}_{k}\in\mathbb{R}^{n}$, $1\leq k\leq N$. The probability that
the absolute error exceeds $t>0$ is at most

(24)

######  Proof 5.6.

Fix $i$, for some $1\leq i\leq n$. The properties in Remark 3.2 allow us to
split off the original diagonal element from the estimator,

$\displaystyle\widehat{a}_{ii}$
$\displaystyle=\frac{1}{N}\sum_{k=1}^{N}\left(\boldsymbol{A}\boldsymbol{w}_{k}\boldsymbol{w}_{k}^{\top}\right)_{ii}=\frac{1}{N}\sum_{k=1}^{N}\sum_{j=1}^{n}{a_{ij}\,(\boldsymbol{w}_{k})_{i}\,(\boldsymbol{w}_{k})_{j}}$
$\displaystyle=a_{ii}+\sum_{k=1}^{N}\sum_{j\neq
i}\underbrace{\frac{a_{ij}}{N}(\boldsymbol{w}_{k})_{i}\,(\boldsymbol{w}_{k})_{j}}_{Z_{kj}^{(i)}}.$

Lemma 5.3 implies that for fixed $i$, the $Z_{kj}^{(i)}$ are independent.
Remark 3.2 implies that they have zero mean, and are bounded by

$\displaystyle-\tfrac{|a_{ij}|}{N}\leq
Z_{kj}^{(i)}\leq\tfrac{|a_{ij}|}{N},\qquad 1\leq k\leq N,\quad 1\leq j\leq
n,\quad j\neq i.$

Hence the absolute error $\widehat{a}_{ii}-a_{ii}$ is a sum of independent
bounded zero-mean random variables, and we can apply Hoeffding’s inequality in
Theorem 2.2

Problem 4

At last, write $\sum_{j\neq
i}a_{ij}^{2}=\|\boldsymbol{a}_{i}\|_{2}^{2}-a_{ii}^{2}$.

Theorem 5.5 implies that the accuracy for estimating a single diagonal element
depends only on the magnitude of the off-diagonal elements in the
corresponding row and column.

We determine the minimal sampling amount required to make the Rademacher Monte
Carlo estimator a componentwise ($\epsilon,\delta)$ diagonal estimator. For
symmetric matrices, this result coincides with the bound in Equation 40 of [3]
but uses a different proof technique.

######  Corollary 5.7.

Let $\boldsymbol{A}\in\mathbb{R}^{n\times n}$ be non-diagonal symmetric, and
let

$\widehat{\boldsymbol{A}}\equiv\frac{1}{N}\sum_{k=1}^{N}\boldsymbol{Aw}_{k}\boldsymbol{w}_{k}^{\top}$

be a Monte Carlo estimator with independent Rademacher vectors
$\boldsymbol{w}_{k}\in\mathbb{R}^{n}$, $1\leq k\leq N$. Pick $\epsilon>0$, and
a diagonal element $a_{ii}\neq 0$ of $\boldsymbol{A}$. For any $0<\delta<1$,
if the sampling amount is at least

$N\geq\left(\frac{\|\boldsymbol{a}_{i}\|_{2}^{2}-a_{ii}^{2}}{a_{ii}^{2}}\right)\frac{2\ln(2/\delta)}{\epsilon^{2}},$

then $|a_{jj}-\widehat{a}_{jj}|\leq\epsilon|a_{jj}|$ holds with probability at
least $1-\delta$.

######  Proof 5.8.

Define the 2-norm offdiagonal column sums

$\mathrm{off}_{i}\equiv(\|\boldsymbol{a}_{i}\|_{2}^{2}-a_{ii}^{2})^{1/2},\qquad
1\leq i\leq n,$

and denote the bound for the failure probability in Theorem 5.5 by

$\delta\equiv 2\exp\left(\frac{-Nt^{2}}{2\,\mathrm{off}_{i}^{2}}\right),$

and solve it for $t$,

$t=\sqrt{\frac{2\,\mathrm{off}_{i}^{2}}{N}\ln(2/\delta)}.$

Restate Theorem 5.5 in terms of the failure probability: With probability at
most $1-\delta$, the absolute error of a specific diagonal element is bounded
above by

$|a_{ii}-\widehat{a}_{ii}|\leq
t=\sqrt{\frac{2\,\mathrm{off}_{i}^{2}}{N}\ln(2/\delta)}.$

Converting this absolute error into a relative error requires
$t\leq\epsilon|a_{ii}|$, which implies

$N\geq\left(\frac{\mathrm{off}_{i}}{|a_{ii}|}\right)^{2}\frac{2\ln(2/\delta)}{\epsilon^{2}}.$

The minimal sampling amount for computing $a_{ii}$ with the Rademacher Monte
Carlo estimator depends on
$\frac{\|\boldsymbol{a}_{i}\|_{2}^{2}-a_{ii}^{2}}{|a_{ii}|}$, which represents
the relative 2-norm deviation of the $i$th column and row of $\boldsymbol{A}$
from diagonality. Thus, the more diagonal the $i$th row and column, the fewer
samples are required for a $(\epsilon,\delta)$ estimator.

###  5.2 Gaussian vectors

We present a componentwise absolute error bound (Theorem 5.9) for Gaussian
Monte Carlo estimators, and a bound on the minimal sampling amount that makes
the Gaussian Monte Carlo estimator a componentwise $(\epsilon,\delta)$
estimators (Corollary 5.11). Our bounds are derived from and identical to
bounds for trace estimators in [7].

######  Theorem 5.9.

Let $\boldsymbol{A}\in\mathbb{R}^{n\times n}$ be non-diagonal symmetric,

$L_{1i}\equiv|a_{ii}|+\|\boldsymbol{a}_{i}\|_{2},\qquad
L_{i2}\equiv|a_{ii}|^{2}+\|\boldsymbol{a}_{i}\|_{2}^{2},\qquad 1\leq i\leq n,$

and let
$\widehat{\boldsymbol{A}}\equiv\frac{1}{N}\sum_{k=1}^{N}\boldsymbol{A}\boldsymbol{z}_{k}\boldsymbol{z}_{k}^{\top}\in\mathbb{R}^{n\times
n}$ be a Monte Carlo estimator with independent Gaussian vectors
$\boldsymbol{z}_{k}\sim\mathcal{N}(\boldsymbol{0},\boldsymbol{I})$ in
$\mathbb{R}^{n}$, $1\leq k\leq N$. If $t>0$, then

(25)

######  Proof 5.10.

Fix $i$ for some $1\leq i\leq n$. Write the diagonal element as an inner
product

$(\boldsymbol{A}\boldsymbol{z}_{k}\boldsymbol{z}_{k}^{\top})_{ii}=\sum_{j=1}^{n}{a_{ij}(\boldsymbol{z}_{k})_{i}(\boldsymbol{z}_{k})_{j}}=\boldsymbol{z}_{k}^{\top}\boldsymbol{B}_{i}\boldsymbol{z}_{k},\qquad
1\leq k\leq N$

involving the symmetric matrix

$\boldsymbol{B}_{i}\equiv\begin{bmatrix}0&&\frac{1}{2}a_{1i}&&0\\\
&&\vdots&&\\\ \frac{1}{2}a_{i1}&\cdots&a_{ii}&\cdots&\frac{1}{2}a_{in}\\\
&&\vdots&&\\\ 0&&\frac{1}{2}a_{ni}&&0\\\ \end{bmatrix}\in{\mathbb{R}}^{n\times
n}\quad\text{with}\quad\operatorname{\mathrm{trace}}(\boldsymbol{B}_{i})=a_{ii}.$

We can think of

$\widehat{a}_{ii}=\left(\frac{1}{N}\sum_{k=1}^{N}{\boldsymbol{A}\boldsymbol{z}_{k}\boldsymbol{z}_{k}^{\top}}\right)_{ii}=\frac{1}{N}\sum_{k=1}^{N}{\left(\boldsymbol{A}\boldsymbol{z}_{k}\boldsymbol{z}_{k}^{\top}\right)_{ii}}=\frac{1}{N}\sum_{k=1}^{N}{\boldsymbol{z}_{k}^{\top}\boldsymbol{B}_{i}\boldsymbol{z}_{k}}$

as a Monte Carlo estimator for
$\operatorname{\mathrm{trace}}(\boldsymbol{B}_{i})$, and apply the bound for
Gaussian trace estimators [7, Theorem 1]

Problem 6

where

$\|\boldsymbol{B}_{i}\|_{F}^{2}=\frac{1}{2}\left(|a_{ii}|^{2}+\|\boldsymbol{a}_{i}\|_{2}^{2}\right)=\frac{L_{i2}}{2},\qquad\|\boldsymbol{B}_{i}\|_{2}=\frac{1}{2}\left(|a_{ii}|+\|\boldsymbol{a}_{i}\|_{2}\right)=\frac{L_{i1}}{2}.$

Theorem 5.9 implies that with Gaussian vectors, the accuracy for estimating a
single diagonal element depends on the magnitude of all elements in the
corresponding row and column. By contrast, the bounds of Theorem 5.5 for
Rademacher vectors depend only on the magnitude of the off-diagonal elements.

We determine the minimal sampling amount required to make the Gaussian Monte
Carlo estimator a componentwise $(\epsilon,\delta)$ estimator.

######  Corollary 5.11.

Let $\boldsymbol{A}\in\mathbb{R}^{n\times n}$ be non-diagonal symmetric,

$\Delta_{1i}\equiv
1+\frac{\|\boldsymbol{a}_{i}\|_{2}}{|a_{ii}|},\qquad\Delta_{i2}\equiv
1+\left(\frac{\|\boldsymbol{a}_{i}\|_{2}}{|a_{ii}|}\right)^{2},\qquad 1\leq
i\leq n,$

and let
$\widehat{\boldsymbol{A}}\equiv\frac{1}{N}\sum_{k=1}^{N}\boldsymbol{A}\boldsymbol{z}_{k}\boldsymbol{z}_{k}^{\top}\in\mathbb{R}^{n\times
n}$ be a Monte Carlo estimator with independent Gaussian vectors
$\boldsymbol{z}_{k}\sim\mathcal{N}(\boldsymbol{0},\boldsymbol{I})$ in
$\mathbb{R}^{n}$, $1\leq k\leq N$. Pick $\epsilon>0$, and a diagonal element
$a_{ii}\neq 0$ of $\boldsymbol{A}$. For any $0<\delta<1$, if the sampling
amount is at least

$N\geq(\Delta_{2i}+\Delta_{1i}\epsilon)\,\frac{2\ln{(2/\delta)}}{\epsilon^{2}}$

then $|\widehat{a}_{ii}-a_{ii}|\leq\epsilon\,|a_{ii}|$ holds with probability
at least $1-\delta$.

######  Proof 5.12.

This follows immediately from the lower bound for $N$ in [7, Theorem 1].

The required sampling amount for computing $a_{ii}$ with the Gaussian Monte
Carlo estimator depends on $\|\boldsymbol{a}_{i}\|_{2}/|a_{ii}|$ which can be
interpreted as the 2-norm derivation of the $i$th column and row of
$\boldsymbol{A}$ from diagonality. The more diagonal the $i$th row and column,
the smaller the sampling amount for the $(\epsilon,\delta)$ estimator.

###  5.3 Normalized Gaussian vectors

We extend and complete the analysis in [4] for a Monte Carlo estimator based
on normalized Gaussian vectors,

(26)
$\widehat{\boldsymbol{A}}\equiv\left(\sum_{k=1}^{N}\boldsymbol{A}\boldsymbol{z}_{k}\boldsymbol{z}_{k}^{\top}\right)\oslash\left(\sum_{k=1}^{N}\boldsymbol{z}_{k}\boldsymbol{z}_{k}^{\top}\right),$

where $\boldsymbol{z}_{k}\in\mathbb{R}^{n}$ are independent random vectors,
and $\oslash$ denotes elementwise division. We derive the distribution of the
componentwise absolute errors (Lemma 5.13), followed by a bound (Theorem
5.16).

We represent the distribution for the absolute errors in the diagonal elements
in terms of a Student $t$-distribution with $N\geq 1$ degrees of freedom [13,
Definition 7.3.3],

(27) $T_{N}\equiv\frac{Z}{\sqrt{U/N}},$

where $Z$ is a Gaussian $\mathcal{N}(0,1)$ random variable, and $U$ an
independent chi-square random variable with $N$ degrees of freedom.

######  Lemma 5.13.

Let $\boldsymbol{A}\in{\mathbb{R}}^{n\times n}$ be symmetric, and (26) be a
Monte Carlo estimator where
$\boldsymbol{z}_{k}\sim\mathcal{N}(\boldsymbol{0},\boldsymbol{I})$ in
$\mathbb{R}^{n}$ are independent Gaussian random vectors, $1\leq k\leq N$. The
absolute errors in the diagonal elements are distributed as

$\widehat{a}_{ii}-a_{ii}\ \sim\
\sqrt{\frac{\|\boldsymbol{a}_{i}\|_{2}^{2}-a_{ii}^{2}}{N}}\,T_{N},\qquad 1\leq
i\leq n.$

######  Proof 5.14.

Due to the normalization in the denominator, we can extract the diagonal
elements of $\boldsymbol{A}$ from the diagonal elements of the Monte Carlo
estimator $\widehat{\boldsymbol{A}}$,

$\widehat{a}_{ii}=\frac{\sum_{k=1}^{N}{\sum_{j=1}^{n}{a_{ij}(\boldsymbol{z}_{k})_{i}(\boldsymbol{z}_{k})_{j}}}}{\sum_{k=1}^{N}{(\boldsymbol{z}_{k})_{i}^{2}}}=a_{ii}+\frac{\sum_{k=1}^{N}{\sum_{j\neq
i}{a_{ij}\,(\boldsymbol{z}_{k})_{i}(\boldsymbol{z}_{k})_{j}}}}{\sum_{k=1}^{N}{(\boldsymbol{z}_{k})_{i}^{2}}},\qquad
1\leq i\leq n.$

Normalize across the $i$th elements of the Gaussian vectors
$\boldsymbol{z}_{k}$ to unit vectors $\boldsymbol{u}_{i}\in\mathbb{R}^{N}$
with elements

$(\boldsymbol{u}_{i})_{k}\equiv\frac{(\boldsymbol{z}_{k})_{i}}{\sqrt{\sum_{k^{\prime}=1}^{N}{(\boldsymbol{z}_{k^{\prime}})_{i}^{2}}}},\qquad
1\leq k\leq N,\quad 1\leq i\leq n.$

Use the denominator to normalize the $i$th component in the $i$th absolute
error,

(28) $\widehat{a}_{ii}-a_{ii}=\frac{\sum_{k=1}^{N}\sum_{j\neq
i}{a_{ij}\,(\boldsymbol{z}_{k})_{j}(\boldsymbol{u}_{i})_{k}}}{\sqrt{\sum_{k^{\prime}=1}^{N}{(\boldsymbol{z}_{k^{\prime}})_{i}^{2}}}},\qquad
1\leq i\leq n.$

The rotational invariance of the standard Gaussian distribution guarantees the
independence of the direction vectors $\boldsymbol{u}_{i}$ and radial
components $\sqrt{\sum_{k=1}^{N}(\boldsymbol{z}_{k})_{i}^{2}}$; see [20,
Exercise 3.3.6]. Hence, the numerator and denominator in (28) are independent.
We can rewrite (28) as

$\widehat{a}_{ii}-a_{ii}=\frac{Z_{i}}{\sqrt{U_{i}/N}},$

where we define the random variables $Z_{i}\equiv\sum_{k=1}^{N}\sum_{j\neq
i}{a_{ij}\,(\boldsymbol{z}_{k})_{j}(\boldsymbol{u}_{i})_{k}}/\sqrt{N}$ and
$U_{i}\equiv\sum_{k=1}^{N}(\boldsymbol{z}_{k})_{i}^{2}$ for $1\leq i\leq n$.
The random variable $U_{i}$ has a chi-square distribution with $N$ degrees of
freedom. The conditional distribution of $Z_{i}$ given $\boldsymbol{u}_{i}$ is
Gaussian (see e.g. Exercise 3.3.3(a) of [20]), with zero mean and variance

$\frac{1}{N}\sum_{k=1}^{N}\sum_{j\neq
i}a_{ij}^{2}(\boldsymbol{u}_{i})_{k}^{2}=\frac{1}{N}\left(\sum_{j\neq
i}a_{ij}^{2}\right)\left(\sum_{k=1}^{N}(\boldsymbol{u}_{i})_{k}^{2}\right)=\frac{1}{N}\sum_{j\neq
i}a_{ij}^{2}=\frac{1}{N}(\|\boldsymbol{a}_{i}\|_{2}^{2}-a_{ii}^{2}).$

Therefore,
$Z_{i}|\boldsymbol{u}_{i}\sim\mathcal{N}(0,\frac{1}{N}(\|\boldsymbol{a}_{i}\|_{2}^{2}-a_{ii}^{2}))$.
However, the conditional distribution is independent of $\boldsymbol{u}_{i}$,
so this is also the unconditional distribution. The claim then follows from
(27).

######  Remark 5.15.

It was observed in [3] that the square of the error
$(\widehat{a}_{ii}-a_{ii})^{2}$ has a scaled $F$-distribution. We note that
the square of a Student $t$-distribution is specifically a scaled
$F$-distribution with one degree of freedom in the numerator. Moreover, for a
single sample $N=1$ the error has a Cauchy distribution, which has undefined
mean and variance.

If $N$ is large, then the $t$-distribution $T_{N}$ can be approximated by a
standard normal distribution. However, $T_{N}$ has wider tails, thus somewhat
weaker tail bounds. Existing tail bounds for the Student $t$-distribution
imply the following concentration inequality for error bounds.

######  Theorem 5.16.

Let $\boldsymbol{A}\in{\mathbb{R}}^{n\times n}$ be symmetric, and (26) be a
Monte Carlo estimator where
$\boldsymbol{z}_{k}\sim\mathcal{N}(\boldsymbol{0},\boldsymbol{I})$ in
$\mathbb{R}^{n}$ are independent Gaussian random vectors, $1\leq k\leq N$. For
any $t>0$,

Problem 7

######  Proof 5.17.

The probability density function of $T_{N}$ is [17]

$f_{N}(x)=c_{N}\left(1+\frac{x^{2}}{N}\right)^{-\frac{N+1}{2}}\qquad\text{where}\quad
c_{N}\equiv\frac{\Gamma((N+1)/2)}{\Gamma(N/2)\sqrt{N\pi}},$

and $1/\pi\leq c_{N}\leq 1/\sqrt{2\pi}$. If $F_{N}(x)$ is the cumulative
distribution function for $T_{N}$ then by [17, Theorem 3.1]

$\displaystyle\prob\left[(\widehat{a}_{ii}-a_{ii})>\sqrt{\frac{\|\boldsymbol{a}_{i}\|_{2}^{2}-a_{ii}^{2}}{N}}x\right]=$
$\displaystyle\>1-F_{N}(x)<\frac{f_{N}(x)}{x}\left(1+\frac{x^{2}}{N}\right)$
$\displaystyle=$
$\displaystyle\>\frac{c_{N}}{x}\left(1+\frac{x^{2}}{N}\right)^{-\frac{N-1}{2}}.$

Take $x=t\,\sqrt{\frac{N}{\|\boldsymbol{a}_{i}\|_{2}^{2}-a_{ii}^{2}}}$ and
bound the upper tail with $c_{N}\leq\sqrt{\frac{1}{2\pi}}$, to obtain

Since $T_{N}$ is symmetric about the origin, the lower tail has the same
bound. Using a union bound gives the desired inequality.

We determine a sampling amount sufficient to make the normalized Gaussian
Monte Carlo estimator a componentwise $(\epsilon,\delta)$ estimator.

######  Corollary 5.18.

Let $\boldsymbol{A}\in\mathbb{R}^{n\times n}$ be non-diagonal symmetric, and
let

$\Psi_{i}\equiv\frac{|a_{ii}|}{(\|\boldsymbol{a}_{i}\|_{2}^{2}-a_{ii}^{2})^{1/2}},\quad
1\leq i\leq n,$

and let $\widehat{\boldsymbol{A}}$ be defined as in (26). Pick $\epsilon>0$
and a diagonal element $a_{ii}\neq 0$ of $\boldsymbol{A}$. For any
$0<\delta<1$, if the sampling amount is positive and at least

$N\geq
1+2\ln\left(\frac{\sqrt{2/\pi}}{\delta\epsilon\Psi_{i}}\right)/\ln(1+\epsilon^{2}\Psi_{i}^{2}),$

then $|\widehat{a}_{ii}-a_{ii}|\leq\epsilon|a_{ii}|$ holds with probability at
least $1-\delta$.

######  Proof 5.19.

In Theorem 5.16, set $t=\epsilon|a_{ii}|$. If the sampling number satisfies
the desired bound, it follows that

$\displaystyle\leq\sqrt{\frac{2\,(\|\boldsymbol{a}_{i}\|_{2}^{2}-a_{ii}^{2})}{\pi
N}}\,\frac{1}{t}\,\left(1+\frac{t^{2}}{\|\boldsymbol{a}_{i}\|_{2}^{2}-a_{ii}^{2}}\right)^{-\frac{N-1}{2}}$
$\displaystyle=\frac{\sqrt{2/\pi}}{\epsilon\Psi_{i}\sqrt{N}}\left(1+\epsilon^{2}\Psi_{i}^{2}\right)^{-\frac{N-1}{2}}$
$\displaystyle\leq\frac{\sqrt{2/\pi}}{\epsilon\Psi_{i}}\left(1+\epsilon^{2}\Psi_{i}^{2}\right)^{-\frac{N-1}{2}}$

where the final inequality holds since $N\geq 1$ by assumption. Set the
failure probability $\delta$ to the right hand side as

$\delta\equiv\frac{\sqrt{2/\pi}}{\epsilon\Psi_{i}}\left(1+\epsilon^{2}\Psi_{i}^{2}\right)^{-\frac{N-1}{2}},$

and solve for $N$.

The larger the value of $\Psi_{i}$ (the same measure of diagonal dominance
that appears in Corollary 5.7), the smaller the sampling amount for the
($\epsilon$, $\delta$) estimator.

##  6 Application: Monte Carlo estimators for a derivative-based global
sensitivity metric

We bound the absolute error (Theorem 6.1) in a Monte Carlo estimator for
global sensitivity analysis, and more specifically for a derivative-based
global sensitivity metric (DGSM) of a function
$f:\mathbb{R}^{n}\rightarrow\mathbb{R}$ whose partial derivatives are square
integrable with respect to a probability density function
$\rho_{\boldsymbol{\mathcal{X}}}(\boldsymbol{x})$. The DGSM is equal to the
diagonal $\mathcal{D}\,(\boldsymbol{C})$ of the matrix

$\displaystyle\boldsymbol{C}=\int_{\mathcal{X}}\nabla f(\boldsymbol{x})[\nabla
f(\boldsymbol{x})]^{\top}\rho_{\boldsymbol{\mathcal{X}}}(\boldsymbol{x})d\boldsymbol{x}.$

The matrix $\boldsymbol{C}$ is well-defined, symmetric positive semidefinite,
and can be interpreted as a second moment matrix of the gradient. We compute
the DGSM with the Monte Carlo estimator

(29)
$\displaystyle\widehat{\boldsymbol{C}}\equiv\frac{1}{N}\sum_{k=1}^{N}\boldsymbol{z}_{k}\boldsymbol{z}_{k}^{\top}\qquad\text{where}\quad\boldsymbol{z}_{k}\equiv\nabla
f(\boldsymbol{x}_{k})$

and $\boldsymbol{x}_{k}$, $1\leq k\leq N$, are independent samples from the
distribution of $\rho_{\boldsymbol{\mathcal{X}}}(\boldsymbol{x})$. Below is a
normwise bound for the error in the DGSM computed by the Monte Carlo estimator
(29). Its derivation is related to the analysis in [12, Section 4].

######  Theorem 6.1.

Let $f:\mathbb{R}^{n}\rightarrow\mathbb{R}$ have square integrable partial
derivatives with respect to the probability density function
$\rho_{\boldsymbol{\mathcal{X}}}(\boldsymbol{x})$, $\|\nabla
f\|_{\infty}\leq\beta$ almost surely, $\widehat{\boldsymbol{C}}$ be the Monte
Carlo estimator in (29), and

$\displaystyle c_{\max}$
$\displaystyle\equiv\|\mathcal{D}\,(\boldsymbol{C})\|_{2}\qquad
S_{1}\equiv\|\mathcal{D}\,(\boldsymbol{C})\left(\beta^{2}\boldsymbol{I}-\mathcal{D}\,(\boldsymbol{C})\right)\|_{2}$
$\displaystyle S_{2}$ $\displaystyle\equiv c_{\max}+\beta^{2},\qquad
d\equiv\frac{\sum_{i=1}^{n}c_{ii}(\beta^{2}-c_{ii})}{S_{1}}.$

If $c_{\max}>0$ and $S_{1}>0$, then

######  Proof 6.2.

Before applying the matrix Bernstein inequality in Theorem 2.3, we need to
verify the assumptions. The Monte Carlo estimate
$\mathcal{D}\,(\widehat{\boldsymbol{C}})$ is an unbiased estimator of the DGSM
$\mathcal{D}\,(\boldsymbol{C})$ whose largest diagonal element is

$c_{\max}=\|\mathcal{D}\,(\boldsymbol{C})\|_{2}=\max_{1\leq i\leq
n}{|c_{ii}|}=\max_{1\leq i\leq n}\mathbb{E}\left[(\nabla
f(\boldsymbol{\mathcal{X}}))_{i}^{2}\right]\leq\beta^{2}.$

The absolute error in the DGSM computed by the Monte Carlo estimator (29) is

$\boldsymbol{Z}=\mathcal{D}\,(\boldsymbol{C})-\mathcal{D}\,(\widehat{\boldsymbol{C}})=\sum_{k=1}^{N}\boldsymbol{S}_{k},\qquad\text{where}\quad\boldsymbol{S}_{k}\equiv\frac{1}{N}\left(\mathcal{D}\,(\boldsymbol{C})-\mathcal{D}\,(\boldsymbol{z}_{k}\boldsymbol{z}_{k}^{\top})\right).$

The summands $\boldsymbol{S}_{k}$ have zero mean and are bounded by

$\displaystyle\|\boldsymbol{S}_{k}\|_{2}$
$\displaystyle\leq\frac{1}{N}\left(\|\mathcal{D}\,(\boldsymbol{C})\|_{2}+\|\mathcal{D}\,(\boldsymbol{z}_{k}\boldsymbol{z}_{k}^{\top})\|_{2}\right)$
$\displaystyle\leq\frac{\|\mathcal{D}\,(\boldsymbol{C})\|_{2}+\beta^{2}}{N}=\frac{c_{\max}+\beta^{2}}{N}=\frac{S_{2}}{N},\qquad
1\leq k\leq N.$

We let $L=S_{2}/N$ so that $\|\boldsymbol{S}_{k}\|_{2}\leq L$. The variance is

(30)
$\operatorname{\mathbb{V}\mathrm{ar}}[\boldsymbol{Z}]=\sum_{k=1}^{N}\mathbb{E}[\boldsymbol{S}_{k}^{2}]=\frac{1}{N^{2}}\sum_{k=1}^{N}\mathbb{E}\left[\left(\mathcal{D}\,(\boldsymbol{C})-\mathcal{D}\,(\boldsymbol{z}_{k}\boldsymbol{z}_{k}^{\top})\right)^{2}\right].$

Linearity of the expectation, the majorization
$\mathcal{D}\,(\boldsymbol{z}_{k}\boldsymbol{z}_{k}^{\top})\preceq\beta^{2}\boldsymbol{I}$,
$\mathbb{E}[\mathcal{D}\,(\boldsymbol{z}_{k}\boldsymbol{z}_{k}^{\top})]=\mathcal{D}\,(\boldsymbol{C})$
and commutativity of diagonal matrices imply for the summands,

$\displaystyle\mathbb{E}\left[\left(\mathcal{D}\,(\boldsymbol{C})-\mathcal{D}\,(\boldsymbol{z}_{k}\boldsymbol{z}_{k}^{\top})\right)^{2}\right]$
$\displaystyle=\mathbb{E}\left[\mathcal{D}\,(\boldsymbol{C})^{2}+(\mathcal{D}\,(\boldsymbol{z}_{k}\boldsymbol{z}_{k}^{\top}))^{2}-2\mathcal{D}\,(\boldsymbol{C})\,\mathcal{D}\,(\boldsymbol{z}_{k}\boldsymbol{z}_{k}^{\top})\right]$
$\displaystyle=\mathbb{E}\left[(\mathcal{D}\,(\boldsymbol{z}_{k}\boldsymbol{z}_{k}^{\top}))^{2}\right]-\mathcal{D}\,(\boldsymbol{C})^{2}$
$\displaystyle\preceq\beta^{2}\mathcal{D}\,(\boldsymbol{C})-\mathcal{D}\,(\boldsymbol{C})^{2},\qquad
1\leq k\leq N.$

Substitute the above into (30)

$\displaystyle\operatorname{\mathbb{V}\mathrm{ar}}[\boldsymbol{Z}]\preceq\boldsymbol{V}\equiv\frac{1}{N}\left(\beta^{2}\mathcal{D}\,(\boldsymbol{C})-\mathcal{D}\,(\boldsymbol{C})^{2}\right),$

and apply Theorem 2.3 with $\nu=\|\boldsymbol{V}\|_{2}=S_{1}/N$ and
$d=\operatorname{\mathrm{intdim}}(\boldsymbol{V})$.

Below is the minimal sampling amount that makes the Monte Carlo estimator in
(29) a normwise $(\epsilon,\delta)$ diagonal estimator.

######  Corollary 6.3.

Let $f:\mathbb{R}^{n}\rightarrow\mathbb{R}$ have square integrable partial
derivatives with respect to the probability density function
$\rho_{\boldsymbol{\mathcal{X}}}(\boldsymbol{x})$, $\|\nabla
f\|_{\infty}\leq\beta$ almost surely, $\widehat{\boldsymbol{C}}$ be the Monte
Carlo estimator in (29), and

$\displaystyle c_{\max}$
$\displaystyle\equiv\|\mathcal{D}\,(\boldsymbol{C})\|_{2},\qquad
S_{1}\equiv\|\mathcal{D}\,(\boldsymbol{C})\left(\beta^{2}\boldsymbol{I}-\mathcal{D}\,(\boldsymbol{C})\right)\|_{2}$
$\displaystyle S_{2}$ $\displaystyle\equiv c_{\max}+\beta^{2},\qquad
d\equiv\frac{\sum_{i=1}^{n}c_{ii}}{S_{1}}.$

Pick $\epsilon>0$. If $c_{\max}>0$ and $S_{1}>0$, then for any $0<\delta<1$,
if the sampling amount is at least

$N\geq\frac{S_{2}}{3\epsilon^{2}}\left(2\epsilon+\frac{6S_{1}}{c_{\max}\,S_{2}}\right)\ln(8d/\delta),$

then
$\|\mathcal{D}\,(\boldsymbol{C})-\mathcal{D}\,(\widehat{\boldsymbol{C}})\|_{2}\leq\epsilon\|\mathcal{D}\,(\boldsymbol{C})\|_{2}$
holds with probability at least $1-\delta$.

######  Proof 6.4.

The proof is similar to that of Corollary 3.6 but is based instead on Theorem
6.1.

###  6.1 Illustrative examples

We determine the constants in Corollary 6.3 for two different functions
$f:\mathbb{R}^{n}\rightarrow\mathbb{R}$, and random variables
$\boldsymbol{\mathcal{X}}\in\mathbb{R}^{n}$ from a uniform distribution over
$\mathcal{X}=[-1,1]^{n}$.

#### Linear Function

Let $f(\boldsymbol{x})=\boldsymbol{h}^{\top}\boldsymbol{x}$ with
$\boldsymbol{h}\in\mathbb{R}^{n}$. Then $\|\nabla
f\|_{\infty}\leq\beta\equiv\|\boldsymbol{h}\|_{\infty}$ almost surely. The
second moment matrix $\boldsymbol{C}=\boldsymbol{hh}^{\top}$ has a largest
diagonal entry $c_{\max}=\|\boldsymbol{h}\|_{\infty}^{2}$. Let
$\boldsymbol{v}\in\mathbb{R}^{n}$ have entries
$v_{i}=h_{i}^{2}\left(\|\boldsymbol{h}\|_{\infty}^{2}-h_{i}^{2}\right)$ for
$1\leq i\leq n$. The constants in Corollary 6.3 are

$S_{1}=\|\boldsymbol{v}\|_{\infty},\quad
S_{2}=2\|\boldsymbol{h}\|_{\infty}^{2},\quad
d=\frac{\sum_{i=1}^{n}v_{i}}{S_{1}}.$

#### Quadratic function

Let
$f(\boldsymbol{x})=\frac{1}{2}\boldsymbol{x}^{\top}\boldsymbol{S}\boldsymbol{x}$
where $\boldsymbol{S}\in\mathbb{R}^{n\times n}$
$\boldsymbol{S}=\boldsymbol{S}^{\top}$ is a symmetric square root of the
positive semidefinite matrix $\boldsymbol{M}=\boldsymbol{S}^{2}$. Then
$\|\nabla f\|_{\infty}\leq\beta\equiv\|\boldsymbol{S}\|_{\infty}$ almost
surely. The second moment matrix $\boldsymbol{C}=\frac{1}{3}\boldsymbol{M}$
has a largest diagonal element
$c_{\max}=\frac{1}{3}\|\mathcal{D}\,(\boldsymbol{M})\|_{\infty}$. Let
$\boldsymbol{v}\in\mathbb{R}^{n}$ be have entries
$v_{i}=\frac{1}{3}m_{ii}\left(\|\boldsymbol{S}\|_{\infty}^{2}-\frac{1}{3}m_{ii}\right)$
for $1\leq i\leq n$. The constants in Corollary 6.3 are

$S_{1}=\|\boldsymbol{v}\|_{\infty},\quad
S_{2}=\frac{1}{3}\,\|\mathcal{D}\,(\boldsymbol{M})\|_{\infty}+\|\boldsymbol{S}\|_{\infty}^{2},\quad
d=\frac{\sum_{i=1}^{n}v_{i}}{S_{1}}.$

##  7 Numerical Experiments

After describing our test matrices (section 7.1), we present four different
types of numerical experiments to illustrate the accuracy of the Monte Carlo
estimators: Rademacher Monte Carlo estimators applied to the test matrices
(section 7.2), accuracy of different Monte Carlo estimators (section 7.3),
effect of the sparsity on the accuracy of Rademacher Monte Carlo estimators
(section 7.4), and accuracy of the DGSM Monte Carlo estimator (section 7.5).

###  7.1 Test Matrices

We perform numerical experiments on three symmetric test matrices from [16] of
dimension $n=100$ that depend on a parameter $\theta$.

  1. 1.

Identity plus rank-1

$\displaystyle\boldsymbol{A}=\boldsymbol{I}+\theta\boldsymbol{e}\boldsymbol{e}^{\top}\qquad\text{where}\quad.01\leq\theta\leq
0.1,$

where $\boldsymbol{e}\in\mathbb{R}^{n}$ is a vector of ones. The constants in
Corollary 3.6 are

$K_{1}=(n-1)\theta^{2},\qquad
K_{2}=(n-1)\theta,\qquad\|\mathcal{D}\,(\boldsymbol{A})\|_{\infty}=1+\theta,$

so that

$\Delta_{1}=\frac{K_{1}}{(1+\theta)^{2}},\qquad\Delta_{2}=\frac{(n-1)\theta}{1+\theta},\qquad
d=n.$

  2. 2.

Rank-1 with decaying elements

$\boldsymbol{A}=\frac{\boldsymbol{xx}^{\top}}{\|\boldsymbol{x}\|_{2}^{2}}\qquad\text{where}\quad\boldsymbol{x}_{j}=e^{-j(1-\theta)},\quad
1\leq j\leq n,\qquad 0.1\leq\theta\leq 1.$

The constants in Corollary 3.6 are

$\displaystyle K_{1}$
$\displaystyle=\left(\frac{x_{1}}{\|\boldsymbol{x}\|_{2}}\right)^{2}\left(1-\left(\frac{x_{1}}{\|\boldsymbol{x}\|_{2}}\right)^{2}\right),\qquad
K_{2}=\frac{x_{1}}{\|\boldsymbol{x}\|_{2}^{2}}\sum_{j>1}x_{j}$

and
$\|\mathcal{D}\,(\boldsymbol{A})\|_{\infty}=\left(\frac{x_{1}}{\|\boldsymbol{x}\|_{2}}\right)^{2}$so
that

$\Delta_{1}=\left(\frac{\|\boldsymbol{x}\|_{2}}{x_{1}}\right)^{2}-1,\qquad\Delta_{2}=\frac{\sum_{j>1}x_{j}}{x_{1}},\qquad
d=\frac{\sum_{i=1}^{n}x_{i}^{2}\left(\|\boldsymbol{x}\|_{2}^{2}-x_{i}^{2}\right)}{x_{1}^{2}(\|\boldsymbol{x}\|_{2}^{2}-x_{1}^{2})}.$

  3. 3.

Tridiagonal Toeplitz matrix

$\boldsymbol{A}=\begin{bmatrix}1&\theta\\\ \theta&1&\ddots\\\
&\ddots&\ddots&\theta\\\ &&\theta&1\end{bmatrix}\qquad\text{where}\quad
0.1\leq\theta\leq 1.$

The constants in Corollary 3.6 are

$\displaystyle K_{1}=2\theta^{2},\qquad
K_{2}=2\theta,\qquad\|\mathcal{D}\,(\boldsymbol{A})\|_{\infty}=1$

so that

$\Delta_{1}=2\theta^{2},\qquad\Delta_{2}=2\theta,\qquad
d=\frac{2(n-1)\theta^{2}}{2\theta^{2}}=(n-1).$

For all the test matrices, the constants $\Delta_{1}$ and $\Delta_{2}$
increase with increasing $\theta$ as the offdiagonal elements become larger in
magnitude relative to the diagonal elements. Therefore, we expect the
Rademacher Monte Carlo estimators to lose accuracy with increasing $\theta$,
as measured by the normwise relative error (NRE) in the computed diagonal
$\mathcal{D}\,(\widehat{\boldsymbol{A}})$,

$\mathrm{NRE}\equiv\frac{\|\mathcal{D}\,(\boldsymbol{A})-\mathcal{D}\,(\widehat{\boldsymbol{A}})\|_{2}}{\|\mathcal{D}\,(\boldsymbol{A})\|_{2}},$

in Figures 1-5.

Figure 1: Rademacher Monte Carlo estimator applied to Test Matrix 1. Big left
panel: NRE for different values of $\theta$ versus sampling amount $N$. Small
panels on the right: NRE (solid black line), and bound (31) (blue dotted line)
versus sampling amount $N$ with failure probability $\delta=10^{-16}$.

###  7.2 Experiment 1: Accuracy of Rademacher Monte Carlo estimator on test
matrices

Figures 1-3 show the NRE of the Rademacher Monte Carlo estimator applied to
the test matrices in section 7.1, and the bounds from the normwise
$(\epsilon,\delta)$ estimators in Corollary 3.6.

The big left panel displays the NRE versus the sampling amount $N$. This NRE
represents the average of the NREs over $10$ different independent runs. The
small panels on the right show the bound $\epsilon$ for the normwise
$(\epsilon,\delta)$ estimators from Corollary 3.6 with failure probability
$\delta=10^{-16}$.

For Corollary 3.6, we solve for $\epsilon$ from the simpler bound

$N\geq\frac{\Delta_{2}}{3\epsilon^{2}}\left(2+6\Delta_{3}\right)\ln(8d/\delta),\qquad\Delta_{3}\equiv\frac{\Delta_{1}}{\Delta_{2}},$

to obtain

(31)
$\epsilon=\sqrt{\frac{\Delta_{2}}{3N}\left(2+6\Delta_{3}\right)\ln(8d/\delta)},\qquad\Delta_{2}=\frac{\|\boldsymbol{A}-\mathcal{D}\,(\boldsymbol{A})\|_{\infty}}{\|\mathcal{D}\,(\boldsymbol{A})\|_{\infty}},\quad\Delta_{3}=\frac{K_{1}}{\|\boldsymbol{A}-\mathcal{D}\,(\boldsymbol{A})\|_{\infty}}.$

The big left panels illustrate that, for a fixed sampling amount $N$, the NRE
for Test Matrices 1 and 3 increases with $\theta$. This is because the
offdiagonals become more dominant as $\theta$ becomes larger.

Figure 2: Rademacher Monte Carlo estimator applied to Test Matrix 2. Big left
panel: NRE for different values of $\theta$ versus sampling amount $N$. Small
panels on the right: NRE (solid black line) and bound (31) (blue dotted line)
versus sampling amount $N$ with failure probability $\delta=10^{-16}$.

Figure 3: Rademacher Monte Carlo estimator applied to Test Matrix 3. Big left
panel: NRE for different values of $\theta$ versus sampling amount $N$. Small
panels on the right: NRE (solid black line) and bound (31) (blue dotted line)
versus sampling amount $N$ with failure probability $\delta=10^{-16}$.

###  7.3 Experiment 2: Different Monte Carlo estimators

We compare the accuracy of the following Monte Carlo estimators on Test Matrix
1 with $\theta=0.01$: Rademacher, Gaussian, sparse Rademacher with $s=3$, and
normalized Gaussian.

For each estimator, Figure 4 shows the mean of the NRE and variance over 100
runs, with the shaded regions representing the $2.5\%$ and $97.5\%$ quantiles.

The normalized Gaussian estimator is about as accurate as the Rademacher
estimator, while the sparse Rademacher with $s=3$ is about as accurate as the
Gaussian estimator. The Gaussian and sparse Rademacher estimators are less
accurate than the Rademacher and normalized Gaussian estimators. The shaded
regions illustrate that, as expected, the sample variance of all estimators
decreases with increasing sampling amount $N$.

Figure 4: Rademacher, Gaussian, sparse Rademacher with $s=3$, and Normalized
Gaussian Monte Carlo estimators applied to Test Matrix 1 with $\theta=0.01$.
Big left panel: NRE mean versus sampling amount $N$ for different estimators.
Small right panels: NRE mean (styled lines), and $2.5\%$ and $97.5\%$
quantiles (shaded regions) versus sampling amount $N$.

###  7.4 Experiment 3: Effect of sparsity in Rademacher vectors

We apply the Rademacher Monte Carlo estimator to Test Matrix 1 with
$\theta=.01$ with four different sparsity levels: $s=1$ (standard Rademacher),
$s=3$ [1], $s=10$, and $s=50$.

For each sampling amount $N$, Figure 5 shows the mean and the variance of the
NRE over 100 runs. It suggests that sparse Rademacher estimators ($s>1$) may
not be able to achieve a single digit of accuracy, unless the sampling amount
is so large as to exceed the matrix dimension.

Figure 5: Sparse Rademacher Monte Carlo estimators with sparsity levels
$s=1,3,10,50$ applied to Test Matrix 1 with $\theta=0.01$. NRE (dotted lines)
and $2.5\%$ and $97.5\%$ quantiles (shaded regions).

###  7.5 Example 4: Bounds for DGSM Monte Carlo estimator

We apply the DGSM Monte Carlo estimator (29) to the diagonal matrix

(32)
$\boldsymbol{S}\equiv\operatorname{\mathrm{diag}}(\boldsymbol{s})\in\mathbb{R}^{n\times
n}\qquad s_{j}\equiv\exp(-10j/n),\qquad 1\leq j\leq n$

from the quadratic function in Section 6.1 for $n=100$, and illustrate the
accuracy of Corollary 6.3.

The left panel of Figure 6 shows the normwise relative error

$\mathrm{NRE}\equiv\frac{\|\mathcal{D}\,(\boldsymbol{C})-\mathcal{D}\,(\boldsymbol{\widehat{C}})\|_{2}}{\|\mathcal{D}\,(\boldsymbol{C})\|_{2}}$

which represents the the average of the NREs over $100$ independent runs.

For Corollary 6.3, we fix the sample size $N$ and solve for $\epsilon$ from
the simpler bound

$N\geq\frac{S_{2}}{3\epsilon^{2}}\left(2+6S_{3}\right)\ln(8d/\delta),\qquad
S_{3}\equiv\frac{S_{1}}{c_{\max}S_{2}},$

to obtain

(33) $\epsilon=\sqrt{\frac{S_{2}}{3N}\left(2+6S_{3}\right)\ln(8d/\delta)}.$

The expressions for $S_{1},S_{2},c_{\max}$ and $d$ for this example have been
derived in subsection 6.1.

Figure 6: DGSM Monte Carlo estimator (29) applied to $100\times 100$ matrix
$\boldsymbol{S}$ in (32). Left panel: NRE mean (solid line) and $2.5\%$ and
$97.5\%$ quantiles (shaded regions) versus sampling amount $N$. Right panel:
NRE and bounds (33) for different failure probabilities $\delta$ versus
sampling amount.

The right panel of Figure 6 illustrates that with less stringent failure
probabilities $\delta$, the relative bounds (33) move closer to the NRE.

###  7.6 Example 5: DGSM on the Circuit model

We apply the Monte Carlo DGSM estimator (29) to the so-called circuit model
from [6]. The quantity of interest being modeled is the midpoint voltage of a
transformerless push-pull circuit, which depends on $n=6$ parameters through a
nonlinear closed-form algebraic expression.

As in [6], we normalize the parameter space to $\mathcal{X}=[-1,1]^{n}$ and
scale the partial derivatives appropriately111MATLAB codes are available in
https://bitbucket.org/paulcon/global-sensitivity-metrics-from-active-
subspaces/src/master/..

Figure 7: DGSM Monte Carlo estimator (29) applied to the Circuit model. NRE
mean (solid line) and $2.5\%$ and $97.5\%$ quantiles (shaded regions) versus
sampling amount $N$.

Figure 7 shows the mean of the NRE and variances over $100$ independent runs
per sampling amount $N$, with the shaded regions representing $2.5\%$ and
$97.5\%$ quantiles. Since the exact expressions for the DGSMs are unavailable,
we use as the exact value a tensor product Gauss-Legendre quadrature-based
approximation with $15$ points per dimension (i.e., $15^{n}$ total points).

##  8 Conclusion and future work

This paper derives probabilistic bounds for the Monte Carlo diagonal
estimators; the bounds come in two flavors, normwise and componentwise bounds
for the absolute and relative errors of the Monte Carlo estimators. There are
several avenues for future work. First, it would be interesting to study the
accuracy of the diagonal estimator for matrix functions when a polynomial, or
rational approximation to the matrix function is used. Second, we are
currently pursuing the extension of the analysis of the Monte Carlo diagonal
estimators to estimators for the selected entries (possibly including
offdiagonal entries) of a matrix.

##  9 Acknowledgements

The authors would like to acknowledge support from the National Science
Foundation through the grants DMS-1745654 (all three authors) and DMS-1845406
(A.K.S.). We would like to acknowledge Alen Alexanderian for helpful
discussions.

## References

  * [1] D. Achlioptas.  Database-friendly random projections: Johnson-Lindenstrauss with binary coins.  J. Comput. System Sci., 66(4):671–687, 2003.
  * [2] H. Avron and S. Toledo.  Randomized algorithms for estimating the trace of an implicit symmetric positive semi-definite matrix.  J. ACM, 58(2):1–34, 2011.
  * [3] R. A. Baston and Y. Nakatsukasa.  Stochastic diagonal estimation: probabilistic bounds and an improved algorithm, 2022.  arXiv:2201.10684.
  * [4] C. Bekas, E. Kokiopoulou, and Y. Saad.  An estimator for the diagonal of a matrix.  Appl. Numer. Math., 57(11-12):1214–1229, 2007.
  * [5] R. Y. Chen, A. Gittens, and J. A. Tropp.  The masked sample covariance estimator: an analysis using matrix concentration inequalities.  Inf. Inference, 1(1):2–20, 2012.
  * [6] P. G. Constantine and P. Diaz.  Global sensitivity metrics from active subspaces.  Reliab. Eng. Syst. Safe., 162:1–13, 2017.
  * [7] A. Cortinovis and D. Kressner.  On randomized trace estimates for indefinite matrices with an application to determinants.  Found. Comput. Math., pages 1–29, 2021.
  * [8] M. F. Hutchinson.  A stochastic estimator of the trace of the influence matrix for Laplacian smoothing splines.  Comm. Statist. Simulation Comput., 18(3):1059–1076, 1989.
  * [9] B. J. Kaperick.  Diagonal estimation with probing methods.  Master’s thesis, Virginia Polytechnic Institute and State University, 2019\.
  * [10] S. Kucherenko and B. Iooss.  Derivative-Based Global Sensitivity Measures, pages 1241–1263.  Springer International Publishing, Cham, 2017.
  * [11] J. H. Laeuchli.  Methods for Estimating The Diagonal of Matrix Functions.  PhD thesis, College of William and Mary, 2016.
  * [12] R. R. Lam, O. Zahm, Y. M. Marzouk, and K. E. Willcox.  Multifidelity dimension reduction via active subspaces.  SIAM J. Sci. Comput., 42(2):A929–A956, 2020.
  * [13] J. L. Larsen and L. M. Morris.  An Introduction to Mathematical Statistics and its Applications.  Pearson Prentice Hall, fourth edition, 2006.
  * [14] P. Li, T. J. Hastie, and K. W. Church.  Very sparse random projections.  In Proceedings of the 12th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, KDD ’06, page 287–296, New York, NY, USA, 2006. Association for Computing Machinery.
  * [15] M. Mitzenmacher and E. Upfal.  Probability and computing.  Cambridge University Press, Cambridge, 2005.  Randomized algorithms and probabilistic analysis.
  * [16] F. Roosta-Khorasani and U. Ascher.  Improved bounds on sample size for implicit matrix trace estimators.  Found. Comput. Math., 15(5):1187–1212, 2015.
  * [17] A. P. Soms.  An asymptotic expansion for the tail area of the t-distribution.  J. Amer. Statist. Assoc., 71(355):728–730, 1976.
  * [18] J. A. Tropp.  An introduction to matrix concentration inequalities.  Found. Trends Mach. Learning, 8(1–2):1–230, 2015.
  * [19] S. Ubaru and Y. Saad.  Applications of trace estimation techniques.  In T. Kozubek, M. Čermák, P. Tichý, R. Blaheta, J. Šístek, D. Lukáš, and J. Jaroš, editors, High Performance Computing in Science and Engineering, pages 19–33, Cham, 2018\. Springer International Publishing.
  * [20] R. Vershynin.  High-dimensional probability: An introduction with applications in data science, volume 47.  Cambridge University Press, 2018.
  * [21] A. J. Wathen.  Preconditioning.  Acta Numer., 24:329–376, 2015.
  * [22] Z. Yao, A. Gholami, S. Shen, M. Mustafa, K. Keutzer, and M. W. Mahoney.  ADAHESSIAN: An adaptive second order optimizer for machine learning.  arXiv preprint arXiv:2006.00719, 2020.