# Allocations when constructing a matrix from columns

**URL:** <https://discourse.julialang.org/t/allocations-when-constructing-a-matrix-from-columns/36795>\
**Category:** Performance\
**Tags:** memory-allocation\
**Created:** [March 31, 2020, 10:59am UTC](https://discourse.julialang.org/t/allocations-when-constructing-a-matrix-from-columns/36795 "2020-03-31T10:59:18Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![mht](https://avatars.discourse-cdn.com/v4/letter/m/ea5d25/32.png) [@mht](https://discourse.julialang.org/u/mht)\
**Post date:** [March 31, 2020, 10:59am UTC](https://discourse.julialang.org/t/allocations-when-constructing-a-matrix-from-columns/36795/1 "2020-03-31T10:59:18Z")

</div>

Hey all!

I’m constructing a matrix from another matrix like this:

```julia
shapematrix(xs) = [xs[:,1].-xs[:,4] xs[:,2].-xs[:,4] xs[:,3].-xs[:,4]]

```

This is somewhat in an inner loop of my program, and I would like to reduce the number of allocations happening. `@time` reports that there are 13 allocations here, when I would think a single one would be sufficient, and it seems no matter how I rewrite the function, the number of allocations is always in the double digits.

Is there a way of writing this that makes Julia be a little more clever about it memory management?

I’m also using Zygote.jl, which doesn’t support mutating `Array`s, so I can’t make a zero array and fill in the entries afterwards.

I’m new-ish to Julia, so maybe there’s something obvious going on here?

---

<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 31, 2020, 1:06pm UTC](https://discourse.julialang.org/t/allocations-when-constructing-a-matrix-from-columns/36795/2 "2020-03-31T13:06:26Z")

</div>

> [@mht](#):
>
> ```julia
> shapematrix(xs) = [xs[:,1].-xs[:,4] xs[:,2].-xs[:,4] xs[:,3].-xs[:,4]]
> 
> ```

Note first that you can use broadcasting to rewrite this:

```julia
shapematrix2(xs) = xs[:, 1:3] .- xs[:, 4]

```

In order to reduce allocations, you can use views:

```julia
shapematrix3(xs) = @views xs[:, 1:3] .- xs[:, 4]

```

Let’s benchmark, and remember to use BenchmarkTools for this, not the `@time` macro, although the allocation results look similar:

```julia
julia> X = rand(1000, 4);

julia> @btime shapematrix($X);
  11.534 μs (14 allocations: 95.05 KiB)

julia> @btime shapematrix2($X);
  5.226 μs (5 allocations: 54.97 KiB)

julia> @btime shapematrix3($X);
  1.895 μs (4 allocations: 23.63 KiB)

```

---

<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 31, 2020, 1:16pm UTC](https://discourse.julialang.org/t/allocations-when-constructing-a-matrix-from-columns/36795/3 "2020-03-31T13:16:57Z")

</div>

> [@mht](#):
>
> I’m also using Zygote.jl, which doesn’t support mutating `Array` s, so I can’t make a zero array and fill in the entries afterwards.

I haven’t used Zygote, but this sounds odd. All of these array operations initialize arrays and then modify them afterwards, even if you don’t see it explicitly in the top-level code.

---

<div class="post-metadata">

**Author:** ![mht](https://avatars.discourse-cdn.com/v4/letter/m/ea5d25/32.png) [@mht](https://discourse.julialang.org/u/mht)\
**Post date:** [March 31, 2020, 1:44pm UTC](https://discourse.julialang.org/t/allocations-when-constructing-a-matrix-from-columns/36795/4 "2020-03-31T13:44:39Z")

</div>

Thanks! I’ve tried to use views before, but for each column and as written in the OP, and this didn’t help anything so I figured that wasn’t the problem.

I’d still like not to use broadcasting excessively, since I find it utterly unreadable, but it’ll do for now.

---

<div class="post-metadata">

**Author:** ![mht](https://avatars.discourse-cdn.com/v4/letter/m/ea5d25/32.png) [@mht](https://discourse.julialang.org/u/mht)\
**Post date:** [March 31, 2020, 1:48pm UTC](https://discourse.julialang.org/t/allocations-when-constructing-a-matrix-from-columns/36795/5 "2020-03-31T13:48:46Z")

</div>

I don’t know anything about Zygotes internals, but an educated guess is that if only internal procedures do the mutation then Zygote has a well defined boundary where it can reason about the operations done to the matrices and calculate gradients based on that.

I know for a fact that the mutation doesn’t work, because if you try you get the error `ERROR: Mutating arrays is not supported` 😄

---

<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 31, 2020, 2:03pm UTC](https://discourse.julialang.org/t/allocations-when-constructing-a-matrix-from-columns/36795/6 "2020-03-31T14:03:33Z")

</div>

> [@mht](#):
>
> I’d still like not to use broadcasting excessively, since I find it utterly unreadable, but it’ll do for now.

Really? It exists to make code easier to read and write, that’s sort of the point. Comparing

```julia
[xs[:,1].-xs[:,4] xs[:,2].-xs[:,4] xs[:,3].-xs[:,4]]
# and 
xs[:, 1:3] .- xs[:, 4]

```

I find the former really hard to read, and the latter very clear.

(Stream of consciousness: “The first line is, hmm a jumble of expressions, oh, it’s horizontal concatenation. And it’s hard to see the pattern, because the indices are jumping, and oh, it’s minus between each term, is it? yes, and 1, 4, 2, something, oh, 1,4,2,4,3,4, ah, 1,2,3 minus the fourth. ok.” The latter is just “Subtract the 4th column from the other three.”)

I’d be interested in hearing what you don’t like about it. Maybe there’s some little mental block that can be cleared?

---

<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 31, 2020, 2:13pm UTC](https://discourse.julialang.org/t/allocations-when-constructing-a-matrix-from-columns/36795/7 "2020-03-31T14:13:57Z")

</div>

In addition to @DNF’s excellent suggestions, if this is really a bottleneck in your code then you could write a micro-optimized mutating version and then define the adjoint manually. See the Zygote manual for details. (I am just pointing out the technical possibility, but I would guess that it is overkill here).

That said, if you are new to Julia then focusing on allocations like this may not be the best allocation (heh) of your time.

---

<div class="post-metadata">

**Author:** ![mht](https://avatars.discourse-cdn.com/v4/letter/m/ea5d25/32.png) [@mht](https://discourse.julialang.org/u/mht)\
**Post date:** [March 31, 2020, 2:20pm UTC](https://discourse.julialang.org/t/allocations-when-constructing-a-matrix-from-columns/36795/8 "2020-03-31T14:20:05Z")

</div>

It’s shorter, no doubt, but it’s also way denser. Originally, I wrote out the 9 terms without any use of `:`, and I think that’s the way I ultimately prefer, but I had to go back and forth a bit, in case there were some performance tricks that would only happen with a specific syntax (and wouldn’t you know!)

Here’s my stream of consciousness:

```julia
[xs[1,1]-xs[1,4] xs[1,2]-xs[1,4] xs[1,3]-xs[1,4] ;
  xs[2,1]-xs[2,4] xs[2,2]-xs[2,4] xs[2,3]-xs[2,4] ;
  xs[3,1]-xs[3,4] xs[3,2]-xs[3,4] xs[3,3]-xs[3,4] ]

```

oh, it’s a 3x3 matrix, and `A[i,j] = ` …uuh lets see … `xs[i,j] - xs[i,4]`.

```julia
[xs[:,1].-xs[:,4] xs[:,2].-xs[:,4] xs[:,3].-xs[:,4]]

```

Okay so it’s an array, and the elements are, oh no wait, `xs[:,1]` is a vector, so we’ll end up with a matrix where … uuh … the vectors here become the columns?, and then we broadcast `xs[:,4]` which is also a vector, so the broadcasting is really just element wise, okay.

```julia
xs[:, 1:3] .- xs[:, 4]

```

Okay so we index the vector, oh with ranges in both places, okay, I guess we’ll get the sub matrix? and then we’ll broadcast sub `xs[:,4]` which is a vector, and subtraction of a vector needs another vector, so probably this will go across the columns? uuh wait `xs` is 3 high and so `xs[:,4]` is 3 long or uuh high, and so we’ll have to go across the columns in `xs`.

* * *

Whether I write the second one, or eventually get to the last one, isn’t really important to me, since I find the first version way clearer. In addition, I would think that spelling it out explicitly would make Julia’s job of generating good code for this easier, so that I don’t have to rewrite all of my functions when I realize that bad code is being generated. Somehow I got it completely backwards this time?

Anyways, thanks for the help! 🙂

---

<div class="post-metadata">

**Author:** ![mht](https://avatars.discourse-cdn.com/v4/letter/m/ea5d25/32.png) [@mht](https://discourse.julialang.org/u/mht)\
**Post date:** [March 31, 2020, 2:25pm UTC](https://discourse.julialang.org/t/allocations-when-constructing-a-matrix-from-columns/36795/9 "2020-03-31T14:25:06Z")

</div>

This single change cut down 1/3 of my allocations for my test case, so this is absolutely a bottleneck.

> and then define the adjoint manually.

I could, but the entire reason I’m using Zygote is so that I don’t have to write out derivatives manually, and so it doesn’t make much sense for me to spend time to learn the ins and outs of Zygote just to write more code, when the whole reason for me using it is to _not_ write code 😄

> That said, if you are new to Julia then focusing on allocations like this may not be the best allocation (heh) of your time.

I appreciate that you’re trying to help, but I got a simulation that runs orders of magnitude too slow than expeted because Julia generates bad code (due to, it seems, me writing bad Julia), and so this is definitely worth my time. I think it’s unfortunate that any time anyone is asking about making things faster people (and not to pick on you in particular!) jumps in screaming premature optimization ☹

---

<div class="post-metadata">

**Author:** ![ericphanson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ericphanson/32/215186_2.png) [@ericphanson](https://discourse.julialang.org/u/ericphanson)\
**Post date:** [March 31, 2020, 2:30pm UTC](https://discourse.julialang.org/t/allocations-when-constructing-a-matrix-from-columns/36795/10 "2020-03-31T14:30:03Z")

</div>

> [@mht](#):
>
> This single change cut down 1/3 of my allocations for my test case, so this is absolutely a bottleneck.

You might get a lot of performance benefit by reusing the memory by passing along the matrix from iteration to iteration and updating it in-place. If you’re matrices are this small, you might also see a lot of benefit from using StaticArrays from StaticArrays.jl (in which case you would not write it in-place). I’m on my phone so I won’t write more here but there’s a lot of examples on these forums of doing these changes and getting big performance improvements.

---

<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 31, 2020, 2:32pm UTC](https://discourse.julialang.org/t/allocations-when-constructing-a-matrix-from-columns/36795/11 "2020-03-31T14:32:34Z")

</div>

> [@mht](#):
>
> the entire reason I’m using Zygote is so that I don’t have to write out derivatives manually

I think you are mistaken about this.

Having AD does not mean that you never, ever have to code derivatives, just that you can focus on which ones you want to code manually (note that you also have a choice here, if you stick to the non-modifying version), and that the chaining will be organized nicely for you.

> [@mht](#):
>
> jumps in screaming premature optimization ☹

I am not sure I understand you. No one is screaming here.

It is very likely that there are alternative solutions depending the context of this function (which we do not have, so it is hard to help), and focusing on optimizing this particular form may not be the best solution. This is natural if you are new to the language, and pointing this out should not be taken personally.

---

<div class="post-metadata">

**Author:** ![mht](https://avatars.discourse-cdn.com/v4/letter/m/ea5d25/32.png) [@mht](https://discourse.julialang.org/u/mht)\
**Post date:** [March 31, 2020, 2:57pm UTC](https://discourse.julialang.org/t/allocations-when-constructing-a-matrix-from-columns/36795/12 "2020-03-31T14:57:00Z")

</div>

This was my first idea too, but I’m using Zygote to differentiate the functions, which complains when I mutate arrays, as written above ☹

---

<div class="post-metadata">

**Author:** ![mht](https://avatars.discourse-cdn.com/v4/letter/m/ea5d25/32.png) [@mht](https://discourse.julialang.org/u/mht)\
**Post date:** [March 31, 2020, 3:05pm UTC](https://discourse.julialang.org/t/allocations-when-constructing-a-matrix-from-columns/36795/13 "2020-03-31T15:05:03Z")

</div>

> Having AD does not mean that you never, ever have to code derivatives, just that you can focus on which ones you want to code manually

I’ve already written some adjoints manually due to exessive allocations in other places; in the limit it looks to me that I’ll have to write my primals carefully to be 100% sure that Julia is generating decent code, _and_ write the adjoints manually to make sure that Zygote understands my primals. At this point, it doesn’t seem to me that using Zygote, or even Julia, makes a lot of sense. Yet, as far as I understand the goals of both projects, this is textbook usage. It’s frustrating, and I don’t understand if I’m just not getting it, or what’s going on.

> I am not sure I understand you. No one is screaming here.

No, but you were suggesting that the this really wasn’t a bottleneck, and that this might not be the best use of my time, which at best, is presumptuous. The reason I’m out here asking for help is, of course, because it _is_ a problem. (Again, I didn’t mean to pick on you, but this is what you always get when asking for help with performance related things pretty much anywhere online, and it’s a pet peeve of mine 😉).

> It is very likely that there are alternative solutions depending the context of this function (which we do not have, so it is hard to help), and focusing on optimizing this particular form may not be the best solution.

This might be, but I intentionally posted a very small sub-problem, since it was a small self-contained part of Julia that I had (and still have!) problems properly understanding.

---

<div class="post-metadata">

**Author:** ![mcabbott](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mcabbott/32/6603_2.png) [@mcabbott](https://discourse.julialang.org/u/mcabbott)\
**Post date:** [March 31, 2020, 3:46pm UTC](https://discourse.julialang.org/t/allocations-when-constructing-a-matrix-from-columns/36795/14 "2020-03-31T15:46:36Z")

</div>

If 3x4 is the size, then the above options are:

```julia
julia> @btime shapematrix(xs) setup=(xs=randn(3,4));
  413.700 ns (13 allocations: 1.23 KiB)

julia> @btime shapematrix3(xs) setup=(xs=randn(3,4));
  64.027 ns (3 allocations: 272 bytes)

```

These are dwarfed by the cost of the gradient… partly because every indexing operation first makes a new `zero(xs)`, then writes into that. Writing even a pretty crude gradient really does pay:

```julia
julia> @btime Zygote.gradient(sum∘shapematrix3, xs) setup=(xs=randn(3,4));
  6.862 μs (122 allocations: 4.70 KiB)

julia> Zygote.@adjoint shapematrix4(xs) = shapematrix4(xs), dys -> (dys .- [0,0,0,1]' .* sum(dys, dims=2),)

julia> @btime Zygote.gradient(sum∘shapematrix4, xs) setup=(xs=randn(3,4));
  773.673 ns (13 allocations: 848 bytes)

```

But for such small arrays, as others have mentioned, you probably want to do something more like this:

```julia
julia> using ForwardDiff, StaticArrays

julia> @btime ForwardDiff.gradient(sum∘shapematrix4, xs) setup=(xs=randn(3,4));
  598.392 ns (5 allocations: 4.02 KiB)

julia> @btime ForwardDiff.gradient(sum∘shapematrix4, xs) setup=(xs=@SArray randn(3,4));
  181.366 ns (2 allocations: 1.34 KiB)

```

You can also do `Zygote.gradient(xs -> sum(Zygote.forwarddiff(shapematrix4, xs)), xs)` but this seems to have more overhead here.

---

<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 31, 2020, 3:56pm UTC](https://discourse.julialang.org/t/allocations-when-constructing-a-matrix-from-columns/36795/15 "2020-03-31T15:56:26Z")

</div>

> [@mht](#):
>
> this is what you always get when asking for help with performance related things pretty much anywhere online, and it’s a pet peeve of mine 😉).

Maybe you’re new here, but this is a board full of almost obsessive microoptimization enthusiasts. You can hardly drop a single line of code here without a hailstorm of optimization advice raining down on you, solicited or unsolicited.

So, if you have a different impression, this must be an off day.

---

<div class="post-metadata">

**Author:** ![mht](https://avatars.discourse-cdn.com/v4/letter/m/ea5d25/32.png) [@mht](https://discourse.julialang.org/u/mht)\
**Post date:** [March 31, 2020, 3:57pm UTC](https://discourse.julialang.org/t/allocations-when-constructing-a-matrix-from-columns/36795/16 "2020-03-31T15:57:40Z")

</div>

The function I’m taking the gradient of is much bigger than my toy example. I’m trying to move things to `SMatrix`, but ran into some problems with Zygote and StaticArrays (which I posted in Zygote’s issue tracker).

Thanks anyways 😄

---

<div class="post-metadata">

**Author:** ![mht](https://avatars.discourse-cdn.com/v4/letter/m/ea5d25/32.png) [@mht](https://discourse.julialang.org/u/mht)\
**Post date:** [March 31, 2020, 3:58pm UTC](https://discourse.julialang.org/t/allocations-when-constructing-a-matrix-from-columns/36795/17 "2020-03-31T15:58:42Z")

</div>

I am indeed new here! 😉

---

<div class="post-metadata">

**Author:** ![mcabbott](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mcabbott/32/6603_2.png) [@mcabbott](https://discourse.julialang.org/u/mcabbott)\
**Post date:** [March 31, 2020, 4:18pm UTC](https://discourse.julialang.org/t/allocations-when-constructing-a-matrix-from-columns/36795/18 "2020-03-31T16:18:10Z")

</div>

OK, but then it’s hard to guess what will help. Care to write a more representative toy problem?

(Likewise the [zygote issue](https://github.com/FluxML/Zygote.jl/issues/570), `gradient(identity, 1)` works pretty well!)

---

<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 31, 2020, 5:10pm UTC](https://discourse.julialang.org/t/allocations-when-constructing-a-matrix-from-columns/36795/19 "2020-03-31T17:10:52Z")

</div>

> [@mht](#):
>
> No, but you were suggesting that the this really wasn’t a bottleneck, and that this might not be the best use of my time, which at best, is presumptuous.

I wasn’t suggesting that this isn’t a bottleneck for you, simply that focusing on the allocations _per se_ (and keeping other things unchanged) may not be the best approach.

Eg if the `xs` small, you can write very clean, idiomatic code using `StaticArrays` (as suggested by others), which should be rather fast, and Zygote will just work fine with it.

The Julia ecosystem has evolved some very efficient techniques for dealing with the issues you are facing, but it is very hard to help without context.

---

<div class="post-metadata">

**Author:** ![Mason](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mason/32/2423_2.png) [@Mason](https://discourse.julialang.org/u/Mason)\
**Post date:** [March 31, 2020, 6:47pm UTC](https://discourse.julialang.org/t/allocations-when-constructing-a-matrix-from-columns/36795/20 "2020-03-31T18:47:01Z")

</div>

I’m not sure if this is a fruitful approach, but if mutation is the only way to get the performance you need, then perhaps the optimal strategy is to use something like [FiniteDiff.jl](https://github.com/JuliaDiff/FiniteDiff.jl) for calculating the pullbacks of the mutating code and then using Zygote’s [custom adjoint machinery](https://fluxml.ai/Zygote.jl/dev/adjoints/) to embed that so that the rest of your code can be handled by Zygote.

I tried to make a toy example, but I don’t have enough information about your actual problem and I also wasn’t too sure about the correct way to mix Zygote.jl with FiniteDiff.jl, but maybe someone here can cook up a good example.

[Next page](https://discourse.julialang.org/t/allocations-when-constructing-a-matrix-from-columns/36795.md?page=2)