Don’t Use Gaussians in Stochastic Trace Estimation

January 28, 2024 by Ethan N. Epperly Leave a comment

Suppose we are interested in estimating the trace $\tr(A) = \sum_{i=1}^n A_{ii}$ of an $n\times n$ matrix $A$ that can be only accessed through matrix–vector products $Ax_1,\ldots,Ax_m$ . The classical method for this purpose is the Girard–Hutchinson estimator

$\hat{\tr} = \frac{1}{m} \left( x_1^\top Ax_1 + \cdots + x_m^\top Ax_m \right),$

where the vectors $x_1,\ldots,x_m$ are independent, identically distributed (iid) random vectors satisfying the isotropy condition

$\expect[x_ix_i^\top] = I.$

Examples of vectors satisfying this condition include

Gaussian: The entries of each $x_i$ are iid standard Gaussian random variables.
Random signs: The entries of each $x_i$ are iid random numbers taking the values $\pm 1$ with equal probabilities (i.e., Rademacher random variables).
Sphere: Each $x_i$ is a uniformly random point on the (Euclidean) sphere of radius $\sqrt{n}$ , $x_i \sim \operatorname{Unif} \{ y \in \mathbb{R}^n : y^\top y = n \}$ .

Stochastic trace estimation has a number of applications: log-determinant computations in machine learning, partition function calculations in statistical physics, generalized cross validation for smoothing splines, and triangle counting in large networks. Several improvements to the basic Girard–Hutchinson estimator have been developed recently. I am partial to XTrace, an improved trace estimator that I developed with my collaborators.

This post is addressed at the question:

Which distribution should be used for the test vectors $x_i$ for stochastic trace estimation?

Since the Girard–Hutchinson estimator is unbiased $\expect[\hat{\tr}] = \tr(A)$ , the variance of $\hat{\tr}$ is equal to the mean-square error. Thus, the lowest variance trace estimate is the most accurate. In my previous post on trace estimation, I discussed formulas for the variance $\Var(\hat{\tr})$ of the Girard–Hutchinson estimator with different choices of test vectors. In that post, I stated the formulas for different choices of test vectors (Gaussian, random signs, sphere) and showed how those formulas could be proven.

In this post, I will take the opportunity to editorialize on which distribution to pick. The thesis of this post is as follows:

The sphere distribution is essentially always preferable to the Gaussian distribution for trace estimation.

To explain why, let’s focus on the case when $A$ is real and symmetric.¹ Let $\lambda_1,\ldots,\lambda_n$ be the eigenvalues of $A$ and define the eigenvalue mean

$\overline{\lambda} = \frac{\lambda_1 + \cdots + \lambda_n}{n}.$

Then the variance of the Girard–Hutchinson estimator with Gaussian vectors $x_i$ is

$\Var(\hat{\tr}_{\rm Gaussian}) = \frac{1}{m} \cdot 2 \sum_{i=1}^n \lambda_i^2.$

For vectors $x_i$ drawn from the sphere, we have

$\Var(\hat{\tr}_{\rm sphere}) = \frac{1}{m} \cdot \frac{n}{n+2} \cdot 2\sum_{i=1}^n (\lambda_i - \overline{\lambda})^2.$

The sphere distribution improves on the Gaussian distribution in two ways. First, the variance of $\Var(\hat{\tr}_{\rm sphere})$ is smaller than $\Var(\hat{\tr}_{\rm Gaussian})$ by a factor of $n/(n+2) < 1$ . This improvement is quite minor. Second, and more importantly, $\Var(\hat{\tr}_{\rm Gaussian})$ is proportional to the sum of $A$ ‘s squared eigenvalues whereas $\Var(\hat{\tr}_{\rm sphere})$ is proportional to the sum of $A$ ‘s squared eigenvalues after having been shifted to be mean-zero!

The difference between Gaussian and sphere test vectors can be large. To see this, consider a $1000\times 1000$ matrix $A$ with eigenvalues uniformly distributed between $0.9$ and $1.1$ with a (Haar orthgonal) random matrix of eigenvectors. For simplicity, since the variance of all Girard–Hutchinson estimates is proportional to $1/m$ , we take $m=1$ . Below show the variance of Girard–Hutchinson estimator for different distributions for the test vector. We see that the sphere distribution leads to a trace estimate which has a variance 300× smaller than the Gaussian distribution. For this example, the sphere and random sign distributions are similar.

