Problems with Graphical Lasso

Question

I'm trying to use the Graphical Lasso algorithm (more specifically the R package glasso) to find an estimated graph representing the connections between a set of nodes by estimating a precision matrix. I have a feature matrix containing the values of multiple features for each of the nodes, and the sample covariance matrix obtained from the product between this matrix and its tranpose is used as the input for the glasso function, along with the l1 regularization coefficient $\lambda$.

However, when testing this with some simple examples I'm not getting the expected results. For example, when using the following feature matrix (5 nodes with 2 features each) as input:

$\begin{bmatrix} 1 & 10 \\ 2 & 10 \\ 3 & 10 \\ 4 & 10 \\ 5 & 10 \end{bmatrix},$

with $\lambda=0.01$, I get (approximately) the following inverse covariance matrix as output:

$P=\begin{bmatrix} 19.8 & -18 & -6.3 & 0 & 2.2 \\ -18 & 34.5 & -12.1 & -4.5 & 0 \\ -6.3 & -12.1 & 36.1 & -12.8 & -3.8 \\ 0 & -4.5 & -12.8 & 34.8 & -15.9 \\ 2.2 & 0 & -3.8 & -15.9 & 20.7 \end{bmatrix}.$

As far as I understand, a value of 0 in this precision matrix indicates that the two corresponding nodes are conditionally independent. As such, it makes sense for entries $p_{1,3}$ and $p_{2,4}$ (along with their symmetrical equivalents) to be zero. However, I don't understand why entry $p_{1,4}$ has a greater value than entry $p_{1,3}$, for example, considering that the feature values in the first node are closer to those of the third node than the fourth. Moreover, I would like have more entries with a value of zero, which means I want a sparser matrix. Therefore, I tried increasing the value of $\lambda$ to 0.03, which results in the following precision matrix:

$P=\begin{bmatrix} 8.3 & -5.8 & -3.2 & -0.6 & 0 \\ -5.8 & 11.9 & -3.6 & -2.5 & -0.1 \\ -3.2 & -3.6 & 13.2 & -3.3 & -2.4 \\ -0.6 & -2.5 & -3.3 & 12.6 & -4.7 \\ 0 & -0.1 & -2.4 & -4.7 & 9.6 \end{bmatrix}.$

Now there are less entries with a value of zero (so the sparsity decreased rather than increasing), and while there was a general decrease in the value of every entry, the larger entries decreased at a higher rate, leading to a more evenly distributed matrix. This is not consistent with the feature selection I'm used to seeing in standard lasso regularization, and looks more like some sort of l2 norm regularization.

I feel like there must be something fundamental that I'm missing completely here. Is this method not supposed to be applied in this way?

Edit: the code I'm using:

D=matrix(c(1, 10, 2, 10, 3, 10, 4, 10, 5 , 10),nrow=5,ncol=2, byrow=TRUE)
m=dim(D)[2]
D = D - mean(D)
covar = (1/m)*D%*%t(D)
D = D/sqrt(max(covar))
covar = (1/m)*D%*%t(D)

a=glasso(s=covar,rho=0.01)

Maybe there's something wrong with the data centering/scaling I'm doing?

Answer 1

(Answering some old questions about glasso in case someone is interested.)

General answer

While glasso is typically used to estimate the inverse covariance matrix $\mathbf\Theta = \mathbf\Sigma^{-1} = \frac{1}{n}(\mathbf{X}^T\mathbf{X})^{-1}$, the question attempts to estimate $\frac{1}{p}(\mathbf{XX}^T)^{-1}$. However, this does not matter for the answer. In fact, the answer is not even specific to glasso, but applies to LASSO regression as well.

There are two (common) misconceptions about (g)lasso in this question:

1. The LASSO penalty selects out less important coefficients.
2. Coefficients monotonically decrease with increased values of $\lambda$ (e.g., here, here and here).

As to (1), people are quick to point to correlation among variables, but the issue is larger than that. LASSO has no way of knowing which estimates are important. A penalty is applied to the sum of all coefficients. It just so happens that coefficients that were smaller to begin with often shrink to $0$ faster.

