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.


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

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


1 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
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.

  • $\begingroup$ 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. $\endgroup$ Nov 19, 2019 at 10:20
  • $\begingroup$ When you divide by $|x|^2$, should you not divide $x^\prime B x$ by this constant as well? $\endgroup$ Nov 19, 2019 at 10:42
  • 1
    $\begingroup$ (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. $\endgroup$
    – whuber
    Nov 19, 2019 at 13:44
  • $\begingroup$ Great, that clears it up. Thank you $\endgroup$ Nov 19, 2019 at 16:43

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.