I have a very large sparse matrix A and need to produce L, U such that A=L*U for the specific type of problem I am solving.
I am aware of LinearAlgebra.lu, where setting the “pivot” flag to false accomplishes this. This function, however, only accepts dense matrices; lu throws an error if I pass it a sparse matrix and also try disabling “pivot”.
The documentation reads:
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]
F further supports the following functions:
\
cond
det
│ Note
│
│ lu(A::SparseMatrixCSC) uses the UMFPACK library that is part of
│ SuiteSparse. As this library only supports sparse matrices with
│ Float64 or ComplexF64 elements, lu converts A into a copy that is
│ of type SparseMatrixCSC{Float64} or SparseMatrixCSC{ComplexF64} as
│ appropriate.
The expression “F.LF.U == (F.Rs . A)[F.p, F.q]” seems quite a bit messier than what I’m looking for.
Of course, I could re-implement LU myself, but I am writing very performance-sensitive code and I don’t think it’s a good idea to try to compete with the built-in linear algebra algorithms.
Another option is to invert Rs, p, and q, but that seems roundabout and inefficient.
What’s the most efficient way to go about this?
Just doing the inversion does not seem too bad.
using LinearAlgebra
using SparseArrays
A = sprand(100, 100, 0.01) + I
F = lu(A)
invperm(p) = (p⁻¹ = zero(p); ((ij) -> p⁻¹[ij[2]] = ij[1]).(enumerate(p)); p⁻¹)
L = inv.(F.Rs) .* F.L[invperm(F.p), :]
U = F.U[:, invperm(F.q)]
L * U ≈ A # => true
Maybe it only working with Float64s has a bigger impact. At least I would like it if Julia was also able to do this with 32-bit sparse matrices.
Edit: Today I used this myself and realized that I was a big dummy. Of course, if you permute the rows and columns of L and U, they will no longer be triagonal.
Realize that unpivoted LU can be numerically unstable in general; are you sure it’s okay in your case? Moreover, in performance-critical sparse-matrix situations the permutation is essential in order to maximize the sparsity of the L and U factors. Or perhaps you happen to know a good permutation for your problem and want to exploit it?
Currently, there is an open issue for this: `lu` support for custom permutation (like in `cholesky`) · Issue #116 · JuliaSparse/SparseArrays.jl · GitHub
Why do you need this? Normally there is a way to work with the PA=LU form for anything you might need to compute…
3 Likes
My use case is that I am implementing algorithm 2.3 in this paper: The Society for Industrial and Applied Mathematics
The LU factorization is an intermediate step in what is more or less a matrix inversion.
If what you said is true, it’s possible that the permutation vector automatically generated by lu is implicitly embedded in this algorithm and I haven’t correctly identified it.
If you just need the LU factorization to solve a system of equations, you can certainly use it with the permutation term. Indeed, you can solve Ax=b with just x = lu(A) \ b and not worry about how it is implemented.
For example, if you look at the algorithm 2.4 in that paper (where you actually use the LU factorizations), it appears that you could implement w = P_F^t U_F^{-1} L_F^{-1} P_F u_F with something like:
w[p_F] = LU_F \ u_F[p_F]
where LU_F is your lu(...) factorization of P_F H_F P_F^t (= L_F U_F in algorithm 2.3) and p_F is the ordering corresponding to their permutation matrix P_F (which you also don’t explicitly construct).
In general, there is a certain amount of translation that has to take place when converting a mathematical algorithm into code. The example I always give my students is that whenever you see x = A^{-1} b, you should read it as “solve Ax=b by the fastest available method”, which typically does not involve computing A^{-1} explicitly.
8 Likes
I wish there was an established/widespread notation for this, eg A \backslash b, so that people would not feel compelled to write down A^{-1} at all when describing an algorithm.
2 Likes
The notation A^{-1} b is part of the full employment act for academic numerical analysts.
8 Likes
This is indeed a much better implementation than what I was thinking of, and this is quite a bit faster than what I had before, thank you!
1 Like