Оптимизация Python: Как ускорить обратную операцию матрицы?

Question

Мой код содержит цикл for с большим количеством итераций.Внутри цикла мне нужно, чтобы матрица умножалась и матрица была обратной (обычно это матрица размером 12 х 12).Мой цикл должен работать 120 000 раз, и в настоящее время я получаю скорость 14 с, что является относительно высоким показателем по сравнению с MATLAB (1 с) и FORTRAN (0,4 с).Ниже приведена функция, которую я пытаюсь оптимизировать:

def fixed_simulator(ref, xg, yg, dt, ndiv, ijk, nst, smx, skx, cdx, smy, sky, cdy):

    gamma = 0.5
    beta = 1.0/6.0

    knx = skx + (gamma/beta/dt)*cdx + (1.0/beta/np.power(dt,2))*smx

    dx1 = np.ones((nst,1), dtype=float)*0.0
    vx1 = np.ones((nst,1), dtype=float)*0.0
    px1 = np.diag(-1.0*smx).reshape(nst,1)*xg[0]
    ax1 = np.matmul(linalg.inv(smx), px1 - np.matmul(cdx, vx1) - np.matmul(skx, dx1))

    # I = np.ones((nst,1), dtype=float)

    dx2 = np.zeros((nst,1), dtype=float)
    vx2 = np.zeros((nst,1), dtype=float)
    px2 = np.zeros((nst,1), dtype=float)
    ax2 = np.zeros((nst,1), dtype=float)

    na1x = (1.0/beta/np.power(dt,2))*smx + (gamma/beta/dt)*cdx
    na2x = (1.0/beta/dt)*smx + (gamma/beta - 1.0)*cdx
    na3x = (1.0/2.0/beta - 1.0)*smx + (gamma*dt/2.0/beta - dt)*cdx

    print(len(xg))

# -----> Below is the loop that's taking long time.  

    for i in range(1,len(xg)):

        px2 = np.diag(-1.0*smx).reshape(nst,1)*xg[i]

        pcx1 = px2 + np.matmul(na1x, dx1) + np.matmul(na2x, vx1) + np.matmul(na3x, ax1)

        dx2 =  np.matmul(np.linalg.inv(smx), pcx1)

        vx2 = (gamma/beta/dt)*(dx2 - dx1) + (1.0 - gamma/beta)*vx1 + dt*(1.0 - gamma/2.0/beta)*ax1

        ax2 = np.matmul(np.linalg.inv(smx), px2 - np.matmul(cdx, vx2) - np.matmul(skx, dx2))

        dx1, vx1, px1, ax1 = dx2, vx2, px2, ax2 

Большую часть времени занимает вычисление обратной и умножающей частей.

Конфигурация Numpy в системе:

blas_mkl_info:
  NOT AVAILABLE
blis_info:
  NOT AVAILABLE
openblas_info:
    library_dirs = ['C:\\projects\\numpy-wheels\\numpy\\build\\openblas']
    libraries = ['openblas']
    language = f77
    define_macros = [('HAVE_CBLAS', None)]
blas_opt_info:
    library_dirs = ['C:\\projects\\numpy-wheels\\numpy\\build\\openblas']
    libraries = ['openblas']
    language = f77
    define_macros = [('HAVE_CBLAS', None)]
lapack_mkl_info:
  NOT AVAILABLE
openblas_lapack_info:
    library_dirs = ['C:\\projects\\numpy-wheels\\numpy\\build\\openblas']
    libraries = ['openblas']
    language = f77
    define_macros = [('HAVE_CBLAS', None)]
lapack_opt_info:
    library_dirs = ['C:\\projects\\numpy-wheels\\numpy\\build\\openblas']
    libraries = ['openblas']
    language = f77
    define_macros = [('HAVE_CBLAS', None)]

