Ускоренная функция оценки для интеграции в scipy

Question

Я пытаюсь передать код из Matlab в SciPy. Вот упрощенная версия кода, который я написал до сих пор: https://gist.github.com/atmo/01b6e007be9ef90e402c. Однако версия Python значительно медленнее, чем Matlab. Я включил профилирующие результаты в суть, и они показывают, что почти 90% времени python тратит оценочную функцию f. Есть ли способ ускорить его оценку, кроме как переписать его на C или Cython?

Answer 1

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

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

В конце концов, я попытался проголосовать до quad, используя np.vectorize.

from scipy.integrate import quad 
from scipy.special import jn as besselj 
from scipy import exp, zeros, linspace 
from scipy.linalg import norm 
import numpy as np 

def complicated_func(lmbd, a, n, k): 
    u,v,w = 5, 3, 2 
    x = a*lmbd 
    fac = exp(2*x) 
    comm = (2*w + x) 
    part1 = ((v**2 + 4*w*(w + 2*x) + 2*x*(x - 1))*fac**5 
     + 2*u*fac**4 
     + (-v**2 - 4*(w*(3*w + 4*x + 1) + x*(x-2)) + 1)*fac**3 
     + (-8*(w + x) + 2)*fac**2 
     + (2*comm*(comm + 1) - 1)*fac) 
    return part1/lmbd *besselj(n+1, lmbd) * besselj(k+1, lmbd) 

def perform_quad(n, k, a): 
    return quad(complicated_func, 0, np.inf, args=(a,n,k))[0] 

def improved_main(): 
    sz = 20 
    amatrix = np.zeros((sz,sz)) 
    ls = -np.linspace(1, 10, 20)/2 
    inds = np.tril_indices(sz) 
    myv3 = np.vectorize(perform_quad) 
    res = myv3(inds[0], inds[1], ls.reshape(-1,1)) 
    results = np.empty(res.shape[0]) 
    for rowind, row in enumerate(res): 
     amatrix[inds] = row 
     symm_matrix = amatrix + amatrix.T - np.diag(amatrix.diagonal()) 
     results[rowind] = norm(symm_matrix) 
    return results 

результаты синхронизации показывают мне увеличение скорости фактора 5 (вы простите меня, если я только запустил его один раз, она занимает достаточно долго, как это):

In [11]: %timeit -n1 -r1 improved_main() 
1 loops, best of 1: 6.92 s per loop 

In [12]: %timeit -n1 -r1 main() 
1 loops, best of 1: 35.9 s per loop 

Там был также если вы заменили v сразу на его квадрат, потому что это единственный раз, когда он используется в этой сложной функции: как ее квадрат.

Там также экстремальное количество повторений в вызовах besselj, но я не вижу, как избежать этого, потому что quad будет определять lmbd, так что вы не можете легко предвычисления этих значений, а затем выполнить поиск.

Если вы прокомментируете improved_main, вы увидите, что количество звонков в complicated_func почти уменьшилось в 2 раза (диагональ еще нужно вычислить). Все остальные увеличения скорости можно отнести к np.vectorize и улучшениям до complicated_func.

У меня нет Matlab в моей системе, поэтому я не могу делать никаких утверждений о его усилении скорости, если вы улучшите сложную функцию там.

Answer 2

Возможно, ваша версия numpy сопоставима с тем, чтобы ускорить работу старших MATLAB. Но новые версии MATLAB выполняют различные формы компиляции «точно в срок», что значительно ускоряет повторные вычисления.

Я предполагаю, что вы можете проглотить код lambda и f и, возможно, сократить время их оценки пополам. Но настоящим убийцей является то, что вы звоните f столько раз.

Для начала я постараюсь сделать предварительные вычисления в f. Например, определите K1=K[1] и используйте K1 в расчетах. Это уменьшит количество вызовов индексирования. Повторяются ли экспоненты? Возможно, замените определение lambda обычным def или соедините его с f.