R Programming – Generate Random Positive Definite Matrix with Zero Constraints

matrixr

How to use R to generate a random symmetric positive definite matrix with zero constraints?

For example, I would like to generate a 4 by 4 random symmetric positive definite matrix $\Omega\in\mathbb{R}^{4\times4}$, and we know $\Omega_{1,2}=\Omega_{2,1}=\Omega_{1,3}=\Omega_{3,1} = 0$. How can I do that in R?

What I had in mind is something like Cholesky decomposition $LL^T=\Omega$, where row $L_i$ and row $L_j$ are orthogonal if $\Omega_{ij}=0$. Possibly solve by the Lagrangian multiplier. But I am not really sure how to implement this. Or if this is possible at all.

Best Answer

Every $d\times d$ symmetric positive (semi)definite matrix $\Sigma$ can be factored as

$$\Sigma = \Lambda^\prime\, Q^\prime \,Q\,\Lambda$$

where $Q$ is an orthonormal matrix and $\Lambda$ is a diagonal matrix with non-negative(positive) entries $\lambda_1, \ldots, \lambda_d.$ ($\Sigma$ is always the covariance matrix of some $d$-variate distribution and $QQ^\prime$ will be its correlation matrix; the $\lambda_i$ are the standard deviations of the marginal distributions.)

Let's interpret this formula. The $(i,j)$ entry $\Sigma_{i,j}$ is the dot product of columns $i$ and $j$ of $Q$, multiplied by $\lambda_i\lambda_j.$ Thus, the zero-constraints on $\Sigma$ are orthogonality constraints on the dot products of the columns of $Q.$

(Notice that all diagonal entries of a positive-definite matrix must be nonzero, so I assume the zero-constraints are all off the diagonal. I also extend any constraint on the $(i,j)$ entry to a constraint on the $(j,i)$ entry, to assure symmetry of the result.)

One (completely general) way to impose such constraints is to generate the columns of $Q$ sequentially. Use any method you please to create a $d\times d$ matrix of initial values. At step $i=1,2,\ldots, d,$ alter column $i$ regressing it on all the columns $1, 2, \ldots, i-1$ of $Q$ that need to be orthogonal to it and retaining the residuals. Normalize those results so their dot product (sum of squares) is unity. That is column $i$ of $Q.$

Having created an instance of $Q,$ randomly generate the diagonal of $\Lambda$ any way you please (as discussed in the closely related answer at https://stats.stackexchange.com/a/215647/919).

The following R function rQ uses iid standard Normal variates for the initial values by default. I have tested it extensively with dimensions $d=1$ through $200,$ checking systematically that the intended constraints hold. I also tested it with Poisson$(0.1)$ variates, which--because they are likely to be zero--generate highly problematic initial solutions.

The principal input to rQ is a logical matrix indicating where the zero-constraints are to be applied. Here is an example with the constraints specified in the question.

set.seed(17)
Q <- matrix(c(FALSE, TRUE, TRUE, FALSE,
              TRUE, FALSE, FALSE, FALSE,
              TRUE, FALSE, FALSE, FALSE,
              FALSE, FALSE, FALSE, FALSE), 4)
Lambda <- rexp(4)
zapsmall(rQ(Q, Lambda))

         [,1]      [,2]      [,3]      [,4]
[1,] 2.646156  0.000000  0.000000  2.249189
[2,] 0.000000  0.079933  0.014089 -0.360013
[3,] 0.000000  0.014089  0.006021 -0.055590
[4,] 2.249189 -0.360013 -0.055590  4.167296

As a convenience, you may pass the diagonal of $\Lambda$ as the second argument to rQ. Its third argument, f, must be a random number generator (or any other function for which f(n) returns a numeric vector of length n).

rQ <- function(Q, Lambda, f=rnorm) {
  normalize <- function(x) {
    v <- zapsmall(c(1, sqrt(sum(x * x))))[2]
    if (v == 0) v <- 1
    x / v
  }
  Q <- Q | t(Q)                    # Force symmetry by applying all constraints
  d <- nrow(Q) 
  if (missing(Lambda)) Lambda <- rep(1, d)
  R <- matrix(f(d^2), d, d)        # An array of column vectors
  for (i in seq_len(d)) {
    j <- which(Q[seq_len(i-1), i]) # Indices of the preceding orthogonal vectors
    R[, i] <- normalize(residuals(.lm.fit(R[, j, drop=FALSE], R[, i])))
  }
  R <- R %*% diag(Lambda)
  crossprod(R)
}
Related Question