Creating a random sparse precision matrix?

Question

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?

How to create an arbitrary covariance matrix is related, but even if we feed a sparse $A$ into Henry's answer, $\Sigma=A^tA$ will be less sparse, and the precision matrix $\Sigma^{-1}$ will be even less sparse. — Stephan Kolassa, Sep 06 '22 at 06:06

score 2 · Accepted Answer · answered Sep 06 '22 at 14:59

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

krkeane · Answer 2 · 2022-09-10T15:27:12.153

... a $K$ dimensional identity matrix, ...

is the inverse $\rho^{-1}$ of the correlation matrix $\rho$ for $K$ uncorrelated variables ($\rho$ could be scaled to a covariance matrix). For multivariate normal distributions, this starting point specifies a maximum entropy distribution. Your insertion of off diagonal elements reduces the entropy of the synthesized distribution.

I address related topics in my thesis in chapters 5 and 6, and an ICPRAM 2019 conference paper. See thesis equations (6.25) and (6.26). The thesis addresses scenarios with larger correlation patterns than pairwise (stock portfolios), but the punch line is roughly

the inverse correlation matrix $\rho^{-1}$ ... is the sum of the outer products of [K-dimensional unit-length vectors with one or two non-zero entries], scaled by the Lagrange parameters required to enforce [univariate and pairwise correlation] constraints.

You need to think about the distribution of number of off-diagonal elements (pairwise correlation constraints) desired, and the distribution of the partial correlation patterns you are trying to represent in the inverse correlation matrix specifying your synthesized distributions.

Creating a random sparse precision matrix?

2 Answers2