Regarding (2), increasing the amount of regularization $\lambda$ typically decreases the size of the coefficients, and ultimately always results in a smaller sum of coefficients, because $$\lim_{\lambda \to \infty} \sum\beta=0.$$

However, shrinkage of the coefficients is not a monotone function of $\lambda$. See: [a], [b], [c]. This happens to be the case if $\mathbf X$ is orthonormal, but it does not hold in general. As a counterexample, see figure 3.15 from Elements of Statistical Learning (p. 75). Another great example of the non-monotonic path taken by LASSO can be seen here.

---

The effect of $\lambda$ in the question

Understanding what happens to the coefficients as you vary $\lambda$ is easier using a regularization path. The glasso package has a built-in function for this called glassopath:

Zooming in shows how $\hat{\theta}_{2, 4}$ has become non-zero by increasing the amount of regularization to $\lambda = 0.03$:

As explained in the first part, this is not unexpected behavior.

Note that neither $\lambda = 0.01$, nor $\lambda = 0.03$ is enough to turn $\mathbf{S} = \hat{\mathbf{\Sigma}}$ into a well-conditioned matrix. A rule of thumb often used is that the condition number should be $$\kappa(\mathbf{S}) \leq 30.$$

We can assess this by plotting the condition number against $\lambda$:

Apparently, a better choice would be $\lambda \gtrapprox 0.1$. We could also try to find the optimal value using cross-validation, as illustrated here: How to determine lambda for graphical lasso?.

---

Code used

require("glasso")
require("plyr")
require("sfsmisc")
D <- matrix(c(1, 10, 2, 10, 3, 10, 4, 10, 5 , 10), nrow = 5,ncol=2, byrow=TRUE)
m <- dim(D)[2]
D <- D - mean(D)
covar <- 1 / m * D %*% t(D)
D <- D / sqrt(max(covar))
covar <- 1 / m * D %*% t(D)
lambdas <- exp(seq(-5, 0, l = 100))
a <- glassopath(s = covar, rho = lambdas, trace = FALSE)$wi
CN <- apply(a, 3, function(x){
  1 / rcond(x, norm = "O")
})
p <- 5
k <- p * (p - 1) / 2
wh <- which(upper.tri(a[, , 1]), arr.ind = TRUE)
a2 <- aaply(a, 3, function(x){
  diag(x) <- 0
  return(x)
})
par(mar = c(5, 4, 0, 0) + 0.1)
plot(NA, xlim = range(lambdas), ylim = range(a2), log = "x", 
     axes = FALSE, xlab = bquote(lambda), ylab = bquote(theta[ij]))
abline(v = c(0.01, 0.03), col = 2, lty = 3)
abline(h = 0, col = "lightgrey")
for(i in 1:k){
  current <- a[wh[i, 1], wh[i, 2], ]
  lines(current ~ lambdas)
}
eaxis(1, n.axp = 1)
axis(1, 0.03, bquote(3 %*% 10^{-2}))
eaxis(2)
plot(NA, xlim = c(2e-2, 4e-2), ylim = c(-1, 1), log = "x", 
     axes = FALSE, xlab = bquote(lambda), ylab = bquote(theta[ij]))
abline(v = c(0.01, 0.03), col = 2, lty = 3)
abline(h = 0, col = "lightgrey")
for(i in 1:k){
  current <- a[wh[i, 1], wh[i, 2], ]
  lines(current ~ lambdas)
}
axis(1, 0.03, bquote(3 %*% 10^{-2}))
eaxis(2)
plot(CN ~ lambdas, log = "x", axes = FALSE, type = "l", lwd = 3, 
     xlab = bquote(lambda), ylab = "Condition number")
eaxis(1, n.axp = 1)
axis(1, 0.03, bquote(3 %*% 10^{-2}))
eaxis(2)
abline(h = 30, col = 4)
abline(v = c(0.01, 0.03), col = 2, lty = 3)
par(family = "serif")
text(0.5, 50, bquote(kappa*(bold(S)) == 30), col = 4)
par(family = "sans", mar = c(5, 4, 4, 3) + 0.1)

---

Reference

Hastie, T., Tibshirani, R., & Friedman, J. (2009). The elements of statistical learning (2nd ed.) https://doi.org/10.1007/978-0-387-84858-7