Hello,
I would like to get you feedback to accelerate my routine glasso for graphical lasso computation (the code follows Friedman et al. 2007 https://arxiv.org/pdf/0708.3517.pdf) where S is the sample covariance matrix and ρ is the scaling factor upfront of the L1 penalty term.
using LinearAlgebra, Statistics, Lasso, SparseArrays, PDMats
function glasso(S, ρ; tol = 1e-3, Niter = 500)
Nx = size(S, 1)
W = S + ρ*I
if Nx == 1
return W
end
W11 = zeros(Nx-1, Nx-1)
b = zeros(Nx-1)
βj = zeros(Nx-1)
Wold = copy(W)
jminus = zeros(Int64, Nx-1)
# Iterate until convergence
for i=1:Niter
for j = 1:Nx
jminus .= setdiff(1:Nx, j)
# Compute a square-root factorization of W11
W11 .= W[jminus, jminus]
W11sqrt = sqrt(W11) #W_11^(1/2)
b .= W11sqrt \ S[jminus,j] # W_11^(-1/2) * s_12
# We need this rescaling by ρ/(Nx-1) due to the definition of the function fit in Lasso.jl,
# that includes a factor 1/(Nx-1) upfront of the log-likelihood term
# fit(LassoPath, X, y, d=Normal(), l=canonicallink(d); ...)
# fits a linear or generalized linear Lasso path given the design
# matrix `X` and response `y`:
# ``\underset{\beta,b_0}{\operatorname{argmin}} -\frac{1}{N} \mathcal{L}(y|X,\beta,b_0) + \lambda\left[(1-\alpha)\frac{1}{2}\|\beta\|_2^2 + \alpha\|\beta\|_1\right]``
# Nx-1 is the number of rows for W11sqrt
βj .= fit(LassoPath, W11sqrt, b, λ=[ρ/(Nx-1)],
standardize=false, intercept=false).coefs
# Update W using that βj is a sparse vector
W[j, jminus] = W11*sparse(βj)
W[jminus, j] = W[j, jminus]
end
# Stop criterion
if norm(W - Wold,2) < tol
break
end
copy!(Wold, W)
end
return W
end
Minimal working example:
Nx = 40
Ne = 100
θX = PDMat(Symmetric(Matrix(2.0*sprand(Nx,Nx,0.3) + 5.0*I)))
ΣX = inv(θX)
πX = Distributions.MvNormal(μX, ΣX)
X = rand(πX, Ne)
S = cov(X')
W = glasso(S, 0.1)