# Outer product & broadcast?

**URL:** <https://discourse.julialang.org/t/outer-product-broadcast/103731>\
**Category:** General Usage\
**Created:** [September 11, 2023, 7:56am UTC](https://discourse.julialang.org/t/outer-product-broadcast/103731 "2023-09-11T07:56:48Z")\
**Posts on this page:** 9\
**Page:** 1

<div class="post-metadata">

**Author:** ![ryofurue](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ryofurue/32/24531_2.png) [@ryofurue](https://discourse.julialang.org/u/ryofurue)\
**Post date:** [September 11, 2023, 7:56am UTC](https://discourse.julialang.org/t/outer-product-broadcast/103731/1 "2023-09-11T07:56:49Z")

</div>

I often need gridded versions of continuous functions. For a function with one argument, the usual broad cast is the most elegant way:

```julia
oned(x) = 2*x # 1D func
xs = 0.0:1.0:10.0 # x grid
oned_gridded = oned.(xs)

```

This method is often found in Julia examples.

Then, what about multiple-argument functions? You can of course use the list comprehension, but it’s a bit more verbose and a little bit more error-prone: `[twod(x,y) for x in xs, y in ys]`.

I then found this thread

> [@Outer broadcasting?](https://discourse.julialang.org/t/outer-broadcasting/22701):
>
> 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 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: 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\> f(a, b, c) = a+b+c f (generic function …

and came up with

```julia
twod(x,y) = "$(x) & $(y)" # 2D func
twod_tuple(tpl) = twod(tpl[1], tpl[2])
xs = 0.0:1.0:10.0 # x grid
ys = 20.0:3.0:50.0 # y grid
twod_gridded = twod_tuple.(Iterators.product(xs,ys))

```

Is this the most succinct solution today? Is there a standard way to convert a multiple-argument function into a single-tuple-argument function?

* * *

By the way, I’ve never been able to remember whether `x in xs, y in ys` or `y in ys, x in xs` is the correct order. That’s one of the reasons why I want to avoid the list comprehension to construct a multidimensional array. For the regular `for` loop,

```julia
for j in axes(arr,2)
   for i in axes(arr,1)
      arr[i,j] = func_of(i,j)

```

is the right order because `i` should change faster. The above is (functionally, at least) equivalent to

```julia
for j in axes(arr,2), i in axes(arr,1)
      arr[i,j] = func_of(i,j)

```

So far so good. But then, you have to flip the order for the list comprehension

```julia
arr = [func_of(i,j) for i in is, j in js]

```

_if I’m not mistaken_.

In contrast, there is no such doubt in `Iterator.product` because it is obviously the standard outer product in linear algebra.

For this reason, I’d probably abolish the regular `for` loop from my code and write it this way:

```julia
for (i,j) in Iterators.product(is,js)
   arr[i,j] = func_of(i,j)

```

After all, a nested `for` loop is an outer product.

Then I hope the above loop can be written like

```julia
for (i,j) in is ⊗ js # outer product

```

---

<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:** [September 11, 2023, 11:11am UTC](https://discourse.julialang.org/t/outer-product-broadcast/103731/2 "2023-09-11T11:11:53Z")

</div>

> [@ryofurue](#):
>
> Is this the most succinct solution today?

You can broadcast over multiple arguments and dimensions (I believe the term actually implies multiple inputs), so you can do

```julia
twod.(x, y') 

```

Notice the `'`. Now the dimensions of length \> 1 will be _broadcast_ (expanded) over the corresponding dimensions in the other argument. So, broadcasting over an Nx1 and 1xM array yields an NxM array. This generalizes to higher dimensions and more inputs.

The broadcasting/dimension expanding behavior is crucial, and distinguishes it from `map`.

---

<div class="post-metadata">

**Author:** ![ryofurue](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ryofurue/32/24531_2.png) [@ryofurue](https://discourse.julialang.org/u/ryofurue)\
**Post date:** [September 11, 2023, 1:30pm UTC](https://discourse.julialang.org/t/outer-product-broadcast/103731/3 "2023-09-11T13:30:05Z")

</div>

> [@DNF](#):
>
> ```julia
> twod.(x, y') 
> 
> ```

That’s indeed what I have been missing! Thanks!

But

> This generalizes to higher dimensions and more inputs.

I’d like to know how.

I’m relatively familiar with matrices and so, the usual linear-algebraic operation **x**  **y** ^T, which in julia is `x y'`, makes perfect sense as an argument to the broadcast function application.

But, beyond 2 dimensions, I’m clueless. In my limited experience with tensors, I’ve seen only this type of notation a\_{i,j,k} = u\_i v\_j w\_k, which is an analogue of the list comprehension.

To apply the `x y'` notation to 3-dimensions, perhaps we would first construct `xy = x y'` and then serialize (1-dimensionalize) `xy` into one dimensional column vector and then `xyz = xy_1d z'` and then somehow recover the first two dimensions. . . .

**[Edit:]** I’ve just been able to implement the idea:

```julia
threed(x,y,z) = "$(x),$(y),$(z)" # 3D func
tuplify(f) = tpl -> f(tpl...)
xs = 1:4 # x grid
ys = 10:10:30 # y grid
zs = 100:100:200 # z grid
gridded1 = (tuplify(threed)).(Iterators.product(xs,ys,zs)) # my original idea.
gridded2 = threed.(reshape(xs, (size(xs,1),1,1) ),
                   reshape(ys, (1, size(ys,1),1)),
                   reshape(zs, (1, 1, size(zs,1))) )

```

(I’ve also improved the way to convert a multi-argument function to a single-tuple function.)  
As in the above, you can use `reshape` to convert the vectors into L × 1 × 1 array, 1 × M × 1 array, and 1 × 1 × N array.

Now, I wonder how to do this succinctly.

---

<div class="post-metadata">

**Author:** ![GunnarFarneback](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gunnarfarneback/32/1827_2.png) [@GunnarFarneback](https://discourse.julialang.org/u/GunnarFarneback)\
**Post date:** [September 11, 2023, 2:25pm UTC](https://discourse.julialang.org/t/outer-product-broadcast/103731/4 "2023-09-11T14:25:19Z")

</div>

Broadcasting automatically extends with trailing singleton dimensions and you can use the colon feature of reshape:

```julia
gridded3 = threed.(xs, reshape(ys, 1, :), reshape(zs, 1, 1, :))

```

Still somewhat too cumbersome for comfort though.

---

<div class="post-metadata">

**Author:** ![bertschi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bertschi/32/33462_2.png) [@bertschi](https://discourse.julialang.org/u/bertschi)\
**Post date:** [September 11, 2023, 3:52pm UTC](https://discourse.julialang.org/t/outer-product-broadcast/103731/5 "2023-09-11T15:52:54Z")

</div>

Your “tuplify” function is already in Base and called `splat`.

---

<div class="post-metadata">

**Author:** ![jar1](https://avatars.discourse-cdn.com/v4/letter/j/c0e974/32.png) [@jar1](https://discourse.julialang.org/u/jar1)\
**Post date:** [September 11, 2023, 4:44pm UTC](https://discourse.julialang.org/t/outer-product-broadcast/103731/6 "2023-09-11T16:44:10Z")

</div>

I like [RectiGrids.jl](https://gitlab.com/aplavin/RectiGrids.jl) for this.

```julia
julia> G = grid(20:30, [:a, :b, :c])
2-dimensional KeyedArray(...) with keys:
↓ 11-element UnitRange{Int64}
→ 3-element Vector{Symbol}
And data, 11×3 RectiGrids.RectiGridArr{Base.OneTo(2), Tuple{Int64, Symbol}, 2, Tuple{Nothing, Nothing}, Tuple{UnitRange{Int64}, Vector{Symbol}}}:
       (:a) (:b) (:c)
 (20) (20, :a) (20, :b) (20, :c)
 (21) (21, :a) (21, :b) (21, :c)
 (22) (22, :a) (22, :b) (22, :c)
 (23) (23, :a) (23, :b) (23, :c)
 (24) (24, :a) (24, :b) (24, :c)
 (25) (25, :a) (25, :b) (25, :c)
 (26) (26, :a) (26, :b) (26, :c)
 (27) (27, :a) (27, :b) (27, :c)
 (28) (28, :a) (28, :b) (28, :c)
 (29) (29, :a) (29, :b) (29, :c)
 (30) (30, :a) (30, :b) (30, :c)

```

---

<div class="post-metadata">

**Author:** ![jishnub](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jishnub/32/33620_2.png) [@jishnub](https://discourse.julialang.org/u/jishnub)\
**Post date:** [September 11, 2023, 5:57pm UTC](https://discourse.julialang.org/t/outer-product-broadcast/103731/7 "2023-09-11T17:57:05Z")

</div>

The syntax may be simplified somewhat by defining a function to perform the reshape:

```julia
julia> ↑(a, ::Val{N}) where {N} = reshape(a, ntuple(_->1, N-1)..., :)
↑ (generic function with 1 method)

julia> xs = 1:2; ys = 10:10:20; zs = 100:100:200;

julia> threed.(xs, ↑(ys, Val(2)), ↑(zs, Val(3)))
2×2×2 Array{String, 3}:
[:, :, 1] =
 "1,10,100" "1,20,100"
 "2,10,100" "2,20,100"

[:, :, 2] =
 "1,10,200" "1,20,200"
 "2,10,200" "2,20,200"

```

---

<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:** [September 11, 2023, 6:32pm UTC](https://discourse.julialang.org/t/outer-product-broadcast/103731/8 "2023-09-11T18:32:57Z")

</div>

You can also mimic numpy’s `newaxis` in Julia, which generalizes to more dimensions:

```julia
const newaxis = [CartesianIndex()]
gridded = threed.(xs, ys[newaxis,:], zs[newaxis,newaxis,:])

```

---

<div class="post-metadata">

**Author:** ![ryofurue](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ryofurue/32/24531_2.png) [@ryofurue](https://discourse.julialang.org/u/ryofurue)\
**Post date:** [September 13, 2023, 7:24am UTC](https://discourse.julialang.org/t/outer-product-broadcast/103731/9 "2023-09-13T07:24:41Z")

</div>

> [@GunnarFarneback](#):
>
> Broadcasting automatically extends with trailing singleton dimensions

Wow! That has made it perfectly clear to me what broadcast is. Before reading the quoted sentence above, it had been quite nebulous to me what broadcast is. Thank you for the great explanation!

> [@GunnarFarneback](#):
>
> ```julia
> gridded3 = threed.(xs, reshape(ys, 1, :), reshape(zs, 1, 1, :))
> 
> ```

Aha! Probably that’s the best method using only core features (broadcast and `reshape`). I didn’t know `reshape` has the same convenience as `view`. (Speaking of that, can we replace `reshape` with `view` in the above? . . . Never mind, I’ll investigate.)

I guess we should advocate this. I’ve been kind of surprised that I keep seeing examples on the Internet like

```julia
gridded = zeros(nx, ny, nz) # Off topic, but I guess the programmer wrote >>
# >> this because Array{Float64}(undef, nx, ny, nz) is cumbersome to write.
for k in 1:nz
  for j in 1:ny
    for i in 1:nx
      gridded[i,j,k] = threed(xs[i], ys[j], zs[k])

```

I’m not saying that explicit looping is necessarily bad. A lot of algorithms are more simply written with indices. But the present example doesn’t belong to that category.