# How to generate a sparse inverse covariance matrix for sampling multivariate Gaussian vectors?

I need to generate a sparse 100x100 precision matrix to sample multivariate Gaussian random vectors using the inverse of it as the covariance matrix. To be a valid precision matrix, the matrix I create should be a positive definite matrix, so I regenerate the matrix until it is positive definite (all its eigenvalues are positive). Here is my R code for this job:

library(pracma)
k = 100
sparsity = .2
while (TRUE) {
# generate the symmetric sparsity mask

# generate the symmetric precision matrix
theta = matrix(rnorm(k^2), k)
theta[lower.tri(theta, diag = TRUE)] = 0
theta = theta + t(theta) + eye(k)

# apply the reqired sparsity

if(sum(eigen(theta)$values > 0) == k) { break } else { print('Theta is not positive definite!') } }  The problem is that this code never ends, which means that that kind of valid precision matrix can never be created. What is the way to achieve this job? ## 3 Answers Have you tried doing a random triangular sparse matrix and then using it as the Cholesky decomposition of your covariance matrix? • Could you explain a little more? – user5054 Jun 4 '15 at 19:54 • Every positive-definite matrix has a Cholesky decomposition that takes the form LL' where L is lower triangular (IIRC the inverse is also true), so you could sample L and compute a positive-definite matrix from it. If L is sparse then LL' is also sparse (make sure L is less sparse then what you want your final matrix to be) – gsmafra Jun 4 '15 at 20:02 In your code, since mask is sparse (only 20% are non zero), the odds are high mask is a singular matrix (i.e., one of its eigenvalues is 0). If this happens, theta is guaranteed to be singular since$rank(theta*mask) \leq rank(mask)$and thus your loop will go on forever There are much easier (and strongly recommended) ways of generating sparse matrices (see http://www.johnmyleswhite.com/notebook/2011/10/31/using-sparse-matrices-in-r/) If you insist on writing your own code, Step 1: Chose$K$sparse vectors ($U_{1}, U_{2} ... U_{K}$) with$K$having length = size of matrix you desire. Make sure that (1) The vectors have norm 1 and (2) the locations of the non-zeros entries in these vectors is never the same for any two$U_{i}$and$U_{j}$for$i \neq j$. Step 2: Chose$K$arbitrary positive numbers$a_{1}, a_{2}, a_{3} ...$Step 3: Your sparse matrix$S = \sum_{i=1}^{i=K} a_{i}U_{i}U_{i}^{T}$. If it is not sparse enough, make the$U\$ vectors sparser. This way you get to regulate the rank too if you want to

• sparsity = 0.2 means that the the nonzero entities in the matrix should be 20% of all entities on average. And sorry, I edited the question, "eigenvectors" was a typo, it should be "eigenvalues", as in the R code. – user5054 Jun 4 '15 at 20:08

I edited my question, it had a wrong statement. Actually, when we add (-min(eig(theta))+q) to the diagonal elements (where q is a small positive number) so that all eigenvalues of theta becomes positive, then the problem is solved. So, below code works:

library(pracma)
k = 100
sparsity = .2
# generate the symmetric sparsity mask