How to get upper triangular Cholesky factor of sparse matrix?

Specific Domains Numerics
1

Hi all. I read the documentation here Linear Algebra · The Julia Language for what is returned when you call cholesky on a SparseMatrixCSC object but I still can’t figure out how to get back the upper triangular factor of the (non permuted) matrix. The docs say to use LtP = F.UP but whatever is returned by the field UP is not being converted into what I want by calling sparse or Matrix on it. Can someone help me figure out how to do this? MWE below:

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];
A = sparse(A)
C = cholesky(A)
sparse(C.UP) # error
Matrix(C.UP) # error
2 
using SparseArrays, LinearAlgebra

A = sparse([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])

C = cholesky(A)

# Option 1: Upper triangular factor in permuted space
U_permuted = sparse(C.L)' 

# Option 2: Upper triangular factor in original space
perm = C.p
invperm_vec = invperm(perm)
L_original = sparse(C.L)[invperm_vec, invperm_vec]
U_original = L_original'

# Verify the factorizations
P = sparse(1:8, perm, ones(8))

# In permuted space: P*A*P' = L_permuted * U_permuted
println("Error in permuted space: ", norm(P*A*P' - sparse(C.L)*U_permuted))

# In original space: A = L_original * U_original
println("Error in original space: ", norm(A - L_original*U_original))

You need to take this with the huge qualification that I’m not even sure that this was covered when I took my maths as an undergrad in the 60s.

1 Like
3

It’s poorly documented, but C.L is the only factor that supports being directly converted to a (sparse) matrix, so you need to do the permutation and transposition yourself, as @technocrat showed. However, you don’t need to permute the second dimension of L— it’s contracted over anyway, so its order is irrelevant as long as it’s consistent between the two factors. Thus, what you want is PtL = sparse(C.L)[invperm(p), :], and, if desired, UP = PtL'.

julia> using SparseArrays, LinearAlgebra

julia> A = sparse([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> C = cholesky(A);

julia> L = sparse(C.L);

julia> PtL = L[invperm(C.p), :];

julia> PtL * PtL'
8×8 SparseMatrixCSC{Float64, Int64} with 30 stored entries:
 10.0   2.0   ⋅   -2.0    ⋅   -4.0    ⋅     ⋅
  2.0   2.0   ⋅   -2.0    ⋅     ⋅     ⋅     ⋅
   ⋅     ⋅   2.0    ⋅     ⋅     ⋅     ⋅     ⋅
 -2.0  -2.0   ⋅    6.0  -2.0  -2.0    ⋅     ⋅
   ⋅     ⋅    ⋅   -2.0  20.0   8.0  -8.0  -8.0
 -4.0    ⋅    ⋅   -2.0   8.0   8.0  -4.0  -4.0
   ⋅     ⋅    ⋅     ⋅   -8.0  -4.0   8.0    ⋅
   ⋅     ⋅    ⋅     ⋅   -8.0  -4.0    ⋅    8.0
2 Likes
4

Many thanks to @technocrat and @danielwe.