Creating a random sparse precision matrix? In my current project, I want to create a random sparse precision matrix $\boldsymbol{P}=\boldsymbol{\Sigma}^{-1}$ (the inverse of a covariance matrix $\boldsymbol{\Sigma}$). My current procedure consists of using a starting with a $K$-dimensional identity matrix, then setting random entries (symmetric across the diagonal) to random values between $0$ and $1$.
However, this does not seem to yield a valid covariance matrix. If I invert $\boldsymbol{P}$ to obtain $\boldsymbol{\Sigma}=\boldsymbol{P}^{-1}$, this covariance matrix has negative values along the diagonal. Clearly, there is more to creating a valid precision matrix.
Hence my question: What am I missing? Can you recommend a procedure to create valid precision matrices?

At https://stats.stackexchange.com/a/215647/919 I explain how to create random $d\times d$ covariance matrices $\Sigma$ by conjugating random (or specified) diagonal matrices $D$ by random orthogonal matrices $P,$
$$\Sigma = P^\prime D P.$$
The precision matrices $V$ will have the same form, but with the reciprocal values on the diagonal, because
$$V = \Sigma^{-1} = (P^\prime D P)^{-1} = P^{-1} D^{-1} (P^\prime)^{-1} = P^\prime D^{-1} P.$$

When you create those orthogonal matrices as products of a controlled small number $n$ of random rotations you can guarantee sparseness, because each product changes the entries in only two columns.


Using the R function rQ shown below, I simulated precision matrices for combinations of $n$ from $1$ through $1000$ and $d$ from $3$ through $100$ and computed the fraction of zeros among the $d(d-1)$ off-diagonal entries. (Note the logarithmic axes for $n$).  The panel headings indicate the dimensions $d.$

This plot demonstrates how you can control the degree of sparseness by selecting an appropriate value of $n$ for the dimension $d.$  The sparseness starts dropping once $n \approx d/3,$ falls exponentially around $n \approx d,$ and has almost disappeared by the time $n = 3d.$  Thus, creating such a precision matrix typically requires on the order of $O(d)$ or fewer matrix multiplications restricted to two rows, a net $O(d^2)$ calculation--asymptotically the same as the number of entries in the matrix.
#
# Create a sparse orthogonal matrix using `n` right multiplications by
# 2 X 2 rotation matrices.
#
rQ <- function(d, n = 0) {
  alpha <- runif(n, 0, 2 * pi)  # Rotation angles
  #
  # Generate random rotation planes.
  #
  i <- sample.int(d, n, replace = TRUE)
  j <- (sample.int(d-1, n, replace = TRUE) + i - 1) %% d + 1
  #
  # Create the result through sequential right multiplication.
  #
  Q <- diag(1, d, d)
  for (k in seq_len(n)) {
    IJ <- c(i[k], j[k])
    c. <- cos(alpha[k])
    s. <- sin(alpha[k])
    Q[, IJ] <- Q[, IJ] %*% matrix(c(c., s., -s., c.), 2) 
  }
  Q
}
#
# Example.
#
sigma <- 1:3                  # Eigenvalues of the covariance matrix
Q <- rQ(length(sigma), 3)     # A random sparse orthogonal matrix
V <- crossprod(Q, Q / sigma)  # The precision matrix

