Совокупный результат умножения матриц - PullRequest
Новая
       75

Совокупный результат умножения матриц

0 голосов
/ 30 апреля 2020

Учитывая список nxn матриц, я хочу вычислить совокупное произведение этих умножений матриц - то есть, учитывая матрицы M0, M1, ..., Mm, я хочу результат R, где R [0] = M1, R [1] = M0 x M1, R [2] = M0 x M1 x M2 и так далее.

Очевидно, что вы можете сделать это с помощью for-loop или хвостовой рекурсии, но я кодирую в python, где это происходит с темпом улитки.

В коде: 

def matrix_mul_cum_sum(M):
   #M[1..m] is an m x n x n  matrix

   if len(M) == 0:
       return []

   result = [M[1]]
   for A in M[1:]:
       result.append(np.mat_mul(result[-1],A))
   return result

Для указанного приложения c n = 4, а m равно 1000 (если это имеет значение). Я надеялся построить для этого формулу numpy, но мне не очень повезло. PyTorch также приемлем, но я не знаю его хорошо. В целом, все, что быстрее, чем версия basi c for-l oop, будет высоко оценено.

Заранее спасибо.

РЕДАКТИРОВАТЬ: Добавление некоторых тестов (@hpaulj упомянул, что np. linalg.multidot здесь очень медленный, и он правильный). 

import numpy as np
np.set_printoptions(suppress=True)
import time
from functools import reduce
n=1000
M = np.random.uniform(0,1,(n,4,4))
attempts = 1000

start = time.time()
for i in range(attempts):
    x1=reduce(np.dot,M)
trd = time.time()-start
trd = np.round(trd/attempts,8)
print('reduce with np.dot: ',trd)


start = time.time()
for i in range(attempts):
    x2=reduce(np.matmul,M)
trm = time.time()-start
trm = np.round(trm/attempts,8)
print('reduce with np.matmul: ',trm)

start = time.time()
for i in range(attempts):
    x3 = M[0]
    for m in M[1:]:
        x3=np.dot(x3,m)
tfd = time.time()-start
tfd = np.round(tfd/attempts,8)
print('for-loop with np.dot:',tfd)

start = time.time()
for i in range(attempts):
    x4 = M[0]
    for m in M[1:]:
        x4=np.matmul(x4,m)
tfm = time.time()-start
tfm = np.round(tfm/attempts,8)
print('for-loop with np.matmul:',tfm)

def tail_rec(x):
    r = x[0]
    helper(x,1,r)
    return r

def helper(x,i,r):
    if i == len(x):
        return
    helper(x,i+1,np.matmul(x[i],r))

start = time.time()
for i in range(attempts):
    x5 = tail_rec(M)
tt = time.time()-start
tt = np.round(tt/attempts,8)
print('tail rec with matmul:',tt)

assert np.allclose(x1,x2,x3,x4,x5)

a = np.array([trd,trm,tfd,tfm,tt])
a = np.round(a/a[:,None],3)
print(a)

output: 

reduce with np.dot:  0.00118456
reduce with np.matmul:  0.00126214
for-loop with np.dot: 0.00122364
for-loop with np.matmul: 0.00142755
tail rec with matmul: 0.00194438
[[1.    1.065 1.033 1.205 1.641]
 [0.939 1.    0.969 1.131 1.541]
 [0.968 1.031 1.    1.167 1.589]
 [0.83  0.884 0.857 1.    1.362]
 [0.609 0.649 0.629 0.734 1.   ]]

Matrix дает относительную производительность методов. запись i, j - это среднее время метода i, разделенное на среднее время метода j ... лучше использовать меньшие числа в столбце.

Похоже, что первые три метода примерно одинаковы.

...