# Surprisingly inaccurate sparse Cholesky factorization

**URL:** <https://discourse.julialang.org/t/surprisingly-inaccurate-sparse-cholesky-factorization/96410>\
**Category:** Numerics\
**Tags:** question, linearalgebra\
**Created:** [March 21, 2023, 5:48pm UTC](https://discourse.julialang.org/t/surprisingly-inaccurate-sparse-cholesky-factorization/96410 "2023-03-21T17:48:31Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![jarapo](https://avatars.discourse-cdn.com/v4/letter/j/6bbea6/32.png) [@jarapo](https://discourse.julialang.org/u/jarapo)\
**Post date:** [March 21, 2023, 5:48pm UTC](https://discourse.julialang.org/t/surprisingly-inaccurate-sparse-cholesky-factorization/96410/1 "2023-03-21T17:48:31Z")

</div>

I am surprised by how inaccurate `cholesky()` can be for sparse matrices; so much so it me wonder if I am not using it correctly.

Define a finite different matrix A:

```julia
L=2;
dvec=-2*ones(2^L);restvec=ones(2^L-1);
A =sparse(-SymTridiagonal(dvec,restvec)+diagm(ones(2^(L))));

```

This gives me

```julia
A
4×4 SparseMatrixCSC{Float64, Int64} with 10 stored entries:
  3.0 -1.0 ⋅ ⋅ 
 -1.0 3.0 -1.0 ⋅ 
   ⋅ -1.0 3.0 -1.0
   ⋅ ⋅ -1.0 3.0

```

Now I obtain the Cholesky factor

```julia
C = cholesky(A);
cholL_sparse = sparse(C.L)
4×4 SparseMatrixCSC{Float64, Int64} with 7 stored entries:
  1.73205 ⋅ ⋅ ⋅ 
 -0.57735 1.63299 ⋅ ⋅ 
   ⋅ ⋅ 1.73205 ⋅ 
   ⋅ -0.612372 -0.57735 1.51383

```

, which is wildly wrong from

```julia
cholL_dense = cholesky(Matrix(A)).L
4×4 LowerTriangular{Float64, Matrix{Float64}}:
  1.73205 ⋅ ⋅ ⋅ 
 -0.57735 1.63299 ⋅ ⋅ 
  0.0 -0.612372 1.62019 ⋅ 
  0.0 0.0 -0.617213 1.61835

```

Indeed

```julia
norm(cholL_sparse*cholL_sparse'-A)
2.0

```

whereas

```julia
norm(cholL_dense*cholL_dense'-A)
4.710277376051325e-16

```

That said, we still get

```julia
cholL_sparse*cholL_sparse' ≈ A[C.p, C.p]
true

```

c.f. [Linear Algebra · The Julia Language](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.cholesky), so something seems to be correct, but I find how it is misleading at best…

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [March 21, 2023, 5:58pm UTC](https://discourse.julialang.org/t/surprisingly-inaccurate-sparse-cholesky-factorization/96410/2 "2023-03-21T17:58:31Z")

</div>

> [@jarapo](#):
>
> which is wildly wrong from

It’s not wrong, it’s a different factorization: the sparse Cholesky factorization is pivoted (i.e. for a permuted A) whereas the dense Choleky factorization is not. The reason for this is that sparse Cholesky uses pivoting to reduce fill-in (i.e. to keep the Cholesky factor as sparse as possible), while in the dense case this is irrelevant.

---

<div class="post-metadata">

**Author:** ![jarapo](https://avatars.discourse-cdn.com/v4/letter/j/6bbea6/32.png) [@jarapo](https://discourse.julialang.org/u/jarapo)\
**Post date:** [March 21, 2023, 6:01pm UTC](https://discourse.julialang.org/t/surprisingly-inaccurate-sparse-cholesky-factorization/96410/3 "2023-03-21T18:01:41Z")

</div>

Thank you. I just read through

> **[Sparse Cholesky decomposition in Julia](http://artadia.blogspot.com/2019/10/sparse-cholesky-decomposition-in-julia.html)**
>
> Julia has nice built-in sparse matrices and algorithms. I spent a whole afternoon trying to extract the sparse Cholesky factorization of a s...

which I found helpful.