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
Iβm pretty sure itβs all from SuiteSparse/UMFPack suitesparse : a suite of sparse matrix software 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.
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?