# Fastest way to index matrices

**URL:** <https://discourse.julialang.org/t/fastest-way-to-index-matrices/21429>\
**Category:** Performance\
**Created:** [March 3, 2019, 6:12pm UTC](https://discourse.julialang.org/t/fastest-way-to-index-matrices/21429 "2019-03-03T18:12:16Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![Donald\_Williams](https://avatars.discourse-cdn.com/v4/letter/d/f07891/32.png) [@Donald\_Williams](https://discourse.julialang.org/u/Donald_Williams)\
**Post date:** [March 3, 2019, 6:12pm UTC](https://discourse.julialang.org/t/fastest-way-to-index-matrices/21429/1 "2019-03-03T18:12:16Z")

</div>

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.

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

```

---

<div class="post-metadata">

**Author:** ![DNF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dnf/32/10191_2.png) [@DNF](https://discourse.julialang.org/u/DNF)\
**Post date:** [March 3, 2019, 8:25pm UTC](https://discourse.julialang.org/t/fastest-way-to-index-matrices/21429/2 "2019-03-03T20:25:21Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Donald\_Williams](https://avatars.discourse-cdn.com/v4/letter/d/f07891/32.png) [@Donald\_Williams](https://discourse.julialang.org/u/Donald_Williams)\
**Post date:** [March 3, 2019, 8:29pm UTC](https://discourse.julialang.org/t/fastest-way-to-index-matrices/21429/3 "2019-03-03T20:29:07Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![DNF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dnf/32/10191_2.png) [@DNF](https://discourse.julialang.org/u/DNF)\
**Post date:** [March 3, 2019, 8:30pm UTC](https://discourse.julialang.org/t/fastest-way-to-index-matrices/21429/4 "2019-03-03T20:30:32Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![Donald\_Williams](https://avatars.discourse-cdn.com/v4/letter/d/f07891/32.png) [@Donald\_Williams](https://discourse.julialang.org/u/Donald_Williams)\
**Post date:** [March 3, 2019, 8:36pm UTC](https://discourse.julialang.org/t/fastest-way-to-index-matrices/21429/5 "2019-03-03T20:36:48Z")

</div>

See here

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

```

---

<div class="post-metadata">

**Author:** ![DNF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dnf/32/10191_2.png) [@DNF](https://discourse.julialang.org/u/DNF)\
**Post date:** [March 3, 2019, 8:47pm UTC](https://discourse.julialang.org/t/fastest-way-to-index-matrices/21429/6 "2019-03-03T20:47:15Z")

</div>

> [@Donald\_Williams](#):
>
> transpose(cor\_mat[i , idx .!= i])

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

---

<div class="post-metadata">

**Author:** ![Donald\_Williams](https://avatars.discourse-cdn.com/v4/letter/d/f07891/32.png) [@Donald\_Williams](https://discourse.julialang.org/u/Donald_Williams)\
**Post date:** [March 3, 2019, 8:49pm UTC](https://discourse.julialang.org/t/fastest-way-to-index-matrices/21429/7 "2019-03-03T20:49:19Z")

</div>

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