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

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?

• What is the difference with this question stackoverflow.com/questions/17295627/… ? – dickoa Jul 9 '13 at 17:31
• Positive-semidefinite matrices have multiple square roots (see the explanation at the end of stats.stackexchange.com/a/71303/919, for instance). You don't necessarily need the one produced by the Cholesky decomposition. Therein lies the heart of the problem: find a method to compute square roots that works even when the matrix is singular. @amoeba The title suggests your interpretation is correct. – whuber Aug 19 '16 at 14:23

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.

• I'm trying to implement Method B in Matlab but I get the following error: Error: "D" was previously used as a variable, conflicting with its use here as the name of a function or command – TEX Feb 26 '20 at 18:26
• Could you help me to understand why? I'm not using any variable D before running the script above. – TEX Feb 26 '20 at 18:27
• It sounds like something is messed up with with either your MATLAB session, (that happens sometimes) or the CVX installation. Try this in a fresh MATLAB session.If that doesn't work, try reinstalling CVX Use CVX 2.2. – Mark L. Stone Feb 26 '20 at 22:33

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")  • +1 but when you say "indefinite" in your first sentence, what exactly do you mean? I checked in Wikipedia and it says that positive semidefinite is not indefinite, i.e. indefinite means that C has both positive and negative eigenvalues. Is that what you mean there? – amoeba Aug 19 '16 at 14:23 • @amoeba Yes, that was a slip. Thanks for noticing. "Indefinite" means the signature of the matrix has both positive and negative signs, whereas "semidefinite" means the signature has only one sign. – whuber Aug 19 '16 at 14:25 • When$C$is not positive semidefinite with small values, this algorithm can flip signs of the resulting covariance matrix because SVD returns values up to a sign. Any comments on fixing this? – Tim Apr 1 '20 at 11:56 • @Tim I don't follow, because this algorithm makes no sense when$C$is not positive semidefinite and I don't understand which "signs of the ... covariance matrix" you might be referring to. – whuber Apr 1 '20 at 12:53 • OK, I missed that this is what you assume. I was thinking of a case where$C$is an inverse of a Hessian and it does not the properties, so like in the question it is not positive semidefinite. – Tim Apr 1 '20 at 13:21 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. • Hi, thx for the links. Respective to the eigen value decomposition, this decomposition does not help, because from there you get complex eigenvalues for square root matrix, but I need reell valued matrix. – Klaus Jul 11 '13 at 23:20 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 • The assumption is that the model is a linear mixed model. And for this case its not relevant to find a correct model for the data, rather the data is given as an example for some calculation. Now there is the possibility that you get a non positive semidefinite matrix as estimation for the covaraince. So what to do from there, if I want to figure out the covariance from the normal distributed population where the data comes from. That the sample is normal distributed is the assuption. – Klaus Jul 10 '13 at 23:23 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

• For users of R.. this might not be a bad "poor man's" version (with less control) of Solution Method B in my answer. – Mark L. Stone Aug 19 '16 at 1:46
• I agree this is not optimal but sometimes it does the trick. – Dr. Mike Aug 19 '16 at 8:48
• @Dr.Mike, I would like to have your input in one of your blog post here in SE – GENIVI-LEARNER Feb 9 '20 at 13:52