# Implementation of graphical lasso with Cholesky factor

**URL:** <https://discourse.julialang.org/t/implementation-of-graphical-lasso-with-cholesky-factor/99016>\
**Category:** Statistics\
**Created:** [May 17, 2023, 8:28pm UTC](https://discourse.julialang.org/t/implementation-of-graphical-lasso-with-cholesky-factor/99016 "2023-05-17T20:28:47Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![mleprovost](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mleprovost/32/7166_2.png) [@mleprovost](https://discourse.julialang.org/u/mleprovost)\
**Post date:** [May 17, 2023, 8:28pm UTC](https://discourse.julialang.org/t/implementation-of-graphical-lasso-with-cholesky-factor/99016/1 "2023-05-17T20:28:47Z")

</div>

Hello,

I am working on an implementation of the graphical Lasso (Glasso) ([https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3019769/pdf/kxm045.pdf](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3019769/pdf/kxm045.pdf)) in pure Julia. In the algorithm, we have to use _a square root_ of sample covariance matrix `W11`. I am thinking of using the Cholesky factor as it will speed-up calculations for the lasso optimization. I am getting error that `W11` is not positive definite. How would you address this?

```julia
using SCS
using Convex
using Distributions
using LinearAlgebra
using Statistics 

function glasso(X, ρ; tol = 1e-2, Niter = 100, withθ = true)
    Nx, Ne = size(X)
    
    S = cov(X')
    
    W = S + ρ*I
    Wold = copy(W)
    W11 = zeros(Nx-1, Nx-1)
    b = zeros(Nx-1)
    B = zeros(Nx-1, Nx-1)
    w12 = zeros(Nx-1)
    
    
    
    # Iterate until convergence
    for i=1:Niter
        for j=1:Nx
            jminus = setdiff(1:Nx, j)
            
            # Compute a Cholesky factorization of W11
            W11 .= W[jminus, jminus]
            W11sqrt .= LowerTriangular(cholesky(Symmetric(W11)).L)
            # W11svd = svd(W11)
            # W11sqrt = W11svd.U * Diagonal(sqrt.(W11svd.S)) * W11svd.Vt #W_11^(1/2)
            
            b .= 0.5*(W11sqrt \ S[jminus,j]) # W_11^(-1/2) * s_12
            
            # Solve the lasso problem
            β_lasso = Variable(Nx-1)
            
            lasso_problem = minimize(sumsquares(W11sqrt*β_lasso - b) + ρ*norm(β_lasso, 1))
            
            
            # Solve the problem by calling solve!
            solve!(lasso_problem, SCS.Optimizer; silent_solver = true)
            
            # TODO check resolution status
            # @assert lasso_problem.status 
            
            w12 .= 2.0*(W11*evaluate(β_lasso))
            
            # Update W
            W[j, jminus] .= copy(w12)
            W[jminus, j] .= copy(w12)
        end
        
        # Stop criterion
        if norm(W - Wold,1) < tol
            break 
        end
        Wold = copy(W)
    end
    
    return W
# if withθ == true
# # Build precision matrix
# θ = zeros(Nx, Nx)
# end
end

```

```julia
Nx = 40
μX = zeros(Nx)
ΣX = PDiagMat(ones(Nx))
πX = Distributions.MvNormal(μX, ΣX)

Ne = 100
X = rand(πX, Ne)

W = glasso(X, 0.01)

```

---

<div class="post-metadata">

**Author:** ![huang\_min](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/huang_min/32/34337_2.png) [@huang\_min](https://discourse.julialang.org/u/huang_min)\
**Post date:** [May 20, 2023, 1:19pm UTC](https://discourse.julialang.org/t/implementation-of-graphical-lasso-with-cholesky-factor/99016/2 "2023-05-20T13:19:22Z")

</div>

[https://discourse.julialang.org/t/glm-jl-logisticregression-errors-matrix-is-not-positive-definite-cholesky-factorization-failed](https://discourse.julialang.org/t/glm-jl-logisticregression-errors-matrix-is-not-positive-definite-cholesky-factorization-failed/28981/13)

---

<div class="post-metadata">

**Author:** ![bertschi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bertschi/32/33462_2.png) [@bertschi](https://discourse.julialang.org/u/bertschi)\
**Post date:** [May 20, 2023, 3:47pm UTC](https://discourse.julialang.org/t/implementation-of-graphical-lasso-with-cholesky-factor/99016/3 "2023-05-20T15:47:02Z")

</div>

In the context of Gaussian processes, the usual fix is to add a small perturbation to the diagonal, i.e., in order to move the eigenvalues away from zero:

```julia
cholesky(Symmetric(A) + 1e-8 * I)

```

Another option is to use

```julia
cholesky(A; check = false)

```

but this might lead to other issues later on (never tried that before).