2016-12-24 3 views
5

Я использую cython в первый раз, чтобы получить некоторую скорость для функции. Функция принимает квадратную матрицу A чисел с плавающей запятой и выводит одно число с плавающей запятой. Функция его вычисления является permanent of a matrixМожно ли оптимизировать этот код cython?

enter image description here

Когда А 30 по 30 мой код занимает около 60 секунд, в настоящее время на моем компьютере.

В приведенном ниже коде реализована формула Balasubramanian-Bax/Franklin-Glynn для постоянной страницы wiki. Я назвал матрицу M.

Одна сложная часть кода - это массив f, который используется для удерживания указателя следующей позиции для перевода в массив d. Массив d содержит значения, равные + -1. Манипуляция f и j в цикле - это просто умный способ быстро обновить код серого.

from __future__ import division 
import numpy as np 
cimport numpy as np 
cimport cython 


DTYPE_int = np.int 
ctypedef np.int_t DTYPE_int_t 
DTYPE_float = np.float64 
ctypedef np.float64_t DTYPE_float_t 

@cython.boundscheck(False) # turn off bounds-checking for entire function 
@cython.wraparound(False) # turn off negative index wrapping for entire function 
def permfunc(np.ndarray [DTYPE_float_t, ndim =2, mode='c'] M): 
    cdef int n = M.shape[0] 
    cdef np.ndarray[DTYPE_float_t, ndim =1, mode='c' ] d = np.ones(n, dtype=DTYPE_float) 
    cdef int j = 0 
    cdef int s = 1 
    cdef np.ndarray [DTYPE_int_t, ndim =1, mode='c'] f = np.arange(n, dtype=DTYPE_int) 
    cdef np.ndarray [DTYPE_float_t, ndim =1, mode='c'] v = M.sum(axis=0) 
    cdef DTYPE_float_t p = 1 
    cdef int i 
    cdef DTYPE_float_t prod 
    for i in range(n): 
     p *= v[i] 
    while (j < n-1): 
     for i in range(n): 
      v[i] -= 2*d[j]*M[j, i] 
     d[j] = -d[j] 
     s = -s 
     prod = 1 
     for i in range(n): 
      prod *= v[i] 
     p += s*prod 
     f[0] = 0 
     f[j] = f[j+1] 
     f[j+1] = j+1 
     j = f[0] 
    return p/2**(n-1) 

Я использовал все простые оптимизации, которые я нашел в учебнике cython. Некоторые аспекты, которые я должен признать, я не совсем понимаю. Например, если я делаю массив d ints, так как значения всегда + -1, код работает примерно на 10% медленнее, поэтому я оставил его как float64.

Есть ли что-нибудь еще, что я могу сделать для ускорения кода?


Это является результатом Cython -a. Как вы можете видеть, все в цикле скомпилировано на C, поэтому основные оптимизации уже сработали.

Result of cython -a

Вот та же функция в NumPy, которая составляет более 100 раз медленнее, чем моя текущая версия Cython.

def npperm(M): 
    n = M.shape[0] 
    d = np.ones(n) 
    j = 0 
    s = 1 
    f = np.arange(n) 
    v = M.sum(axis=0) 
    p = np.prod(v) 
    while (j < n-1): 
     v -= 2*d[j]*M[j] 
     d[j] = -d[j] 
     s = -s 
     prod = np.prod(v) 
     p += s*prod 
     f[0] = 0 
     f[j] = f[j+1] 
     f[j+1] = j+1 
     j = f[0] 
    return p/2**(n-1) 

Timings обновленный

Вот тайминги (с использованием IPython) моей версии Cython, версия Numpy и улучшение romeric к коду Cython. Я установил семя для воспроизводимости.

from scipy.stats import ortho_group 
import pyximport; pyximport.install() 
import permlib # This loads in the functions from permlib.pyx 
import numpy as np; np.random.seed(7) 
M = ortho_group.rvs(23) #Creates a random orthogonal matrix 
%timeit permlib.npperm(M) # The numpy version 
1 loop, best of 3: 44.5 s per loop 
%timeit permlib.permfunc(M) # The cython version 
1 loop, best of 3: 273 ms per loop 
%timeit permlib.permfunc_modified(M) #romeric's improvement 
10 loops, best of 3: 198 ms per loop 
M = ortho_group.rvs(28) 
%timeit permlib.permfunc(M) # The cython version run on a 28x28 matrix 
1 loop, best of 3: 15.8 s per loop 
%timeit permlib.permfunc_modified(M) # romeric's improvement run on a 28x28 matrix 
1 loop, best of 3: 12.4 s per loop 

Может код Cython быть ускорен на всех?

Я использую GCC и CPU является FX AMD 8350.

+1

Да вы могли бы спросить это на [codereview.se]. – usr2564301

+2

@RadLexus Спасибо. Однако, похоже, вопросы цитона встречаются редко. Было 30 когда-либо! – eleanora

+5

@eleanora: с такими рассуждениями, что это число остается низким. – Mat