Which Distribution Should You Use: Signs vs. Sphere

The main point of this post is to argue against using the Gaussian distribution. But which distribution should you use: Random signs? The sphere distribution? The answer, for most applications, is one of those two, but exactly which depends on the properties of the matrix $A$ .

The variance of the Girard–Hutchinson estimator with the random signs estimator is

$\Var(\hat{\tr}_{\rm signs}) = 2 \sum_{i\ne j} A_{ij}^2.$

Thus, $\Var(\hat{\tr}_{\rm signs})$ depends on the size of the off-diagonal entries of $A$ ; $\Var(\hat{\tr}_{\rm signs})$ does not depend on the diagonal of $A$ at all! For matrices with small off-diagonal entries (such as diagonally dominant matrices), the random signs distribution is often the best.

However, for other problems, the sphere distribution is preferable to random signs. The sphere distribution is rotation-invariant, so $\Var(\hat{\tr}_{\rm sphere})$ is independent of the eigenvectors of the (symmetric) matrix $A$ , depending only on $A$ ‘s eigenvalues. By contrast, the variance of the Girard–Hutchinson estimator with the random signs distribution can significantly depend on the eigenvectors of the matrix $A$ . For a given set of eigenvalues and the worst-case choice of eigenvectors, $\Var(\hat{\tr}_{\rm sphere})$ will always be smaller than $\Var(\hat{\tr}_{\rm signs})$ . In fact, $\Var(\hat{\tr}_{\rm sphere})$ is the minimum variance distribution for Girard–Hutchinson trace estimation for a matrix with fixed eigenvalues and worst-case eigenvectors; see this section of my previous post for details.

In my experience, random signs and the sphere distribution are both perfectly adequate for trace estimation and either is a sensible default if you’re developing software. The Gaussian distribution on the other hand… don’t use it unless you have a good reason to.

Don’t Solve the Normal Equations

July 26, 2022 by Ethan N. Epperly 1 Comment

The (ordinary) linear least squares problem is as follows: given an $m\times n$ matrix $A$ and a vector $b$ of length $m$ , find the vector $x$ such that $Ax$ is as close to $b$ as possible, when measured using the two-norm $\| \cdot \|$ . That is, we seek to

(1) $\begin{equation*} \mbox{find } x \in\mathbb{R}^n \mbox{ such that }\| b - Ax \|^2 = \sum_{i=1}^m \left(b_i - \sum_{j=1}^n A_{ij} x_j \right)^2 \mbox{ is minimized}. \end{equation*}$

From this equation, the name “least squares” is self-explanatory: we seek $x$ which minimizes the sum of the squared discrepancies between the entries of $b$ and $Ax$ .

The least squares problem is ubiquitous in science, engineering, mathematics, and statistics. If we think of each row $a_i$ of $A$ as an input and its corresponding entry $b_i$ of $b$ as an output, then the solution $x$ to the least squares model gives the coefficients of a linear model for the input–output relationship. Given a new previously unseen input $a_{\rm new}$ , our model predicts the output $b_{\rm new}$ is approximately $b_{\rm new} \approx a_{\rm new}^\top x = \sum_{i=1}^n x_i (a_{\rm new})_i$ . The vector $x$ consists of coefficients for this linear model. The least squares solution satisfies the property that the average squared difference between the output $b_i$ and the prediction $a_i^\top x$ is as small as it could possibly be for all choices of coefficient vectors $x$ .

How do we solve the least squares problem? A classical solution approach, ubiquitous in textbooks, is to solve a system of linear equations known as the normal equations. The normal equations associated with the least squares problem (1) are given by

(2) $\begin{equation*} A^\top A \,x = A^\top b. \end{equation*}$

This system of equations always has a solution. If $A$ has full column-rank, then $A^\top A$ is invertible and the unique least squares solution to (1) is given by $(A^\top A)^{-1} A^\top b$ . We assume that $A$ has full column-rankQ for the rest of this discussion. To solve the normal equations in software, we compute $A^\top A$ and $A^\top b$ and solve (2) using a linear solver like MATLAB’s “\”.¹ (As is generally true in matrix computations, it is almost never a good idea to explicitly form the inverse of the matrix $A^\top A$ , or indeed any matrix.) We also can solve the normal equations using an iterative method like (preconditioned) conjugate gradient.

