NumPy arrays iterate over the left-most axis first. Thus if B has shape (2,3,4), then B[0] has shape (3,4) and B[1] has shape (3,4). In this sense, you could think of B as 2 arrays of shape (3,4). You can sort of see the two arrays in the repr of B:

In [233]: B = np.arange(2*3*4).reshape((2,3,4))
array([[[ 0,  1,  2,  3],
        [ 4,  5,  6,  7],       <-- first (3,4) array 
        [ 8,  9, 10, 11]],

       [[12, 13, 14, 15],
        [16, 17, 18, 19],      <-- second (3,4) array 
        [20, 21, 22, 23]]])

You can also think of B as containing four 2x3 arrays by iterating over the last index first:

for i in range(4):
    print(B[:,:,i])

# [[ 0  4  8]
#  [12 16 20]]
# [[ 1  5  9]
#  [13 17 21]]
# [[ 2  6 10]
#  [14 18 22]]
# [[ 3  7 11]
#  [15 19 23]]

but you could just as easily think of B as three 2x4 arrays:

for i in range(3):
    print(B[:,i,:])

# [[ 0  1  2  3]
#  [12 13 14 15]]
# [[ 4  5  6  7]
#  [16 17 18 19]]
# [[ 8  9 10 11]
#  [20 21 22 23]]

NumPy arrays are completely flexible this way. But as far as the repr of B is concerned, what you see corresponds to two (3x4) arrays since B iterates over the left-most axis first.

for arr in B:
    print(arr)

# [[ 0  1  2  3]
#  [ 4  5  6  7]
#  [ 8  9 10 11]]
# [[12 13 14 15]
#  [16 17 18 19]
#  [20 21 22 23]]

arr[:][123][:] is processed piece by piece, not as a whole.

arr[:]  # just a copy of `arr`; it is still 3d
arr[123]  # select the 123 item along the first dimension
# oops, 1st dim is only 31, hence the error.

arr[:, 123, :] is processed as whole expression. Select one 'item' along the middle axis, and return everything along the other two.

arr[12][123] would work, because that first select one 2d array from the 1st axis. Now [123] works on that 285 length dimension, returning a 1d array. Repeated indexing [][].. works just often enough to confuse new programmrs, but usually it is not the right expression.

This computation is not slow because of Numpy, but because of your hardware, and actually on most modern hardware. The main reason it is slow is because of random RAM accesses resulting in a latency-bound computation.

The input array is big and thus cannot be stored in CPU cache but in RAM. Modern RAM can have a quite high throughput but each fetch require a pretty big latency (about 80ns per random fetch on recent x86 processors). While the throughput of new RAM devices tends to significantly improve over time, the it is barely the case for latency. Fetching one (8-bytes) double-precision floating-point number at a time sequentially would result in a throughput of size_of_double/RAM_latency = 8/80e-9 ≈ 95 MiB/s. This is a fraction of what modern mainstream RAM devices can do (dozens of GiB/s)

To solve this problem, modern processors can fetch several memory block at a time and try to predict RAM accesses and load data ahead of time using prefetching units and speculation. This works well for predictible access patterns (especially contiguous loads) but not at all with a random access pattern like in your code. Still modern processors succeed to fetch multiple memory blocks in parallel with on a sequential code, but not enough so this kind of code can be fast enough (the throughput of your code is about 400 MiB/s). Not to mention mainstream x86 processors load systematically cache-lines of 64-bytes from RAM devices while you only need 8-bytes.

One solution is to parallelize this operation. However, this is not very efficient because of cache-lines (you will barely get more than 10% of the maximal throughput).

An alternative solution is to transpose your input data so that so that the fetched memory blocks can be more contiguous. Here is an example:

transposedLookup = np.array(lookup.T)
%timeit transposedLookup[m[:,0], m[:,1],:].T

Note that the first transposition will be rather slow (mainly because it is not yet optimized by Numpy, but also because of the access pattern) and requires twice the amount of RAM. You can use this solution to speed up the transposition. It is also possible to transpose the input matrix in-place if it is cubic. It would be even better if you could generate the array directly in its transposed form.

Also note that the second transposition is fast because it is done lazily, but even an eagerly transposition is still several times faster than the original code.

Here are the timings on my machine:

Original code: 14.1 ms
Eager Numpy transposition: 2.6 ms
Lazy Numpy transposition: 0.6 ms

EDIT: note that the same thing apply for m.

Maybe you figured out the answer but anyway here is my answer:

