Broadcast product of arrays

New to Julia
1

I would like to call a function with all combination of values from several arrays, and get the corresponding array of values. Is a convenient general way to do that? Broadcasting comes first to my mind, but I wouldn’t call it convenient for this case, e.g.:

array_a = 1:3
array_b = 1:5
array_c = 1:7
f.(array_a, reshape(array_b, 1, length(array_b)), reshape(array_c, 1, 1, length(array_c)))

array_a = 1:3
array_b = [1 2; 3 4]
array_c = 1:7
f.(array_a, reshape(array_b, 1, size(array_b)...), reshape(array_c, 1, 1, 1, length(array_c)))

First it is very verbose and repetitive, and second it is error-prone especially if some of those arrays have the same size. I failed to find a better way, is there any?

2

Consider using Iterators.product:

[f(t...) for t = Iterators.product(array_a, array_b, array_c)]
2 Likes
3

Oh, really! I didn’t even look at that method thinking that Iterators are “single-dimensional” only. I guess that will work fine for me.
If only it could also transparently be applied to cases like f.(scalar, scalar, arr, arr, scalar, a=b, c=d)…

4

I wrote a package for doing such things, in which your second example looks like this:

using TensorCast
@cast res[α,β,β′,γ] := f(array_a[α], array_b[β,β′], array_c[γ])

It should be no problem to give f other scalar arguments.

2 Likes
5

Interesting package indeed, it looks very general!

6

You could just use a comprehension [f(a,b) for a in A, b in B], depending on what shape you want for the result array.

5 Likes
7

A comprehension is the way to go, or maybe Iterators.product. But note also that you’re making reshape unnecessarily verbose:

f.(a, reshape(b, 1, :), reshape(c, 1, 1, :))

You could also use transpose(b) or even b' for real data.

8

Also nice suggestion, thank you - didn’t know about reshape with colon.

However, one needs to be very careful with transpose in Julia: at first I used a and b' for two-dimensional broadcast as you say, and surprisingly got wrong results. After digging in it turned out that transpose is not the same transpose I got used to in other languages, but also performs conjugation. Yes, my arrays were complex in that case, but after that I’m certainly not going to use transpose for broadcasting. Broadcasting syntax should not depend on the array element type, otherwise you can easily get errors like this when some eltype changes.

And anyway, transpose solves only a single simple case of multidimensional broadcasting, it is not of much help for > 2 arrays, or when they are > 1 dimensional (see 1st post for example of both). And it doesn’t work at all with non-numeric types.

9

Yes, I would not use transpose or ' in this case.

Be aware, though, that transpose does not do complex conjugation, while ' does (' is the same as adjoint). Both of them work recursively on their elements, which is important to be aware of, if the elements of your array can themselves be transposed.

10

Ah, ok. Actually approach to transpose in julia feels like an inconvenience for me, and very simple stuff like “swap rows and columns in e.g. a matplotlib plot” cannot be written as ax' like in Python (ax.T there of course), but only as permutedims(ax, 2, 1). However after reading the discussion I see that the majority is happy with the current state, so it will stay like this.

11

Are you sure? You should be able to use ' or transpose in most cases. Are you talking about higher-dimensional arrays? You would have the same problem in Matlab, btw. And in Matlab ' also complex conjugates.

12

Yes, sure - I tried that :slight_smile: Minimal example:

arr = [:a :b; :c :d]
a'  # MethodError: no method matching adjoint(::Symbol)
arr |> transpose  # MethodError: no method matching transpose(::Symbol)

Same goes for any non-numerical array.

13

Wow, that caught me by surprise. To be honest, I personally think that should work. The reason it’s not is that there is no method for transpose(::Symbol). Might be worth opening a discussion thread on that topic specifically.

But how does this look in python?

14

In numpy transpose arr.T just swaps rows and cols of an array, without any regards to its internal structure (eltype), and without modifying elements of course, unlike arr' in julia.

And yes, there is not transpose(::Symbol), as well as no transpose(::anything else but for numerics). And as I understood the corresponding discussion, it is by design.

15

I mean, what does the array look like? I have only ever seen numerical arrays in numpy.

16

Numpy supports any kind of objects in elements - e.g. after fig, ax = plt.subplots(2, 2) ax is a 2x2 numpy array of matplotlib axes.

17

OK, didn’t know that. Anyway, I think it’s strange that transpose doesn’t work on non-numeric matrices. I would ask about it specifically if I were you, unless you’ve seen somewhere that it is definitely intensional.

Anyway, note that for matrices, you don’t need to specify dimension to swap, since there are only two:

julia> arr = [:a :b; :c :d]
2×2 Array{Symbol,2}:
 :a  :b
 :c  :d

julia> permutedims(arr)
2×2 Array{Symbol,2}:
 :a  :c
 :b  :d
18

Cool, didn’t know about that as well (permutedims without explicit dims)!
And as for transpose, it is definitely intensional:

  • permutedims(m::AbstractMatrix) is now short for permutedims(m, (2,1)) , and is now a more convenient way of making a “shallow transpose” of a 2D array. This is the recommended approach for manipulating arrays of data, rather than the recursively defined, linear-algebra function transpose . Similarly, permutedims(v::AbstractVector) will create a row matrix (#24839).
19

Note that for a matrix, you don’t need to specify dimensions, just permutedims(ax) is fine. See ?permutedims.

20

Note that transpose(mat) is lazy, while permutedims(mat) makes a copy. There is also PermutedDimsArray(mat,(2,1)) which is lazy and not recursive.

(It would probably make sense to have PermutedDimsArray(mat) too, for matrices, but currently this does not exist.)

1 Like