Find $\delta$ such that sparse covariance matrix is positive definite

Question

Good day, I was looking through some papers to help with my project assignment that wants me to implements 2 lasso approaches. I am having trouble simulating the samples from a MVN distribution.

$\Theta = B + δI_p \in R_{p\times p}$, where $I_p$ is the identity matrix, each off-diagonal entry in $B$ (symmetric matrix) is generated independently and equals $0.5$ with probability $0.1$ or $0$ with probability $0.9$. $\delta > 0 $is chosen such that $\Theta$ is positive definite. Finally, the matrix is standardized to have unit diagonals (transforming from covariance matrix to correlation matrix). The sparsity pattern in $\Theta$ corresponds to the true edge set $E=\{ (j,l):c_{jl}\neq0, 1 \leq j,l \leq p,j\neq l\}$.

I had a look at some papers that generated this matrix but I had no luck in finding the code. I assume I have to use scio and glasso packages for this from what I gathered and then apply mvrnorn to generate the samples.

EDIT:

Having attempted this, I think I managed to generate the matrix $B$ as described.

set.seed(123)
m = matrix(NA, ncol = 100, nrow = 100)
m[lower.tri(m)] = 0.5*rbinom(100*(100-1)/2,1,0.1)
m[upper.tri(m)] = 0.5*rbinom(100*(100-1)/2,1,0.1)
diag(m) = 0
theta=nearPD(m, keepDiag=FALSE)

I also used the nearPD function to get the $\Theta$ matrix. But I am not sure how I can do this without using a prebuilt function. i.e. use the format of the question $\Theta = B + δI_p \in R_{p\times p}$ and find a suitable $\delta$.

How do I find this $\delta$?

Answer 1

The recipe isn't sufficiently definite to determine $\delta,$ but it does constrain $\delta$ within a well-defined interval.

---

Apply the definition of positive-definiteness: a matrix $A$ is positive-definite if and only if for all nonzero vectors $x,$ $x^\prime A x \gt 0.$ Thus, setting $A=B+\delta I_p,$ we seek any values of $\delta$ for which

$$0 \lt x^\prime(B+\delta I_p) x = x^\prime B x + \delta\left(x^\prime x\right).\tag{*}$$

for all $x \ne 0.$ This simplifies a little because both sides are quadratic forms in $x,$ whence the inequality is preserved upon dividing them both by $|x|^2,$ yielding

$$0 \lt x^\prime B x + \delta$$

for all $x$ with $|x|=1.$

Now the set of all such vectors is the $p-1$-sphere $S^{p-1}$ (a compact manifold) and the function $x\to x^\prime B x$ for $x\in S^{p-1}$ is continuous. Therefore this function attains a minimum, finite value. Indeed, the theory of symmetric quadratic forms shows that $B$ can be diagonalized: there is an orthogonal matrix $Q$ for which

$$x^\prime B x = (Qx)^\prime A (Qx)$$

and $A$ is a diagonal matrix with eigenvalues $\lambda_1 \le \ldots \le \lambda_p$ on the diagonal. That is, in terms of the new variables $(y_1,y_2,\ldots, y_p) = y=Qx,$ the form is the weighted sum of squares

$$y^\prime A y = \lambda_1 y_1^2 + \lambda_2 y_2^2 + \cdots + \lambda_p y_p^2.$$

Because $|y|=|x|=1,$

$$\lambda_1 \le y^\prime A y = (Qx)^\prime A (Qx) = x^\prime B x \le \lambda_p$$

for all $x$ with $|x|=1.$

Therefore, you require that $\delta+\lambda_i \gt 0$ for all $i.$ The answer is now clear:

The solutions to $(*)$ are all $\delta \gt - \lambda_1 = -\min\{\lambda_1,\ldots,\lambda_p\}.$

---

In your code example, the matrix m you create is extremely unlikely to be symmetric, because its upper and lower triangles are filled with independent random values. Corresponding values will agree with probability $0.1^2+0.9^2,$ causing all $\binom{100}{2}$ to agree (and therefore make m symmetric) with probability

$$(0.1^2 + 0.9^2)^{\binom{100}{2}} \lt 10^{-982}.$$

You need a different approach. To create a matrix as described in the quotation, fill one triangle and copy it into the other. Here's one way using base R functions:

p <- 100
set.seed(123)
B <- matrix(sample(c(0,0.5), p^2, replace=TRUE, prob=c(0.9,0.1)), p)
i <- lower.tri(B)
B[i] <- t(B)[i]
diag(B) <- rep(0, p)

The calculation of the smallest possible value of $\delta$ is immediate:

(delta <- -min(eigen(B, symmetric=TRUE, only.values=TRUE)$values))

[1] 2.838141

It is up to you to choose $\delta$ among the (positive) real numbers larger than this.