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?

Thanks

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?

Thanks

0 Likes

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
```

0 Likes

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]
```

3 Likes

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

0 Likes