4
$\begingroup$

I am using glasso function from glasso package, as follow:

obj <- glasso(var(X), rho = 0.09, zero = info, approx = TRUE)

Regardless of rho value, all of entries in obj$w, estimated covariance matrix, are zero. Do you have any idea why this happens?

For your information, the dimension of var(X) is 1990 x 1990 and the number of rows in info is 1959841.

EDIT: You can download X and info variables as RData file from this link: https://www.dropbox.com/s/t9s4iw6ulbys72o/varX.RData?dl=0

$\endgroup$
6
  • $\begingroup$ I don't know anything about this, but a quick perusal of cran.r-project.org/web/packages/glasso/glasso.pdf , shows that your 'info' is being used to specify entries in the estimated (inverse) covariance matrix which are constrained to be zero. Are you constraining all entries to be zero, or at least enough of them that all zeros is the only solution it can find? What are the contents of 'info'? $\endgroup$ Commented Apr 23, 2016 at 16:01
  • $\begingroup$ @MarkL.Stone Indeed I force the algorithm to make most of the entries zero and only about 5300 entries are non-zero which I would like to estimate. $\endgroup$ Commented Apr 23, 2016 at 16:05
  • $\begingroup$ First of all, per the documentation, each element (k,j) is constrained to be zero if (j,k) is constrained to be zero. Does that still leave you with 5300 elements not constrained to be zero? Second of all, however many unconstrained elements you think that leaves you, I presume the estimated matrix must be positive semi-definite (PSD). The elements constrained to be zero can indirectly constrain other elements of the matrix due to the PSD constraint. Are you constraining any diagonal elements to be zero? I think you will find the specification of elements constrained to be zero to be key. $\endgroup$ Commented Apr 23, 2016 at 16:18
  • $\begingroup$ @MarkL.Stone Thanks. The weird thing is that even by running glasso(var(X), rho = 0.01, zero = info[-c(1000:1950000),], approx=TRUE), the estimated cov matrix is still zero!! $\endgroup$ Commented Apr 23, 2016 at 16:37
  • $\begingroup$ What happens when you set zero = NULL (i.e., not specify zero)? Sorry, I don't know what zeros you have specified when you do zero = info[-c(1000:1950000),]. Are you leaving out your first 999 constraints (whatever they are)? $\endgroup$ Commented Apr 23, 2016 at 16:44

3 Answers 3

1
$\begingroup$

Both other answers are wrong and misunderstand the Meinshausen–Bühlmann approach. The method is implemented correctly, but does not return what you might expect. It is also not necessarily faster.


A quick answer

If you use approx = TRUE, proceed as follows:

A <- glasso(S, rho = 0, approx = TRUE)$wi
Pcor <- 1/2 * (A + t(A))

That's it. You now have a matrix of approximate partial correlations.

The covariance matrix is not estimated by this method, which is why $w is all zero.

Comparison to an exact solution

Let's use some example data and $\lambda = 0$ to see that this does indeed do what we want.

Example data:

set.seed(2024)
n  <- 20
p  <- 3
X  <- matrix(rnorm(n * p), nrow = n, ncol = p)
Xc <- scale(X, center = TRUE, scale = FALSE)
S  <- 1/n * crossprod(Xc)

Meinshausen–Bühlmann approach:

library("glasso")
MB <- glasso(S, rho = 0, approx = TRUE)
(Pcor <- 1/2 * (MB$wi + t(MB$wi)))
#            [,1]       [,2]       [,3]
# [1,]  0.0000000  0.2150574 -0.2212086
# [2,]  0.2150574  0.0000000 -0.1252797
# [3,] -0.2212086 -0.1252797  0.0000000

Exact solution:

P <- solve(S)
P <- -P
diag(P) <- -diag(P)
(Pcor <- cov2cor(P))
#            [,1]       [,2]       [,3]
# [1,]  1.0000000  0.2057225 -0.2164845
# [2,]  0.2057225  1.0000000 -0.1247561
# [3,] -0.2164845 -0.1247561  1.0000000

Some background

Inverting a large matrix can be slow, even with the coordinate descent algorithm used in glasso. The option approx = TRUE uses an alternative based on neighborhood selection (Meinshausen & Bühlmann, 2006). Importantly, it does not estimate the empirical covariance matrix $\mathbf{S}$, nor the precision matrix $\mathbf\Theta$. Instead, it relies on the observation that partial correlations are closely related to regression coefficients (1), and attempts to estimate non-zero off-diagonal elements directly using (LASSO-penalized) regression.

You can demonstrate that this is what happens by implementing your own version:

library("glmnet")
A <- matrix(0, nrow = p, ncol = p)
for(i in 1:p){
  y_current <- X[, i]
  X_current <- X[, -i]
  LASSO <- glmnet(X_current, y_current, alpha = 1, lambda = 0)
  A[-i, i] <- coef(LASSO)[-1]
}
(Pcor <- 1/2 * (A + t(A)))
#            [,1]       [,2]       [,3]
# [1,]  0.0000000  0.2150574 -0.2212086
# [2,]  0.2150574  0.0000000 -0.1252797
# [3,] -0.2212086 -0.1252797  0.0000000

(Compare this to the output of glasso(..., approx = TRUE) above.)