cProfile Результат

         2157895 function calls in 2.519 seconds

   Ordered by: cumulative time

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
        1    1.474    1.474    2.519    2.519 C:\Users\Naseef\OneDrive\04AllPhDPrograms\mhps\mhps\fixed.py:154(fixed_simulator)
   839163    0.556    0.000    0.556    0.000 {built-in method numpy.core.multiarray.matmul}
   119881    0.105    0.000    0.439    0.000 C:\Users\Naseef\AppData\Local\Programs\Python\Python36-32\lib\site-packages\numpy\lib\twodim_base.py:197(diag)
   119881    0.083    0.000    0.256    0.000 C:\Users\Naseef\AppData\Local\Programs\Python\Python36-32\lib\site-packages\numpy\core\fromnumeric.py:1294(diagonal)
   239762    0.049    0.000    0.107    0.000 C:\Users\Naseef\AppData\Local\Programs\Python\Python36-32\lib\site-packages\numpy\core\numeric.py:504(asanyarray)
   119881    0.103    0.000    0.103    0.000 {method 'diagonal' of 'numpy.ndarray' objects}
   239763    0.058    0.000    0.058    0.000 {built-in method numpy.core.multiarray.array}
   119881    0.049    0.000    0.049    0.000 {method 'reshape' of 'numpy.ndarray' objects}
   119881    0.022    0.000    0.022    0.000 {built-in method builtins.isinstance}
   239764    0.019    0.000    0.019    0.000 {built-in method builtins.len}
        1    0.000    0.000    0.000    0.000 C:\Users\Naseef\AppData\Local\Programs\Python\Python36-32\lib\site-packages\numpy\linalg\linalg.py:468(inv)
        2    0.000    0.000    0.000    0.000 C:\Users\Naseef\AppData\Local\Programs\Python\Python36-32\lib\site-packages\numpy\core\numeric.py:156(ones)
        1    0.000    0.000    0.000    0.000 C:\Users\Naseef\AppData\Local\Programs\Python\Python36-32\lib\site-packages\numpy\linalg\linalg.py:141(_commonType)
        2    0.000    0.000    0.000    0.000 {built-in method numpy.core.multiarray.empty}
        1    0.000    0.000    0.000    0.000 {built-in method builtins.print}
        2    0.000    0.000    0.000    0.000 C:\Users\Naseef\AppData\Local\Programs\Python\Python36-32\lib\site-packages\progressbar\utils.py:28(write)
        1    0.000    0.000    0.000    0.000 C:\Users\Naseef\AppData\Local\Programs\Python\Python36-32\lib\site-packages\numpy\linalg\linalg.py:108(_makearray)
        2    0.000    0.000    0.000    0.000 {built-in method numpy.core.multiarray.copyto}
        4    0.000    0.000    0.000    0.000 {built-in method numpy.core.multiarray.zeros}
        2    0.000    0.000    0.000    0.000 C:\Users\Naseef\AppData\Local\Programs\Python\Python36-32\lib\site-packages\progressbar\bar.py:547(update)
        1    0.000    0.000    0.000    0.000 C:\Users\Naseef\AppData\Local\Programs\Python\Python36-32\lib\site-packages\numpy\linalg\linalg.py:126(_realType)
        1    0.000    0.000    0.000    0.000 {method 'astype' of 'numpy.ndarray' objects}
        1    0.000    0.000    0.000    0.000 C:\Users\Naseef\AppData\Local\Programs\Python\Python36-32\lib\site-packages\numpy\core\numeric.py:433(asarray)
        2    0.000    0.000    0.000    0.000 {method 'write' of '_io.StringIO' objects}
        2    0.000    0.000    0.000    0.000 C:\Users\Naseef\AppData\Local\Programs\Python\Python36-32\lib\site-packages\numpy\linalg\linalg.py:113(isComplexType)
        1    0.000    0.000    0.000    0.000 C:\Users\Naseef\AppData\Local\Programs\Python\Python36-32\lib\site-packages\numpy\linalg\linalg.py:103(get_linalg_error_extobj)
        1    0.000    0.000    0.000    0.000 C:\Users\Naseef\AppData\Local\Programs\Python\Python36-32\lib\site-packages\numpy\linalg\linalg.py:200(_assertRankAtLeast2)
        1    0.000    0.000    0.000    0.000 C:\Users\Naseef\AppData\Local\Programs\Python\Python36-32\lib\site-packages\numpy\linalg\linalg.py:211(_assertNdSquareness)
        2    0.000    0.000    0.000    0.000 {built-in method time.perf_counter}
        3    0.000    0.000    0.000    0.000 {built-in method builtins.issubclass}
        1    0.000    0.000    0.000    0.000 {method '__array_prepare__' of 'numpy.ndarray' objects}
        1    0.000    0.000    0.000    0.000 {method 'disable' of '_lsprof.Profiler' objects}
        1    0.000    0.000    0.000    0.000 {method 'get' of 'dict' objects}
        1    0.000    0.000    0.000    0.000 {built-in method builtins.getattr}

