General options to avoid nested for-loops for 3d matrix operations.

Question

Following my research about nested loops I'd like to put up up this general question: I have to run various functions on 3d-matrices (multi band image files) with a shape of e.g. (2500, 100, 200) or (bands, rows, cols).

The general (slow) procedure how I approach processing functions looks like this:

import numpy as np

import time

x = np.random.random_sample([2500, 100, 200])

#  1. & 2. derivative example

def derivate(matrix):

    d1 = np.empty(shape=np.shape(matrix))

    d2 = np.empty(shape=np.shape(matrix))

    for row in range(matrix.shape[1]):

        for col in range(matrix.shape[2]):

            d1[:, row, col] = np.gradient(matrix[:, row, col])

            d2[:, row, col] = np.gradient(d1[:, row, col])

    return d1, d2

t1 = time.perf_counter()

y, z = derivate(x)

t2 = time.perf_counter()

print(t2 - t1)

>>> 4.38 seconds

Since many functions in this style already exist in my project environment I am looking for the easiest solution to drastically speed up those nested for-loop functions.

I have seen that many questions here address nested for-loops but I didn't discover any general and transferable solution. I would very much appreciate if many available approaches were presented, cheers!

Answer 1

You can just use numpy's built in derivate_vec method like this

import numpy as np

def derivate_vec(matrix):

    d1 = np.gradient(matrix, axis=0)

    d2 = np.gradient(d1, axis=0)

    return d1, d2