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)
}
