Why does inversion of a covariance matrix yield partial correlations between random variables? I heard that partial correlations between random variables can be found by inverting the covariance matrix and taking appropriate cells from such resulting precision matrix (this fact is mentioned in http://en.wikipedia.org/wiki/Partial_correlation, but without a proof).
Why is this the case?
 A: Note that the sign of the answer actually depends on how you define partial correlation. There is a difference between regressing $X_i$ and $X_j$ on the other $n - 1$ variables separately vs. regressing $X_i$ and $X_j$ on the other $n - 2$ variables together. Under the second definition, let the correlation between residuals $\epsilon_i$ and $\epsilon_j$ be $\rho$. Then the partial correlation of the two (regressing $\epsilon_i$ on $\epsilon_j$ and vice versa) is $-\rho$. 
This explains the confusion in the comments above, as well as on Wikipedia. The second definition is used universally from what I can tell, so there should be a negative sign.
I originally posted an edit to the other answer, but made a mistake - sorry about that!
A: When a multivariate random variable $(X_1,X_2,\ldots,X_n)$ has a nondegenerate covariance matrix $\mathbb{C} = (\gamma_{ij}) = (\text{Cov}(X_i,X_j))$, the set of all real linear combinations of the $X_i$ forms an $n$-dimensional real vector space with basis $E=(X_1,X_2,\ldots, X_n)$ and a non-degenerate inner product given by
$$\langle X_i,X_j \rangle = \gamma_{ij}\ .$$
Its dual basis with respect to this inner product, $E^{*} = (X_1^{*},X_2^{*}, \ldots, X_n^{*})$, is uniquely defined by the relationships
$$\langle X_i^{*}, X_j \rangle = \delta_{ij}\ ,$$
the Kronecker delta (equal to $1$ when $i=j$ and $0$ otherwise).
The dual basis is of interest here because the partial correlation of $X_i$ and $X_j$ is obtained as the correlation between the part of $X_i$ that is left after projecting it into the space spanned by all the other vectors (let's simply call it its "residual", $X_{i\circ}$) and the comparable part of $X_j$, its residual $X_{j\circ}$. Yet $X_i^{*}$ is a vector that is orthogonal to all vectors besides $X_i$ and has positive inner product with $X_i$ whence $X_{i\circ}$ must be some non-negative multiple of $X_i^{*}$, and likewise for $X_j$.  Let us therefore write
$$X_{i\circ} = \lambda_i X_i^{*},\ X_{j\circ} = \lambda_j X_j^{*}$$
for positive real numbers $\lambda_i$ and $\lambda_j$.
The partial correlation is the normalized dot product of the residuals, which is unchanged by rescaling:
$$\rho_{ij\circ} = \frac{\langle X_{i\circ}, X_{j\circ} \rangle}{\sqrt{\langle X_{i\circ}, X_{i\circ} \rangle\langle X_{j\circ}, X_{j\circ} \rangle}} =  \frac{\lambda_i\lambda_j\langle X_{i}^{*}, X_{j}^{*} \rangle}{\sqrt{\lambda_i^2\langle X_{i}^{*}, X_{i}^{*} \rangle\lambda_j^2\langle X_{j}^{*}, X_{j}^{*} \rangle}} = \frac{\langle X_{i}^{*}, X_{j}^{*} \rangle}{\sqrt{\langle X_{i}^{*}, X_{i}^{*} \rangle\langle X_{j}^{*}, X_{j}^{*} \rangle}}\ .$$
(In either case the partial correlation will be zero whenever the residuals are orthogonal, whether or not they are nonzero.)
We need to find the inner products of dual basis elements.  To this end, expand the dual basis elements in terms of the original basis $E$:
$$X_i^{*} = \sum_{j=1}^n \beta_{ij} X_j\ .$$
Then by definition
$$\delta_{ik} = \langle X_i^{*}, X_k \rangle = \sum_{j=1}^n \beta_{ij}\langle X_j, X_k \rangle = \sum_{j=1}^n \beta_{ij}\gamma_{jk}\ .$$
In matrix notation with $\mathbb{I} = (\delta_{ij})$ the identity matrix and $\mathbb{B} = (\beta_{ij})$ the change-of-basis matrix, this states
$$\mathbb{I} = \mathbb{BC}\ .$$
That is, $\mathbb{B} = \mathbb{C}^{-1}$, which is exactly what the Wikipedia article is asserting.  The previous formula for the partial correlation gives
$$\rho_{ij\cdot} =  \frac{\beta_{ij}}{\sqrt{\beta_{ii} \beta_{jj}}} = \frac{\mathbb{C}^{-1}_{ij}}{\sqrt{\mathbb{C}^{-1}_{ii} \mathbb{C}^{-1}_{jj}}}\ .$$
A: For another perspective, this will examine the left inverse of a finite data matrix $A$.  We can consider the data to be a sample rather than a theoretical distribution.  While any distribution -- even continuous -- will have a covariance matrix, you can't generally talk about a data matrix unless you get into infinite vectors and/or special inner products.   
So we have a finite sample in an n-by-m data matrix $A$.  Let each column be one random variable. Then it's $n$ samples and $m$ random variables. Let $A$'s columns (the random variables) be linearly independent (this is independence in the linear algebra sense, not as in independent random variables).
Let $A$ be mean-centered already. Then,
$$
C = \frac{1}{n}A^TA
$$
is our covariance matrix.  It's invertible since $A$'s columns are linearly independent.
And we'll use later that $C^{-1} = n(A^TA)^{-1}$
The left inverse of $A$ is
$B = (A^TA)^{-1}A^T$.
And we have
$BA = I_{m-by-m}$.
What do we know about $B$?


*

*It's m-by-n.  There's a row of $B$ corresponding to each column of $A$.

*Because $BA = I$, we know the inner product of the $i$th row of $B$ with the $i$th column in $A$ equals 1 (diagonal of $I$).

*An inner product of the $i$th row of $B$ with a $j$th ($i \neq j$) column of $A$ is 0 (off-diagonal of $I$).

*The right-most term in the expression for $B$ is $A^T$.  Therefore $B$'s rows are in the rowspace of $A^T$, the column space of $A$.

*by (4) and the fact that $A$'s columns are mean-centered, $B$'s rows must also be mean-centered.


Let $x_i$ be the $i$th column of $A$.
The only vectors that have a non-zero inner product with the $x_i$, zero inner product with all other $x_j$, and are linear combinations of the columns of $A$, are vectors parallel to the residual of $x_i$ after projecting it into the space spanned by all the other $x_j$.
Call these residuals $r_{i}$. And call the projection (the linear regression result) $p_i$. So the $i$th row of $B$ must be parallel to $r_i$ (6).
Now we know its direction, but what about magnitude?  Let $b_i$ be the $i$th row of $B$. 
$$
\begin{align}
1 & = b_i \cdot x_i &&\text{by (2)} \\
    & = b_i \cdot (p_i + r_i)  &&\text{$x_i$ is the sum of its projection and residual}\\
    & = (b_i \cdot p_i) + (b_i \cdot r_i) &&\text{linearity of dot product} \\
    & = 0 + (b_i \cdot r_i) &&\text{by (3), and that $p_i$ is a linear combination of the $x_j$s ($j \neq i$)} \\
    & = (c_i r_i) \cdot r_i &&\text{for some constant $c_i$, by (6)} \\
\end{align}
$$
Therefore, $c_i = \dfrac{1}{r_i \cdot r_i} = \dfrac{1}{\|r_i\|^2}$, so $b_i = \dfrac{r_i}{\|r_i\|^2}$.
We now know what each row of $B$ looks like.  Notice
$BB^T = ((A^TA)^{-1}A^T)(A((A^TA)^{-1})^T) = (A^TA)^{-1} = \frac{1}{n}C^{-1}$
We can look at any $i,j$th element
$C^{-1}_{ij} = n(BB^T)_{ij} = n (b_i \cdot b_j) = n\dfrac{r_i \cdot r_j}{\|r_i\|^2\|r_j\|^2}$
The $(r_i \cdot r_j)$ part of that should tell you we're getting close to covariances and correlations of these residuals. Conveniently, the diagonal elements look like
$C^{-1}_{ii} =  n\dfrac{r_i \cdot r_i}{\|r_i\|^2\|r_i\|^2} = n\dfrac{1}{\|r_i\|^2}$.
This quantity is exactly 1 over the variance of the residual $r_i$, $\dfrac{\|r_i\|^2}{n}$ (the $n$ makes it a variance instead of a squared vector magnitude).
Then to get partial correlations you just need to combine the elements of $C^{-1}$ in the way others have shown.


*

*Gilbert Strang lecture on left inverses

*Gilbert Strang lecture on projection, residuals
