# How do I do a plain old "A=LU" decomposition of sparse matrices?

**URL:** <https://discourse.julialang.org/t/how-do-i-do-a-plain-old-a-lu-decomposition-of-sparse-matrices/49619>\
**Category:** Numerics\
**Tags:** linearalgebra\
**Created:** [November 5, 2020, 9:18am UTC](https://discourse.julialang.org/t/how-do-i-do-a-plain-old-a-lu-decomposition-of-sparse-matrices/49619 "2020-11-05T09:18:44Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![casey\_watkins](https://avatars.discourse-cdn.com/v4/letter/c/e95f7d/32.png) [@casey\_watkins](https://discourse.julialang.org/u/casey_watkins)\
**Post date:** [November 5, 2020, 9:18am UTC](https://discourse.julialang.org/t/how-do-i-do-a-plain-old-a-lu-decomposition-of-sparse-matrices/49619/1 "2020-11-05T09:18:44Z")

</div>

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:

```julia
 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.L_F.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?

---

<div class="post-metadata">

**Author:** ![martenlienen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/martenlienen/32/18572_2.png) [@martenlienen](https://discourse.julialang.org/u/martenlienen)\
**Post date:** [November 5, 2020, 9:56am UTC](https://discourse.julialang.org/t/how-do-i-do-a-plain-old-a-lu-decomposition-of-sparse-matrices/49619/2 "2020-11-05T09:56:01Z")

</div>

Just doing the inversion does not seem too bad.

```julia
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 `Float64`s 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.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [November 11, 2020, 11:49pm UTC](https://discourse.julialang.org/t/how-do-i-do-a-plain-old-a-lu-decomposition-of-sparse-matrices/49619/3 "2020-11-11T23:49:25Z")

</div>

> [@casey\_watkins](#):
>
> 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.

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](https://github.com/JuliaLang/julia/issues/36720)

Why do you need this? Normally there is a way to work with the PA=LU form for anything you might need to compute…

---

<div class="post-metadata">

**Author:** ![casey\_watkins](https://avatars.discourse-cdn.com/v4/letter/c/e95f7d/32.png) [@casey\_watkins](https://discourse.julialang.org/u/casey_watkins)\
**Post date:** [November 12, 2020, 8:04am UTC](https://discourse.julialang.org/t/how-do-i-do-a-plain-old-a-lu-decomposition-of-sparse-matrices/49619/4 "2020-11-12T08:04:05Z")

</div>

My use case is that I am implementing algorithm 2.3 in this paper: [The Society for Industrial and Applied Mathematics](https://epubs.siam.org/doi/abs/10.1137/100804644)

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.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [November 12, 2020, 1:26pm UTC](https://discourse.julialang.org/t/how-do-i-do-a-plain-old-a-lu-decomposition-of-sparse-matrices/49619/5 "2020-11-12T13:26:25Z")

</div>

> [@casey\_watkins](#):
>
> The LU factorization is an intermediate step in what is more or less a matrix inversion.

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:

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

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [November 12, 2020, 2:12pm UTC](https://discourse.julialang.org/t/how-do-i-do-a-plain-old-a-lu-decomposition-of-sparse-matrices/49619/6 "2020-11-12T14:12:12Z")

</div>

> [@stevengj](#):
>
> whenever you see x=A−1bx = A^{-1} b , you should read it as “solve Ax=bAx=b by the fastest available method”

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.

---

<div class="post-metadata">

**Author:** ![ctkelley](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ctkelley/32/10684_2.png) [@ctkelley](https://discourse.julialang.org/u/ctkelley)\
**Post date:** [November 12, 2020, 2:27pm UTC](https://discourse.julialang.org/t/how-do-i-do-a-plain-old-a-lu-decomposition-of-sparse-matrices/49619/7 "2020-11-12T14:27:38Z")

</div>

The notation A^{-1} b is part of the full employment act for academic numerical analysts.

---

<div class="post-metadata">

**Author:** ![casey\_watkins](https://avatars.discourse-cdn.com/v4/letter/c/e95f7d/32.png) [@casey\_watkins](https://discourse.julialang.org/u/casey_watkins)\
**Post date:** [November 19, 2020, 1:16am UTC](https://discourse.julialang.org/t/how-do-i-do-a-plain-old-a-lu-decomposition-of-sparse-matrices/49619/8 "2020-11-19T01:16:29Z")

</div>

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!