Don’t Use Gaussians in Stochastic Trace Estimation

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 GirardHutchinson 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

Stochastic trace estimation has a number of applications: log-determinant computations in machine learningpartition function calculations in statistical physicsgeneralized 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.1The same principles hold in the general case, but the variance formulas are more delicate to state. See my previous post for the formulas. 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.

DistributionVariance (divided by \tr(A)^2)
Gaussian2.0\times 10^{-3}
Sphere6.7\times 10^{-6}
Random signs6.7\times 10^{-6}

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

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 “\”.1Even better, we could us a Cholesky decomposition since the matrix A^\top A is positive definite. (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.2One can represent rational numbers on a computer as fractions of integers and operations can be done exactly. However, this is prone to gross inefficiencies as the number of digits in the rational numbers can grow to be very large, making the storage and time to solve linear algebra problems with rationals dramatically more expensive. For these reasons, the vast majority of numerical computations use floating point numbers which store only a finite number of digits for any given real number. In this model, except for extremely rare circumstances, rounding errors during arithmetic operations are a fact of life. 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.3This is a simplification in multiple ways. First, computers store numbers in binary and thus, rather than storing 16 decimal digits after the decimal point, they store 52 binary digits. This amounts to roughly 16 decimal digits. Secondly, there are different formats for storing real numbers as floating point on a computer with different amounts of stored digits. The widely used IEEE double precision format has about 16 decimal digits of accuracy; the IEEE single precision format has roughly 8. 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).4In fact, there are multiple definitions of the condition number depending on the norm which one uses the measure the sizes of vectors. Since we use the 2-norm, the appropriate 2-norm condition number \kappa(M) is the ratio \kappa(M) = \sigma_{\rm max}(M)/\sigma_{\rm min}(M) of the largest and smallest singular values of 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.5In MATLAB, the least squares problem can be solved with QR factorization by calling “A\b”. Without going into details, the upshot is that the least squares solution by a well-designed6One way of computing the QR factorization is by Gram–Schmidt orthogonalization, but the accuracy properties of this are poor too. A gold-standard way of computing the QR factorization by means of Householder reflectors, which has excellent accuracy properties. 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 behavior7More precisely, the rule of thumb is like \|\hat{x} - x\|/\|x\| \lessapprox \kappa(A)\times 10^{-16} \times(1+ \kappa(A)\| b - Ax \|/(\|A\|\|b\|)). So even if we solve the least squares problem with QR factorization, we still get a squared condition number in our error bound, but this condition number squared is multiplied by the residual \|b-Ax\|, which is small if the least squares fit is good. The least squares solution is usually only interesting when the residual is small, thus justifying dropping it in the rule of thumb.

(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,8E.g., it can work for sparse matrices while the normal equations often do not, it has superior accuracy to Gaussian elimination with partial pivoting even for solving linear systems, the “Q” matrix in the QR factorization can be represented implicitly as a product of easy-to-compute-with Householder reflectors which is much more efficient whenm\gg n, etc. 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

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.1Everything carries over essentially unchanged for complex-valued matrices B with U and V being unitary matrices and every (\cdot)^\top being replaced by (\cdot)^* for (\cdot)^* the Hermitian transpose. 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:

  1. 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.2Relatedly, the two-norm condition number \kappa(B) := \sigma_{\rm max}(B) / \sigma_{\rm min}(B) of B^\top B is twice that of B. The dynamic range of our numerical computations has been cut in half!
  2. 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.
  3. 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.
  4. 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.3For instance, one cannot prove Weyl’s perturbation theorem for singular values by considering B^\top B and applying Weyl’s perturbation theorem for symmetric eigenvalues.

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 dilation4As Stewart and Sun detail in Section 4 of Chapter 1 of their monograph Matrix Perturbation Theory, the connections between the Hermitian dilation and the SVD go back to the discovery of the SVD itself, as it is used in Jordan’s construction of the SVD in 1874. (The SVD was also independently discovered by Beltrami the year previous.) Stewart and Sun refer to this matrix as the Jordan-Wiedlant matrix associated with B, as they attribute the widespread use of the matrix today to the work of Wiedlant. We shall stick to the term Hermitian dilation to refer to this matrix. 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.5This follows by noting \operatorname{rank}(\mathcal{H}(B)) = 2\operatorname{rank}(B) and thus \pm \sigma_i for i = 1,2,\ldots,\operatorname{rank}(B) account for all the nonzero eigenvalues of \mathcal{H}(B). 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.6The linearity of the Hermitian dilation gives direct extensions of most results about the symmetric eigenvalues to singular values; see Exercise 22.

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.7This is relevant if, say, we want to find the small singular values by inverse iteration. Positive definite linear systems are easier to solve by either direct or iterative methods. 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.