# Array indexing/slicing is slow

**URL:** <https://discourse.julialang.org/t/array-indexing-slicing-is-slow/97633>\
**Category:** General Usage\
**Tags:** array\
**Created:** [April 18, 2023, 9:39pm UTC](https://discourse.julialang.org/t/array-indexing-slicing-is-slow/97633 "2023-04-18T21:39:37Z")\
**Posts on this page:** 10\
**Page:** 1

<div class="post-metadata">

**Author:** ![danicaratelli](https://avatars.discourse-cdn.com/v4/letter/d/90db22/32.png) [@danicaratelli](https://discourse.julialang.org/u/danicaratelli)\
**Post date:** [April 18, 2023, 9:39pm UTC](https://discourse.julialang.org/t/array-indexing-slicing-is-slow/97633/1 "2023-04-18T21:39:37Z")

</div>

I am doing very simple calculations in loops. It looks like the array extraction is the most expensive part of this. Is there an obvious way to speed it up? Here is a MWE of what I am doing. Any help would be greatly appreciated!

`Is1` and `Is2` are matrices of indices, and `Pis1` and `Pis2` are matrices of probabilities. I compute the following

```julia
Threads.@threads for iA = 1:Nstates            
    @simd for iy = 1:nY
        @simd for iP = 1:nP
            @simd for ip = 1:np
                # relevant continuation value
                I1, PI1, I2, PI2 = Is1[iy, iP, ip, iA], Pis1[iy, iP, ip, iA], Is2[iy, iP, ip, iA], Pis2[iy, iP, ip, iA]
                Xcont = compute_cont(X, I1, PI1, I2, PI2 );
                # then I use Xcont in other computations but the above line is what is slow
            end
        end
    end
end

function compute_cont(X, I1, I2, PI1, PI2)
    I1_next, I2_next = min(nY, I1+1), min(nP, I2+1);
    
    output = @views X[:, :, I1, I2, :] * PI1 * PI2 + X[:, :, I1_next, I2, :] * (1 - PI1) * PI2 
                    + X[:, :, I1, I2_next, :] * PI1 * (1 - PI2) + X[:, :, I1_next, I2_next, :] * (1 - PI1) * (1 - PI2);
    
    return output
end

```

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [April 18, 2023, 10:07pm UTC](https://discourse.julialang.org/t/array-indexing-slicing-is-slow/97633/2 "2023-04-18T22:07:54Z")

</div>

> [@danicaratelli](#):
>
> ```julia
> output = @views X[:, :, I1, I2, :] * PI1 * PI2 + X[:, :, I1_next, I2, :] * (1 - PI1) * PI2 
> + X[:, :, I1, I2_next, :] * PI1 * (1 - PI2) + X[:, :, I1_next, I2_next, :] * (1 - PI1) * (1 - PI2);
> 
> ```

Every time you execute this line, you are allocating a dozen temporary arrays. Read about [fusing array operations in Julia](https://julialang.org/blog/2017/01/moredots/), and write into a [pre-allocated](https://docs.julialang.org/en/v1/manual/performance-tips/#Pre-allocating-outputs) `Xcont` array.

Also, put all of your [loops into functions](https://docs.julialang.org/en/v1/manual/performance-tips/#Performance-critical-code-should-be-inside-a-function).

(Also think about your dimension ordering — ideally, order your arrays so that you are doing contiguous slices and other contiguous access, i.e. [use the first dimensions for slices and inner loops](https://docs.julialang.org/en/v1/manual/performance-tips/#man-performance-column-major).)

---

<div class="post-metadata">

**Author:** ![danicaratelli](https://avatars.discourse-cdn.com/v4/letter/d/90db22/32.png) [@danicaratelli](https://discourse.julialang.org/u/danicaratelli)\
**Post date:** [April 18, 2023, 11:11pm UTC](https://discourse.julialang.org/t/array-indexing-slicing-is-slow/97633/3 "2023-04-18T23:11:44Z")

</div>

Thank you! Just a couple of points:

- what is pre-allocating `Xcont` gives an out-of-memory error?

- I try re-defining the `compute_cont` function as follows (which I think should not reallocate the temporary arrays) but it actually takes **much** longer… Maybe I am “fusing” the arrays incorrectly.

```julia
function compute_cont(X, I1, I2, PI1, PI2)
    I1_next, I2_next = min(nY, I1+1), min(nP, I2+1);
    output = zeros(size(X, 1), size(X, 2), size(X, 5));
    for i=1:size(X, 1)
        for j=1:size(X, 2)
            for k=1:size(X, 5)
                output = @views X[i, j, I1, I2, k] * PI1 * PI2 + X[i, j, I1_next, I2, k] * (1 - PI1) * PI2 
                                + X[i, j, I1, I2_next, k] * PI1 * (1 - PI2) + X[i, j, I1_next, I2_next, k] * (1 - PI1) * (1 - PI2);
            end
        end
    end
    
    return output
end

```

---

<div class="post-metadata">

**Author:** ![Jeff\_Emanuel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jeff_emanuel/32/15440_2.png) [@Jeff\_Emanuel](https://discourse.julialang.org/u/Jeff_Emanuel)\
**Post date:** [April 18, 2023, 11:37pm UTC](https://discourse.julialang.org/t/array-indexing-slicing-is-slow/97633/4 "2023-04-18T23:37:58Z")

</div>

> [@danicaratelli](#):
>
> ` output =`

```julia
  output[i, j, k] = 

```

---

<div class="post-metadata">

**Author:** ![danicaratelli](https://avatars.discourse-cdn.com/v4/letter/d/90db22/32.png) [@danicaratelli](https://discourse.julialang.org/u/danicaratelli)\
**Post date:** [April 18, 2023, 11:43pm UTC](https://discourse.julialang.org/t/array-indexing-slicing-is-slow/97633/5 "2023-04-18T23:43:28Z")

</div>

Yes, sorry: I have that in the code (just a mistake reporting it here).

---

<div class="post-metadata">

**Author:** ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)\
**Post date:** [April 19, 2023, 12:16am UTC](https://discourse.julialang.org/t/array-indexing-slicing-is-slow/97633/6 "2023-04-19T00:16:04Z")

</div>

> [@danicaratelli](#):
>
> `output = zeros(size(X...`

Allocate this outside this function, and pass it as a parameter. Allocate just one per thread.

(Also invert the order of the loops, run over i first, then j, then k)

---

<div class="post-metadata">

**Author:** ![roflmaostc](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/roflmaostc/32/30123_2.png) [@roflmaostc](https://discourse.julialang.org/u/roflmaostc)\
**Post date:** [April 19, 2023, 8:43am UTC](https://discourse.julialang.org/t/array-indexing-slicing-is-slow/97633/7 "2023-04-19T08:43:14Z")

</div>

can you provide a minimal working example (MWE)?

Like providing the correct datatypes and shapes for `X, I1, ...`

---

<div class="post-metadata">

**Author:** ![danicaratelli](https://avatars.discourse-cdn.com/v4/letter/d/90db22/32.png) [@danicaratelli](https://discourse.julialang.org/u/danicaratelli)\
**Post date:** [April 19, 2023, 1:11pm UTC](https://discourse.julialang.org/t/array-indexing-slicing-is-slow/97633/8 "2023-04-19T13:11:47Z")

</div>

Hello,

Yes: this should do it. Thank you in advance for your help!

```julia
Is1 = rand(1:3,3,3,3,5);
Pis1 = rand(3,3,3,5);
Is2 = rand(1:3,3,3,3,5);
Pis2 = rand(3,3,3,5);

X = rand(3,3,3,3,5);

function compute_cont(X, I1, I2, PI1, PI2)
    I1_next, I2_next = min(3, I1+1), min(3, I2+1);
    
    output = @views X[:, :, I1, I2, :] * PI1 * PI2 + X[:, :, I1_next, I2, :] * (1 - PI1) * PI2 + X[:, :, I1, I2_next, :] * PI1 * (1 - PI2) + X[:, :, I1_next, I2_next, :] * (1 - PI1) * (1 - PI2);
    
    return output
end

for i = 1:3            
    for j = 1:3
        for k = 1:3
            for w = 1:5
                # relevant continuation value
                I1, PI1, I2, PI2 = Is1[i,j,k,w], Pis1[i,j,k,w], Is2[i,j,k,w], Pis2[i,j,k,w];
                Xcont = compute_cont(X, I1, I2, PI1, PI2);  
                # some other operations occur on Xcont              
            end
        end
    end
end

```

---

<div class="post-metadata">

**Author:** ![Dan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dan/32/42581_2.png) [@Dan](https://discourse.julialang.org/u/Dan)\
**Post date:** [April 19, 2023, 2:05pm UTC](https://discourse.julialang.org/t/array-indexing-slicing-is-slow/97633/9 "2023-04-19T14:05:24Z")

</div>

Some rewrite + minimal optimization to better conform to must-read: [Performance Tips · The Julia Language](https://docs.julialang.org/en/v1/manual/performance-tips/#man-performance-tips)

```julia
Is1 = rand(1:3,3,3,3,5);
Pis1 = rand(3,3,3,5);
Is2 = rand(1:3,3,3,3,5);
Pis2 = rand(3,3,3,5);

X = rand(3,3,3,3,5);

function compute_cont!(Xcont, X, I1, I2, PI1, PI2)
    I1_next, I2_next = min(3, I1+1), min(3, I2+1);
    
    @. Xcont = @views X[:, :, I1, I2, :] * PI1 * PI2 + X[:, :, I1_next, I2, :] * (1 - PI1) * PI2 + X[:, :, I1, I2_next, :] * PI1 * (1 - PI2) + X[:, :, I1_next, I2_next, :] * (1 - PI1) * (1 - PI2);
    
    return Xcont
end

function do_thing(X, Is1, Pis1, Is2,Pis2)
    Xcont = similar(@view X[:,:,1,1,:])
    @inbounds for I in CartesianIndices(Is1)
        # relevant continuation value
        I1, PI1, I2, PI2 = Is1[I], Pis1[I], Is2[I], Pis2[I];
        compute_cont!(Xcont, X, I1, I2, PI1, PI2);  
        # some other operations occur on Xcont              
    end
    return Xcont
end

using BenchmarkTools
@btime do_thing(X, Is1, Pis1, Is2, Pis2);

```

---

<div class="post-metadata">

**Author:** ![mihalybaci](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mihalybaci/32/13528_2.png) [@mihalybaci](https://discourse.julialang.org/u/mihalybaci)\
**Post date:** [April 19, 2023, 2:07pm UTC](https://discourse.julialang.org/t/array-indexing-slicing-is-slow/97633/10 "2023-04-19T14:07:27Z")

</div>

Echoing what @stevengj said, make sure to put the loops intto a function and fuse the loops. Simply changing

```julia
output = @views X....

```

to

```julia
output = @. @views X....

```

took the runtime on my computer from 90 ms to 21 ms, so a nearly 5x speed up by adding two characters.

(@Dan beat me to it and his in-place version is probably even faster)