ответ

1

Этот ответ основан на ранее опубликованном коде @romeric. Я исправил код и упростил его, и добавил директиву компилятора cdivision.

@cython.boundscheck(False) 
@cython.wraparound(False) 
@cython.cdivision(True) 
def permfunc_modified_2(np.ndarray [double, ndim =2, mode='c'] M): 
    cdef: 
     int n = M.shape[0], s=1, i, j 
     int *f = <int*>malloc(n*sizeof(int)) 
     double *d = <double*>malloc(n*sizeof(double)) 
     double *v = <double*>malloc(n*sizeof(double)) 
     double p = 1, prod 

    for i in range(n): 
     v[i] = 0. 
     for j in range(n): 
      v[i] += M[j,i] 
     p *= v[i] 
     f[i] = i 
     d[i] = 1 
    j = 0 
    while (j < n-1): 
     prod = 1. 
     for i in range(n): 
      v[i] -= 2.*d[j]*M[j, i] 
      prod *= v[i] 
     d[j] = -d[j] 
     s = -s    
     p += s*prod 
     f[0] = 0 
     f[j] = f[j+1] 
     f[j+1] = j+1 
     j = f[0] 

    free(d) 
    free(f) 
    free(v) 
    return p/pow(2.,(n-1)) 

оригинальный код @romeric не инициализировать v, так что вы получите разные результаты иногда. Кроме того, я объединил две петли до while и две петли внутри while соответственно.

Наконец, сравнение

In [1]: from scipy.stats import ortho_group 
In [2]: import permlib 
In [3]: import numpy as np; np.random.seed(7) 
In [4]: M = ortho_group.rvs(5) 
In [5]: np.equal(permlib.permfunc(M), permlib.permfunc_modified_2(M)) 
Out[5]: True 
In [6]: %timeit permfunc(M) 
10000 loops, best of 3: 20.5 µs per loop 
In [7]: %timeit permlib.permfunc_modified_2(M) 
1000000 loops, best of 3: 1.21 µs per loop 
In [8]: M = ortho_group.rvs(15) 
In [9]: np.equal(permlib.permfunc(M), permlib.permfunc_modified_2(M)) 
Out[9]: True 
In [10]: %timeit permlib.permfunc(M) 
1000 loops, best of 3: 1.03 ms per loop 
In [11]: %timeit permlib.permfunc_modified_2(M) 
1000 loops, best of 3: 432 µs per loop 
In [12]: M = ortho_group.rvs(28) 
In [13]: np.equal(permlib.permfunc(M), permlib.permfunc_modified_2(M)) 
Out[13]: True 
In [14]: %timeit permlib.permfunc(M) 
1 loop, best of 3: 14 s per loop 
In [15]: %timeit permlib.permfunc_modified_2(M) 
1 loop, best of 3: 5.73 s per loop 
+0

Блестящий! Я не мог понять, почему код иногда дал неправильный ответ. Спасибо. Я думаю, что умножение на 2 можно вывести из цикла, верно? – eleanora

+0

Да, вы можете написать '' 'd [i] = 2.''' и' '' v [i] - = d [j] * M [j, i] '' ', но это не изменится исполнение кода – sebacastroh

0

Ну, один очевидный оптимизации для установки D [я] до -2 и +2 и избежать умножения на 2. Я подозреваю, это не будет иметь никакого значения, но все же.

Другой должен убедиться, что компилятор C++, который компилирует полученный код, имеет все оптимизаторы, включенные (особенно векторизация).

Цикл, который вычисляет новый v [i] s, может быть распараллелен с помощью Cython's support of OpenMP. При 30 итерациях это также может не иметь значения.

+0

Спасибо за идеи. Вы правы относительно коэффициента 2. Код на языке Cython - C (не C++). Я попробую сыграть с флагами компилятора и посмотреть, не имеет значения. Я думаю, мне нужно будет распараллелить основной цикл while, поскольку вы предлагаете, что было бы выполнимо, если бы я использовал другой способ вычисления кодов серого. – eleanora

3

Существует не так много, что вы можете сделать с вашей функцией cython, так как она уже хорошо оптимизирована. Тем не менее, вы по-прежнему сможете получить умеренное ускорение, полностью избегая вызовов до numpy.

import numpy as np 
cimport numpy as np 
cimport cython 
from libc.stdlib cimport malloc, free 
from libc.math cimport pow 

cdef inline double sum_axis(double *v, double *M, int n): 
    cdef: 
     int i, j 
    for i in range(n): 
     for j in range(n): 
      v[i] += M[j*n+i] 


