Let's say I have 3 variables, A, B, and C and I want to generate data where A and B are correlated at r=x, A and C are correlated at r=y, and B and C are correlated at r=z.

  1. Is there any algorithm that will tell me, given specified values for x, y, and z, if any set of variances for A, B and C will yield a positive definite covariance matrix?
  2. Given the existence of such an algorithm, is there another algorithm that, in cases where a positive definite covariance matrix is possible, will give me a set of variances that achieve a positive definite matrix?
  • 1
    $\begingroup$ I'm not sure I quite understand your question. Is there a reason you don't like or can't do an eigendecomposition? Without loss of generality you can assume all variances are one. $\endgroup$ – cardinal Feb 22 '12 at 15:53
  • $\begingroup$ No reason but for ignorance! Care to post an answer that elaborates on how to do an eigendecomposition (ideally in R)? $\endgroup$ – Mike Lawrence Feb 22 '12 at 16:17
  • $\begingroup$ Sure. Give me a bit and I'll be glad to do so. $\endgroup$ – cardinal Feb 22 '12 at 16:20
  • $\begingroup$ I think the essence of the answer can be captured in statistical language: correlations are covariances of standardized versions of the variables. Therefore you can arbitrarily rescale (and recenter) the variables without changing the correlation matrix. In conjunction with any algorithm to test for symmetry and positive (semi)definiteness, that fully answers question 1 and question 2. Eigendecompositions are useful (for constructing variables with a given correlation matrix) but are really not essential to understanding the situation. $\endgroup$ – whuber Feb 22 '12 at 17:57
  • $\begingroup$ @whuber: The eigendecomposition helps me understand the situation since it tells me exactly what a covariance matrix must look like both operationally and intuitively: It must be one that takes an arbitrary vector, rotates it into a new coordinate space, then fiddles with its coordinates in such a way that it stays strictly inside the orthant it started in, and then rotates the new version back to the original space. In the rotated space a vector will lie exactly on the boundary between orthants after the "fiddling" if and only if it started on the same boundary. $\endgroup$ – cardinal Feb 22 '12 at 19:23

To follow up on @cardinal's comment: your $x$, $y$, and $z$ define a $(3 \times 3)$ correlation matrix $R$. Since a correlation matrix also is a possible covariance matrix (of standardized variables), it has to be positive definite. This is the case if all eigenvalues are $> 0$. If $R$ is indeed positive definite, then all vectors $\boldsymbol{s}$ of variances (i.e., numbers $> 0$) will turn $\boldsymbol{R}$ into a positive definite covariance matrix $\boldsymbol{\Sigma} = \boldsymbol{D}_{s}^{1/2} \boldsymbol{R} \boldsymbol{D}_{s}^{1/2}$, where $\boldsymbol{D}_{s}^{1/2}$ is the square root of the diagonal matrix made from $\boldsymbol{s}$.

So just construct $R$ from $x, y, z$, and check if the eigenvalues are all $> 0$. If so, you're good, and you can transform any set of data to have a corresponding covariance matrix with arbitrary variances:

x <- 0.5
y <- 0.3                            # changing this to -0.6 makes it not pos.def.
z <- 0.4
R <- matrix(numeric(3*3), nrow=3)   # will be the correlation matrix
diag(R) <- 1                        # set diagonal to 1
R[upper.tri(R)] <- c(x, y, z)       # fill in x, y, z to upper right
R[lower.tri(R)] <- c(x, y, z)       # fill in x, y, z to lower left
eigen(R)$values                     # get eigenvalues to check if pos.def.


[1] 1.8055810 0.7124457 0.4819732

So our $\boldsymbol{R}$ here is positive definite. Now construct the corresponding covariance matrix from arbitrary variances.

vars  <- c(4, 16, 9)                # the variances
Sigma <- diag(sqrt(vars)) %*% R %*% diag(sqrt(vars))

Generate some data matrix $\boldsymbol{X}$ that we will transform to later have exactly that covariance matrix.

library(mvtnorm)                    # for rmvnorm()
N  <- 100                           # number of simulated observations
mu <- c(1, 2, 3)                    # some arbitrary centroid
X  <- round(rmvnorm(n=N, mean=mu, sigma=Sigma))

To do that, we first orthonormalize matrix $\boldsymbol{X}$, giving matrix $\boldsymbol{Y}$ with covariance matrix $\boldsymbol{I}$ (identity).

orthGS <- function(X) {             # implement Gram-Schmidt algorithm
    Id <- diag(nrow(X))
    for(i in 2:ncol(X)) {
        A <- X[ , 1:(i-1), drop=FALSE]
        Q <- qr.Q(qr(A))
        P <- tcrossprod(Q)
        X[ , i] <- (Id-P) %*% X[ , i]
    scale(X, center=FALSE, scale=sqrt(colSums(X^2)))

Xctr <- scale(X, center=TRUE, scale=FALSE)  # centered version of X
Y    <- orthGS(Xctr)                        # Y is orthonormal

Transform matrix $\boldsymbol{Y}$ to have covariance matrix $\boldsymbol{\Sigma}$ and centroid $\boldsymbol{\mu}$.

Edit: what's going on here: Do a spectral decomposition $\boldsymbol{\Sigma} = \boldsymbol{G} \boldsymbol{D} \boldsymbol{G}^{t}$, where $\boldsymbol{G}$ is the matrix of normalized eigenvectors of $\boldsymbol{\Sigma}$, and $\boldsymbol{D}$ is the corresponding matrix of eigenvalues. Now matrix $\boldsymbol{G} \boldsymbol{D}^{1/2} \boldsymbol{Y}$ has covariance matrix $\boldsymbol{G} \boldsymbol{D}^{1/2} Cov(\boldsymbol{Y}) \boldsymbol{D}^{1/2} \boldsymbol{G}^{t} = \boldsymbol{G} \boldsymbol{D} \boldsymbol{G}^{t} = \boldsymbol{\Sigma}$, as $Cov(\boldsymbol{Y}) = \boldsymbol{I}$.

eig    <- eigen(Sigma)
A      <- eig$vectors %*% sqrt(diag(eig$values))
XX1ctr <- t(A %*% t(Y)) * sqrt(nrow(Y))
XX1    <- sweep(XX1ctr, 2, mu, "+")         # move centroid to mu

Check that the correlation matrix is really $\boldsymbol{R}$.

> all.equal(cor(XX1), R)
[1] TRUE

For other purposes, the question might now be: How do I find a positive definite matrix that is "very similar" to a pre-specified one that is not positive definite. That I don't know.

Edit: corrected some square roots

  • 1
    $\begingroup$ (+1) This answer hits all the main points. $\endgroup$ – cardinal Feb 22 '12 at 17:33

If you're not given a complete collection of correlations, but only a partial set, then you can construct a semidefinite programming feasibility problem that can be used to determine whether or not there is a correlation matrix that has the specified correlations. This can be solved by various SDP solvers such as SeDuMi and SDPT3.


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.