1
$\begingroup$

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?

$\endgroup$
1
$\begingroup$

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

$\endgroup$
  • 1
    $\begingroup$ Could you explain a little more? $\endgroup$ – user5054 Jun 4 '15 at 19:54
  • $\begingroup$ 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) $\endgroup$ – gsmafra Jun 4 '15 at 20:02
1
$\begingroup$

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

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – user5054 Jun 4 '15 at 20:08
1
$\begingroup$

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)
$\endgroup$

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.