# Randomized trace technique

I've met the following randomized trace technique in M. Seeger, “Low rank updates for the Cholesky decomposition,” University of California at Berkeley, Tech. Rep, 2007.

$$\operatorname{tr}(\mathbf{A}) = {E[\mathbf{x}^T \mathbf{A} \mathbf{x}]}$$

where $\mathbf{x} \sim N(\mathbf{0},\mathbf{I})$.

As a person without deep mathematics background, I wonder how this equality can be achieved. Moreover, how can we interpret $\mathbf{x}^T \mathbf{A} \mathbf{x}$, for example geometrically? Where should I look in order to understand the meaning of taking the inner product of a vector and its range value? Why is the mean equal to the sum of the eigenvalues? Besides theoretical property, what is its practical importance?

I've written a MATLAB code snippet to see whether it works

#% tr(A) == E[x'Ax], x ~ N(0,I)

N = 100000;
n = 3;
x = randn([n N]); % samples
A = magic(n); % any n by n matrix A

y = zeros(1, N);
for i = 1:N
y(i) = x(:,i)' * A * x(:,i);
end
mean(y)
trace(A)


The trace is 15 where the approximation is 14.9696.

## 3 Answers


Bare-hands approach

Let $\A = (a_{ij})$ be an arbitrary $n \times n$ matrix. By definition $\tr(\A) = \sum_{i=1}^n a_{ii}$. Then, $$\tr(\A) = \sum_{i=1}^n a_{ii} = \sum_{i=1}^n a_{ii} \e \x_i^2 = \sum_{i=1}^n a_{ii} \e \x_i^2 + \sum_{i\neq j} a_{ij} \e \x_i \x_j ,$$ and so we are done.

In case that's not quite obvious, note that the right-hand side, by linearity of expectation, is $$\sum_{i=1}^n a_{ii} \e \x_i^2 + \sum_{i\neq j} a_{ij} \e \x_i \x_j = \e\Big(\sum_{i=1}^n \sum_{j=1}^n a_{ij} \x_i \x_j \Big) = \e(\x^T \A \x)$$

Proof via trace properties

There is another way to write this that is suggestive, but relies, conceptually on slightly more advanced tools. We need that both expectation and the trace operator are linear and that, for any two matrices $\A$ and $\newcommand{\B}{\mathbf{B}}\B$ of appropriate dimensions, $\tr(\A\B) = \tr(\B\A)$. Then, since $\x^T \A \x = \tr(\x^T \A \x)$, we have $$\e(\x^T \A \x) = \e( \tr(\x^T \A \x) ) = \e( \tr(\A \x \x^T) ) = \tr( \e( \A \x \x^T ) ) = \tr( \A \e \x \x^T ),$$ and so, $$\e(\x^T \A \x) = \tr(\A \mathbf{I}) = \tr(\A) .$$

Quadratic forms, inner products and ellipsoids

If $\A$ is positive definite, then an inner product on $\mathbf{R}^n$ can be defined via $\langle \x, \mathbf{y} \rangle_{\A} = \x^T \A \mathbf{y}$ and $\mathcal{E}_{\A} = \{\x: \x^T \A \x = 1\}$ defines an ellipsoid in $\mathbf{R}^n$ centered at the origin.

• $\mathbf{x}_i$ is the $i$th coordinate of the vector $\mathbf{x}$. The others are simply typos. Sorry about that. I was trying to follow your notation as closely as possible. I would normally use $\mathbf{X} = (X_i)$ with $X_i$ as the coordinates of the random variable $\mathbf{X}$. But, I didn't want to (potentially) confuse. – cardinal Jun 16 '11 at 19:45
• Actually, it is consistent within the answer. I just wanted to make sure that the subscripted variables are the elements of the vector. Now, it is clear. – petrichor Jun 16 '11 at 21:07
• Well, it is consistent (now) because I edited it! :) Thanks for pointing out the typos. I'll try to add a little more about the geometry at some point over the next couple of days. – cardinal Jun 16 '11 at 22:07

If $A$ is symmetric positive definite, then $A = U^tDU$ with $U$ orthonormal, and $D$ diagonal with eigenvalues on the diagonal. Since $x$ has identity covariance matrix, and $U$ is orthonormal, $Ux$ is also has an identity covariance matrix. Hence writing $y = Ux$, we have $E[x^TAx] = E[y^tDy]$. Since the expectation operator is linear, this is just $\sum_{i=0}^n \lambda_i E[y_i^2]$. Each $y_i$ is chi-square with 1 degree of freedom, so has expected value 1. Hence the expectation is the sum of eigenvalues.

Geometrically, symmetric positive definite matrices $A$ are in 1-1 correspondence with ellipsoids -- given by the equation $x^TAx = 1$. The lengths of the ellipsoid's axes are given by $1/\sqrt\lambda_i$ where $\lambda_i$ are the eigenvalues.

When $A = C^{-1}$ where $C$ is the covariance matrix, this is the square of the Mahalanobis distance.

Let me address the "what is its practical importance" part of the question. There are many situations in which we have the ability to compute matrix vector products $Ax$ efficiently even if we don't have a stored copy of the matrix $A$ or don't have enough storage to save a copy of $A$. For example, $A$ might be of size 100,000 by 100,000 and fully dense- it would require 80 gigabytes of RAM to store such a matrix in double preciison floating point format.

Randomized algorithms like this can be used to estimate the trace of $A$ or (using a related algorithm) individual diagonal entries of $A$.

Some applications of this technique to large scale geophysical inversion problems are discussed in

J. K. MacCarthy, B. Borchers, and R. C. Aster. Efficient stochastic estimation of the model resolution matrix diagonal a nd generalized cross validation for large geophysical inverse problems. Journal of Geophysical Research, 116, B10304, 2011. Link to the paper

• +1 I met with randomized algorithms this semester and fascinated with them. Let me add another nice article. Nathan Halko, Per-Gunnar Martinsson, Joel A. Tropp, "Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions", 2010, arxiv.org/abs/0909.4061 – petrichor Jan 13 '12 at 17:26