# Why re-use factorization(A) rather than inv(A)

**URL:** <https://discourse.julialang.org/t/why-re-use-factorization-a-rather-than-inv-a/113234>\
**Category:** Numerics\
**Tags:** performance, linearalgebra\
**Created:** [April 19, 2024, 7:16am UTC](https://discourse.julialang.org/t/why-re-use-factorization-a-rather-than-inv-a/113234 "2024-04-19T07:16:56Z")\
**Posts on this page:** 1\
**Showing post:** 5

<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:** [April 19, 2024, 12:25pm UTC](https://discourse.julialang.org/t/why-re-use-factorization-a-rather-than-inv-a/113234/5 "2024-04-19T12:25:00Z")

</div>

> [@photor](#):
>
> Note that the matrix inversion needs to be done only once because the matrix A never changes. In the following you only need to do matrix-vector multiplication, which has a better prefactor than doing left division every time, as you said.

The constant factor of matrix–vector is only slightly better than left-division by a factorization, so (as I explained above) you have to be solving a _lot_ of right-hand sides (proportional to the size of the matrix) before computing the matrix inversion is worth it.

For example, with a 1000 \times 1000 matrix:

```julia
julia> using BenchmarkTools, LinearAlgebra

julia> LinearAlgebra.BLAS.set_num_threads(1) # single-threaded benchmarking is simpler

julia> A = rand(1000,1000); b = rand(1000); x = similar(b);

julia> F = @btime lu($A);
  18.455 ms (6 allocations: 7.64 MiB)

julia> A⁻¹ = @btime inv($A);
  54.006 ms (10 allocations: 8.13 MiB)

julia> @btime ldiv!($x, $F, $b);
  203.052 μs (0 allocations: 0 bytes)

julia> @btime mul!($x, $A⁻¹, $b);
  159.604 μs (0 allocations: 0 bytes)

```

So, on my machine, you’d have to be solving around **800 right-hand sides (one by one)** for the cost of a 1000x1000 matrix inversion to break even.

If you are solving lots of right hand sides at once, i.e. you are computing X = A^{-1} B where B is a matrix, then you need a lot more to hit the break-even point. On my machine the break-even comes at about **2700 simultaneous right-hand sides** for a 1000x1000 matrix:

```julia
julia> n = 2700; B = rand(1000,n); X = similar(B);

julia> @btime ldiv!($B, $F, $B);
  131.404 ms (4 allocations: 20.60 MiB)

julia> @btime mul!($X, $A⁻¹, $B);
  95.540 ms (0 allocations: 0 bytes)

```

(so that the performance loss of `ldiv!` vs `mul!` cancels the performance loss of `inv` vs `lu`. Together, they both take about 150ms in total.)

PS. Obviously, if the matrix is tiny, e.g. 2x2 or 3x3, then the constant factors are completely different, especially if you are using StaticArrays.jl.

---

_[View the full topic](https://discourse.julialang.org/t/why-re-use-factorization-a-rather-than-inv-a/113234)._