Answer 1

Есть два шага, которые я обычно делаю, прежде чем пытаться ускорить простой код.

1. Профилируйте код, чтобы найти то, что занимает больше всего времени
2. Создайте пару тестоввызов кода для проверки правильности работы оптимизированной версии

Контрольные примеры должны быть простыми и быстрыми, но при этом отражать реальные данные.Вы предоставили профилирование, но не тестировали, так что дальнейшие действия будут непроверенными.Просматривая тестовые случаи с наибольшим временем выполнения, становится ясно, что время выполнения исходит из цикла и в основном из матричных операций в нем.

119881 0.049 0.000 0.049 0.000 {метод 'reshape' объектов 'numpy.ndarray'} 239763 0,058 0,000 0,058 0,000 {встроенный метод numpy.core.multiarray.array} 119881 0,103 0,000 0,103 0,000 {метод 'diagonal' из объектов 'numpy.ndarray'} 239762 0,049 0,000 0,107 0,000 ... \ core \ numeric.py: 504 (asanyarray) 119881 0,083 0,000 0,256 0,000 ... \ core \ fromnumeric.py: 1294 (диагональ) 119881 0,105 0,000 0,439 0,000 ... \ lib \ twodim_base.py: 197 (диагональ) 839163 0,556 0,000 0,556 0,000 {встроенный метод numpy.core.multiarray.matmul}

Первая странность в том, что np.linalg.inv(smx) не отображается в медленных операциях.Я думаю, что вы неправильно поняли совет комментаторов и полностью удалили его из основного цикла.Он все еще должен быть в основном цикле, но вызываться только один раз.

for i in range(1,len(xg)):
    ....
    smxinv = np.linalg.inv(smx) ## Calculate inverse once per loop
    dx2 =  np.matmul(smxinv, pcx1)
...
ax2 = np.matmul(smxinv, px2 - np.matmul(cdx, vx2) - np.matmul(skx, dx2))
...

Самая медленная операция - matmul.Это не так уж удивительно - оно вызывается семь раз в основном цикле.Каждый вызов имеет уникальные аргументы, поэтому я не вижу простого способа ускорить его.Далее diag и diagonal.Они создают диагональный массив, который имеет в основном нулевые записи, поэтому перемещение создания за пределы цикла и только обновление ненулевых записей должно обеспечить ускорение.

##  Pre allocate px2 array (may not have a large effect)
px2 = np.diag(1).reshape(nst,1)
px2i = where(px2) ## Setup index of non-zero entries

for i in range(1,len(xg)):
    px2[px2i] = -smx*xg[i]  ## This should be equivalent
    ...

Это также удаляет вызов для изменения формы.Вы также можете предварительно вычислить некоторые константы и избежать нескольких вычислений в каждом цикле, но это, вероятно, не окажет большого влияния на общее время выполнения.

Каждый из этих шагов необходимо выполнить в тестовом примере, чтобы убедиться, чтоони не изменяют поведение функции, а затем профилируются, чтобы увидеть, насколько улучшено (если есть) улучшение.Хотя я ожидаю, что вы можете получить его менее чем за четыре или пять секунд, Python не сможет сравниться с производительностью скомпилированных языков.