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

Question

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
    mask = rand(k)
    mask = mask * (mask < sparsity)
    mask[lower.tri(mask, diag = TRUE)] = 0
    mask = mask + t(mask) + eye(k)
    mask[mask > 0] = 1

    # 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
    theta = theta * mask

    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?

Answer 1

Have you tried doing a random triangular sparse matrix and then using it as the Cholesky decomposition of your covariance matrix?

Answer 2

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

Answer 3

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
mask = rand(k)
mask = mask * (mask < sparsity)
mask[lower.tri(mask, diag = TRUE)] = 0
mask = mask + t(mask) + eye(k)
mask[mask > 0] = 1

# 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
theta = theta * mask

# force it to be positive definite
theta = theta - (min(eig(theta))-.1) * eye(k)