Extracting Lower Triangular Factor of Sparse Cholesky

General Usage
1

This question pertains to issue #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 only works on diagonal matricies.

# 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> 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> 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 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)
2

Pivoting?

3

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> 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
1 Like
4

This is still not working in 1.1.

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
5

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:

Matrix(sparse(C.L))

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

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.