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

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$$?

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. • I tried to implement this, but when I apply the is.positive.definite() function, it returns false. Also, I cannot seem to invert the matrix which again shows that it is not positive definite. Commented Nov 19, 2019 at 10:20 • When you divide by$|x|^2$, should you not divide$x^\prime B x$by this constant as well? Commented Nov 19, 2019 at 10:42 • (1) There is no such thing as is.positive.definite in base R, so I cannot comment on it. (2) Yes,$x^\prime B x$is divided by$|x|^2.$When$|x|=1,$that leaves us with$x^\prime B x.$(3) If you choose$\delta$to equal the smallest possible value, you guarantee non-invertibility, because the matrix is only positive semidefinite. That's why the inequality statements in this answer use strict inequality$\lt$rather than partial inequality$\le.\$ Choosing a larger value assures invertibility.
– whuber
Commented Nov 19, 2019 at 13:44
• Great, that clears it up. Thank you Commented Nov 19, 2019 at 16:43