Outer product & broadcast?

General Usage
1

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

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

and came up with

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,

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

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

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:

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

for (i,j) in is ⊗ js # outer product
1 Like
2

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

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.

4 Likes
3

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:

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.

4

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

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

Still somewhat too cumbersome for comfort though.

1 Like
5

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

1 Like
6

I like RectiGrids.jl for this.

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)
7

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

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"
8

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

const newaxis = [CartesianIndex()]
gridded = threed.(xs, ys[newaxis,:], zs[newaxis,newaxis,:])
2 Likes
9

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!

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

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.