Матрица покинула деление сложенных массивов, используя NumPy - PullRequest
Новая
       41

Матрица покинула деление сложенных массивов, используя NumPy

0 голосов
/ 10 июля 2019

Я работаю над программой для решения уравнений Блоха (или, точнее, Блоха МакКоннелла) в python.Таким образом, уравнение, которое нужно решить:

Bloch McConnell equation

, где A равно NxN matrix, A + его псевдообратное и M0 и B - векторы размера N.

Особенность в том, что я хочу решить уравнениедля нескольких смещений (и, следовательно, нескольких матриц A ) одновременно.Новые размеры:

  • A: MxNxN
  • b: Nx1
  • M0: MxNx1

Стандартная версия программы(используя цикл по 1-му измерению размера M) работает нормально, но я застрял в одной точке в «параллельной версии».

В данный момент мой код выглядит так: 

def bmcsim(A, b, M0, timestep):
    ex = myexpm(A*timestep) # returns stacked array of size MxNxN
    M = np.zeros_like(M0)
    for n in range(ex.shape[0]):
        A_tmp = A[n,:,:]
        A_b = np.linalg.lstsq(A_tmp ,b, rcond=None)[0]
        M[n,:,:] = np.abs(np.real(np.dot(ex[n,:,:], M0[n,:,:] + A_b) - A_b))      
    return M

и я бы хотел избавиться от этой for n in range(ex.shape[0]) петли.К сожалению, np.linalg.lstsq не работает для стековых массивов, не так ли?В myexpm используется np.apply_along_axis для другой проблемы: 

def myexpm(A):
    vals,vects = np.linalg.eig(A)
    tmp = np.einsum('ijk,ikl->ijl', vects, np.apply_along_axis(np.diag, -1, np.exp(vals)))
    return np.einsum('ijk,ikl->ijl', tmp, np.linalg.inv(vects))

Однако это работает только для входных данных 1D.Есть ли что-то подобное, что я могу использовать с np.linalg.lstsq?np.dot в bmcsim будет заменено на np.einsum, как в myexpm Я полагаю, или есть лучшие способы?

Спасибо за вашу помощь!

Обновление: Я только что понял, что могу заменить np.linalg.lstsq(A,b) на np.linalg.solve(A.T.dot(A), A.T.dot(b)) и сумел избавиться от цикла следующим образом: 

def bmcsim2(A, b, M0, timestep):
    ex = myexpm(A*timestep)
    b_stack = np.repeat(b[np.newaxis, :, :], offsets.size, axis=0)
    tmp_left = np.einsum('kji,ikl->ijl', np.transpose(A), A)
    tmp_right = np.einsum('kji,ikl->ijl', np.transpose(A), b_stack)
    A_b_stack = np.linalg.solve(tmp_left , tmp_right )

    return np.abs(np.real(np.einsum('ijk,ikl->ijl',ex, M0+A_b_stack ) - A_b_stack ))

Это примерно в 3 раза быстрее, но все же немного сложнее,Я надеюсь, что есть лучший способ (короче / проще), может быть, даже быстрее?!

...