# How to optimise a linear solve in a hot loop?

**URL:** <https://discourse.julialang.org/t/how-to-optimise-a-linear-solve-in-a-hot-loop/120835>\
**Category:** Performance\
**Tags:** linearalgebra, linearsolve\
**Created:** [October 3, 2024, 12:55am UTC](https://discourse.julialang.org/t/how-to-optimise-a-linear-solve-in-a-hot-loop/120835 "2024-10-03T00:55:56Z")\
**Posts on this page:** 1\
**Showing post:** 6

<div class="post-metadata">

**Author:** ![RoyiAvital](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/royiavital/32/571_2.png) [@RoyiAvital](https://discourse.julialang.org/u/RoyiAvital)\
**Post date:** [October 3, 2024, 6:15am UTC](https://discourse.julialang.org/t/how-to-optimise-a-linear-solve-in-a-hot-loop/120835/6 "2024-10-03T06:15:52Z")

</div>

Unless you have a connection between each update, as @Oscar_Smith suggested, the only thing you can optimize is the actual solution.

1. Use the fastest solver. In case you are on `x64` and the arrays are dense, try using [`MKL.jl`](https://github.com/JuliaLinearAlgebra/MKL.jl). On Apple Silicon [`AppleAccelerate.jl`](https://github.com/JuliaLinearAlgebra/AppleAccelerate.jl) can do some wonders as well.
2. Choose manually the optimized decomposition for your matrix by doing the decomposition steps manually. For inspiration, have a look at [MATLAB’s `mldivide()`](https://www.mathworks.com/help/matlab/ref/double.mldivide.html).
3. Reduce allocations of the decomposition and solution by using [`FastLapackInterface.jl`](https://github.com/DynareJulia/FastLapackInterface.jl).

It won’t reduce much compared to the actual factorization and solving, but it is a starting point.

Even if the connection in structure of \boldsymbol{A} between iterations can not be exploited for iterative decomposition, if the solution of the previous iteration is similar to the next one, then you may consider iterative method with warmup.

Another idea, **in case your system is guaranteed to have a good condition number** , is to take advantage of the small number of columns relative to the number of rows and solve the [Normal Equations](https://en.wikipedia.org/wiki/Ordinary_least_squares#Normal_equations). Solving \boldsymbol{C} \boldsymbol{x} = \boldsymbol{d} where \boldsymbol{C} = \boldsymbol{A}^{T} \boldsymbol{A} and \boldsymbol{d} = \boldsymbol{A}^{T} \boldsymbol{b}.  
Now the system is a square system with 100^2 rows.  
You will be able to use Cholesky (Pivoting, no check, See [Cholesky Decomposition of Low Rank Positive Semidefinite Matrix](https://discourse.julialang.org/t/70397)) for decomposition (As \boldsymbol{C} is SPSD to the least) and it will be as fast as you can get. Still use 1, 3 from above (2 is covered by Cholesky).  
**It will square your condition number** , so use it wisely.

---

_[View the full topic](https://discourse.julialang.org/t/how-to-optimise-a-linear-solve-in-a-hot-loop/120835)._