Arithmetics on Array of Arrays

General Usage
1

I have an array of vectors, represented as static arrays (SArray), and want to perform arithmetic operations like addition (vec + vec), scaling (vec * vec) and so on. Obvious ways I tried turn out to behave very unintuitively:

using StaticArrays

a = SVector.([(1, 2, 3), (4, 5, 6)])
b = SVector(1, 2, 3)
a .+ 1  # OK: SArray{Tuple{3},Int64,1,3}[[2, 3, 4], [5, 6, 7]]
a .* 2  # OK: SArray{Tuple{3},Int64,1,3}[[2, 4, 6], [8, 10, 12]]
1 ./ a  # Wrong: LinearAlgebra.Transpose{Float64,SArray{Tuple{3},Float64,1,3}}[[0.0714286 0.142857 0.214286], [0.0519481 0.0649351 0.0779221]]
a .+ b  # DimensionMismatch("arrays could not be broadcast to a common size")
a .+ [b]  # OK: SArray{Tuple{3},Int64,1,3}[[2, 4, 6], [5, 7, 9]]
a .* [b]  # MethodError: no method matching *(::SArray{Tuple{3},Int64,1,3}, ::SArray{Tuple{3},Int64,1,3})
a ./ [b]  # Wrong: SArray{Tuple{3,3},Float64,2,9}[[0.0714286 0.142857 0.214286; 0.142857 0.285714 0.428571; 0.214286 0.428571 0.642857], [0.285714 0.571429 0.857143; 0.357143 0.714286 1.07143; 0.428571 0.857143 1.28571]]

After reading and thinking about how it works, I can understand how these results appeared. But it still a) is far from intuitive for such a common operation, and b) I don’t see any clear way to actually perform elementwise vector arithmetics here.
Am I missing something?

2

Hi. If you complain about “unituitive” behavior, which is clearly somewhat subjective, you should at least specify what result would be intuitive for you, instead of just stating “wrong”.

Nonetheless, let me guess what might be confusing for you.

1 ./ a : Maybe unituitive but certainly not wrong. What’s happening here is that Julia is taking linear algebra seriously and for v a vector computes it’s Moore–Penrose (pseudo)inverse for 1 / v. The dot broadcasts, so you get a vector of two pseudo inverses.

a ./ [b] : Basically the same explanation as above, as it calculates a[1] * (1 / b) and a[2] * (1 / b) if you will, where 1 / b is, again, the pseudo inverse.

a .+ b : This is equivalent to broadcast(+, a, b). The problem here is that both a and b are (incompatible) vectors. What you probably want is broadcast(+, a, Ref(b)) where b is treated as a scalar. The equivalent short form is a .+ Ref(b).

a .* [b] : Julia, again, takes linear algebra seriously. Hence, column vector * column vector isn’t defined. What you (probably) want is row vector * column vector, i.e. v' * v instead of v * v. So transpose.(a) .* [b] might give you the expected result.

Hope that helps.

2 Likes
3

Sorry, I wrote OK/wrong/unintuitive based on the expectations of plain elementwise operations. Probably I didn’t state it clearly enough in the post.

As I say, after reading the docs and playing around I understand what my original lines of code actually do. Your suggestion on Ref(b) indeed works, althrough I’m not quite sure it looks better than [b]. transpose.(a) .* [b] computes dot product - which also can be useful, so thanks for pointing at this. However, I’m looking for elementwise multiplication/division - any idea how to do them?

On a more general note, it’s quite strange to see stuff from linear algebra like matrix multiplication and inverses suddenly appearing in operations on arrays: they don’t always behave like matrices, and can even have any number of dimensions.

4

The issue is that by “elementwise” you actually mean “elementswise elementwise”. If I understand correctly you want to multiple every component of b to every component of every vector in a (here is the “elementwise elementwise”). One straightforward way to get this would for example be

[x .* b for x in a]

Another way is broadcast((x,y)->x .* y, a, Ref(b)) or equivalently ((x,y)->x .* y).(a, Ref(b)).

5

Thank you, this works! And the first way with comprehension even is obvious. However, it looks like often the other two solutions are required, and they certainly look less intuitive. Comprehension will not work e.g. when a and b are both arrays (of the same size), or when we want to keep the type of outer array if it’s not plain Array.

Still, it’s surprisingly different for plain numbers (where a ./ b would work for any dimensions, including scalars) and for 3-vectors for example (where one needs ((x,y)->x .* y).(a, Ref(b)) or ((x,y)->x .* y).(a, b) or ((x,y)->x .* y).(Ref(a), b)).

6

I think it’s awesome :slight_smile: In fact, not being able to multiply a matrix by a vector in the mathematical sense easily is one of the things that always annoyed me in other programming languages (in the context of scientific computing).

7

I don’t quite get the “easily” part - in python you either use @ instead of *, or just use matrix type instead of array.

8

You don’t need a comprehension in this case. If a and b are both arrays (i.e. vectors) of same size you can just do element wise multiplication a .* b.

I see your point, but you can, for example, do @SVector [x .* b for x in a].

9

IIUC, you don’t need a comprehension in this case. If a and b are both arrays (i.e. vectors) of same size you can just do element wise multiplication a .* b .

It doesn’t look like this: a .* a gives MethodError: no method matching *(::SArray{Tuple{3},Int64,1,3}, ::SArray{Tuple{3},Int64,1,3}).

You mean @SVector [x .* b for x in a] ?

Basically yes, but it does require to explicitly state the type of outer container, which doesn’t help writing generic code. And it’s even done differently for SVector (macro) or many others like AxisArray (function, but with extra parameters).

10

I was talking about the b .* b case. Again, your a isn’t a vector but a vector of vectors. Hence, a .* a basically means [a[i] * a[i] for i in eachindex(a)] which fails, because vector * vector isn’t defined in math (see above).

11

Ok, I see. So, if I understand correctly, there are no plans to make elementwise operations less painful and more consistent, at least in the future? I mean b = 2; a .* b scales each vector by 2; b = [2, 3]; a .* b scales them by corresponding consecutive numbers from b; but if I need to scale each component separately, it requires completely different syntax, which further depends on whether the scaling coefficient is the same for each vector or not.

12

Personaly, I don’t agree with your “less painful” and “more consistent” assessment in any case. When you say “elementwise operations”, as I have repeatedly argued in this thread, you actually mean “elementwise elementwise operations”, that is you are talking about scaling every single component of multiple vectors in an array. You don’t need any special syntax to “scale each component (of a simple vector) separately”:

julia> a = rand(1:10, 3)
3-element Array{Int64,1}:
 5
 7
 9

julia> b = collect(1:3)
3-element Array{Int64,1}:
 1
 2
 3

julia> a .* b
3-element Array{Int64,1}:
  5
 14
 27

Most importantly though, I’m not one of Julia’s core developers. Therefore I can’t make any reliable statements about the future of Julia - which I also did not do at any point in this thread.

1 Like
13

I don’t agree with … “more consistent” assessment

For plain numbers as array elements: a .* b works no matter if a and/or b are arrays.
For 2d vectors (pairs of numbers) as array elements: ((x,y)->x .* y).(a, Ref(b)) or ((x,y)->x .* y).(a, b) or ((x,y)->x .* y).(Ref(a), b) is needed.

It is especially confusing, if the types are Vector{Point} or something like that, where Point is a 2d/3d static vector as defined in e.g. GeometryTypes package. At first I was sure it’s a bug in that package, and only when digging deeper could understand the issue.