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.
 A: 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.
A: NB The stated result does not depend on any assumption of normality or even independence of the coordinates of $\newcommand{\x}{\mathbf{x}}\newcommand{\e}{\mathbb{E}}\newcommand{\tr}{\mathbf{tr}}\newcommand{\A}{\mathbf{A}}\x$. It does not depend on $\A$ being positive definite either. Indeed, suppose only that the coordinates of $\x$ have zero mean, variance of one and are uncorrelated (but not necessarily independent); that is, $\e \x_i = 0$, $\e \x_i^2 = 1$, and $\e \x_i \x_j = 0$ for all $i \neq j$.
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.
A: 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