The purpose of the article is to advocate against the use of the normal equations for solving the least squares problems, at least in most cases. So what’s wrong with the normal equations? The problem is not that the normal equations aren’t mathematically correct. Instead, the problem is that the normal equations often lead to poor accuracy for the least squares solution using computer arithmetic.

Most of the time when using computers, we store real numbers as floating point numbers.² At a coarse level, the right model to have in your head is that real numbers on a computer are stored in scientific notation with only 16 decimal digits after the decimal point.³ When two numbers are added, subtracted, multiplied, and divided, the answer is computed and then rounded to 16 decimal digits; any extra digits of information are thrown away. Thus, the result of our arithmetic on a computer is the true answer to the arithmetic problem plus a small rounding error. These rounding errors are small individually, but solving an even modestly sized linear algebra problem requires thousands of such operations. Making sure many small errors don’t pile up into a big error is part of the subtle art of numerical computation.

To make a gross simplification, if one solves a system of linear equations $Mx = c$ on a computer using a well-designed piece of software, one obtains an approximate solution $\hat{x}$ which is, after accounting for the accumulation of rounding errors, close to $x$ . But just how close the computed solution $\hat{x}$ and the true solution $x$ are depends on how “nice” the matrix $M$ is. The “niceness” of a matrix $M$ is quantified by a quantity known as the condition number of $M$ , which we denote $\kappa(M)$ .⁴ As a rough rule of thumb, the relative error between $x$ and $\hat{x}$ is roughly bounded as

(3) $\begin{equation*} \frac{\| \hat{x} - x \|}{\|x\|} \lessapprox \kappa(M)\times 10^{-16}. \end{equation*}$

The “ $10^{-16}$ corresponds to the fact we have roughly 16 decimal digits of accuracy in double precision floating point arithmetic. Thus, if the condition number of $M$ is roughly $10^{10}$ , then we should expect around $6$ digits of accuracy in our computed solution.

The accuracy of the least squares problem is governed by its own condition number $\kappa(A)$ . We would hope that we can solve the least squares problem with an accuracy like the rule-of-thumb error bound (3) we had for linear systems of equations, namely a bound like $\|\hat{x} - x\|/\|x\| \lessapprox \kappa(A)\times 10^{-16}$ . But this is not the kind of accuracy we get for the least squares problem when we solve it using the normal equations. Instead, we get accuracy like

(4) $\begin{equation*} \frac{\| \hat{x} - x \|}{\|x\|} \lessapprox \left(\kappa(A)\right)^2\times 10^{-16}. \end{equation*}$

By solving the normal equations we effectively square the condition number! Perhaps this is not surprising as the normal equations also more-or-less square the matrix $A$ by computing $A^\top A$ . This squared condition number drastically effects the accuracy of the computed solution. If the condition number of $A$ is $10^{8}$ , then the normal equations give us absolute nonsense for $\hat{x}$ ; we expect to get no digits of the answer $x$ correct. Contrast this to above, where we were able to get $6$ correct digits in the solution to $Mx = c$ despite the condition number of $M$ being $100$ times larger than $A$ !

All of this would be just a sad fact of life for the least squares problem if the normal equations and their poor accuracy properties were the best we could do for the least squares problem. But we can do better! One can solve linear least squares problems by computing a so-called QR factorization of the matrix $A$ .⁵ Without going into details, the upshot is that the least squares solution by a well-designed⁶ QR factorization requires a similar amount of time to solving the normal equations and has dramatically improved accuracy properties, achieving the desirable rule-of-thumb behavior⁷

(5) $\begin{equation*} \frac{\| \hat{x} - x \|}{\|x\|} \lessapprox \kappa(A)\times 10^{-16}. \end{equation*}$

