What Is the Inverse Covariance Matrix?

In statistics and optimization, we often multiply a vector $\mathbf{x}$ by the inverse of its covariance matrix $\boldsymbol{\Sigma}$:

$\boldsymbol{\Sigma}^{-1} \mathbf{x}. \tag{1}$

For example, you see this in the probability density function of the multivariate normal distribution or in the analytic solution to mean–variance optimization. What does this operation do? I’ve written about this before, when I discussed the covariance matrix $\boldsymbol{\Sigma}$ as high-dimensional variance, but that was cursory at best. The goal of this post is to explore this operation in more detail.

If $\boldsymbol{\Sigma}$ is a $1 \times 1$ matrix, then clearly $\mathbf{x}$ is a $1$-vector, and this just normalizes the scalar value:

$\frac{x}{\sigma^2}. \tag{2}$

And if $\boldsymbol{\Sigma}$ is a diagonal $P \times P$ matrix, then clearly $\mathbf{x}$ is a $P$-vector, and this just normalizes each component of $\mathbf{x}$:

$\boldsymbol{\Sigma}^{-1} \mathbf{x} = \begin{bmatrix} x_1 / \sigma_1^2 \\ x_2 / \sigma_2^2 \\ \vdots \\ x_P / \sigma_P^2 \end{bmatrix}. \tag{3}$

The interesting case is when $\boldsymbol{\Sigma}$ is a non-diagonal $P \times P$ matrix. In that case, we can visualize the operation via the eigendecomposition, sometimes called the spectral decomposition. Recall that the eigendecomposition of a square matrix is

$\boldsymbol{\Sigma} = \mathbf{Q} \boldsymbol{\Lambda} \mathbf{Q}^{\top}, \tag{4}$

where $\mathbf{Q}$ is a square matrix whose $n$-th column is the $n$-th eigenvector of $\boldsymbol{\Sigma}$ and where $\boldsymbol{\Lambda}$ is a diagonal matrix of associated eigenvalues $\{\lambda_1, \dots, \lambda_P \}$. Since $\mathbf{Q}$ is an orthonormal basis and $\boldsymbol{\Lambda}$ is a diagonal matrix, we can think of it as: rotate our data into a vector space in which the transformed data is orthogonal, scale in this space, and rotate back to the original vector space. If this kind of geometric interpretation makes no sense, please see my post on the singular value decomposition (SVD).

Now let’s see this in equations. Raising this matrix to a power can be simplified as:

$\boldsymbol{\Sigma}^n = \underbrace{\left( \mathbf{Q} \boldsymbol{\Lambda} \mathbf{Q}^{\top} \right) \left( \mathbf{Q} \boldsymbol{\Lambda} \mathbf{Q}^{\top} \right) \dots \left( \mathbf{Q} \boldsymbol{\Lambda} \mathbf{Q}^{\top} \right)}_{\text{$k$ times}} = \mathbf{Q} \underbrace{\boldsymbol{\Lambda} \boldsymbol{\Lambda} \dots \boldsymbol{\Lambda}}_{\text{$k$ times}} \mathbf{Q}^{\top} = \mathbf{Q} \boldsymbol{\Lambda}^k \mathbf{Q}^{\top}. \tag{5}$

This works because $\mathbf{Q}$ is orthonormal and thus

$\mathbf{Q} \mathbf{Q}^{\top} = \mathbf{Q}^{\top} \mathbf{Q} = \mathbf{I}, \tag{6}$

where $\mathbf{I}$ is the identity matrix. And this means that $\boldsymbol{\Lambda}^k$ is just a diagonal matrix where each element is raised to the power $k$.

It follows that the inverse matrix $\boldsymbol{\Sigma}^{-1}$ is the original spectral decomposition but with $\boldsymbol{\Lambda}^{-1}$ instead of $\boldsymbol{\Lambda}$. And we can interpret multiplying a vector $\mathbf{x}$ by this inverse as follows:

$\begin{aligned} \mathbf{z} &:= \mathbf{Q}^{\top} \mathbf{x} && \text{rotate into orthogonal vector space} \\ \tilde{\mathbf{z}} &:= \boldsymbol{\Lambda}^{-1} \mathbf{z} && \text{scale by the eigenvalues} \\ \tilde{\mathbf{x}} &:= \mathbf{Q} \tilde{\mathbf{z}} && \text{rotate back into original vector space} \end{aligned} \tag{7}$

What will this operation actually do? We can visualize it in Figure $1$. In my mind, the geometric fact that jumps out the most is that when we divide by eigenvalues, we convert our ellipsoid into another ellipsoid, stretching along the dimension that originally was the smaller of the two. So for example, if $\boldsymbol{\Sigma}$ was the empirical covariance matrix of $\mathbf{x}$, we would normalize a la Equation $2$ but in $P$ dimensions. Variables of $\mathbf{x}$ with high variance get downweighted.

