2
$\begingroup$

Given a symmetric, positive definite, covariance matrix $S$, how can I calculate a matrix $C$ such that

$CSC^T$ is the Identity matrix $I$?

$\endgroup$
  • 2
    $\begingroup$ Any invertible matrix for which $AA^T = S$ with $C=A^{-1}$ will work, so for example, Cholesky decomposition is one way - it's a simple algorithm for finding a lower triangular $A$ (or equivalently, in many implementations, an upper triangular $A^T$). That is, it finds $S = LL^T$, whence $L^{-1}S(L^{-1})^T = I$. But a number of other choices are available. $\endgroup$ – Glen_b Aug 16 '13 at 0:57
  • 2
    $\begingroup$ To appreciate @Glen_b's last point, observe that there are many orthogonal matrices $\mathbb{O}$: by definition, these preserve the identity under conjugation, $\mathbb{OIO}'=\mathbb{I}$. If $\mathbb{CSC}'=\mathbb{I}$, then clearly $\mathbb{I=OIO'=OCSC'O'=(OC)S(OC)'}$, whence $\mathbb{OC}$ works just as well. Conversely, $\mathbb{CSC'=I=BSB'}$ implies $\mathbb{B}$ is invertible and $\mathbb{B^{-1}C}$ is orthogonal. (The set of $n$ by $n$ orthogonal matrices is naturally a submanifold of $\mathbb{R}^{n^2}$ of dimension $n(n-1)/2$: that's the sense in which there are quite a lot of them.) $\endgroup$ – whuber Aug 16 '13 at 1:57
  • $\begingroup$ I expected this question would already have an answer for me to point to but I'm not seeing anything exactly on topic here. If nobody points a duplicate out or posts an answer, I may expand on my comment a little and turn it into one so that the question is answered at least once (and that expansion would include an explanation about how it's usually the case that if $S = X'X$ we do something else). $\endgroup$ – Glen_b Aug 16 '13 at 2:03
3
$\begingroup$

Since I haven't found a duplicate and nobody else has pointed one out, I have expanded my comment into an answer and added an example to show you how it works.

Any invertible matrix $A$ for which $AA^T=S$ will work by taking $C=A^{−1}$ .

So for example, Cholesky decomposition is one common way to do it - it's a simple algorithm for finding a lower triangular $A$ (or equivalently, in many implementations, an upper triangular $A^T$). That is, it finds $S=LL^T$, whence $L^{−1}S(L^{−1})^T=I$.

Example:

Consider the matrix $$S=\left(\begin{matrix} 4 & 1 \\ 1 & 2.5 \end{matrix}\right) $$

We write $S = LL^T$, which implies that

$$s_{11} = l_{11}^2$$

That is, $l_{11}=2$ (we'll take positive square roots each time to make our diagonals on $L$ positive).

It also implies that $s_{21} = l_{21}l_{11}$, so $l_{21} = s_{21}/l_{11}$, i.e. $l_{21} = 1/2$. Then we see that $s_{22} = l_{21}^2+l_{22}^2$, so $l_{22}^2 = s_{22} - l_{21}^2 = 2.5 - 0.5^2 = 2.25$, in which case $l_{22} = 1.5$, so

$$L=\left(\begin{matrix} 2 & 0 \\ 0.5 & 1.5 \end{matrix}\right) $$

You can confirm that

$$C = L^{-1}=\left(\begin{matrix} \frac{1}{2} & 0 \\ -\frac{1}{6} & \frac{2}{3} \end{matrix}\right) $$

and that $CSC^T = I$.

In practice we don't invert $L$, however, but solve a linear system, it's more stable.

Another possibility is to use Singular Value Decomposition to find a matrix to do the diagonalizing, and then scale by the singular values.

But a number of other choices are available. (As whuber points out, a very large number.)

In practice, frequently we know $S$ is of the form $X^TX$ for some $X$. In that case, it's usually better to work on $X$ than $S$, via the QR decomposition.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.