Better way of replacing Matrix values with another Array values

General Usage Performance
,
1

Dear Julia community,

I want to change some values of an existing Matrix with some values that are on another array, but I find that this part of my code is the bottleneck of what I am doing. Is there a better way to do this?

This would be a MWE.

function store_values!(A::Matrix{T}, B::Array{T, 4}) where {T}
	@inbounds A[:, :] = [
	    B[1,1,1,1] B[1,1,2,2] B[1,1,3,3] B[1,1,1,2] B[1,1,2,3] B[1,1,3,1];
	    B[2,2,1,1] B[2,2,2,2] B[2,2,2,3] B[2,2,1,2] B[2,2,2,3] B[2,2,3,1];
	    B[3,3,1,1] B[3,3,2,2] B[3,3,3,3] B[3,3,1,2] B[3,3,2,3] B[3,3,3,1];
	    B[1,2,1,1] B[1,2,2,2] B[1,2,3,3] B[1,2,1,2] B[1,2,2,3] B[1,2,3,1];
	    B[2,3,1,1] B[2,3,2,2] B[2,3,3,3] B[2,3,1,2] B[2,3,2,3] B[2,3,3,1];
	    B[3,1,1,1] B[3,1,2,2] B[3,1,3,3] B[3,1,1,2] B[3,1,2,3] B[3,1,3,1]
	]
	
	return nothing
end #store_values!

using BenchmarkTools
A = zeros(6, 6)
B = rand(3,3,3,3)
@btime store_values!(A, B) # 1.034 μs (38 allocations: 1.22 KiB)

I have read the docs, but did not find something useful. Please let me know if I am missing something.

Thanks!
Eduardo

2

Maybe not necessarily the prettiest solution, but performance-wise about as fast as it gets:

function store_values2!(A::Matrix{T}, B::Array{T, 4}) where {T}
	@inbounds begin
        A[1,:] .= B[1,1,1,1], B[1,1,2,2], B[1,1,3,3], B[1,1,1,2], B[1,1,2,3], B[1,1,3,1]
	    A[2,:] .= B[2,2,1,1], B[2,2,2,2], B[2,2,2,3], B[2,2,1,2], B[2,2,2,3], B[2,2,3,1]
	    A[3,:] .= B[3,3,1,1], B[3,3,2,2], B[3,3,3,3], B[3,3,1,2], B[3,3,2,3], B[3,3,3,1]
	    A[4,:] .= B[1,2,1,1], B[1,2,2,2], B[1,2,3,3], B[1,2,1,2], B[1,2,2,3], B[1,2,3,1]
	    A[5,:] .= B[2,3,1,1], B[2,3,2,2], B[2,3,3,3], B[2,3,1,2], B[2,3,2,3], B[2,3,3,1]
	    A[6,:] .= B[3,1,1,1], B[3,1,2,2], B[3,1,3,3], B[3,1,1,2], B[3,1,2,3], B[3,1,3,1]
    end
	
	return nothing
end #store_values!

using BenchmarkTools
A = zeros(6, 6)
B = rand(3,3,3,3)
store_values!(A, B)

store_values2!(A2, B)
println(A2 == A) #true
A2 = zeros(6,6)
@btime store_values!($A, $B) #785.106 ns (38 allocations: 1.22 KiB)
@btime store_values2!($A, $B) #39.556 ns (0 allocations: 0 bytes)

Note that the key here is to avoid allocating a new matrix in the function and instead actually inserting everything element-wise. As long as you do that the result should be fast.

3 Likes
3

Writing those out seems like an invitation to have hard-to-find typos. You can make a loop, but why the strange selection of indices of B? (If the last pair were (1,3) then it would be sort-of the upper triangle.)

julia> function store!(A::Array, B::Array)
         size(B) == (3,3,3,3) || throw("wrong size B!")
         size(A) == (6,6) || throw("wrong size A!")
         inds = ((1,1), (2,2), (3,3), (1,2), (2,3), (3,1))
         for j in 1:6
          for i in 1:6
            @inbounds A[i,j] = B[inds[i]..., inds[j]...]
          end
         end
         A
       end;

julia> @btime store_values!($A, $B);
  586.806 ns (38 allocations: 1.22 KiB)

julia> @btime store_values2!($A, $B);
  19.140 ns (0 allocations: 0 bytes)

julia> @btime store!($A, $B);
  13.179 ns (0 allocations: 0 bytes)
6 Likes
4

On my Win10 laptop, Julia 1.6.1, it is more like:

@btime store_values!($A, $B)  # 661 ns (38 allocations: 1.22 KiB)
@btime store_values2!($A, $B) #  57 ns (0 allocations: 0 bytes)
@btime store!($A, $B)         #  18 ns (0 allocations: 0 bytes)
1 Like
5

Thank you! I will keep in mind this.

I was working version to have a loop and less possibility of typos, but I like more your approach and it performs way better too :slight_smile:. Thank you!

And about the indices of B, It is a way to represent a 4th order tensor with symmetries in a 6x6 Matrix (you can see more: Voigt notation).

1 Like
6

While you now have a code thats really fast and does the job, there is probably a library that does what you want, so you do not have to worry much about errors in the implementation. In case you need to upscale your problem in the future in might also be a good idea not to explicitly write out matrix indices.
So in that case, maybe the Tensors.jl package does what you need as well, although I have not tried it myself:
https://juliahub.com/docs/Tensors/F7rKl/1.4.3/man/other_operators/#Tensors.tovoigt

2 Likes
7

Actually, I am using that package for other operations and it is amazing, but in my case after the operations with Tensors.jl I need the output in Voigt’s notation and I found out that using the function tovoigt!(...) was the bottleneck of my code, that is why I was trying to find other ways to do this.

tovoigt!(...) # 246.420 ns (0 allocations: 0 bytes)

Maybe in the future (I am still a newbie), I try to make a general function with the comments that are here and make a ‘change’ request.

2 Likes