Consider this:

a = array(1:60, c(3,4,5))

# Extract (*,2,2)^th elements
a[cbind(1:3,2,2)]
                                
# ^^ returns a vector of length 3, c(16, 17, 18)
# Note this is asymptotically inefficient -- cbind(1:k,2,2) is a kx3 matrix
# So it's not making use of slicing
 
# Extract (1, y_i, z_i)^th elements
N = 100000
y = sample(1:4, N, replace = TRUE)
z = sample(1:5, N, replace = TRUE)
a[cbind(1, y, z)]

# ^^ returns a vector of length N

How do I efficiently extract the (*, y_i, z_i) elements, with the result as a Nx3 matrix (i.e. 2d array)?

Note this question is similar but the only answer proceeds elementwise + so will be slow.