@cython.boundscheck(False) 
@cython.wraparound(False) 
def permfunc_modified(np.ndarray [double, ndim =2, mode='c'] M): 
    cdef: 
     int n = M.shape[0], j=0, s=1, i 
     int *f = <int*>malloc(n*sizeof(int)) 
     double *d = <double*>malloc(n*sizeof(double)) 
     double *v = <double*>malloc(n*sizeof(double)) 
     double p = 1, prod 

    sum_axis(v,&M[0,0],n) 

    for i in range(n): 
     p *= v[i] 
     f[i] = i 
     d[i] = 1 

    while (j < n-1): 
     for i in range(n): 
      v[i] -= 2.*d[j]*M[j, i] 
     d[j] = -d[j] 
     s = -s 
     prod = 1 
     for i in range(n): 
      prod *= v[i] 
     p += s*prod 
     f[0] = 0 
     f[j] = f[j+1] 
     f[j+1] = j+1 
     j = f[0] 

    free(d) 
    free(f) 
    free(v) 
    return p/pow(2.,(n-1)) 

Вот существенные проверки и тайминги:

In [1]: n = 12 
In [2]: M = np.random.rand(n,n) 
In [3]: np.allclose(permfunc_modified(M),permfunc(M)) 
True 
In [4]: n = 28 
In [5]: M = np.random.rand(n,n) 
In [6]: %timeit permfunc(M) # your version 
1 loop, best of 3: 28.9 s per loop 
In [7]: %timeit permfunc_modified(M) # modified version posted above 
1 loop, best of 3: 21.4 s per loop 

EDIT Позволяет выполнить некоторые основные SSE векторизации путем разматывания внутренней prod петли, то есть изменить цикл в приведенной выше коде на следующее

# define t1, t2 and t3 earlier as doubles 
t1,t2,t3=1.,1.,1. 
for i in range(0,n-1,2): 
    t1 *= v[i] 
    t2 *= v[i+1] 
# define k earlier as int 
for k in range(i+2,n): 
    t3 *= v[k] 
p += s*(t1*t2*t3) 

и теперь время

In [8]: %timeit permfunc_modified_vec(M) # vectorised 
1 loop, best of 3: 14.0 s per loop 

Так почти 2X ускорения по сравнению с первоначальным оптимизированным кодом Cython, не плохо.

+0

Это замечательно. Спасибо. Два вопроса. a) Может ли это помочь, если бы я использовал одиночные прецизионные поплавки, а не удваивал и б) Есть ли способ помочь компилятору векторизовать код? – eleanora

+0

Изменение двойного на одиночное может или не поможет в зависимости от архитектуры и компилятора. Что касается векторизации, проверьте мой обновленный ответ – romeric

+0

@eleanora Oh. При разворачивании вы должны убедиться, что размер цикла для цикла является кратным размеру шага (в этом случае 2). Я обновил ответ сейчас. Кроме того, ваш код имеет проблему переполнения для varaible 'p', поскольку он очень быстро взрывается. В результате «SSE» и скалярные коды приведут к разным результатам для более крупных «M». – romeric

3

отказа от ответственности: Я ядро ​​DEV инструмента нижеперечисленного

В качестве альтернативы Cython, вы можете дать попробовать на pythran. Один аннотации к исходному коду Numpy:

#pythran export npperm(float[:,:]) 
import numpy as np 
def npperm(M): 
    n = M.shape[0] 
    d = np.ones(n) 
    j = 0 
    s = 1 
    f = np.arange(n) 
    v = M.sum(axis=0) 
    p = np.prod(v) 
    while (j < n-1): 
     v -= 2*d[j]*M[j] 
     d[j] = -d[j] 
     s = -s 
     prod = np.prod(v) 
     p += s*prod 
     f[0] = 0 
     f[j] = f[j+1] 
     f[j+1] = j+1 
     j = f[0] 
    return p/2**(n-1) 

Собран с:

> pythran perm.py 

дает прирост скорости, похожий на тот, с Cython:

> # numpy version 
> python -mtimeit -r3 -n1 -s 'from scipy.stats import ortho_group; from perm import npperm ; import numpy as np; np.random.seed(7);M = ortho_group.rvs(23)' 'npperm(M)' 
1 loops, best of 3: 21.7 sec per loop 
> # pythran version 
> pythran perm.py 
> python -mtimeit -r3 -n1 -s 'from scipy.stats import ortho_group; from perm import npperm ; import numpy as np; np.random.seed(7);M = ortho_group.rvs(23)' 'npperm(M)' 
1 loops, best of 3: 171 msec per loop 

Без необходимости реализовав sum_axis (это уже сделано pythran).

Более интересный, pythran способен распознавать несколько распараллеливаемых (в смысле генерации SSE/AVX встроенных функций) шаблонов, с помощью флага опции:

> pythran perm.py -DUSE_BOOST_SIMD -march=native 
> python -mtimeit -r3 -n10 -s 'from scipy.stats import ortho_group; from perm import npperm ; import numpy as np; np.random.seed(7);M = ortho_group.rvs(23)' 'npperm(M)' 
10 loops, best of 3: 93.2 msec per loop 

, что делает окончательное x232 ускорение с Numpy версии, ускорение, сравнимое с развернутой версией cython, без большой ручной настройки.

+0

Это очень интересное спасибо! Я проведу ваш код правильно через несколько дней и отчитаюсь. – eleanora

Смежные вопросы