The asymmetry in the solution is a consequence of the fact that regressing y ~ x $\neq$ x ~ y. Some papers try using an OR rule or AND rule to decide which elements should be non-zero, like in Banerjee et al. (2008), but the symmetrization trick shown here provides a much more reasonable approximation in my opinion.

When should you use Meinshausen–Bühlmann (approx = TRUE)?

The main advantage over the block-coordinate descent algorithm in glasso is that each element $\theta_{ij}$ is estimated separately, so this is very easy to parallelize. However, Witten et al. (2011) showed that there is a simple way to check whether $\mathbf\Theta$ is block-diagonal for a given value of $\lambda$. If it is, then we can separately estimate each block $k \in K$ of $$\mathbf\Theta = \begin{pmatrix} \mathbf\Theta_1 & & & \\ & \mathbf\Theta_2 & & \\ & & \ddots & \\ & & & \mathbf\Theta_K\end{pmatrix}.$$

The authors conclude:

Strikingly, when the graphical lasso is performed with a large tuning parameter value, then the algorithms proposed in this article lead to such massive speed improvements that computing the exact graphical lasso solution is much faster even than computing the Meinshausen and Buhlmann (2006) approximate solution.

From glasso 1.7 on, if you expect $\mathbf\Theta$ to be sparse, approx = FALSE can be orders of magnitude faster.


References

Meinshausen, N., & Bühlmann, P. (2006). High-dimensional graphs and variable selection with the Lasso. Annals of Statistics, 34(3), 1436–1462. https://doi.org/10.1214/009053606000000281

Banerjee, O., El Ghaoui, L., & d’Aspremont, A. (2008). Model Selection Through Sparse Maximum Likelihood Estimation for Multivariate Gaussian or Binary Data. J. Mach. Learn. Res., 9, 485–516. https://www.jmlr.org/papers/volume9/banerjee08a/banerjee08a.pdf

Witten, D. M., Friedman, J. H., & Simon, N. (2011). New Insights and Faster Computations for the Graphical Lasso. Journal of Computational and Graphical Statistics, 20(4), 892–900. http://www.jstor.org/stable/23248939

$\endgroup$
0
$\begingroup$

I was having the same problem and discovered it was related to the parameter approx=TRUE, which lets glasso use an approximation method used to speed up the inversion of the covariance matrix, so the problem may be related to the original covariance matrix s not fulfilling the assumptions used for this method.

In my case my matrix wasn't that big (800x800), so I was able to fix the problem just by setting the approx parameter to FALSE and being a little patient, but I know this may not be feasible for larger matrices, so I hope someone could share a better solution for this problem.

$\endgroup$
1
  • 1
    $\begingroup$ This answer is incorrect, because approx = TRUE does not estimate the inverse covariance matrix at all. The approximation is for deciding which off-diagonal elements are non-zero, and has a relation to the matrix of partial correlations. $\endgroup$ Commented Aug 7 at 7:43
-1
$\begingroup$

I did some experimentation with the glasso function and I find that it does not return the covariance matrix w, and the inverse covariance matrix wi lacks the diagonal terms and is not symmetric.

Setting up ...

library(glasso)

n <- 20
k <- 3
x <- matrix(rnorm(n*k),ncol=k)

Since $n > k,$ the covariance matrix can be inverted directly.

> var(x)
            [,1]        [,2]       [,3]
[1,]  1.37818760 -0.03676559 -0.0446064
[2,] -0.03676559  1.15692561 -0.4592093
[3,] -0.04460640 -0.45920931  1.1981595
> solve(var(x))
           [,1]       [,2]       [,3]
[1,] 0.72803016 0.03997524 0.04242491
[2,] 0.03997524 1.02163778 0.39304343
[3,] 0.04242491 0.39304343 0.98683155

With approx=FALSE the glasso function uses the lasso method and gets the right answer.

> glasso(var(x),rho=0,approx=FALSE)[c('w','wi')]
$w
            [,1]        [,2]        [,3]
[1,]  1.37818760 -0.03676557 -0.04460639
[2,] -0.03676557  1.15692561 -0.45920931
[3,] -0.04460639 -0.45920931  1.19815952

$wi
           [,1]       [,2]       [,3]
[1,] 0.72803010 0.03997522 0.04242489
[2,] 0.03997386 1.02163778 0.39304343
[3,] 0.04242438 0.39304343 0.98683155

Setting approx=TRUE causes the function to use the Meinshausen-Buhlmann method, and that seems to be poorly implemented.

> glasso(var(x),rho=0,approx=TRUE)[c('w','wi')]
$w
     [,1] [,2] [,3]
[1,]    0    0    0
[2,]    0    0    0
[3,]    0    0    0

$wi
            [,1]        [,2]        [,3]
[1,]  0.00000000 -0.03912856 -0.04299102
[2,] -0.05490688  0.00000000 -0.39828827
[3,] -0.05827283 -0.38471896  0.00000000

Now that I look at it, the wi returned by setting approx=TRUE looks like the covariance, not inverse covariance, matrix. This is a serious error in glasso function. Or maybe it's an error in the documentation. Either way, it's serious.

$\endgroup$
1
  • $\begingroup$ (-1). I am downvoting this answer because it claims to have found a "serious error" while the method is working exactly as intended. $\endgroup$ Commented Aug 7 at 7:49

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.