# Outer broadcasting?

**URL:** <https://discourse.julialang.org/t/outer-broadcasting/22701>\
**Category:** General Usage\
**Created:** [April 3, 2019, 3:03pm UTC](https://discourse.julialang.org/t/outer-broadcasting/22701 "2019-04-03T15:03:44Z")\
**Posts on this page:** 12\
**Page:** 1

<div class="post-metadata">

**Author:** ![tim](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tim/32/6848_2.png) [@tim](https://discourse.julialang.org/u/tim)\
**Post date:** [April 3, 2019, 3:03pm UTC](https://discourse.julialang.org/t/outer-broadcasting/22701/1 "2019-04-03T15:03:44Z")

</div>

Hi,

I would like to do something like the [outer product](https://en.wikipedia.org/wiki/Outer_product) for general operations. Basically, I f I have a function

```julia
f(x1, x2, [...], xn) = [...]

```

and I provide it with n arrays `a1` … `an` as input, then I would like to apply `f` on all combination of input elements, such that the resulting array has the shape:

```julia
outer(f, a1, a2, [...], an) |> size == (size(a1)..., size(a2)..., [...], size(an)...)

```

_How can I do that in a clean and performant way?_

_Example:_

```julia
julia> f(a, b, c) = a+b+c
f (generic function with 1 method)

julia> x = 1:300; y = 1:400; z = 1:500;

# performant but super hard to read
julia> f.(x, reshape(y, (ones(Int64, ndims(x))..., size(y)...)), reshape(z, (ones(Int64, ndims(x) + ndims(y))..., size(z)...))) |> size == (size(x)..., size(y)..., size(z)...)
true

julia> @benchmark f.(x, reshape(y, (ones(Int64, ndims(x))..., size(y)...)), reshape(z, (ones(Int64, ndims(x) + ndims(y))..., size(z)...)))
BenchmarkTools.Trial: 
  memory estimate: 457.76 MiB
  allocs estimate: 13
  --------------
  minimum time: 141.622 ms (0.29% GC)
  median time: 150.833 ms (9.09% GC)
  mean time: 158.361 ms (10.34% GC)
  maximum time: 193.617 ms (22.94% GC)
  --------------
  samples: 32
  evals/sample: 1

# easy to read but slow
julia> broadcast( u -> f(u...), Base.Iterators.product(x, y, z)) |> size == (size(x)..., size(y)..., size(z)...)
true

julia> @benchmark broadcast( u -> f(u...), Base.Iterators.product(x, y, z))
BenchmarkTools.Trial: 
  memory estimate: 1.79 GiB
  allocs estimate: 6
  --------------
  minimum time: 731.471 ms (1.30% GC)
  median time: 799.173 ms (13.66% GC)
  mean time: 792.291 ms (11.79% GC)
  maximum time: 820.479 ms (11.84% GC)
  --------------
  samples: 7
  evals/sample: 1

```

The first version seems performant, the second one seems clean. Is there a standard function that I am overlooking?

---

<div class="post-metadata">

**Author:** ![yha](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yha/32/3502_2.png) [@yha](https://discourse.julialang.org/u/yha)\
**Post date:** [April 3, 2019, 3:34pm UTC](https://discourse.julialang.org/t/outer-broadcasting/22701/2 "2019-04-03T15:34:26Z")

</div>

`Iterators.product(arrays...)` returns an outer product iterator like you’re looking for:

```julia
julia> a = 1:4;
julia> b = [(x,y) for x=5:7, y=8:9]
3×2 Array{Tuple{Int64,Int64},2}:
 (5, 8) (5, 9)
 (6, 8) (6, 9)
 (7, 8) (7, 9)

julia> collect(Iterators.product(a,b))
4×3×2 Array{Tuple{Int64,Tuple{Int64,Int64}},3}:
[:, :, 1] =
 (1, (5, 8)) (1, (6, 8)) (1, (7, 8))
 (2, (5, 8)) (2, (6, 8)) (2, (7, 8))
 (3, (5, 8)) (3, (6, 8)) (3, (7, 8))
 (4, (5, 8)) (4, (6, 8)) (4, (7, 8))

[:, :, 2] =
 (1, (5, 9)) (1, (6, 9)) (1, (7, 9))
 (2, (5, 9)) (2, (6, 9)) (2, (7, 9))
 (3, (5, 9)) (3, (6, 9)) (3, (7, 9))
 (4, (5, 9)) (4, (6, 9)) (4, (7, 9))

```

So you can do

```julia
outer(f,arrays...) = [f(x...) for x in Iterators.product(arrays...)]

```

---

<div class="post-metadata">

**Author:** ![tim](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tim/32/6848_2.png) [@tim](https://discourse.julialang.org/u/tim)\
**Post date:** [April 4, 2019, 1:15pm UTC](https://discourse.julialang.org/t/outer-broadcasting/22701/3 "2019-04-04T13:15:53Z")

</div>

Thanks a lot, @yha. You basically suggested the same as my second solution, just with a liste comprehension, right? And seemingly your idea is way faster than mine with the braodcasting, so I am going to use yours 😄

But do you know why the list comprehension is so much faster than the broadcasting? I didn’t expect that.

```julia
julia> f(a, b, c) = a+b+c
f (generic function with 1 method)

julia> x = 1:300; y = 1:400; z = 1:500;

julia> using BenchmarkTools

julia> @benchmark broadcast( u -> f(u...), Base.Iterators.product(x, y, z))
BenchmarkTools.Trial: 
  memory estimate: 1.79 GiB
  allocs estimate: 6
  --------------
  minimum time: 733.547 ms (7.83% GC)
  median time: 760.383 ms (12.34% GC)
  mean time: 765.740 ms (12.62% GC)
  maximum time: 797.509 ms (15.32% GC)
  --------------
  samples: 7
  evals/sample: 1

julia> @benchmark [f(q...) for q in Iterators.product(x, y, z)]
BenchmarkTools.Trial: 
  memory estimate: 457.76 MiB
  allocs estimate: 5
  --------------
  minimum time: 192.790 ms (6.51% GC)
  median time: 199.686 ms (6.35% GC)
  mean time: 202.157 ms (7.98% GC)
  maximum time: 238.016 ms (19.89% GC)
  --------------
  samples: 25
  evals/sample: 1

```

---

<div class="post-metadata">

**Author:** ![rdeits](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rdeits/32/286_2.png) [@rdeits](https://discourse.julialang.org/u/rdeits)\
**Post date:** [April 4, 2019, 1:29pm UTC](https://discourse.julialang.org/t/outer-broadcasting/22701/4 "2019-04-04T13:29:49Z")

</div>

Another solution, which is faster than both of the options presented so far, is to use broadcasting directly, without any need for `Iterators.product`:

```julia
julia> @btime f.($x, reshape($y, (1, :)), reshape($z, (1, 1, :)));
  114.924 ms (3 allocations: 457.76 MiB)

```

---

<div class="post-metadata">

**Author:** ![tim](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tim/32/6848_2.png) [@tim](https://discourse.julialang.org/u/tim)\
**Post date:** [April 4, 2019, 1:52pm UTC](https://discourse.julialang.org/t/outer-broadcasting/22701/5 "2019-04-04T13:52:56Z")

</div>

Thanks @rdeits. Actually, my first version is the generalization of your suggestion. And on my computer they run both with the same speed, I assume your cpu is faster 😄

But again, the moment one generalizes this idea, it becomes rather hard to read.

```julia
julia> f(a, b, c) = a+b+c
f (generic function with 1 method)

julia> x = 1:300; y = 1:400; z = 1:500;

julia> using BenchmarkTools

julia> @benchmark f.(x, reshape(y, (ones(Int64, ndims(x))..., size(y)...)), reshape(z, (ones(Int64, ndims(x) + ndims(y))..., size(z)...)))
BenchmarkTools.Trial: 
  memory estimate: 457.76 MiB
  allocs estimate: 13
  --------------
  minimum time: 146.591 ms (8.67% GC)
  median time: 162.362 ms (8.05% GC)
  mean time: 193.845 ms (10.20% GC)
  maximum time: 367.354 ms (17.39% GC)
  --------------
  samples: 26
  evals/sample: 1

julia> @benchmark f.(x, reshape(y, (1, :)), reshape(z, (1, 1, :)))
BenchmarkTools.Trial: 
  memory estimate: 457.76 MiB
  allocs estimate: 7
  --------------
  minimum time: 147.229 ms (8.59% GC)
  median time: 150.185 ms (8.76% GC)
  mean time: 155.161 ms (10.41% GC)
  maximum time: 182.251 ms (26.30% GC)
  --------------
  samples: 33
  evals/sample: 1

```

---

<div class="post-metadata">

**Author:** ![rdeits](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rdeits/32/286_2.png) [@rdeits](https://discourse.julialang.org/u/rdeits)\
**Post date:** [April 4, 2019, 2:01pm UTC](https://discourse.julialang.org/t/outer-broadcasting/22701/6 "2019-04-04T14:01:25Z")

</div>

Whoops, sorry, I missed that 🙂

In that case, this seems like a great time for some metaprogramming! How does this look?

```julia
julia> @generated function outer(f, args::Tuple{Vararg{Any, N}}) where {N}
         expr = Expr(:call, :broadcast, :f)
         for i in 1:N
           if i == 1
             push!(expr.args, :(args[1]))
           else
             push!(expr.args, Expr(:call, :reshape, :(args[$i]), Expr(:tuple, [1 for _ in 1:i - 1]..., :(:))))
           end
         end
         expr
       end
outer (generic function with 1 method)

julia> @btime outer($f, ($x, $y, $z));
  108.493 ms (3 allocations: 457.76 MiB)

```

---

<div class="post-metadata">

**Author:** ![mbauman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mbauman/32/31082_2.png) [@mbauman](https://discourse.julialang.org/u/mbauman)\
**Post date:** [April 4, 2019, 2:22pm UTC](https://discourse.julialang.org/t/outer-broadcasting/22701/7 "2019-04-04T14:22:08Z")

</div>

The list comprehension over `Iterators.product` is faster than broadcasting because broadcasting requires index-ability. Since `Iterators.product` is a simple iterable, `broadcast` makes a copy of it with `collect` before doing its thing.

This is something I do think it’d be nice to have a better solution for — we’ve been brainstorming an API for it on the back-burner. Suggestions have included things like `f.(orthogonalize(x, y, z)...)` or `f.(x, ^y, ^z)`.

---

<div class="post-metadata">

**Author:** ![yha](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yha/32/3502_2.png) [@yha](https://discourse.julialang.org/u/yha)\
**Post date:** [April 4, 2019, 3:37pm UTC](https://discourse.julialang.org/t/outer-broadcasting/22701/8 "2019-04-04T15:37:53Z")

</div>

> [@tim](#):
>
> Thanks a lot, @yha. You basically suggested the same as my second solution, just with a liste comprehension, right? And seemingly your idea is way faster than mine with the braodcasting, so I am going to use yours 😄

Yes, apparently I skimmed your post too quickly, because I missed the fact that you were already using `Iterators.product` in your second solution.

---

<div class="post-metadata">

**Author:** ![yha](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yha/32/3502_2.png) [@yha](https://discourse.julialang.org/u/yha)\
**Post date:** [April 4, 2019, 3:42pm UTC](https://discourse.julialang.org/t/outer-broadcasting/22701/9 "2019-04-04T15:42:05Z")

</div>

Can’t `Iterators.product` itself be made indexable when its inputs are? It already has shape when its input have shape – an unusual property for iterators outside the Julia world – so it doesn’t seem like a big leap.  
Though maybe the fact that it’s in `Iterators` would make it less discoverable, since it’s somewhat surprising that “iterators” have shape and/or are indexable.

---

<div class="post-metadata">

**Author:** ![mbauman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mbauman/32/31082_2.png) [@mbauman](https://discourse.julialang.org/u/mbauman)\
**Post date:** [April 4, 2019, 4:30pm UTC](https://discourse.julialang.org/t/outer-broadcasting/22701/10 "2019-04-04T16:30:01Z")

</div>

> [@yha](#):
>
> Can’t `Iterators.product` itself be made indexable when its inputs are?

Yup, it sure could. It’d be as simple as defining `axes`, `getindex`, and `Broadcast.broadcastable` (to behave as itself) for `Base.Iterators.ProductIterator{Tuple{Vararg{AbstractArray}}}`.

---

<div class="post-metadata">

**Author:** ![tim](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tim/32/6848_2.png) [@tim](https://discourse.julialang.org/u/tim)\
**Post date:** [April 5, 2019, 6:42am UTC](https://discourse.julialang.org/t/outer-broadcasting/22701/11 "2019-04-05T06:42:29Z")

</div>

Then maybe I will suggest that. I generally like to use broadcasting so julia can collapse some of the loops when concatenating multiple functions.

I also like the idea of `orthogonolize` or the similar. As you said you have been brainstorming: Is there an issue on that?

---

<div class="post-metadata">

**Author:** ![mbauman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mbauman/32/31082_2.png) [@mbauman](https://discourse.julialang.org/u/mbauman)\
**Post date:** [April 5, 2019, 2:18pm UTC](https://discourse.julialang.org/t/outer-broadcasting/22701/12 "2019-04-05T14:18:19Z")

</div>

> [@tim](#):
>
> As you said you have been brainstorming: Is there an issue on that?

It’s a bit of a tangent to the actual issue, but we’ve been talking about it in [What to do about nonscalar indexing? · Issue #30845 · JuliaLang/julia · GitHub](https://github.com/JuliaLang/julia/issues/30845).