I have not described how the QR factorization is accurately computed nor how to use the QR factorization to solve least squares problems nor even what the QR factorization is. All of these topics are explained excellently by the standard textbooks in this area, as well as by publicly available resources like Wikipedia. There’s much more that can be said about the many benefits of solving the least squares problem with the QR factorization,⁸ but in the interest of brevity let me just say this: TL;DR when presented in the wild with a least squares problem, the solution method one should default to is one based on a well-implemented QR factorization, not solving the normal equations.

Suppose for whatever reason we don’t have a high quality QR factorization algorithm at our disposal. Must we then resort to the normal equations? Even in this case, there is a way we can reduce the problem of solving a least squares problems to a linear system of equations without squaring the condition number! (For those interested, to do this, we recognize the normal equations as a Schur complement of a somewhat larger system of linear equations and then solve that. See Eq. (7) in this post for more discussion of this approach.)

The title of this post Don’t Solve the Normal Equations is deliberately overstated. There are times when solving the normal equations is appropriate. If $A$ is well-conditioned with a small condition number, squaring the condition number might not be that bad. If the matrix $A$ is too large to store in memory, one might want to solve the least squares problem using the normal equations and the conjugate gradient method.

However, the dramatically reduced accuracy of solving the normal equations should disqualify the approach from being the de-facto way of solving least squares problems. Unless you have good reason to think otherwise, when you see $A^\top A$ , solve a different way.

The Better Way to Convert an SVD into a Symmetric Eigenvalue Problem

March 18, 2021 by Ethan N. Epperly 6 Comments

A singular value decomposition of an $m\times n$ matrix $B$ is a factorization of the form $B = U\Sigma V^\top$ , where $U$ and $V$ are square, orthogonal matrices and $\Sigma$ is a diagonal matrix with $(i,i)$ th entry $\sigma_i \ge 0$ .¹ The diagonal entries of $\Sigma$ are referred to as the singular values of $B$ and are conventionally ordered $\sigma_{\rm max} = \sigma_1 \ge \sigma_2 \ge \cdots \ge \sigma_{\min(m,n)} = \sigma_{\rm min}$ . The columns of the matrices $U$ and $V$ are referred to as the right- and left- singular vectors of $B$ and satisfy the relations $Bv_i = \sigma_i u_i$ and $B^\top u_i = \sigma_i v_i$ .

One can obtain the singular values and right and left singular vectors of $B$ from the eigenvalues and eigenvectors of $B^\top B$ and $BB^\top$ . This follows from the calculations $B^\top B = V\Sigma^2 V^\top$ and $B^\top B = U\Sigma^2 U^\top$ . In other words, the nonzero singular values of $B$ are the square roots of the nonzero eigenvalues of $B^\top B$ and $BB^\top$ . If one merely solves one of these problems, computing $\Sigma$ along with $U$ or $V$ , one can obtain the other matrix $V$ or $U$ by computing $U = BV \Sigma^{-1}$ or $V = B^\top U \Sigma^{-1}$ . (These formulas are valid for invertible square matrices $B$ , but similar formulas hold for singular or rectangular $B$ to compute the singular vectors with nonzero singular values.)

This approach is often unundesirable for several reasons. Here are a few I’m aware of:

Accuracy: Roughly speaking, in double-precision arithmetic, accurate stable numerical methods can resolve differences on the order of 16 orders of magnitude. This means an accurately computed SVD of $B$ can resolve the roughly 16 orders of magnitude of decaying singular values, with singular values smaller than that difficult to compute accurately. By computing $B^\top B$ , we square all of our singular values, so resolving 16 orders of magnitude of the eigenvalues of $B^\top B$ means we only resolve 8 orders of magnitude of the singular values of $B$ .² The dynamic range of our numerical computations has been cut in half!
Loss of orthogonality: While $U = BV \Sigma^{-1}$ and $V = B^\top U \Sigma^{-1}$ are valid formulas in exact arithmetic, they fair poorly when implemented numerically. Specifically, the numerically computed values $U_{\rm numerical}$ and $V_{\rm numerical}$ may not be orthogonal matrices with, for example, $U_{\rm numerical}^\top U_{\rm numerical}$ not even close to the identity matrix. One can, of course, orthogonalize the computed $U$ or $V$ , but this doesn’t fix the underlying problem that $U$ or $V$ have not been computed accurately.
Loss of structure: If $B$ possesses additional structure (e.g. sparsity), this structure may be lost or reduced by computing the product $B^\top B$ .
Nonlinearity: Even if we’re not actually computing the SVD numerically but doing analysis with pencil and paper, finding the SVD of $B$ from $B^\top B$ has the disadvantage of performing a nonlinear transformation on $B$ . This prevents us from utilizing additive perturbation theorems for sums of symmetric matrices in our analysis.³

