# Solve a system of linear equations many times without allocating memory

**URL:** <https://discourse.julialang.org/t/solve-a-system-of-linear-equations-many-times-without-allocating-memory/36228>\
**Category:** Numerics\
**Created:** [March 19, 2020, 9:16pm UTC](https://discourse.julialang.org/t/solve-a-system-of-linear-equations-many-times-without-allocating-memory/36228 "2020-03-19T21:16:16Z")\
**Posts on this page:** 19\
**Page:** 1

<div class="post-metadata">

**Author:** ![MartinOtter](https://avatars.discourse-cdn.com/v4/letter/m/87869e/32.png) [@MartinOtter](https://discourse.julialang.org/u/MartinOtter)\
**Post date:** [March 19, 2020, 9:16pm UTC](https://discourse.julialang.org/t/solve-a-system-of-linear-equations-many-times-without-allocating-memory/36228/1 "2020-03-19T21:16:16Z")

</div>

I have the problem that during simulation a linear equation system A\*x=b is solved many times. The desired solution is to allocate memory for all arrays involved in the operation once before simulation starts and during simulation only update the array elements and then compute a new solution. I did not yet find a way to implement this in Julia. Example:

```julia
module Test_lu

using LinearAlgebra

nx = 10
A = rand(nx,nx)
b = rand(nx)

for i=1:5 # Solve liner equation system 5-times
    for j=1:nx
        for k=1:nx
            A[j,k] = A[j,k] + i # New values for A
        end
        b[j] = b[j] + i # New values for b
    end
    @time w = lu!(A)
    @time ldiv!(w,b)
end

end

# Gives output:
  0.000010 seconds (2 allocations: 192 bytes)
  0.000006 seconds
  0.000003 seconds (2 allocations: 192 bytes)
  0.000001 seconds
  0.000003 seconds (2 allocations: 192 bytes)
  0.000002 seconds
  0.000003 seconds (2 allocations: 192 bytes)
  0.000004 seconds
  0.000006 seconds (2 allocations: 192 bytes)
  0.000004 seconds

```

Its clear that lu! always allocates memory because it returns a new vector (here: w), probably the pivot vector. Is there a way to avoid this (so provide a pre-allocated w vector)?  
Or is there another way to solve a linear equation system many times (with changed elements of A,b), without allocating memory in every iteration?

Probably, I could directly use LAPACK functions and then it should be possible to formulate it, but then I loose the possibility to change the type of the matrices (e.g. use Measurements.jl).

---

<div class="post-metadata">

**Author:** ![jlchan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jlchan/32/10958_2.png) [@jlchan](https://discourse.julialang.org/u/jlchan)\
**Post date:** [March 19, 2020, 9:24pm UTC](https://discourse.julialang.org/t/solve-a-system-of-linear-equations-many-times-without-allocating-memory/36228/2 "2020-03-19T21:24:04Z")

</div>

2 allocations doesn’t seem like it would contribute significantly to runtime. Are there specific reasons you wish to avoid all allocations?

---

<div class="post-metadata">

**Author:** ![MartinOtter](https://avatars.discourse-cdn.com/v4/letter/m/87869e/32.png) [@MartinOtter](https://discourse.julialang.org/u/MartinOtter)\
**Post date:** [March 19, 2020, 9:40pm UTC](https://discourse.julialang.org/t/solve-a-system-of-linear-equations-many-times-without-allocating-memory/36228/3 "2020-03-19T21:40:39Z")

</div>

Julia is very nice to use, but the drawback is that it allocates memory at many places. Some time ago, I systematically went through some part of my code called during simulation to get rid of unnecessary memory allocations in order to have a fair comparision with available (non-Julia) tools in my area. At the end, the speed-up (just by removing unnecessary memory allocations) was about 10. I do not ask currently, whether one particular place really counts for simulation speed-up or not, but just want to get rid of memory allocations whenever possible. Note, all this would be uncritical, if Julia would have some O(1) caching mechanism, so that previously freed memory can be very efficiently re-allocated again. But this seems to be not the case.

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [March 19, 2020, 9:49pm UTC](https://discourse.julialang.org/t/solve-a-system-of-linear-equations-many-times-without-allocating-memory/36228/4 "2020-03-19T21:49:50Z")

</div>

Do you need to compute the lu factorization? If not, you could just use `ldiv!(X,b)`

---

<div class="post-metadata">

**Author:** ![MartinOtter](https://avatars.discourse-cdn.com/v4/letter/m/87869e/32.png) [@MartinOtter](https://discourse.julialang.org/u/MartinOtter)\
**Post date:** [March 19, 2020, 10:02pm UTC](https://discourse.julialang.org/t/solve-a-system-of-linear-equations-many-times-without-allocating-memory/36228/5 "2020-03-19T22:02:36Z")

</div>

I do not need the lu factorization, just the solution of the equation system. I just tried your proposal but got errors (I am using Julia 1.3.1):

```julia
nx = 10
A = rand(nx,nx)
b = rand(nx)
ldiv!(A,b) 
# ERROR: MethodError: no method matching ldiv!(::Array{Float64,2}, ::Array{Float64,1})

```

Is this a feature of a newer Julia version?

---

<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:** [March 19, 2020, 10:18pm UTC](https://discourse.julialang.org/t/solve-a-system-of-linear-equations-many-times-without-allocating-memory/36228/6 "2020-03-19T22:18:20Z")

</div>

It is `ldiv!(Y, A, B)` I believe.

---

<div class="post-metadata">

**Author:** ![DNF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dnf/32/10191_2.png) [@DNF](https://discourse.julialang.org/u/DNF)\
**Post date:** [March 19, 2020, 11:02pm UTC](https://discourse.julialang.org/t/solve-a-system-of-linear-equations-many-times-without-allocating-memory/36228/7 "2020-03-19T23:02:45Z")

</div>

You should definitely use BenchmarkTools for this. You will never see zero allocations using the `@time` macro.

The timing estimates are probably also wildly inaccurate. Microbenchmarks are not a good use case for `@time`. Use BenchmarkTools, and remember to interpolate.

---

<div class="post-metadata">

**Author:** ![MartinOtter](https://avatars.discourse-cdn.com/v4/letter/m/87869e/32.png) [@MartinOtter](https://discourse.julialang.org/u/MartinOtter)\
**Post date:** [March 20, 2020, 6:23am UTC](https://discourse.julialang.org/t/solve-a-system-of-linear-equations-many-times-without-allocating-memory/36228/8 "2020-03-20T06:23:01Z")

</div>

> [@PetrKryslUCSD](#):
>
> It is `ldiv!(Y, A, B)` I believe.

I tried:

```julia
nx = 10
A = rand(nx,nx)
b = rand(nx)
x = zeros(nx)
@time ldiv!(x,A,b)
#ERROR: MethodError: no method matching ldiv!(::Array{Float64,1}, ::Array{Float64,2}, ::Array{Float64,1})

```

> [@DNF](#):
>
> You should definitely use BenchmarkTools for this. You will never see zero allocations using the `@time` macro.
> 
> The timing estimates are probably also wildly inaccurate. Microbenchmarks are not a good use case for `@time` . Use BenchmarkTools, and remember to interpolate

lu! returns a struct that contains a reference to the overwritten A and a (new) pivot vector. So, it is clear that lu! always allocates memory.

---

<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:** [March 20, 2020, 6:50am UTC](https://discourse.julialang.org/t/solve-a-system-of-linear-equations-many-times-without-allocating-memory/36228/9 "2020-03-20T06:50:37Z")

</div>

> [@MartinOtter](#):
>
> So, it is clear that lu! always allocates memory.

Not if you disable pivoting. You may be looking for

```julia
Alu = lu!(A, Val(false))
ldiv!(Alu, b)

```

Note, however, that you may worried about something that is unlikely to be relevant: if `A` is large, then the computation itself will dominate timing, if `A` is small, you should using StaticArrays which again do not allocate.

Also, please do read the docstrings (`?lu!`, `?ldiv!`) instead of programming by trial and error. Most of these things are documented.

---

<div class="post-metadata">

**Author:** ![MartinOtter](https://avatars.discourse-cdn.com/v4/letter/m/87869e/32.png) [@MartinOtter](https://discourse.julialang.org/u/MartinOtter)\
**Post date:** [March 20, 2020, 7:48am UTC](https://discourse.julialang.org/t/solve-a-system-of-linear-equations-many-times-without-allocating-memory/36228/10 "2020-03-20T07:48:20Z")

</div>

> [@MartinOtter](#):
>
> > [@MartinOtter](#):
> >
> > So, it is clear that lu! always allocates memory.
> 
> Not if you disable pivoting. You may be looking for
> 
> ```julia
> Alu = lu!(A, Val(false))
> ldiv!(Alu, b)
> 
> ```
> 
> Note, however, that you may worried about something that is unlikely to be relevant: if `A` is large, then the computation itself will dominate timing, if `A` is small, you should using StaticArrays which again do not allocate.
> 
> Also, please do read the docstrings ( `?lu!` , `?ldiv!` ) instead of programming by trial and error. Most of these things are documented.

Disabling pivoting is no alternativ, because a solution of a regular system may fail.

I am contributing to a generic modeling and simulation system (Modia/Modia3D/ModiaMath) and depending on user models a few or many, small or large linear equation systems will appear.

My summary of the remarks is: _In Julia there is currently no generic library function available to solve repeately a (small or large) linear equation system without allocating memory in every iteration. For the time being I will just use lu! and ldiv! and at some time in the future provide a special library function that either uses LAPACK functions directly (for Float64) and otherwise uses lu! and ldiv!._

---

<div class="post-metadata">

**Author:** ![DNF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dnf/32/10191_2.png) [@DNF](https://discourse.julialang.org/u/DNF)\
**Post date:** [March 20, 2020, 8:04am UTC](https://discourse.julialang.org/t/solve-a-system-of-linear-equations-many-times-without-allocating-memory/36228/11 "2020-03-20T08:04:14Z")

</div>

> [@MartinOtter](#):
>
> lu! returns a struct that contains a reference to the overwritten A and a (new) pivot vector. So, it is clear that lu! always allocates memory.

This doesn’t change the fact that you should always use BenchmarkTools for microbenchmarks:

```julia
julia> @time ldiv!(w, b);
  0.000021 seconds (4 allocations: 160 bytes)

julia> @btime ldiv!($w, $b);
  1.530 μs (0 allocations: 0 bytes)

julia> 0.000021/1.53e-6
13.72549019607843

```

Or perhaps better, not sure what `ldiv!` does (edited this a couple of times):

```julia
for i=1:5 # Solve liner equation system 5-times
    for j=1:nx
        for k=1:nx
            A[j,k] = A[j,k] + i # New values for A
        end
        b[j] = b[j] + i # New values for b
    end
    @btime lu!(A_) setup=(A_=copy(A))
    @btime ldiv!(w, b_) setup=(w=lu!(copy(A)); b_=copy(b))
end

1.089 μs (2 allocations: 192 bytes)
356.595 ns (0 allocations: 0 bytes)
1.116 μs (2 allocations: 192 bytes)
356.665 ns (0 allocations: 0 bytes)
1.112 μs (2 allocations: 192 bytes)
356.832 ns (0 allocations: 0 bytes)
1.107 μs (2 allocations: 192 bytes)
356.900 ns (0 allocations: 0 bytes)
1.109 μs (2 allocations: 192 bytes)
356.928 ns (0 allocations: 0 bytes)

```

You have to do be a bit careful with mutating functions (not sure if I did this correctly, but to me it appears to be a 10-60x difference in measurement result.)

---

<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:** [March 20, 2020, 9:26am UTC](https://discourse.julialang.org/t/solve-a-system-of-linear-equations-many-times-without-allocating-memory/36228/12 "2020-03-20T09:26:11Z")

</div>

> [@MartinOtter](#):
>
> in the future provide a special library function

I think it would make sense to have `lu!(Alu, A)` with pivoting. Consider making a PR to LinearAlgebra.

When using Julia, keep in mind that library functions do not cover every conceivable use case. Some are omitted by design, but some are waiting to be implemented. We are talking about a few lines of code anyway.

---

<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:** [March 20, 2020, 3:21pm UTC](https://discourse.julialang.org/t/solve-a-system-of-linear-equations-many-times-without-allocating-memory/36228/13 "2020-03-20T15:21:14Z")

</div>

```julia
a = rand(6, 6); b = rand(6); x = rand(6)
a = lu!(a) 
ldiv!(x, a, b) 

```

---

<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:** [March 20, 2020, 4:58pm UTC](https://discourse.julialang.org/t/solve-a-system-of-linear-equations-many-times-without-allocating-memory/36228/14 "2020-03-20T16:58:41Z")

</div>

> [@Tamas\_Papp](#):
>
> I think it would make sense to have `lu!(Alu, A)` with pivoting. Consider making a PR to LinearAlgebra.

IIRC, the issue is that LAPACK’s factorization routines typically require some additional workspace vectors to be supplied, in addition to storage for the matrix itself, so these have to be allocated in `lu!`. To be completely allocation-free, the caller would need to pass an additional workspace parameter.

(Of course, you could always call the LAPACK routines directly from your Julia code and do whatever low-level hackery you want.)

---

<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:** [March 21, 2020, 4:31pm UTC](https://discourse.julialang.org/t/solve-a-system-of-linear-equations-many-times-without-allocating-memory/36228/15 "2020-03-21T16:31:53Z")

</div>

> [@stevengj](#):
>
> To be completely allocation-free, the caller would need to pass an additional workspace parameter.

I think that should be feasible. Alternatively, something like

```julia
Alu = lu(A)
lu!(Alu, A) # use previously allocated results if types are compatible

```

could also work and could provide a nice general interface for similar methods.

---

<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:** [March 21, 2020, 4:32pm UTC](https://discourse.julialang.org/t/solve-a-system-of-linear-equations-many-times-without-allocating-memory/36228/16 "2020-03-21T16:32:42Z")

</div>

That said, I still think that worrying about allocations is not something one should consider unless they are willing to pay quite a bit in code complexity for very marginal improvements. Here is a function for those who want to experiment:

```julia
using LinearAlgebra, BenchmarkTools

function f(n, m)
    A = rand(n, n)
    b = rand(n)
    z = 0.0
    for _ in 1:m
        z += (A \ b)[1]
        A .+= 1
        b .+= 1
    end
    z
end

@benchmark f(20, 20) # try n = 10, 20, 30, 50

```

where the outer loop was included to make profiling easier. Around 1–2% of the time is spent on GC even for small matrices. LAPACK time dominates everything. This becomes much more prominent for larger matrices.

---

<div class="post-metadata">

**Author:** ![MartinOtter](https://avatars.discourse-cdn.com/v4/letter/m/87869e/32.png) [@MartinOtter](https://discourse.julialang.org/u/MartinOtter)\
**Post date:** [April 7, 2020, 10:48am UTC](https://discourse.julialang.org/t/solve-a-system-of-linear-equations-many-times-without-allocating-memory/36228/17 "2020-04-07T10:48:30Z")

</div>

> [@Tamas\_Papp](#):
>
> > [@stevengj](#):
> >
> > To be completely allocation-free, the caller would need to pass an additional workspace parameter.
> 
> I think that should be feasible. Alternatively, something like
> 
> ```julia
> Alu = lu(A)
> lu!(Alu, A) # use previously allocated results if types are compatible
> 
> ```
> 
> could also work and could provide a nice general interface for similar methods.

Yes, this is a nice general interface (especially, since the caller does not need to know the type of the internally allocated memory).

---

<div class="post-metadata">

**Author:** ![MartinOtter](https://avatars.discourse-cdn.com/v4/letter/m/87869e/32.png) [@MartinOtter](https://discourse.julialang.org/u/MartinOtter)\
**Post date:** [April 7, 2020, 11:08am UTC](https://discourse.julialang.org/t/solve-a-system-of-linear-equations-many-times-without-allocating-memory/36228/18 "2020-04-07T11:08:58Z")

</div>

> [@Tamas\_Papp](#):
>
> That said, I still think that worrying about allocations is not something one should consider unless they are willing to pay quite a bit in code complexity for very marginal improvements. …

When you have many calls, as in simulation or optimization, permanently allocating memory (even small amounts), might be not acceptable:

- Offline simulation: Assume every call allocates 1kbyte and this is done every milli-second, your main memory has 16Gbyte → Simulating more as 4.4 hours (16e6\*0.001/60/60) will fill the complete main memory. Of course, garbage collection will “somehow” take care of it, but this will take time and this is hard to measure.

- Reatime simulation: If hardware is involved, especially for a safety critical operation, it is completely forbidden to allocate memory during operation. Hard realtime simulation might not be possible currently with Julia, but I guess this can be achieved at some point in the future.

---

<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:** [April 7, 2020, 11:28am UTC](https://discourse.julialang.org/t/solve-a-system-of-linear-equations-many-times-without-allocating-memory/36228/19 "2020-04-07T11:28:34Z")

</div>

> [@MartinOtter](#):
>
> Of course, garbage collection will “somehow” take care of it

Garbage collection will take care of it by (surprise) collecting unused memory. Yes, it takes time, but modern GC can be very efficient. Sometimes more efficient than manual memory management.

I am not saying that one should never worry about allocations — they can be a key optimization in hot loops. That said, writing code with nontrivial logic and control flow with absolutely 0 allocations is rarely worth it in Julia, especially if one wants to take advantage of facilities like AD. The effort expended on the last 0.1% of time spent on GC can be usually allocated (😉) much better.

> [@MartinOtter](#):
>
> Hard realtime simulation might not be possible currently with Julia, but I guess this can be achieved at some point in the future.

Hard realtime is something you design a language for from the very beginning. Allocations are only part of the story. It is a very constraining requirement for which specialized tools exist (and consequently it can be very expensive). I don’t think it is planned for Julia in any sense.