# How to pre-allocate matrices to avoid constructing new matrices during the for-loop?

**URL:** <https://discourse.julialang.org/t/how-to-pre-allocate-matrices-to-avoid-constructing-new-matrices-during-the-for-loop/7866>\
**Category:** Performance\
**Tags:** memory-allocation\
**Created:** [December 20, 2017, 12:09am UTC](https://discourse.julialang.org/t/how-to-pre-allocate-matrices-to-avoid-constructing-new-matrices-during-the-for-loop/7866 "2017-12-20T00:09:54Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![hongyan](https://avatars.discourse-cdn.com/v4/letter/h/4da419/32.png) [@hongyan](https://discourse.julialang.org/u/hongyan)\
**Post date:** [December 20, 2017, 12:09am UTC](https://discourse.julialang.org/t/how-to-pre-allocate-matrices-to-avoid-constructing-new-matrices-during-the-for-loop/7866/1 "2017-12-20T00:09:54Z")

</div>

Hi guys, I have some codes as follows:

```julia
function sim(M,R,S,q,numIters,numGrids)
    for i = 1:numIters
        tmp = R + reshape(q' * S, numGrids, numGrids);
        q = (M - tmp) \ ((M + tmp) * q);
    end
    return q;
end

n = 100;
numIters = 10000;
M = rand(n,n);
R = rand(n,n);
S = rand(n,n*n);
q = rand(n);
@time sim(M,R,S,q,numIters,n);

```

If you run the codes, it will show **23.508989 seconds (220.00 k allocations: 3.765 GB, 1.33% gc time)**

Clearly there is a lot of memory allocations. The main reasons are that I need to construct the intermediate matrix **tmp** in each for loop. I know that I should pre-allocate the matrix to avoid this issue, but how? I’m very new in Julia and really appreciate any comments and suggestions!

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [December 20, 2017, 7:49am UTC](https://discourse.julialang.org/t/how-to-pre-allocate-matrices-to-avoid-constructing-new-matrices-during-the-for-loop/7866/2 "2017-12-20T07:49:35Z")

</div>

There is an issue with your code. Do you need to initialise tmp?

Then you have to look for inplace `reshape` and `\`. A starter is (to me this shaves of 1s out of 6):

```julia
function sim!(M,R,S,q,numIters,numGrids,tmp)
    tmp .= tmp .* 0
    for i = 1:numIters
        tmp .= R .+ reshape(q' * S, numGrids, numGrids);
        q .= (M .- tmp) \ ((M .+ tmp) * q);
    end
end

n = 100;
numIters = 10000;
M = rand(n,n);
R = rand(n,n);
S = rand(n,n*n);
q = rand(n);

tmp = rand(n,n);
@time sim!(M,R,S,q,numIters,n,tmp);

```

---

<div class="post-metadata">

**Author:** ![Evizero](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/evizero/32/10118_2.png) [@Evizero](https://discourse.julialang.org/u/Evizero)\
**Post date:** [December 20, 2017, 9:25am UTC](https://discourse.julialang.org/t/how-to-pre-allocate-matrices-to-avoid-constructing-new-matrices-during-the-for-loop/7866/4 "2017-12-20T09:25:10Z")

</div>

> [@rveltz](#):
>
> Then you have to look for inplace reshape

reshape shouldn’t be an allocation bottleneck. The array it returns shares the same underlying data. If it can’t it returns a `ReshapedArray` which is just a view.

```julia
julia> @btime reshape(1:100000, 2, :);
  25.831 ns (0 allocations: 0 bytes)

julia> @btime reshape($(collect(1:100000)), 2, :);
  37.775 ns (2 allocations: 96 bytes)

```

---

<div class="post-metadata">

**Author:** ![Christopher\_Fisher](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/christopher_fisher/32/26132_2.png) [@Christopher\_Fisher](https://discourse.julialang.org/u/Christopher_Fisher)\
**Post date:** [December 20, 2017, 9:43am UTC](https://discourse.julialang.org/t/how-to-pre-allocate-matrices-to-avoid-constructing-new-matrices-during-the-for-loop/7866/5 "2017-12-20T09:43:44Z")

</div>

You can also pre-allocate temp before the for loop with `temp = fill(0.0,size(M))`, unless you need it outside the function for some reason.

---

<div class="post-metadata">

**Author:** ![mohamed82008](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mohamed82008/32/18171_2.png) [@mohamed82008](https://discourse.julialang.org/u/mohamed82008)\
**Post date:** [December 20, 2017, 9:50am UTC](https://discourse.julialang.org/t/how-to-pre-allocate-matrices-to-avoid-constructing-new-matrices-during-the-for-loop/7866/6 "2017-12-20T09:50:41Z")

</div>

The following code cuts down allocations:

```julia
function sim2(M,R,S,q,numIters,numGrids)
     tmp1 = zeros(size(q,2), size(S,2))
     tmp2 = zeros(numGrids, numGrids)
     tmp3 = zeros(M)
     tmp4 = zeros(M)
     tmp5 = zeros(size(M,1), size(q,2))
     for i = 1:numIters
         At_mul_B!(tmp1, q, S)
         tmp2 .= R .+ reshape(tmp1, numGrids, numGrids)
         tmp3 .= M .- tmp2
         tmp4 .= M .+ tmp2
         A_mul_B!(tmp5, tmp4, q)
         q = tmp3 \ tmp5
     end
     return q
end
#258.348244 seconds (100.01 k allocations: 782.930 MiB, 0.03% gc time)

```

The bottleneck is the system solve `q = tmp3 \ tmp5` because it solves the system and allocates for the output every iteration. If you pre-factorize tmp3, you can use the inplace version `A_ldiv_B!`. The following code is significantly faster but it re-arranges the lines. It is obviously a cheat because it only factorizes once, but if you can cheat like this in your actual code then go for it!

```julia
function sim3(M,R,S,q,numIters,numGrids)
   tmp1 = zeros(size(q,2), size(S,2))
   tmp2 = zeros(numGrids, numGrids)
   tmp3 = zeros(M)
   tmp4 = zeros(M)
   tmp5 = zeros(size(M,1), size(q,2))

   At_mul_B!(tmp1, q, S)
   tmp2 .= R .+ reshape(tmp1, numGrids, numGrids)
   tmp3 .= M .- tmp2
   tmp4 .= M .+ tmp2
   fact = lufact(tmp4)
   for i = 1:numIters-1
       A_ldiv_B!(q, fact, tmp5)
       At_mul_B!(tmp1, q, S)
       tmp2 .= R .+ reshape(tmp1, numGrids, numGrids)
       tmp3 .= M .- tmp2
       tmp4 .= M .+ tmp2
       A_mul_B!(tmp5, tmp4, q)
   end
   A_ldiv_B!(q, fact, tmp5)

   return q
end
# 10.253734 seconds (30.02 k allocations: 1.452 MiB)

```

Try these on your machine, and check the timings.

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [December 20, 2017, 10:11am UTC](https://discourse.julialang.org/t/how-to-pre-allocate-matrices-to-avoid-constructing-new-matrices-during-the-for-loop/7866/7 "2017-12-20T10:11:40Z")

</div>

You can use an in-place version of an iterative solver when solving the linear system, e.g.,

```julia
using IterativeSolvers

function sim2(M,R,S,q,numIters,numGrids)
     tmp1 = zeros(size(q,2), size(S,2))
     tmp2 = zeros(numGrids, numGrids)
     tmp3 = zeros(M)
     tmp4 = zeros(M)
     tmp5 = zeros(size(M,1), size(q,2))
     for i = 1:numIters
         At_mul_B!(tmp1, q, S)
         tmp2 .= R .+ reshape(tmp1, numGrids, numGrids)
         tmp3 .= M .- tmp2
         tmp4 .= M .+ tmp2
         A_mul_B!(tmp5, tmp4, q)
         cg!(q, tmp3, tmp5) # Solve for q in-place using iterative solver ConjugateGradients
     end
     return q
end

```

If `q` is not expected to change much between iterations, the last `q` will be a very good initial guess for the iterative solver and the solution will be found very fast.

---

<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:** [December 20, 2017, 11:17am UTC](https://discourse.julialang.org/t/how-to-pre-allocate-matrices-to-avoid-constructing-new-matrices-during-the-for-loop/7866/8 "2017-12-20T11:17:00Z")

</div>

You can also cut down on the allocations caused by the solver by using lufact!, without cheating. This is not entirely allocation-free but uses much fewer. Although I agree with @baggepinnen that iterative solvers are almost surely better here, depending on numerical values (moderately small updates).

```julia
function sim_lu(M,R,S,q,numIters,numGrids)
     tmp1 = zeros(size(q,2), size(S,2))
     tmp2 = zeros(numGrids, numGrids)
     tmp3 = zeros(M)
     tmp4 = zeros(M)
     tmp5 = zeros(size(M,1), size(q,2))
     for i = 1:numIters
         At_mul_B!(tmp1, q, S)
         tmp2 .= R .+ reshape(tmp1, numGrids, numGrids)
         tmp3 .= M .- tmp2
         tmp4 .= M .+ tmp2
         A_mul_B!(tmp5, tmp4, q)
         fact = lufact!(tmp3)
         A_ldiv_B!(q,fact,tmp5)
      end
     return q
end

```