Fastest way to index matrices

General Usage Performance
1

I am wondering best practices for indexing matrices in Julia. The following seems less than ideal, but struggling to find a good way to do, e.g., cor_mat[ - j, -j], as most examples are cor_mat[ j , j ] online. Note that I have both p, the number of variables, and samples, the number of prospective samples. This is a test of one part of a function that will need need at least 5,000 samples.

using LinearAlgebra

function timing(p , n, samples)
     
     b = zeros(1, p)
     bb = Array{Float64}(b)
     cor_mat =  Matrix{Float64}(I,p,p)
     zstar = Array{Float64}(undef, n, p)
     x = ones(n)
     m = x * bb
     idx = collect(1:1:p)
    

     for j in 1:samples
          for i in 1:p
     mm = transpose(m[:,1])  + transpose(cor_mat[1 , idx .!= 1 ]) *
                          inv( cor_mat[idx .!= 1 , idx .!= 1 ]) *
                          transpose(zstar[:,idx .!= 1] - m[:,idx .!= 1])
   
     ss = cor_mat[1, 1] - transpose(cor_mat[1 , idx .!= 1 ])  *
                          inv( cor_mat[idx .!= 1 , idx .!= 1 ]) *
                          cor_mat[1 , idx .!= 1 ]       
end
end
end


@time timing(20, 200, 5000)
15.637497 seconds (10.35 M allocations: 14.218 GiB, 8.29% gc time)
2

Do you think you can make your MWE (minimal working example) more minimal?

It’s a bit difficult to know what to say about the code, since it is doing so much redundant stuff. You can remove most of the code, and the double loop just re-does the same calculation over and over, without changing anything.

Also, strange things, such as instantiating b, then copying it to bb and discarding b; and then you only need bb for creating m, which is just a zero matrix.

It’s hard to know what’s important and what isn’t.

1 Like
3

I am initializing because in the larger functions it needs to be initialized as well, although it would only be initialized one in that case. Let me clean it up a bit.

But all of it will be needed in the end, as that zero matrix will not be zeroes on the next iteration in the larger function.

4

Perhaps you can move the initialization stuff outside the core function, and make sure that the core function only does what is essential and non-reducible?

1 Like
5

See here

   p = 20
   n = 20

     b = zeros(1, p)
     cor_mat =  Matrix{Float64}(I,p,p)
     zstar = Array{Float64}(undef, n, p)
     x = ones(n)
     m = x * b
     idx = collect(1:1:p)


 @time for j in 1:5000

          for i in 1:p
     mm = transpose(m[:,i])  + transpose(cor_mat[i , idx .!= i ]) *
                          inv( cor_mat[idx .!= i , idx .!= i ]) *
                          transpose(zstar[:,idx .!= i] - m[:,idx .!= i])

     ss = cor_mat[i, i] - transpose(cor_mat[i , idx .!= i ])  *
                          inv( cor_mat[idx .!= i , idx .!= i ]) *
                          cor_mat[i , idx .!= i ]
end
end

13.469376 seconds (9.38 M allocations: 6.043 GiB, 12.41% gc time)
6

You’re calculating transpose(cor_mat[i , idx .!= i ]), inv( cor_mat[idx .!= i , idx .!= i ]), etc. twice in each loop. Is this necessary?

Also, you should probably use @views, and in general it is faster and more accurate to do A \ M rather than inv(A) * M.

7

Indeed. I had a hard time defining things in the loop, as my primary language is R.