Hi, I would like to perform the following element operations to a matrix W (m rows and n columns):
W(i, j) := a*W(i, j) - b*d(j)y(i),
where a and b are scalars, d(j) are the elements of a vector d with n elements, and y(i) are the elements of a vector y with m elements.
In Julia, this can be easily programmed using the following code:
W .= a.*W .- b.*(y*d')
Although the dots prevent some memory allocation, matrix product y*d' allocates memory. (This operation is to be performed multiple time inside a cycle.)
Is there any way to do the whole operation with simple operations and without memory allocation (if there was no matrix product this could be easy done with broadcasting)? Or the only way out is to do two nested for cycles?
I can pre-allocate memory for the matrix product, but I would like to avoid that.
Sure, you just need one more dot:
W = rand(3,4); a,b = rand(2); d = rand(4); y = rand(3);
W2 = similar(W);
f1!(W2, W,a,b,d,y) = W2 .= a.*W .- b.*(y*d')
f2!(W2, W,a,b,d,y) = W2 .= a.*W .- b.*(y .* d')
using Einsum
f3!(W2, W,a,b,d,y) = @einsum W2[i,j] = a * W[i,j] - b * y[i] * d[j]
julia> @btime f1!($W2, $W,$a,$b,$d,$y);
92.883 ns (1 allocation: 176 bytes)
julia> @btime f2!($W2, $W,$a,$b,$d,$y);
37.475 ns (0 allocations: 0 bytes)
julia> @btime f3!($W2, $W,$a,$b,$d,$y);
28.712 ns (0 allocations: 0 bytes)
4 Likes
Wow, beautiful! Thank you very much.
Could one enforce matrix multiplication if there were two squares matrices with the same size? I guess in this case the .* would do element-wise multiplication.
Not sure I follow. On matrices A * B and A .* B are quite different, it’s just for vectors that v * w' and v .* w' happen to agree. But there is mul!, if that’s what you want?
It can be done using mul!, but I will have to use two line of code. Matrix multiplication first and the element-wise operations after. Just thought it would look nicer in a single line of code just with operators.
Thank you very much.