In 3d array the slicing syntax is denoted by [matrices, rows, columns]

import numpy as np

# Below we are creating a 3D matrix similar to your matrix where there are 
# 5 matrices each containing 3 rows and 4 columns

routine_matrix = np.arange(5*3*4).reshape((5,3,4))
print(routine_matrix) 

routine_matrix:

[[[ 0  1  2  3]
  [ 4  5  6  7]
  [ 8  9 10 11]]

 [[12 13 14 15]
  [16 17 18 19]
  [20 21 22 23]]

 [[24 25 26 27]
  [28 29 30 31]
  [32 33 34 35]]

 [[36 37 38 39]
  [40 41 42 43]
  [44 45 46 47]]

 [[48 49 50 51]
  [52 53 54 55]
  [56 57 58 59]]]

Now, As per your question you are interested in 1, 13, 25, 37, 49 that is first element of 0th row and 1st column of each internal matrix

So, in order to achieve that we do

print(routine_matrix[:, 0, 1])

Understand slicing here:

• ':' says that choose all the matrices
• '0' says that choose only the 0th row of all the matrices

i.e: [0 1 2 3], [12, 13, 14, 15], [24, 25, 26, 27], [36, 37, 38, 39], [48, 49, 50, 51]

• '1' says that of all the rows selected above just choose the first column (Array index starts from 0 in python)

i.e: [1, 13, 25, 37, 49]

Hence the output:

[ 1 13 25 37 49]

In your case it will be ['aa' 'bb' 'cc' 'dd' 'ee']

arr = np.array([[[1,2,3], [4,5,6]], [[7,8,9],[10,11,12]]])

# array([[[ 1,  2,  3],
#         [ 4,  5,  6]],

#        [[ 7,  8,  9],
#         [10, 11, 12]]])

arr.shape # ------> (2,2,3)
# Think of them as axis
# lets create the first 2 axis of (`2`, ...)

#         |(1)
#         |
#         |
#         |         
#         |---------(0)


# now lets create second 2 axis of (2, `2`, ..)

#            (1,1)
#         |
#         |---(1,0)
#         |
#         |
#         |
#         |         |(0,1)
#         |---------|---(0,0)

# now lets create the last 3 axis of (2, 2, `3`)

#           /``(1,1,0) = 10
#          |-- (1,1,1) = 11
#         | \__(1,1,2) = 12
#         |
#         |  /``(1,0,0) = 7
#         |--|--(1,0,1) = 8
#         |  \__(1,0,2) = 9
#         |
#         |
#         |         /``(0,1,0) = 4
#         |         |--(0,1,1) = 5
#         |         \__(0,1,2) = 6
#         |         |
#         |         |
#         |---------|---/``(0,0,0) = 1
#                       |--(0,0,1) = 2
#                       \__(0,0,2) = 3

# now suppose you ask
arr[0, :, :] # give me the first axis alon with all it's component

#         |
#         |         /``(0,1,0) = 4
#         |         |--(0,1,1) = 5
#         |         \__(0,1,2) = 6
#         |         |
#         |         |
#         |---------|---/``(0,0,0) = 1
#                       |--(0,0,1) = 2
#                       \__(0,0,2) = 3

# So it will print 

# array([[1, 2, 3],
#        [4, 5, 6]])

arr[:, 0, :] # you ask take all the first axis 1ut give me only the first axis of the first axis and all its components

#           
#         
#         
#         
#         |  /``(1,0,0) = 7
#         |--|--(1,0,1) = 8
#         |  \__(1,0,2) = 9
#         |
#         |
#         |         
#         |         
#         |         
#         |         
#         |         
#         |---------|---/``(0,0,0) = 1
#                       |--(0,0,1) = 2
#                       \__(0,0,2) = 3

# so you get the output

# array([[1, 2, 3],
#        [7, 8, 9]])

# like wise you ask
print(arr[0:2, : , 2])
# so you are saying give (0,1) first axis, all of its children and only 3rd (index starts at 0 so 2 means 3) children
# 0:2 means 0 to 2 `excluding` 2; 0:5 means 0,1,2,3,4

#           
#          |
#         | \__(1,1,2) = 12
#         |
#         |  
#         |--
#         |  \__(1,0,2) = 9
#         |
#         |
#         |        
#         |         
#         |         \__(0,1,2) = 6
#         |         |
#         |         |
#         |---------|---/
#                       |
#                       \__(0,0,2) = 3

# so you get

# array([[ 3,  6],
#        [ 9, 12]])