# Why lu() factorization is not reducing the execution time?

**URL:** <https://discourse.julialang.org/t/why-lu-factorization-is-not-reducing-the-execution-time/68664>\
**Category:** Performance\
**Tags:** question, linearalgebra\
**Created:** [September 24, 2021, 1:23am UTC](https://discourse.julialang.org/t/why-lu-factorization-is-not-reducing-the-execution-time/68664 "2021-09-24T01:23:23Z")\
**Posts on this page:** 11\
**Page:** 1

<div class="post-metadata">

**Author:** ![Amro](https://avatars.discourse-cdn.com/v4/letter/a/4491bb/32.png) [@Amro](https://discourse.julialang.org/u/Amro)\
**Post date:** [September 24, 2021, 1:23am UTC](https://discourse.julialang.org/t/why-lu-factorization-is-not-reducing-the-execution-time/68664/1 "2021-09-24T01:23:23Z")

</div>

I am using lu() to reduce the computation time but it turned out not as expected. Any idea why is that?

```julia
julia> A = [1 3 3 2;2 6 9 7; -1 -3 3 4]
3×4 Matrix{Int64}:
  1 3 3 2
  2 6 9 7
 -1 -3 3 4

julia> lu(A,check=false)
Failed factorization of type LU{Float64, Matrix{Float64}}

julia> b = [1,2,3]
3-element Vector{Int64}:
 1
 2
 3

julia> @btime A\b
  3.388 μs (37 allocations: 70.50 KiB)
4-element Vector{Float64}:
 -0.09523809523809525
 -0.28571428571428575
  0.21904761904761905
  0.3142857142857143

julia> @btime U\(L\b[p])
  4.029 μs (39 allocations: 71.08 KiB)
4-element Vector{Float64}:
 -3.333333333333332
  1.999999999999999
 -1.6666666666666656
  1.6666666666666674

```

---

<div class="post-metadata">

**Author:** ![jling](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jling/32/212909_2.png) [@jling](https://discourse.julialang.org/u/jling)\
**Post date:** [September 24, 2021, 1:38am UTC](https://discourse.julialang.org/t/why-lu-factorization-is-not-reducing-the-execution-time/68664/2 "2021-09-24T01:38:25Z")

</div>

first of all you’re benchmarking in global scope with non-const variables, second, the results are not the same?..

---

<div class="post-metadata">

**Author:** ![jling](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jling/32/212909_2.png) [@jling](https://discourse.julialang.org/u/jling)\
**Post date:** [September 24, 2021, 1:43am UTC](https://discourse.julialang.org/t/why-lu-factorization-is-not-reducing-the-execution-time/68664/4 "2021-09-24T01:43:00Z")

</div>

what is U,L,p

---

<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:** [September 24, 2021, 3:20am UTC](https://discourse.julialang.org/t/why-lu-factorization-is-not-reducing-the-execution-time/68664/6 "2021-09-24T03:20:35Z")

</div>

First, the right way to use an LU factorization is to call `\` with the factorization object, not to individually apply the `L` and `U` factors yourself.:

```julia
julia> A = randn(1000,1000); b = randn(1000);

julia> LU = lu(A);

julia> @btime $A \ $b;
  9.904 ms (4 allocations: 7.64 MiB)

julia> @btime $LU \ $b;
  241.441 μs (1 allocation: 7.94 KiB)

```

Second, your matrix is too small:

> [@Amro](#):
>
> `b = [1,2,3]`

For a 3x3 matrix, the overhead of LAPACK etcetera is enormous. [StaticArrays](https://github.com/JuliaArrays/StaticArrays.jl) will be much faster for such small problems if the size is fixed:

```julia
julia> A = randn(3,3); b = randn(3);

julia> LU = lu(A);

julia> sA = SMatrix{3,3}(A); sb = SVector{3}(b);

julia> @btime $A \ $b;
  417.249 ns (3 allocations: 384 bytes)

julia> @btime $LU \ $b;
  167.488 ns (1 allocation: 112 bytes)

julia> @btime $sA \ $sb;
  6.537 ns (0 allocations: 0 bytes)

```

---

<div class="post-metadata">

**Author:** ![stillyslalom](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stillyslalom/32/45687_2.png) [@stillyslalom](https://discourse.julialang.org/u/stillyslalom)\
**Post date:** [September 24, 2021, 2:17pm UTC](https://discourse.julialang.org/t/why-lu-factorization-is-not-reducing-the-execution-time/68664/10 "2021-09-24T14:17:34Z")

</div>

Your matrix is rank-deficient, so the decomposition fails - unless you’re tied to a particular decomposition, you can instead use `factorize(A)`, which dispatches to a polyalgorithm according to the properties of `A` (see this chart, or the [docs](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.factorize)). For example,

```julia
julia> f = factorize(A)
QRPivoted{Float64, Matrix{Float64}}
Q factor:
3×3 LinearAlgebra.QRPackedQ{Float64, Matrix{Float64}}:
 -0.301511 -0.275241 -0.912871
 -0.904534 -0.220193 0.365148
 -0.301511 0.935819 -0.182574
R factor:
3×4 Matrix{Float64}:
 -9.94987 -5.4272 -8.14081 -1.80907
  0.0 -4.95434 1.65145 -1.65145
  0.0 0.0 1.50772e-15 2.36998e-16
permutation:
4-element Vector{Int64}:
 3
 2
 4
 1

julia> f\b
4-element Vector{Float64}:
 -0.09523809523809522
 -0.2857142857142857
  0.2190476190476191
  0.3142857142857143

```

Julia’s speed requires that the compiler can figure out the machine representation (the type) of all the variables participating in a calculation. With non-constant global variables, you can redefine their types in the middle of a calculation, or between function calls, which makes it impossible for the compiler to optimize functions that depend on such variables. There’s little point to benchmarking a function if it’s hobbled by global variables - we know it’ll be slow.

---

<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:** [September 24, 2021, 3:15pm UTC](https://discourse.julialang.org/t/why-lu-factorization-is-not-reducing-the-execution-time/68664/12 "2021-09-24T15:15:12Z")

</div>

> [@Amro](#):
>
> Can u tell me what should I do in my MWE to over come the error?

The error message (“matrix is not square”) tells you exactly what is wrong — you can’t do an LU factorization of a non-square (3x4) matrix.

(In real applications, usually you know what kind of problem you are solving in advance, so you don’t need to resort to multiple factorization branches at runtime.)

---

<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:** [September 24, 2021, 4:28pm UTC](https://discourse.julialang.org/t/why-lu-factorization-is-not-reducing-the-execution-time/68664/14 "2021-09-24T16:28:46Z")

</div>

> [@Amro](#):
>
> What does it mean when we type `A \ b` as in below. it is not `inv(A)*b` , right?

See the [`\` documentation](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#Base.:%5C-Tuple%7BAbstractMatrix%7BT%7D%20where%20T,%20AbstractVecOrMat%7BT%7D%20where%20T%7D) (type `?\` in the REPL):

> For rectangular `A` the result is the minimum-norm least squares solution computed by a pivoted QR factorization of `A` and a rank estimate of `A` based on the R factor.

In your case, you have a 3x4 matrix `A`, so Ax=b is an _under_-determined system of equations with infinitely many solutions, and as explained above it returns the solution with the smallest possible `norm`.

`inv(A)` does not exist for a non-square matrix. However, `A \ b` _is_ equivalent (numerically, not computationally) to `pinv(A) * b`, where `pinv` is the [pseudo-inverse](https://en.wikipedia.org/wiki/Moore%E2%80%93Penrose_inverse). Computationally, one rarely needs to compute a matrix inverse (or pseudo-inverse) explicitly.

> [@Amro](#):
>
> what is the function in julia to check if matrix is square rather than typing `if size(A,1) == size(A,2)` ?

Well, you could type `==(size(A)...)` if you want something more compact. There is also [`LinearAlgebra.checksquare`](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.checksquare), but that throws an error if the matrix is non-square (it’s intended for use in functions that _require_ a square matrix.

What is your application that you don’t know ahead of time if you want an exact solution or a least-square solution (whether your problem is over/underdetermined or square/non-singular)?

You can do `qr(A, val(true)) \ b` if you know that you _always_ want to use the QR minimum-norm/least-square solution (which gives the exact solution if `A` is square and invertible).

---

<div class="post-metadata">

**Author:** ![apo383](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/apo383/32/11272_2.png) [@apo383](https://discourse.julialang.org/u/apo383)\
**Post date:** [September 24, 2021, 5:10pm UTC](https://discourse.julialang.org/t/why-lu-factorization-is-not-reducing-the-execution-time/68664/16 "2021-09-24T17:10:59Z")

</div>

Please take a look at the [documentation](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#Base.:%5C-Tuple%7BAbstractMatrix%7BT%7D%20where%20T,%20AbstractVecOrMat%7BT%7D%20where%20T%7D), which should answer most of your questions.

It is generally faster and numerically better to do A\B rather than inv(A)\*B. More info is available in any modern linear algebra class or book that talks about computations. If the book uses Matlab then they will explicitly use \. Edit: fix backslash

---

<div class="post-metadata">

**Author:** ![Amro](https://avatars.discourse-cdn.com/v4/letter/a/4491bb/32.png) [@Amro](https://discourse.julialang.org/u/Amro)\
**Post date:** [September 24, 2021, 5:11pm UTC](https://discourse.julialang.org/t/why-lu-factorization-is-not-reducing-the-execution-time/68664/17 "2021-09-24T17:11:58Z")

</div>

Thank you for your guide.

---

<div class="post-metadata">

**Author:** ![jling](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jling/32/212909_2.png) [@jling](https://discourse.julialang.org/u/jling)\
**Post date:** [September 24, 2021, 5:22pm UTC](https://discourse.julialang.org/t/why-lu-factorization-is-not-reducing-the-execution-time/68664/18 "2021-09-24T17:22:01Z")

</div>

> [@Amro](#):
>
> Using sparse matrix is only helpful when the matrix is large

when the matrix is sparse…

---

<div class="post-metadata">

**Author:** ![jling](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jling/32/212909_2.png) [@jling](https://discourse.julialang.org/u/jling)\
**Post date:** [September 24, 2021, 5:24pm UTC](https://discourse.julialang.org/t/why-lu-factorization-is-not-reducing-the-execution-time/68664/19 "2021-09-24T17:24:39Z")

</div>

> [@Amro](#):
>
> LU factorization should always gives faster execution than not factorization

yes because doing A and B is always gonna be slower than just doing A, in this case, `lu()` is needed for `inv()`:

> <https://github.com/JuliaLang/julia/blob/53da724417802166a5fd7fd4ef05bd7c85d4e1c5/stdlib/LinearAlgebra/src/dense.jl#L879>