Умножение N матриц - символьный расчет - PullRequest
Новая
       78

Умножение N матриц - символьный расчет

0 голосов
/ 20 января 2019

Какой самый эффективный (самый быстрый) способ умножить 20 одинаковых матриц 6х6 (М)? 

N = 20
w = sy.Symbol("w");v = sy.Symbol("v");p = sy.Symbol("p");q = sy.Symbol("q");c = 1;n = 1;nc = 1
M = np.array([[w*p*q,w*q,0,0,0,0], 
              [0,0,v,0,0,0], 
              [0,0,0,nc,0,c], 
              [0,0,0,0,v,0], 
              [w,w,v,nc,0,c],
              [0,0,0,n,0,1]])
Mi = np.array([[w*p*q,w*q,0,0,0,0],
               [0,0,v,0,0,0],
               [0,0,0,nc,0,c],
               [0,0,0,0,v,0], 
               [w,w,v,nc,0,c],
               [0,0,0,n,0,1]])
for l in range(N-1):
    M = np.dot(M, Mi)
difZ = sy.diff(Z2,w)
expr = w*(np.divide(difZ,Z2))
Z_lamda = sy.lambdify([w,v,p,q], expr, "numpy")

1 Ответ

0 голосов
/ 20 января 2019

Для вашего особого случая использования я бы рекомендовал использовать numpy.linalg.matrix_power (который не был упомянут в связанном вопросе).

Времена

Вот код установки, который я использовал: 

import numpy as np
import sympy as sy
sy.init_printing(pretty_print=False)

N = 20
w = sy.Symbol("w");v = sy.Symbol("v");p = sy.Symbol("p");q = sy.Symbol("q");c = 1;n = 1;nc = 1
M = np.array([[w*p*q,w*q,0,0,0,0], 
              [0,0,v,0,0,0], 
              [0,0,0,nc,0,c], 
              [0,0,0,0,v,0], 
              [w,w,v,nc,0,c],
              [0,0,0,n,0,1]])
Mi = M.copy()

и вот некоторые моменты времени, сравнивающие ваш оригинальный итеративный dot подход к matrix_power: 

%%timeit
M = Mi.copy()
for _ in range(N-1):
    M = np.dot(M, Mi)
# 527 ms ± 14.9 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

%%timeit
np.linalg.matrix_power(Mi, N)
# 6.63 ms ± 96.6 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

Так что matrix_power примерно в 80 раз быстрее.

Дополнительный бонус: matrix_power лучше работает с массивами выражений Sympy

По какой-то причине matrix_power работает лучшес Sympy, чем итерационный метод dot.Выражения в результирующем массиве будут более упрощенными с меньшим количеством терминов.Вот как можно подсчитать термины в массиве результатов: 

import numpy as np
import sympy as sy

def countterms(arr):
    return np.sum([len(e.args) for e in arr.flat])

N = 20
w = sy.Symbol("w");v = sy.Symbol("v");p = sy.Symbol("p");q = sy.Symbol("q");c = 1;n = 1;nc = 1
M = np.array([[w*p*q,w*q,0,0,0,0], 
              [0,0,v,0,0,0], 
              [0,0,0,nc,0,c], 
              [0,0,0,0,v,0], 
              [w,w,v,nc,0,c],
              [0,0,0,n,0,1]])
Mi = M.copy()

for _ in range(N-1):
    M = np.dot(M, Mi)

Mpow = np.linalg.matrix_power(Mi, N)

print("%d terms total in looped dot result\n" % countterms(M))
print("%d terms total in matrix_power result\n" % countterms(Mpow))

Вывод: 

650 terms total in looped dot result

216 terms total in matrix_power result

В частности, print(Mpow) работает намного, намного быстрее, чем print(M).

...