# Performance gotcha in linear algebra lu()

**URL:** <https://discourse.julialang.org/t/performance-gotcha-in-linear-algebra-lu/18128>\
**Category:** General Usage\
**Tags:** performance, linearalgebra\
**Created:** [November 29, 2018, 5:14am UTC](https://discourse.julialang.org/t/performance-gotcha-in-linear-algebra-lu/18128 "2018-11-29T05:14:59Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![PetrKryslUCSD](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petrkryslucsd/32/215825_2.png) [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)\
**Post date:** [November 29, 2018, 5:14am UTC](https://discourse.julialang.org/t/performance-gotcha-in-linear-algebra-lu/18128/1 "2018-11-29T05:14:59Z")

</div>

Curiosity killed the cat. I know. But I was curious: how would a naïve implementation of the LU factorization fare compared to what is built in and supported by BLAS 1, 2, 3? The result surprised me as I was able to beat the built-in `lu!()`! Not by much, but still. (EDIT: Note well, the beaten version was the generic LU factorization, not the blas-supported one.)

![image](https://global.discourse-cdn.com/julialang/original/3X/9/a/9ab26f1690d0abaa8f186a2b6b9d7bbb7c17d1e9.png)

The complete code is [here](https://gist.github.com/PetrKryslUCSD/203d1d57f009c0cbf704a2d3fbc97ae6). (Note: both versions of the factorization are without pivoting.)

My configuration:

```julia
julia> versioninfo(verbose = true)                      
Julia Version 1.0.1                                     
Commit 0d713926f8 (2018-09-29 19:05 UTC)                
Platform Info:                                          
  OS: Windows (x86_64-w64-mingw32)                      
      Microsoft Windows [Version 10.0.17134.345]        
  CPU: Intel(R) Core(TM) i7-6650U CPU @ 2.20GHz:        
              speed user nice sys idle irq                                   
       #1 2208 MHz 3902015 0 3664234 33590234 406906 ticks                              
       #2 2208 MHz 3563296 0 2253140 35339640 36734 ticks                              
       #3 2208 MHz 4818578 0 2935328 33402171 40687 ticks                              
       #4 2208 MHz 5637328 0 2455843 33062906 31609 ticks                              
                                                                                                                       
  Memory: 15.927024841308594 GB (7806.58984375 MB free)                                                                
  Uptime: 41156.1347257 sec                                                                                            
  Load Avg: 0.0 0.0 0.0                                                                                             
  WORD_SIZE: 64                                                                                                        
  LIBM: libopenlibm                                                                                                    
  LLVM: libLLVM-6.0.0 (ORCJIT, skylake)       

```

I’m sure this is not the definitive investigation of the matter, but you may find it of interest.

---

<div class="post-metadata">

**Author:** ![antoine-levitt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/antoine-levitt/32/4008_2.png) [@antoine-levitt](https://discourse.julialang.org/u/antoine-levitt)\
**Post date:** [November 29, 2018, 10:32am UTC](https://discourse.julialang.org/t/performance-gotcha-in-linear-algebra-lu/18128/2 "2018-11-29T10:32:07Z")

</div>

I can reproduce, I get almost 2x speedup on N=1000. I’m very confused by why BLAS isn’t more efficient.

---

<div class="post-metadata">

**Author:** ![foobar\_lv2](https://avatars.discourse-cdn.com/v4/letter/f/ee59a6/32.png) [@foobar\_lv2](https://discourse.julialang.org/u/foobar_lv2)\
**Post date:** [November 29, 2018, 1:02pm UTC](https://discourse.julialang.org/t/performance-gotcha-in-linear-algebra-lu/18128/3 "2018-11-29T13:02:44Z")

</div>

That’s because `lu(a, Val(false))` does not call into BLAS and instead uses the fallback `generic_lufact!`. Nevertheless, your naive code benchmarks faster than the `LinearAlgebra` generic fallback, so it would be interesting to see why and whether the stdlib version can be improved.

Tip: Don’t benchmark the inplace versions, use `copy`. Otherwise you will perform the next lu on the result of the last lu, with the possibility of hitting slow subnormals in the benchmark loop. The cost of creating a copy is negligible compared to the cost of the lu factorization.

---

<div class="post-metadata">

**Author:** ![juthohaegeman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juthohaegeman/32/8620_2.png) [@juthohaegeman](https://discourse.julialang.org/u/juthohaegeman)\
**Post date:** [November 29, 2018, 1:40pm UTC](https://discourse.julialang.org/t/performance-gotcha-in-linear-algebra-lu/18128/4 "2018-11-29T13:40:45Z")

</div>

> [@foobar\_lv2](#):
>
> Tip: Don’t benchmark the inplace versions, use `copy` .

Or use: `@btime lu!(b) setup=(b=copy($a))`

---

<div class="post-metadata">

**Author:** ![foobar\_lv2](https://avatars.discourse-cdn.com/v4/letter/f/ee59a6/32.png) [@foobar\_lv2](https://discourse.julialang.org/u/foobar_lv2)\
**Post date:** [November 29, 2018, 1:54pm UTC](https://discourse.julialang.org/t/performance-gotcha-in-linear-algebra-lu/18128/5 "2018-11-29T13:54:23Z")

</div>

Nope, `setup` does not reliably work for this.

```julia
julia> a=Ref(3.0);
julia> function f(a)
       a[] = -sqrt(a[])
       nothing
       end
julia> @btime f(a) setup=(a[]=1.0)
ERROR: DomainError with -1.0:

```

Need to do

```julia
julia> @btime f(a) setup=(a[]=1.0) evals=1
  49.000 ns (0 allocations: 0 bytes)
julia> @btime f($a) setup=(a[]=1.0) evals=1
  30.000 ns (0 allocations: 0 bytes)

```

But if you use `@benchmarkable` to set `setup` and `evals`, and then `tune!`, like `BenchBasemarks.jl`, then the evals is overwritten:

```julia
julia> bm = @benchmarkable f($a) setup=(a[]=1.0) evals=1
julia> @benchmark bm
BenchmarkTools.Trial:
[...]
julia> tune!(bm)
ERROR: DomainError with -1.0:

```

That’s a bug / bad design decision imo, and an issue is already open on basebenchmarks and benchmarktools.

---

<div class="post-metadata">

**Author:** ![PetrKryslUCSD](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petrkryslucsd/32/215825_2.png) [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)\
**Post date:** [November 29, 2018, 4:03pm UTC](https://discourse.julialang.org/t/performance-gotcha-in-linear-algebra-lu/18128/6 "2018-11-29T16:03:08Z")

</div>

> [@juthohaegeman](#):
>
> @btime lu!(b) setup=(b=copy($a))

With the suggestion of @juthohaegeman (concerning the setup of the benchmark test), and in order to test the blas-supported solver (running `@btime lu!(b, Val(true)) setup=(b=copy($a))`, i. e. with pivoting enabled), I get the following graph:

 ![image](https://global.discourse-cdn.com/julialang/original/3X/0/b/0b7dde7e6fd20e3ddb4901f919944e37fc5b43eb.png)

I do have some explanations now.

1. The built-in `generic_lufact!` leaves some performance on the table in the [loop](https://github.com/JuliaLang/julia/blob/d789231e9985537686052db9b2314c0d51656308/stdlib/LinearAlgebra/src/lu.jl#L126)  
 ![image](https://global.discourse-cdn.com/julialang/original/3X/6/5/65900343dd94dec22e8eb14400fab8518d3950c7.png)

The performance of the naïve implementation is matched (for larger number of equations) if this loop is replaced as

```julia
            # Update the rest
            for j = k+1:n
                nAkj = -A[k,j]
                for i = k+1:m
                    A[i,j] += A[i,k]*nAkj
                end
            end

```

1. When pivoting is turned off, `generic_lufact!` still allocates the pivot-position vector, and goes through the motions of working with it. This costs some performance for small matrix sizes. I suppose that is the price of doing business with the `LU` factorization object which expects the pivot vector during construction.

**I guess there is a lesson here: For larger matrices, even if one is sure that no pivoting is necessary, it should not be turned off, because the call is routed to `generic_lufact!` with an attendant order-of-magnitude slow down. For small matrices it may pay to avoid the creation of the factorization object and use the naïve factorization.**

---

<div class="post-metadata">

**Author:** ![antoine-levitt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/antoine-levitt/32/4008_2.png) [@antoine-levitt](https://discourse.julialang.org/u/antoine-levitt)\
**Post date:** [November 29, 2018, 5:55pm UTC](https://discourse.julialang.org/t/performance-gotcha-in-linear-algebra-lu/18128/7 "2018-11-29T17:55:20Z")

</div>

It should probably be noted in the docstring that the routine without pivoting is not BLAS, as this is a non trivial performance gotcha. OpenBLAS is not very fast for small matrices, it would be interesting to see what MKL does here.

---

<div class="post-metadata">

**Author:** ![foobar\_lv2](https://avatars.discourse-cdn.com/v4/letter/f/ee59a6/32.png) [@foobar\_lv2](https://discourse.julialang.org/u/foobar_lv2)\
**Post date:** [November 29, 2018, 9:02pm UTC](https://discourse.julialang.org/t/performance-gotcha-in-linear-algebra-lu/18128/8 "2018-11-29T21:02:09Z")

</div>

The other kind of gotcha in loops is also really important: We know that the write to `A[i,j]`does not alias the read of `A[k,j]`, but the compiler doesn’t know that. At least that’s what I think is going on.

Do you plan to make a PR?

---

<div class="post-metadata">

**Author:** ![PetrKryslUCSD](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petrkryslucsd/32/215825_2.png) [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)\
**Post date:** [November 30, 2018, 12:13am UTC](https://discourse.julialang.org/t/performance-gotcha-in-linear-algebra-lu/18128/9 "2018-11-30T00:13:06Z")

</div>

I have to admit I don’t know yet what to propose. One thing is the documentation, but there’s also the possibility of producing some code that would allow for additional optimization of the choice of the method based on the size of the problem.

---

<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 30, 2018, 4:49am UTC](https://discourse.julialang.org/t/performance-gotcha-in-linear-algebra-lu/18128/10 "2018-11-30T04:49:03Z")

</div>

The `pivot=false` option is basically only useful for pedagogy (comparing to hand calculations or toy code for tiny matrices), and should never be used in real problems because it can make the calculation numerically unstable. Improving its performance would serve no practical purpose.

I agree that the documentation should be clearer on the meaning of the `pivot` argument.

---

<div class="post-metadata">

**Author:** ![PetrKryslUCSD](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petrkryslucsd/32/215825_2.png) [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)\
**Post date:** [November 30, 2018, 4:30pm UTC](https://discourse.julialang.org/t/performance-gotcha-in-linear-algebra-lu/18128/11 "2018-11-30T16:30:59Z")

</div>

> [@stevengj](#):
>
> The `pivot=false` option is basically only useful for pedagogy (comparing to hand calculations or toy code for tiny matrices), and should never be used in real problems because it can make the calculation numerically unstable. Improving its performance would serve no practical purpose.

That is a fair point. But we shouldn’t miss the important part of the story: the performance of the BLAS-supported LU factorization is erratic for smaller matrix sizes, and in fact it is outperformed by the generic Julia version (both WITH PIVOTING).

With just a small change to the code (the “better” implementation of the generic LU factorization), the performance can be improved further.

![image](https://global.discourse-cdn.com/julialang/original/3X/b/7/b7eda50c4a101ca766f96b8d5f3a8633cbea2756.png)

The break-even point between `lu!` and `generic_lufact!` is around 200 equations. NB: For 10 equations or less, the static-array solver implementation may provide further improvements in speed with the generic version of the factorization. (I don’t know if the static array can be passed to the BLAS-supported solver.)

I think it might be of interest to compare with the MKL solver. If anyone has access to it, would you please run the code in [Testing LU · GitHub](https://gist.github.com/PetrKryslUCSD/758d19f03ee2f643bf1f76fe9f5c40ac) and post the generated graph?

EDIT: Additional results for complex matrices. Same computer as above.  
 ![image](https://global.discourse-cdn.com/julialang/original/3X/4/9/49517f2aad7e6a19a5383131b467a8a44450b184.png)

---

<div class="post-metadata">

**Author:** ![PetrKryslUCSD](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petrkryslucsd/32/215825_2.png) [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)\
**Post date:** [November 30, 2018, 4:49pm UTC](https://discourse.julialang.org/t/performance-gotcha-in-linear-algebra-lu/18128/12 "2018-11-30T16:49:42Z")

</div>

Do you have access to MKL-Julia? Could you run the test [https://gist.github.com/PetrKryslUCSD/758d19f03ee2f643bf1f76fe9f5c40ac](https://gist.github.com/PetrKryslUCSD/758d19f03ee2f643bf1f76fe9f5c40ac)?

---

<div class="post-metadata">

**Author:** ![PetrKryslUCSD](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petrkryslucsd/32/215825_2.png) [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)\
**Post date:** [November 30, 2018, 5:54pm UTC](https://discourse.julialang.org/t/performance-gotcha-in-linear-algebra-lu/18128/13 "2018-11-30T17:54:40Z")

</div>

Result from a Linux machine:  
 ![image](https://global.discourse-cdn.com/julialang/original/3X/4/b/4b2a3ec942e763582f5938d13bd96c58c70ee819.png)

```julia
julia> versioninfo(verbose = true)
Julia Version 1.1.0-DEV.509
Commit 1db6047018 (2018-10-21 02:53 UTC)
Platform Info:
  OS: Linux (x86_64-pc-linux-gnu)
      "CentOS release 6.9 (Final)"
  uname: Linux 2.6.32-358.el6.x86_64 #1 SMP Fri Feb 22 00:31:26 UTC 2013 x86_64 x86_64
  CPU: AMD Opteron(tm) Processor 6380 :
                 speed user nice sys idle irq
       #1-64 1400 MHz 4259689828 s 22696 s 423913693 s 79618358095 s 42188 s

  Memory: 252.2855682373047 GB (198149.6875 MB free)
  Uptime: 1.3177519e7 sec
  Load Avg: 1.37744140625 1.42236328125 1.591796875
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-6.0.1 (ORCJIT, bdver1)

```

The “better” generic LU is now slower than the original!? No idea why.

---

<div class="post-metadata">

**Author:** ![antoine-levitt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/antoine-levitt/32/4008_2.png) [@antoine-levitt](https://discourse.julialang.org/u/antoine-levitt)\
**Post date:** [November 30, 2018, 6:45pm UTC](https://discourse.julialang.org/t/performance-gotcha-in-linear-algebra-lu/18128/14 "2018-11-30T18:45:48Z")

</div>

Unfortunately I don’t at the moment.

Also be careful that benchmarking on small matrices is tricky (cache and branch predictions effects), although that probably doesn’t explain why BLAS is so slow.

---

<div class="post-metadata">

**Author:** ![dpo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpo/32/3335_2.png) [@dpo](https://discourse.julialang.org/u/dpo)\
**Post date:** [November 30, 2018, 6:56pm UTC](https://discourse.julialang.org/t/performance-gotcha-in-linear-algebra-lu/18128/15 "2018-11-30T18:56:58Z")

</div>

I’m curious whether you see similar improvements for the Cholesky decomposition, where pivoting isn’t needed.

---

<div class="post-metadata">

**Author:** ![PetrKryslUCSD](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petrkryslucsd/32/215825_2.png) [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)\
**Post date:** [November 30, 2018, 6:58pm UTC](https://discourse.julialang.org/t/performance-gotcha-in-linear-algebra-lu/18128/16 "2018-11-30T18:58:43Z")

</div>

Correct, that might be an interesting experiment as well.

---

<div class="post-metadata">

**Author:** ![Juan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juan/32/7657_2.png) [@Juan](https://discourse.julialang.org/u/Juan)\
**Post date:** [November 30, 2018, 10:32pm UTC](https://discourse.julialang.org/t/performance-gotcha-in-linear-algebra-lu/18128/17 "2018-11-30T22:32:59Z")

</div>

If I had Julia with MKL built… Will all packages automatically benefit from it or just the packages that were explicitly designed to use it?

---

<div class="post-metadata">

**Author:** ![kristoffer.carlsson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kristoffer.carlsson/32/22_2.png) [@kristoffer.carlsson](https://discourse.julialang.org/u/kristoffer.carlsson)\
**Post date:** [November 30, 2018, 10:36pm UTC](https://discourse.julialang.org/t/performance-gotcha-in-linear-algebra-lu/18128/18 "2018-11-30T22:36:31Z")

</div>

Any function that ends up calling BLAS would use MKL instead.

---

<div class="post-metadata">

**Author:** ![Ajaychat3](https://avatars.discourse-cdn.com/v4/letter/a/ecd19e/32.png) [@Ajaychat3](https://discourse.julialang.org/u/Ajaychat3)\
**Post date:** [December 1, 2018, 6:17am UTC](https://discourse.julialang.org/t/performance-gotcha-in-linear-algebra-lu/18128/19 "2018-12-01T06:17:51Z")

</div>

If I install Mkl without compiling Julia with it, would I be able to use mkl for compilation of my code? If so how?

---

<div class="post-metadata">

**Author:** ![tkoolen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tkoolen/32/1603_2.png) [@tkoolen](https://discourse.julialang.org/u/tkoolen)\
**Post date:** [December 1, 2018, 7:37am UTC](https://discourse.julialang.org/t/performance-gotcha-in-linear-algebra-lu/18128/20 "2018-12-01T07:37:46Z")

</div>

I suppose you can `ccall` specific functions from the dynamic libs, but why would you want to?

[Next page](https://discourse.julialang.org/t/performance-gotcha-in-linear-algebra-lu/18128.md?page=2)