Может ли этот код cython быть оптимизирован?

Вы мало что можете сделать со своей функцией 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, неплохо.

отказ от ответственности: я основной разработчик инструмента, упомянутого ниже.

В качестве альтернативы на Cython, вы можете дать Pythran попробовать. Одной аннотации в исходный код включает в себя:

#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 (Об этом заботится Пифран).

Что еще интереснее, Пифран способен распознавать несколько векторизуемых (в смысле порождающих SSE / AVX intrinsics) паттерны, через флаг опции:

> 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, без большой ручной настройки.

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

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

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