Broadcasting in matrix multiplication

General Usage
1

Consider a 1x2 matrix A and a 1d array of 2x2 matrices B:

A = [1 2]
B = [rand(2,2) for _ in 1:3]

Obviously, B[1] is a 2x2 matrix, B[2] and B[3] as well.

Now, I would like to get an array of matrix products [A*B[1], A*B[2], ... A*B[end]]. How can I get it? I am trying the dot in

julia> R = A.*B
3×2 Array{Array{Float64,2},2}:
 [0.460514 0.415326; 0.0842543 0.33096]  [0.921028 0.830652; 0.168509 0.66192]
 [0.908421 0.342417; 0.103733 0.155627]  [1.81684 0.684833; 0.207466 0.311255]
 [0.828274 0.451481; 0.053429 0.810234]  [1.65655 0.902963; 0.106858 1.62047]

but obviously it is not the expected result

julia> R[1]
2×2 Array{Float64,2}:
 0.460514   0.415326
 0.0842543  0.33096

Instead, R[1] (and for other indices as well) is expected to be a 1x2 matrix.

2 Likes
2

You want Ref(A).*B. Ref is a single element container, so broadcast will broadcast over it, not A.

6 Likes
3

Wow, that was superfast :slight_smile: Thanks. Now that I have a solution, I still have to understand it. This Ref stuff is completely new to me. Thanks.

4

It’s conceptually the same as [A].*B You are basically just creating a very low overhead container, so that A is treated as an element rather than a collection.

3 Likes
5

You can also just use (A,) (a one element tuple). I personally like it better than Ref because I thought that was mostly meant for ccall usage.

3 Likes
6

You could also use a comprehension, which is clearer, IMO, or even a generator if you don’t want to materialize.

julia> R = [A*b for b in B]
3-element Array{Array{Float64,2},1}:
 [1.1064518582166465 1.4257048365761784]
 [0.7462570773129014 0.2360058367114488]
 [0.2639640481443961 1.0567511644538756]

julia> R = (A*b for b in B)
Base.Generator{Array{Array{Float64,2},1},var"#11#12"}(var"#11#12"(), [[0.28439915229561796 0.024986306638530964; 0.41102635296051426 0.7003592649688237], [0.6502316144726392 0.15084252545701005; 0.04801273142013107 0.042581655627219384], [0.0710881212956167 0.14551768207448657; 0.0964379634243897 0.45561674118969453]])
3 Likes
7

Excellent, thanks. Indeed, it is clearer.

8

It is indeed clearer. In the case A and each element in B are square matrices, how could you do it inplace?

1 Like
9

B .= [A*b for b in B] works, but does do some allocations

for 2x2 this seems to work

using LinearAlgebra
for b in B
       mul!(b,A,b)
end

but when it gets bigger it attempts calling BLAS which doesn’t like the output being an alias to the input

1 Like
10

I suspect that this may actually be worse than B = [A*b for b in B], because either way it allocates on the right-hand side, but then it additionally copies over to the left hand side instead of just rebinding the variable, which is basically free.

1 Like
11

True, but it probably depends on the relative sizes of B and A and memory pressure etc. One big allocation for a new B, might or might not be relevant compared to the numerous “small” allocations for A.
Also the inplace version might still be relevant if for example B is a view into a larger array. But now I’m probably reaching :slight_smile:
Either way, the request was for an inplace version and I was hoping someone would suggest a clever way to fix the allocations. (The fact that the for loop with the mul! gave different behaviour for different sizes also thought me something new today).

1 Like
12

Actually, my (untested) thesis is that the dotted assignment is strictly worse than the undotted one, because the dotted one does everything that the undotted one does, plus extra copying at the end. I think.

1 Like