There are times where these problems are insignificant and this approach is sensible: we shall return to this point in a bit. However, these problems should disqualify this approach from being the de facto way we reduce SVD computation to a symmetric eigenvalue problem. This is especially true since we have a better way.

The better way is by constructing the so-called Hermitian dilation⁴ of $B$ , which is defined to be the matrix

(1) $\begin{equation*} \mathcal{H}(B) = \begin{bmatrix} 0 & B \\ B^\top & 0 \end{bmatrix}. \end{equation*}$

One can show that the nonzero eigenvalues of $\mathcal{H}(B)$ are precisely plus-or-minus the singular values of $B$ . More specifically, we have

(2) $\begin{equation*} \mathcal{H}(B) \begin{bmatrix} u_i \\ \pm v_i \end{bmatrix} = \pm \sigma_i \begin{bmatrix} u_i \\ \pm v_i \end{bmatrix}. \end{equation*}$

All of the remaining eigenvalues of $\mathcal{H}(B)$ not of this form are zero.⁵ Thus, the singular value decomposition of $B$ is entirely encoded in the eigenvalue decomposition of $\mathcal{H}(B)$ .

This approach of using the Hermitian dilation $\mathcal{H}(B)$ to compute the SVD of $B$ fixes all the issues identified with the “ $B^\top B$ ” approach. We are able to accurately resolve a full 16 orders of magnitude of singular values. The computed singular vectors are accurate and numerically orthogonal provided we use an accurate method for the symmetric eigenvalue problem. The Hermitian dilation $\mathcal{H}(B)$ preserves important structural characteristics in $B$ like sparsity. For purposes of theoretical analysis, the mapping $B \mapsto \mathcal{H}(B)$ is linear.⁶

Often one can work with the Hermitian dilation only implicitly: the matrix $\mathcal{H}(B)$ need not actually be stored in memory with all its extra zeros. The programmer designs and implements an algorithm with $\mathcal{H}(B)$ in mind, but deals with the matrix $B$ directly for their computations. In a pinch, however, forming $\mathcal{H}(B)$ directly in software and utilizing symmetric eigenvalue routines directly is often not too much less efficient than a dedicated SVD routine and can cut down on programmer effort significantly.

As with all things in life, there’s no free lunch here. There are a couple of downsides to the Hermitian dilation approach. First, $\mathcal{H}(B)$ is, except for the trivial case $B = 0$ , an indefinite symmetric matrix. By constast, $B^\top B$ and $BB^\top$ are positive semidefinite, which can be helpful in some contexts.⁷ Further, if $n\ll m$ (respectively, $m \ll n$ ), then $B^\top B$ (respectively, $BB^\top$ ) is tiny compared to $\mathcal{H}(B)$ , so it might be considerably cheaper to compute an eigenvalue decomposition of $B^\top B$ (or $BB^\top$ ) than $\mathcal{H}(B)$ .

Despite the somewhat salacious title of this article, the $B^\top B$ and Hermitian dilation approaches both have their role, and the purpose of this article is not to say the $B^\top B$ approach should be thrown in the dustbin. However, in my experience, I frequently hear the $B^\top B$ approach stated as the definitive way of converting an SVD into an eigenvalue problem, with the Hermitian dilation approach not even mentioned. This, in my opinion, is backwards. For accuracy reasons alone, the Hermitian dilation should be the go-to tool for turning SVDs into symmetric eigenvalue problems, with the $B^\top B$ approach only used when the problem is known to have singular values which don’t span many orders of magnitude or $B$ is tall and skinny and the computational cost savings of the $B^\top B$ approach are vital.

Distribution	Variance (divided by $\tr(A)^2$ )
Gaussian	$2.0\times 10^{-3}$
Sphere	$6.7\times 10^{-6}$
Random signs	$6.7\times 10^{-6}$