Generate normally distributed random numbers with non positive-definite covariance matrix

Question

I estimated the sample covariance matrix $C$ of a sample and get a symmetric matrix. With $C$, I would like to create $n$-variate normal distributed r.n. but therefore I need the Cholesky decomposition of $C$. What should I do if $C$ is not positive definite?

Answer 1

Solution Method A:

1. If C is not symmetric, then symmetrize it. D <-- $0.5(C + C^T)$
2. Add a multiple of the Identity matrix to the symmetrized C sufficient to make it positive definite with whatever margin, m, is desired, i.e., such that smallest eigenvalue of new matrix has minimum eigenvalue = m. Specifically, D <-- $D + (m - min(eigenvalue(D)))I$, where I is the identity matrix. D contains the desired positive definite covariance matrix.

In MATLAB, the code would be

D = 0.5 * (C + C');
D =  D + (m - min(eig(CD)) * eye(size(D));

Solution Method B: Formulate and solve a Convex SDP (Semidefinite Program) to find the nearest matrix D to C according to the frobenius norm of their difference, such that D is positive definite, having specified minimum eigenvalue m.

Using CVX under MATLAB, the code would be:

n = size(C,1);
cvx_begin
variable D(n,n)
minimize(norm(D-C,'fro'))
D -m *eye(n) == semidefinite(n)
cvx_end

Comparison of Solution Methods: Apart from symmetrizing the initial matrix, solution method A adjusts (increases) only the diagonal elements by some common amount, and leaves the off-diagonal elements unchanged. Solution method B finds the nearest (to the original matrix) positive definite matrix having the specified minimum eigenvalue, in the sense of minimum frobenius norm of the difference of the positive definite matrix D and the original matrix C, which is based on the sums of squared differences of all elements of D - C, to include the off-diagonal elements. So by adjusting off-diagonal elements, it may reduce the amount by which diagonal elements need to be increased, and diagoanl elements are not necessarily all increased by the same amount.

Answer 2

The question concerns how to generate random variates from a multivariate Normal distribution with a (possibly) singular covariance matrix $\mathbb{C}$. This answer explains one way that will work for any covariance matrix. It provides an R implementation that tests its accuracy.

---

Algebraic analysis of the covariance matrix

Because $\mathbb{C}$ is a covariance matrix, it necessarily is symmetric and positive-semidefinite. To complete the background information, let $\mu$ be the vector of desired means.

Because $\mathbb{C}$ is symmetric, its Singular Value Decomposition (SVD) and its eigendecomposition will automatically have the form

$$\mathbb{C} = \mathbb{V\, D^2\, V^\prime}$$

for some orthogonal matrix $\mathbb{V}$ and diagonal matrix $\mathbb{D}^2$. In general the diagonal elements of $\mathbb{D}^2$ are nonnegative (implying they all have real square roots: choose the positive ones to form the diagonal matrix $\mathbb{D}$). The information we have about $\mathbb{C}$ says that one or more of those diagonal elements are zero--but that won't affect any of the subsequent operations, nor will it prevent the SVD from being computed.

Generating multivariate random values

Let $X$ have a standard multivariate Normal distribution: each component has zero mean, unit variance, and all covariances are zero: its covariance matrix is the identity $\mathbb{I}$. Then the random variable $Y=\mathbb{VD}X$ has covariance matrix

$$\operatorname{Cov}(Y) = \mathbb{E}(Y Y^\prime) = \mathbb{E}(\mathbb{V D}X\, X^\prime \mathbb{D^\prime V^\prime}) = \mathbb{V D}\mathbb{E}(X X^\prime)\mathbb{D V^\prime} = \mathbb{V D I D V^\prime} = \mathbb{V D^2 V^\prime} = \mathbb{C}.$$

Consequently the random variable $\mu + \mathbb{Y}$ has a multivariate Normal distribution with mean $\mu$ and covariance matrix $\mathbb{C}$.

Computation and Example code

The following R code generates a covariance matrix of given dimensions and rank, analyzes it with the SVD (or, in commented-out code, with an eigendecomposition), uses that analysis to generate a specified number of realizations of $Y$ (with mean vector $0$), and then compares the covariance matrix of those data to the intended covariance matrix both numerically and graphically. As shown, it generates $10,000$ realizations where the dimension of $Y$ is $100$ and the rank of $C$ is $50$. The output is

        rank           L2 
5.000000e+01 8.846689e-05 

That is, the rank of the data is also $50$ and the covariance matrix as estimated from the data is within distance $8\times 10^{-5}$ of $C$--which is close. As a more detailed check, the coefficients of $C$ are plotted against those of its estimate. They all lie close to the line of equality:

The code exactly parallels the preceding analysis and so should be self-explanatory (even to non-R users, who might emulate it in their favorite application environment). One thing it reveals is the need for caution when using floating-point algorithms: the entries of $\mathbb{D}^2$ can easily be negative (but tiny) due to imprecision. Such entries need to be zeroed out before computing the square root to find $\mathbb{D}$ itself.

n <- 100         # Dimension
rank <- 50
n.values <- 1e4  # Number of random vectors to generate
set.seed(17)
#
# Create an indefinite covariance matrix.
#
r <- min(rank, n)+1
X <- matrix(rnorm(r*n), r)
C <- cov(X)
#
# Analyze C preparatory to generating random values.
# `zapsmall` removes zeros that, due to floating point imprecision, might
# have been rendered as tiny negative values.
#
s <- svd(C)
V <- s$v
D <- sqrt(zapsmall(diag(s$d)))
# s <- eigen(C)
# V <- s$vectors
# D <- sqrt(zapsmall(diag(s$values)))
#
# Generate random values.
#
X <- (V %*% D) %*% matrix(rnorm(n*n.values), n)
#
# Verify their covariance has the desired rank and is close to `C`.
#
s <- svd(Sigma <- cov(t(X)))
(c(rank=sum(zapsmall(s$d) > 0), L2=sqrt(mean(Sigma - C)^2)))

plot(as.vector(C), as.vector(Sigma), col="#00000040",
     xlab="Intended Covariances",
     ylab="Estimated Covariances")
abline(c(0,1), col="Gray")

Answer 3

One way would be to compute the matrix from an eigenvalue decomposition. Now I'll admit I don't know too much of the Math behind these processes but from my research it seems fruitful to look at this help file:

http://stat.ethz.ch/R-manual/R-patched/library/Matrix/html/chol.html

and some other related commands in R.

Also, check out 'nearPD' in the Matrix package.

Sorry I couldn't be of more help but I hope my searching around can help push you in the right direction.

Answer 4

I would begin by thinking about the model you are estimating.

If a covariance matrix is not positive semi-definite, it may indicate that you have a colinearity problem in your variables which would indicate a problem with the model and should not necessarily be solved by numerical methods.

If the matrix is not positive semidefinite for numerical reasons, then there some solutions which can be read about here

Answer 5

You can get the results from the nearPD function in the Matrix package in R. This will give you a real valued matrix back.

library(Matrix)
A <- matrix(1, 3,3); A[1,3] <- A[3,1] <- 0
n.A <- nearPD(A, corr=T, do2eigen=FALSE)
n.A$mat

# 3 x 3 Matrix of class "dpoMatrix"
#           [,1]      [,2]      [,3]
# [1,] 1.0000000 0.7606899 0.1572981
# [2,] 0.7606899 1.0000000 0.7606899
# [3,] 0.1572981 0.7606899 1.0000000