# Find a design matrix giving a prescribed covariance matrix

Given a symmetrical positive-definite matrix $L$, how do I find a matrix $X$ such that $L = X^T X$?

• This isn't really on-topic, but the Cholesky decomposition is typically the most convenient. Commented May 29, 2015 at 8:33
• It is not possible to recover design matrix $X$ from the product $X^TX$, without any additional information, since for any orthogonal matrix $O$, the product $(OX)^T(OX)=X^TX$ does not change. Commented May 29, 2015 at 9:59
• @Dougal What is the reason you think this question is not on-topic here?
– whuber
Commented May 29, 2015 at 13:37
• @whuber I somehow missed the title of the question in the review page, and thought it was purely a linear algebra question. With the edited title it's clearer that it's appropriate here. Commented May 29, 2015 at 17:03

Essentially, you can reverse the principal components process.

Let's start by making some uncorrelated, centered and normalized data, you can create the columns in your matrix by drawing from any distribution as long as the draws are independent. Here I draw from a centered and standardized uniform distribution $X \sim U(-.5, .5) / \sqrt{12}$, but only because it is convenient and illustrates the process:

> X <- matrix(runif(120000, -.5, .5)/ (1/sqrt(12)),
nrow=40000)
> t(X) %*% X / 400000
[,1]          [,2]          [,3]
[1,]  0.0996386570 -0.0002036308 -0.0006817355
[2,] -0.0002036308  0.1000222332  0.0015659048
[3,] -0.0006817355  0.0015659048  0.1000916918


That is:

$$\frac{1}{N} X^t X = I$$

Now let's write down a symmetric, positive definite matrix that we want to be our final correlation matrix:

S <- matrix(c(  1, -.5,   0,
-.5,   1, -.5,
0, -.5,   1), nrow=3, byrow=TRUE
)


Use the sqrtm function to compute the matrix square root:

> rtS <- sqrtm(S)
> t(rtS) %*% rtS
[,1] [,2]          [,3]
[1,]  1.000000e+00 -0.5  5.551115e-17
[2,] -5.000000e-01  1.0 -5.000000e-01
[3,]  5.551115e-17 -0.5  1.000000e+00


Then form the matrix product $X\sqrt{S}$. This has the property that:

$$(X\sqrt{S})^t X\sqrt{S} = \sqrt{S}^t X^t X \sqrt{S} = N \sqrt{S}^t \sqrt{S} = NS$$

so

$$\frac{1}{N} (X\sqrt{S})^t (X\sqrt{S}) = S$$

so $X \sqrt{S}$ is the data matrix you want

> t(X %*% sqrtm(S)) %*% (X %*% sqrtm(S)) / 40000
[,1]       [,2]         [,3]
[1,]  0.998555833 -0.4987640 -0.009608847
[2,] -0.498764047  0.9921826 -0.483641573
[3,] -0.009608847 -0.4836416  0.993164671