# Repeated vector vector multiplication to obtain matrix

**URL:** <https://discourse.julialang.org/t/repeated-vector-vector-multiplication-to-obtain-matrix/42980>\
**Category:** Performance\
**Tags:** question\
**Created:** [July 13, 2020, 12:36pm UTC](https://discourse.julialang.org/t/repeated-vector-vector-multiplication-to-obtain-matrix/42980 "2020-07-13T12:36:59Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![morkip](https://avatars.discourse-cdn.com/v4/letter/m/35a633/32.png) [@morkip](https://discourse.julialang.org/u/morkip)\
**Post date:** [July 13, 2020, 12:36pm UTC](https://discourse.julialang.org/t/repeated-vector-vector-multiplication-to-obtain-matrix/42980/1 "2020-07-13T12:36:59Z")

</div>

In my application I am running many iterations of a multiplication like this `mat = vec1 .* vec2'`.

The current implementation is like this

```julia
dP_dw = zeros(Float64, (number_of_states, settings.n_vis * 2, settings.n_hid))
for a = 1:number_of_states
    dP_dw[a, :, :] = states[a, :] .* tanh_term[a, :]'
end

```

I was wondering if there was a more ‘julianic’ / faster way of doing this instead of for loops using the broadcast function, but I cannot put it together as customizing it seems tough.

---

<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:** [July 13, 2020, 12:57pm UTC](https://discourse.julialang.org/t/repeated-vector-vector-multiplication-to-obtain-matrix/42980/2 "2020-07-13T12:57:24Z")

</div>

Here’s a quick attempt:

```julia
function outers1(states, tanh_term) 
    number_of_states = size(states, 1)
    size(states, 1) == number_of_states || error()
    dP_dw = zeros(Float64, (number_of_states, size(states,2), size(tanh_term,2)))
    @inbounds for a = 1:number_of_states
        @views dP_dw[a, :, :] .= states[a, :] .* tanh_term[a, :]'
    end
    dP_dw
end

using Einsum, Test
outers2(states, tanh_term) = @einsum dP_dw[a,s,t] := states[a, s] * tanh_term[a, t]
outers3(states, tanh_term) = @vielsum dP_dw[a,s,t] := states[a, s] * tanh_term[a, t]

outers4(states, tanh_term) = states .* reshape(tanh_term, size(states, 1), 1, :)

N = 50; states, tanh_term = randn(N,N), randn(N,N);
@test outers1(states, tanh_term) ≈ outers2(states, tanh_term) ≈ outers4(states, tanh_term)

#== results ==#

julia> @btime outers0($states, $tanh_term); # as in question
  364.486 μs (202 allocations: 1.96 MiB)

julia> @btime outers1($states, $tanh_term); # with @inbounds, @views, and .=
  189.673 μs (2 allocations: 976.64 KiB)

julia> @btime outers2($states, $tanh_term); 
  62.388 μs (2 allocations: 976.64 KiB)

julia> @btime outers3($states, $tanh_term);
  52.462 μs (62 allocations: 984.86 KiB)

julia> @btime outers4($states, $tanh_term);                   
  63.387 μs (4 allocations: 976.75 KiB) 

```

If you look at `@macroexpand1 @einsum dP_dw[a,s,t] := states[a, s] * tanh_term[a, t]`, the main change is that it orders the loops with `a` innermost.

---

<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:** [July 13, 2020, 1:31pm UTC](https://discourse.julialang.org/t/repeated-vector-vector-multiplication-to-obtain-matrix/42980/3 "2020-07-13T13:31:28Z")

</div>

> [@morkip](#):
>
> In my application I am running many iterations of a multiplication like this `mat = vec1 .* vec2'` .  
> […]  
> I was wondering if there was a more ‘julianic’ / faster way of doing this instead of for loops using the broadcast function, but I cannot put it together as customizing it seems tough.

First reconsider whether you should construct such matrices at all. You are constructing rank-1 matrices M = uv^\*, which requires \Theta(n^2) storage and time for n-component vectors, and naively multiplying such a matrix by a vector (Mx) requires \Theta(n^2) time. In contrast, storing u and v separately requires \Theta(n) storage and you can multiply (uv^\*)x = u(v^\*x) in \Theta(n) time.

If you must compute the matrices explicitly, then

1. Realize that `dP_dw[a, :, :] = states[a, :] .* tanh_term[a, :]'` first allocates two vectrors (`states[a, :]` and `tanh_term[a, :]`) since [slicing makes a copy](https://docs.julialang.org/en/v1/manual/performance-tips/index.html#Consider-using-views-for-slices-1) in Julia (unless you use `@views`), then allocates a matrix, then writes that matrix into `dP_dw`. That is a lot of allocations.

2. You can partly eliminate the allocations by using `.=` and `@views` as in @mcabbot’s post, but I think you may need Julia 1.5 for this to elide allocations of view objects?

3. If you have small fixed-size arrays whose size is known at compile-time, e.g. 3-component or 4-component arrays, you should use the StaticArrays package (which will completely unroll your loops and eliminate all allocations).

(I’m assuming that your critical code is actually in a function, of course — [don’t use global](https://docs.julialang.org/en/v1/manual/performance-tips/index.html#Avoid-global-variables-1) loops.)

---

<div class="post-metadata">

**Author:** ![morkip](https://avatars.discourse-cdn.com/v4/letter/m/35a633/32.png) [@morkip](https://discourse.julialang.org/u/morkip)\
**Post date:** [July 13, 2020, 1:51pm UTC](https://discourse.julialang.org/t/repeated-vector-vector-multiplication-to-obtain-matrix/42980/4 "2020-07-13T13:51:31Z")

</div>

Thanks for the effort, I appreciate it! Impressive speed up!

---

<div class="post-metadata">

**Author:** ![morkip](https://avatars.discourse-cdn.com/v4/letter/m/35a633/32.png) [@morkip](https://discourse.julialang.org/u/morkip)\
**Post date:** [July 13, 2020, 2:09pm UTC](https://discourse.julialang.org/t/repeated-vector-vector-multiplication-to-obtain-matrix/42980/5 "2020-07-13T14:09:32Z")

</div>

Sadly I do need them, as these matrices represent derivatives in a neural network and I need every possible combination of the different neurons (fully connected layers), but it’s a good point.

`a` will usually be as large as 10000, while the second dimension `s` in @mcabbott post will be as large as 40 and `t` can go up to 20 - 30. StaticArrays won’t help sadly.

`states` is a Boolean Array. Does Einsum / Vielsum in some way or another incorporate this?

I’ll do a bit of benchmarking myself and see if I can make any further improvements. Thanks a lot again to both of you! (and no worries, everything’s in a function)