Looking for documentation related to sparse matrix factorization


This is surely a dumb question, but I can’t seem to find documentation describing LU factorization for sparse matrices. Could anyone direct me to the docs?




I’m pretty sure it’s all from SuiteSparse/UMFPack http://faculty.cse.tamu.edu/davis/suitesparse.html which uses a graph theoretic algorithm to determine optimal pivoting, which depends on the specific entries to avoid ill-conditioning.

In other words, it’s complicated and best treated as a black-box unless it’s your area of research.

1 Like


Maybe I’m not being clear; I’m looking for the documentation describing what arguments the lu() function takes, what it returns etc when given a sparse input. For example, the documentation here does not seem correct; try:

using LinearAlgebra
using SparseArrays

A = spzeros(2,2)
A[1:2, 1:2] = [1 1; -1 1]
F = lu(A)
print(F.L * F.U == A[F.p, :]) # false


For sparse matrices, output of lu is slightly different. The following is from help mode of Julia REPL (notice the very last line):

  lu(A::SparseMatrixCSC; check = true) -> F::UmfpackLU

  Compute the LU factorization of a sparse matrix A.

  For sparse A with real or complex element type, the return type of F is UmfpackLU{Tv, Ti}, with Tv = Float64 or ComplexF64 respectively
  and Ti is an integer type (Int32 or Int64).

  When check = true, an error is thrown if the decomposition fails. When check = false, responsibility for checking the decomposition's
  validity (via issuccess) lies with the user.

  The individual components of the factorization F can be accessed by indexing:

Component Description                    
––––––––– –––––––––––––––––––––––––––––––
L         L (lower triangular) part of LU
U         U (upper triangular) part of LU
p         right permutation Vector       
q         left permutation Vector        
Rs        Vector of scaling factors      
:         (L,U,p,q,Rs) components        

  The relation between F and A is

  F.L*F.U == (F.Rs .* A)[F.p, F.q]


Thanks! This is exactly what I was looking for. Is there any reason why this isn’t part of the Julia docs?