# Arithmetics on Array of Arrays

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?

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

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.

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

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

I think itâ€™s awesome 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).

I donâ€™t quite get the â€śeasilyâ€ť part - in python you either use `@` instead of `*`, or just use `matrix` type instead of `array`.

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]`.

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

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

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.

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

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.