Learning Graph Structures, Graphical Lasso and Its Applications - Part 7: Julia implementation of Glasso
Published:
The Julia code here is for illustration purpose only and tries to be a literal translation of the algorithm steps in previous sections for easy understanding. It tries to be a minimalistic walkthrough of the algorithm rather than a numerically stable implementation, which would become topics for following up articles.
Efficiency and numerical stability have not yet been taken care of, which however will be the topic for the next article, which aims to venture into the big data world and tries to solve the inverse covariance matrix estimation problem efficiently when faced with huge number of variables. Detailed performance calibration and comparison will be documented there too.
#Import lasso solver: https://github.com/JuliaStats/Lasso.jl
using Lasso
function glasso(S, lambda, maxiter=10, tol=1e-5)
if lambda == 0
return S, inv(S)
else
W = deepcopy(S)
_, p = size(W)
W_0 = zeros(p, p)
#Loop over the max number of iterations:
for iter in 1:maxiter
#Iterate over each dimension
for dim in 1:p
#Single out W11 and S12, input to Lasso
W11 = W[1:end.!=dim, 1:end.!=dim]
S12 = S[dim, 1:end.!=dim]
f = Lasso.fit(LassoPath, sqrt(W11), inv(sqrt(W11))*S12)
coefs = f.coefs[1:end,end]
#Update current W
W_0[dim, 1:end.!=dim] = coefs
W_0[1:end.!=dim, dim] = coefs
W_0[dim, dim] = S[dim, dim]
end
W = deepcopy(W_0)
end
#T is the inverse covariance matrix
T = Array{Float64}(undef, p, p)
#Update T using W
for dim in 1:p
coefs = inv(W[1:end.!=dim, 1:end.!=dim])*W[1:end.!=dim, dim]
T[dim, dim] = 1.0/(W[dim, dim] - W[dim, 1:end.!=dim]'coefs)
T[1:end.!=dim, dim] = -T[dim, dim]*coefs
T[dim, 1:end.!=dim] = -T[dim, dim]*coefs
end
return W, T
end
end
#Testing with data
println("Running glasso")
#S = [4 1.2 1.6; 1.2 9 2.4; 1.6 2.4 16]
#W, T = glasso(S, 0.0001)
#println(W)
#println(T)Checking our results with those of the glasso package in R and the one in Python scikit-learn, we see that the results are identical.
Stay tuned for the next article, which focuses on optimizing the algorithm for big data and building an application on it.