# Extracting Lower Triangular Factor of Sparse Cholesky

**URL:** <https://discourse.julialang.org/t/extracting-lower-triangular-factor-of-sparse-cholesky/23095>\
**Category:** General Usage\
**Tags:** linearalgebra\
**Created:** [April 12, 2019, 8:52pm UTC](https://discourse.julialang.org/t/extracting-lower-triangular-factor-of-sparse-cholesky/23095 "2019-04-12T20:52:21Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![suchitm](https://avatars.discourse-cdn.com/v4/letter/s/7c8e57/32.png) [@suchitm](https://discourse.julialang.org/u/suchitm)\
**Post date:** [April 12, 2019, 8:52pm UTC](https://discourse.julialang.org/t/extracting-lower-triangular-factor-of-sparse-cholesky/23095/1 "2019-04-12T20:52:21Z")

</div>

This question pertains to [issue #26147](https://github.com/JuliaLang/julia/issues/26147). I am having the same issues with the sparse cholesky factorization as the user who posted the last comment. The results for the lower triangular matrix from a sparse cholesky are not the same as the results from its dense counterpart. The solution in [this issue](https://discourse.julialang.org/t/get-cholesky-factors-of-sparse-input-matrix-in-julia-v0-7-1-0/13987) only works on diagonal matricies.

```julia
# setup
using LinearAlgebra
using SparseArrays
A = sparse(Float64[10 1 1 1; 1 10 0 0; 1 0 10 0; 1 0 0 10])
B = Matrix(A)
a_chol = cholesky(A)
b_chol = cholesky(B)
c = randn(4)

```

The lower part of the cholesky of A is different than that of B (dense A).

```julia
julia> a_chol.L |> sparse |> Matrix
4×4 Array{Float64,2}:
 3.16228 0.0 0.0 0.0
 0.0 3.16228 0.0 0.0
 0.0 0.0 3.16228 0.0
 0.316228 0.316228 0.316228 3.11448

julia> b_chol.L
4×4 LowerTriangular{Float64,Array{Float64,2}}:
 3.16228 ⋅ ⋅ ⋅
 0.316228 3.14643 ⋅ ⋅
 0.316228 -0.0317821 3.14627 ⋅
 0.316228 -0.0317821 -0.0321048 3.1461

```

While the lower parts are not equivalent, solving a linear system leads to almost identical results.

```julia
julia> a_chol \ c
4-element Array{Float64,1}:
  0.09736352798735841
 -0.09335725756564012
  0.13720111810643743
  0.0318010099509956  

julia> b_chol \ c
4-element Array{Float64,1}:
  0.09736352798735844
 -0.0933572575656401
  0.13720111810643745
  0.0318010099509956

```

Version Info:

```julia
Julia Version 1.0.0
Commit 5d4eaca0c9 (2018-08-08 20:58 UTC)
Platform Info:
  OS: macOS (x86_64-apple-darwin14.5.0)
  CPU: Intel(R) Core(TM) i5-8259U CPU @ 2.30GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-6.0.0 (ORCJIT, skylake)

```

---

<div class="post-metadata">

**Author:** ![PetrKryslUCSD](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petrkryslucsd/32/215825_2.png) [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)\
**Post date:** [April 12, 2019, 9:33pm UTC](https://discourse.julialang.org/t/extracting-lower-triangular-factor-of-sparse-cholesky/23095/2 "2019-04-12T21:33:47Z")

</div>

Pivoting?

---

<div class="post-metadata">

**Author:** ![suchitm](https://avatars.discourse-cdn.com/v4/letter/s/7c8e57/32.png) [@suchitm](https://discourse.julialang.org/u/suchitm)\
**Post date:** [April 12, 2019, 9:53pm UTC](https://discourse.julialang.org/t/extracting-lower-triangular-factor-of-sparse-cholesky/23095/3 "2019-04-12T21:53:32Z")

</div>

The suggested way to pivot in issue #26147 is `Lp = sparse(C.L)[C.p, :]`. That doesn’t even yield a lower triangular matrix in this case.

```julia

julia> Matrix(sparse(a_chol.L)[a_chol.p, :])
4×4 Array{Float64,2}:
 0.316228 0.316228 0.316228 3.11448
 0.0 0.0 3.16228 0.0    
 0.0 3.16228 0.0 0.0    
 3.16228 0.0 0.0 0.0  

```

---

<div class="post-metadata">

**Author:** ![Andres\_Legarra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/andres_legarra/32/8786_2.png) [@Andres\_Legarra](https://discourse.julialang.org/u/Andres_Legarra)\
**Post date:** [October 4, 2019, 6:56pm UTC](https://discourse.julialang.org/t/extracting-lower-triangular-factor-of-sparse-cholesky/23095/4 "2019-10-04T18:56:22Z")

</div>

This is still not working in 1.1.

```julia
julia> A=[ 10.0 2.0 0.0 -2.0 0.0 -4.0 0.0 0.0;
         2.0 2.0 0.0 -2.0 0.0 0.0 0.0 0.0;
         0.0 0.0 2.0 0.0 0.0 0.0 0.0 0.0;
        -2.0 -2.0 0.0 6.0 -2.0 -2.0 0.0 0.0;
         0.0 0.0 0.0 -2.0 20.0 8.0 -8.0 -8.0;
        -4.0 0.0 0.0 -2.0 8.0 8.0 -4.0 -4.0;
         0.0 0.0 0.0 0.0 -8.0 -4.0 8.0 0.0;
         0.0 0.0 0.0 0.0 -8.0 -4.0 0.0 8.0];

julia> Asp=sparse(A);

julia> Matrix(sparse(C.L)[C.p,:])
8×8 Array{Float64,2}:
  0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.41421
  0.0 0.0 0.0 -2.0 1.41421 -1.41421 1.41421 0.0    
  0.0 0.0 0.0 0.0 1.41421 0.0 0.0 0.0    
  0.0 0.0 -1.0 -1.0 -1.41421 1.41421 0.0 0.0    
  0.0 2.82843 0.0 0.0 0.0 0.0 0.0 0.0    
 -1.41421 -1.41421 8.88178e-16 2.0 0.0 0.0 0.0 0.0    
  2.82843 0.0 0.0 0.0 0.0 0.0 0.0 0.0    
 -2.82843 -2.82843 2.0 0.0 0.0 0.0 0.0 0.0   

```

---

<div class="post-metadata">

**Author:** ![Andres\_Legarra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/andres_legarra/32/8786_2.png) [@Andres\_Legarra](https://discourse.julialang.org/u/Andres_Legarra)\
**Post date:** [October 7, 2019, 7:35am UTC](https://discourse.julialang.org/t/extracting-lower-triangular-factor-of-sparse-cholesky/23095/5 "2019-10-07T07:35:20Z")

</div>

After reading help carefully, I realised that @suchitm and I are confused. The Sparse Cholesky gives a decomposition **LL’** of a _permuted_ matrix **PAP’** , not of **A** itself. The factor **L** :

```julia
Matrix(sparse(C.L)) 

```

is lower triangular but is the decomposition of **PAP’**. To extract the decomposition of **A** we need

```julia
sparse(C.L)[invperm(C.p),:]

```

(which multiplied by its transpose gives **A** ) but this **P’L** which is not triangular because it can’t be - a strictly triangular decomposition would use **P=I** and would be little sparse in most cases. More explanations [here](https://artadia.blogspot.com/2019/10/sparse-cholesky-decomposition-in-julia.html).