Figure 1. Step-by-step multiplication of random vectors $\mathbf{x}$ with the components of the spectral decomposition of their inverse covariance matrix $\boldsymbol{\Sigma}^{-1}$. (Left) Random samples $\mathbf{x}$ from a multivariate normal distribution. The red and black arrows are the axes of an ellipsoid, defined by the square root of the eigenvalues $\{\lambda_1, \lambda_2\}$ of $\boldsymbol{\Sigma}^{-1}$. (Middle left) The random vectors and ellipsoid after multiplication by $\mathbf{Q}^{\top}$, a rotation matrix. (Middle right) The random vectors and ellipsoid after further rescaling by $\{\lambda_1^{-1}, \lambda_2^{-1}\}$. (Right) The random vectors and ellipsoid, rotated back into the original vector space by the rotation matrix $\mathbf{Q}$.

This might be surprising at first. Perhaps you expected this operation to “whiten” the data in the sense of being the multivariate analog to dividing by the standard deviation. But it makes sense that it is not, because as we saw above, if $\boldsymbol{\Sigma}$ is a digaonal matrix, then $\boldsymbol{\Sigma}^{-1}$ is like dividing by variances, not standard deviations.

This suggests that to whiten the data, we simply need to multiply $\mathbf{x}$ by $\boldsymbol{\Sigma}$ raised to the power $k = -1/2$. See Figure $2$. We can see that after the initial rotation, dividing by the square root of the eigenvalues converts our ellipsoid into a sphere. The final rotation simply rotates the whitened data back into “data space”.

Figure 2. Figure $1$ but with $\boldsymbol{\Sigma}^{-1/2}$. Thus, the ellipsoid is transformed into a circle, since each axis is divided by its the square root of its eigenvalue, i.e. by its singular value.

This should make sense if you understand the SVD and that the singular values are the square roots of the eigenvalues. To see this, think about the SVD of our data matrix:

$\mathbf{X} = \mathbf{U S V}^{\top}. \tag{8}$

Here, $\mathbf{X}$ has shape $N \times P$, where $N$ is the number of observations and $P$ is again the dimensionality of our data. Intuitively, the singular values of a matrix $\mathbf{X}$ are in the units of $\mathbf{X}$, since $\mathbf{U}$ and $\mathbf{V}$ are just rotations. And what are the eigenvalues? They’re the square of the singular values, since:

$\boldsymbol{\Sigma} := \mathbf{X}^{\top} \mathbf{X} = (\mathbf{U S V}^{\top})^{\top} (\mathbf{U S V}^{\top}) = \mathbf{V S}^2 \mathbf{V}^{\top}. \tag{9}$

Geometrically, we stretch by $\mathbf{X}^{\top}$ and then by $\mathbf{X}$—so we’re double stretching. See Suraj Rampure’s notes for a bit more detail. But the high-level point is that $\boldsymbol{\Sigma}^{-1} \mathbf{x}$ is the multivariate analog of normalizing by the variance, while $\boldsymbol{\Sigma}^{-1/2} \mathbf{x}$ is multivariate analog of normalizing by the standard deviation.

Finally, this means we can generalize to any power of $k$ of $\boldsymbol{\Sigma}$. When $k$ is large, our data is stretched along it’s current direction. When $k$ is negative, the data is stretched in the opposite direciton. And when $k \approx 0$, then the data is unchanged (Figure $3$).

Figure 3. Figures $1$ and $2$ but with other powers $k$.

This view of the inverse covariance matrix is helpful in interpreting other operations. For example, consider this derivation, with hopefully obvious notation, where now $\mathbf{u}$ and $\mathbf{v}$ denote arbitrary vectors unrelated to the SVD:

$\begin{aligned} \mathbf{u}^{\top} \mathbf{D}^{-1} \mathbf{v} &= \mathbf{u}^{\top} \mathbf{Q} \boldsymbol{\Lambda}^{-1} \mathbf{Q}^{\top} \mathbf{v} \\ &= \mathbf{u}_Q^{\top} \boldsymbol{\Lambda}^{-1} \mathbf{v}_Q \\ &= \mathbf{u}_Q^{\top} \boldsymbol{\Lambda}^{-1/2} \boldsymbol{\Lambda}^{-1/2} \mathbf{v}_Q \\ &= \mathbf{u}_{Q\Lambda}^{\top} \mathbf{v}_{Q\Lambda}. \end{aligned} \tag{10}$

This shows that $\mathbf{u}^{\top} \mathbf{D}^{-1} \mathbf{v}$ is just a standard Euclidean inner product $\mathbf{u}_{Q\Lambda}^{\top} \mathbf{v}_{Q\Lambda}$ after we rotate and whiten both vectors. The inverse matrix $\mathbf{D}^{-1}$ tells us how to transform each vector into a space where the standard dot product is the correct notion of distance.

For a more complicated example, consider a common modeling assumption, which is that a variable $\mathbf{x}_n \in \mathbb{R}^{P}$ is a linear function of a low-dimensional factor $\mathbf{f}_n \in \mathbb{R}^K$ where $K \ll P$. We can write this in matrix notation as

$\mathbf{X} = \mathbf{F} \mathbf{B} + \mathbf{E}, \tag{11}$

where $\mathbf{X} \in \mathbb{R}^{N \times P}$, $\mathbf{F} \in \mathbb{R}^{N \times K}$, $\mathbf{B} \in \mathbb{R}^{K \times P}$, and $\mathbf{E} \in \mathbb{R}^{N \times P}$. Under some common assumptions, we can show that the covariance matrix of $\mathbf{x}$ is

$\boldsymbol{\Sigma}_x := \text{cov}[\mathbf{x}, \mathbf{x}] = \mathbf{B}^{\top} \boldsymbol{\Sigma}_f \mathbf{B} + \boldsymbol{\Sigma}_{e}. \tag{12}$

For more details, we see factor analysis. In this scenario, the Woodbury identity tells us how to invert our low-rank covariance matrix $\boldsymbol{\Sigma}_x$:

$\boldsymbol{\Sigma}_x^{-1} = \boldsymbol{\Sigma}_{e}^{-1} - \boldsymbol{\Sigma}_{e}^{-1} \mathbf{B}^{\top} \left(\boldsymbol{\Sigma}_f^{-1} + \mathbf{B} \boldsymbol{\Sigma}_{e}^{-1} \mathbf{B}^{\top} \right)^{-1} \mathbf{B} \boldsymbol{\Sigma}_{e}^{-1}. \tag{13}$

Now here’s the question: what does it mean to normalize an observation $\mathbf{x}_n$ by this quantity? To understand, let’s rewrite this as

$\begin{aligned} \boldsymbol{\Sigma}_x^{-1} \mathbf{x}_n &= \boldsymbol{\Sigma}_{e}^{-1} \mathbf{x}_n - \boldsymbol{\Sigma}_{e}^{-1} \mathbf{B}^{\top} \left(\boldsymbol{\Sigma}_f^{-1} + \mathbf{B} \boldsymbol{\Sigma}_{e}^{-1} \mathbf{B}^{\top} \right)^{-1} \mathbf{B} \boldsymbol{\Sigma}_{e}^{-1} \mathbf{x}_n \\ &= \boldsymbol{\Sigma}_{e}^{-1} \left( \mathbf{x}_n - \mathbf{B}^{\top}_n \hat{\mathbf{f}}_n \right) \\ &= \boldsymbol{\Sigma}_{e}^{-1} \tilde{\mathbf{x}}_n, \end{aligned} \tag{14}$

where $\hat{\mathbf{f}}_n$ is the generalized least squares estimator of the factor $\mathbf{f}_n$ using weights $\boldsymbol{\Sigma}_e$ and Ridge penality $\boldsymbol{\Sigma}_f$,

$\hat{\mathbf{f}}_n := \left(\boldsymbol{\Sigma}_f^{-1} + \mathbf{B} \boldsymbol{\Sigma}_{e}^{-1} \mathbf{B}^{\top} \right)^{-1} \mathbf{B} \boldsymbol{\Sigma}_{e}^{-1} \mathbf{x}_n, \tag{15}$

and where $\tilde{\mathbf{x}}_n$ is thus our observation residualized by its estimated factor. This is why, if you squint your eyes, parts of the Woodbury identity look like the normal equation for ordinary least squares. It’s a projection. In words, we can say: we estimate the factor component of $\mathbf{x}_n$, subtract that factor component out, and then whiten the residual against the idiosyncratic variances.

We can concretize this a bit with an example. In mean–variance optimization, one’s optimal portfolio position $\mathbf{w}^{\star}$ without costs or constraints is essentially one’s risk-adjusted alpha:

$\mathbf{w}^{\star} \propto \boldsymbol{\Sigma}^{-1} \mathbf{a}. \tag{16}$

In finance, people commonly assume a factor structure for $\boldsymbol{\Sigma}$. See risk factor modeling for details. For example, we might want to explicitly model that two companies in the same industry “load” onto a shared industry factor and model their correlation through that factor rather than by directly estimating $\boldsymbol{\Sigma}$. In this model, our optimal position is

$\mathbf{w}^{\star} \propto \boldsymbol{\Sigma}_e^{-1} \left( \mathbf{a} - \mathbf{B}^{\top} \hat{\mathbf{f}} \right). \tag{17}$

Given everything above, we can interpret this optimal position as our factor-residualized alpha, adjusted by its idiosyncratic risk. In words, we only care about the alpha that isn’t correlated to other alpha through a common factor, and then we rescale this residualized alpha by its idiosyncractic variance, such that, all things being equal, we take larger positions in less volatile names.