這是一個嵌套回圈，其中內部索引依賴于外部索引，具有以下引數：
f = rand(1,70299)
nech=24*30*24
N=length(f);
xh=(1:nech)/24;
在 MATLAB 中：
sf2(1:nech)=0.;
sf2vel(1:nech)=0.;
count(1:nech)=0.;
for i=1:nech
for j=1:N-i-1
sf2(i)=sf2(i) (f(j i)-f(j))^2;
count(i)=count(i) 1;
end
sf2(i)=sf2(i)/count(i);
end
在 Python 中：
def structFunPython(f,N,nech):
sf2 = np.zeros(N)
count = np.zeros(N)
for i in range(nech):
indN = np.arange(1,N-i-1)
for j in indN:
sf2[i] = np.power((f[i j]-f[j]),2)
count[i] = 1
sf2[i] = sf2[i]/count[i]
return sf2
使用 cython：
import cython
cimport numpy as np
import numpy as np
def structFun(np.ndarray f,N,nech):
cdef np.ndarray sf2 = np.zeros(N), count = np.zeros(N),
for i in range(nech):
indN = np.arange(1,N-i-1)
for j in indN:
sf2[i] = np.power((f[i j]-f[j]),2)
count[i] = 1
sf2[i] = sf2[i]/count[i]
return sf2
Executions times:
Matlab: 7.8377 sec
Python: 3651.35 sec
Cython: 3336.21 sec
I have a hard time believing Python and Cython (especially Cython) are that slow for the same computation, so I think I must have made an error in my Python/Cython loops, but I can't see where.
uj5u.com熱心網友回復：
免責宣告：正如@norok2 在評論中指出的那樣，N*(N-1)/2由于使用pdist. 對您而言N = 70299，這意味著陣列中有大約 18.5 GB 的雙精度。其他索引陣列將具有相似的大小。因此，除非您的某些用例較小N，否則此答案中的矢量化方法僅在您有大量記憶體時才可行。
正如其他人所指出的，僅將代碼從一種語言翻譯成另一種語言不會導致兩種語言的最佳代碼。并且單獨使用 Cython 并不能保證加速，就像單獨使用 NumPy 不能保證加速一樣。
norok2 的回答很好地向您展示了如何使用numba或類似的東西來編譯您的數字代碼。這可以為您提供與 MATLAB 中的性能非常相似的東西，因為 MATLAB 有自己的即時 (JIT) 編譯器。還有回旋余地來優化您的編譯代碼，因為多個實作最終可能會產生截然不同的性能。
無論如何，我想說的是，您還可以通過使用 NumPy 和 SciPy 中的高級功能來加速您的代碼。特別是，您想要計算一組 1d 點之間的成對平方距離。這就是scipy.spatial.distance.pdist可以為您做的（使用'sqeuclidean'平方歐幾里得范數）。好處是它只計算每個成對距離一次（這對 CPU 和記憶體性能很有幫助），但缺點是選擇你想要總結的貢獻有點麻煩。
無論如何，這是與您的 Python 實作相對應的代碼（帶有內部回圈使用np.arange(1, N-i)而不是的修復np.arange(1, N-i-1)）：
from scipy.spatial.distance import pdist
def pdisty(f, nech):
offset = f.size - nech
res = np.zeros_like(f, shape=nech)
dists = pdist(f[:, None], metric='sqeuclidean')
counts = np.arange(offset, f.size)[::-1]
inds = np.repeat(np.arange(res.size), counts)
np.add.at(res, inds, dists[:inds.size])
res /= counts
return res
這里發生的是
- 我們計算每對唯一陣列值的成對距離并將其存盤到
dists - 我們計算每個點所涉及的對數（這是我們最后必須標準化的），將其存盤到
counts - 找出
dists它對應的每個值的一維索引（這是困難的部分），將其存盤在inds - 用于
np.add.at累積對適當輸出指數的每個貢獻 - 用計數標準化。
以下是一些時間N = 1000，norok2 的答案func2()中的相應函式在哪里：
>>> %timeit structFunPython(f, f.size - 1)
... %timeit func2(f, f.size - 1)
... %timeit pdisty(f, f.size - 1)
1.48 s ± 89.6 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
274 ms ± 2.71 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
36.7 ms ± 1.21 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
上述解決方案要快得多，但當然它仍然比完全編譯的解決方案慢。如果您對附加依賴項或在系統上安裝 llvm 有疑問，這可能是一個合理的折衷方案。底線是代碼應該適應你試圖優化它的語言。
為了完整起見，這里是我用于比較的實作（我稍微更改了簽名，因為N可以從輸入陣列計算，并且我修復了一些 fencepost 錯誤）：
def structFunPython(f, nech):
"""Very slightly modified from the question"""
N = f.size
sf2 = np.zeros(nech)
count = np.zeros(nech)
for i in range(nech):
indN = np.arange(1,N-i)
for j in indN:
sf2[i] = np.power((f[i j]-f[i]),2)
count[i] = 1
sf2[i] = sf2[i]/count[i]
return sf2
def func2(f_arr, nech):
"""Very slightly modified from norok2's answer
See https://stackoverflow.com/a/71704834/5067311
"""
n = f_arr.size
sf2 = np.zeros(nech)
for i in range(nech):
for j in range(1, n - i):
sf2[i] = (f_arr[i j] - f_arr[i]) ** 2
sf2[i] /= (n - i - 1)
return sf2
通過這些定義，所有三個函式在機器精度內給出相同的結果：
rng = np.random.default_rng()
N = 1000
nech = N - 2
f = rng.random(1000)
assert np.allclose(structFunPython(f, nech), func2(f, nech))
assert np.allclose(structFunPython(f, nech), pdisty(f, nech))
uj5u.com熱心網友回復：
我將 MATLAB 代碼重寫為等效 Python 代碼的方式可能是（注意在 MATLAB 中從 1 開始，在 Python 中從 0 開始的索引......因為我不知道在沒有背景關系的情況下應該如何調整它，所以我采用了最簡單的方法）：
import numpy as np
def func(f_arr, nech, n):
sf2 = np.zeros(nech)
count = np.zeros(nech)
for i in range(nech):
for j in range(n - i):
sf2[i] = (f_arr[i j] - f_arr[i]) ** 2
count[i] = 1
sf2[i] /= count[i]
return sf2
請注意，這count[i] = 1是無用的，因為 的最終值count[i]是預先知道的，實際上整體count是無用的，例如：
import numpy as np
def func2(f_arr, nech, n):
sf2 = np.zeros(nech)
for i in range(nech):
for j in range(n - i):
sf2[i] = (f_arr[i j] - f_arr[i]) ** 2
sf2[i] /= (n - i)
return sf2
加速
這是 Numba 加速的手動案例。這就像添加/使用 Numba@njit裝飾器一樣簡單：
import numba as nb
func_nb = nb.njit(func)
func2_nb = nb.njit(func2)
現在，、func和都執行相同的計算：func_nbfunc2func2_nb
nech = n = 6
f_arr = np.arange(n)
print(func(f_arr, nech, n))
# [9.16666667 6. 3.5 1.66666667 0.5 0. ]
print(func_nb(f_arr, nech, n))
# [9.16666667 6. 3.5 1.66666667 0.5 0. ]
print(func2(f_arr, nech, n))
# [9.16666667 6. 3.5 1.66666667 0.5 0. ]
print(func2_nb(f_arr, nech, n))
# [9.16666667 6. 3.5 1.66666667 0.5 0. ]
如果你真的需要堅持使用 Cython，這里有一個基于func2：
%%cython -c-O3 -c-march=native -a
#cython: language_level=3, boundscheck=False, wraparound=False, initializedcheck=False, cdivision=True, infer_types=True
import numpy as np
import cython as cy
cpdef func2_cy(f_arr, nech, n):
sf2 = np.zeros(nech)
_func2_cy(f_arr.astype(np.float_), sf2, nech, n)
return sf2
cdef _func2_cy(double[:] f_arr, double[:] sf2, cy.int nech, cy.int n):
for i in range(nech):
for j in range(1, n - i):
sf2[i] = sf2[i] (f_arr[i j] - f_arr[i]) ** 2
sf2[i] = sf2[i] / (n - i)
與 Numba 相比，這寫起來要復雜得多，但實作了相似的性能。訣竅是擁有一個_func2_cy與 Python 幾乎沒有互動的函式（閱讀：它以 C 速度運行）。
結果再次與以下相同func2：
nech = n = 6
f_arr = np.arange(n)
print(func2_cy(f_arr, nech, n))
# [9.16666667 6. 3.5 1.66666667 0.5 0. ]
計時
With some little toy benchmarking we get a feeling of the speed-ups, including a comparable function written as you did, and the vectorized solution proposed in Andras Deak's very fine answer (but fixing the indices to match the above):
def func_OP(f, nech, n):
sf2 = np.zeros(n)
count = np.zeros(n)
for i in range(nech):
indN = np.arange(n - i) # <-- indexing fixed
for j in indN:
sf2[i] = np.power((f[i j]-f[i]),2)
count[i] = 1
sf2[i] = sf2[i] / count[i]
return sf2
func_OP_nb = nb.njit(func_OP)
def func_pdisty(f, nech, n):
res = np.zeros(nech)
dists = scipy.spatial.distance.pdist(f[:, None], metric='sqeuclidean')
counts = np.arange(n - 1, n - nech - 1, -1)
inds = np.repeat(np.arange(res.size), counts)
np.add.at(res, inds, dists[:inds.size])
res /= (counts 1)
return res
nech = n = 6
f_arr = np.arange(n)
print(func_OP(f_arr, nech, n))
# [9.16666667 6. 3.5 1.66666667 0.5 0. ]
print(func_pdisty(f_arr, nech, n))
# [9.16666667 6. 3.5 1.66666667 0.5 0. ]
nech = n = 1000
f_arr = np.arange(n)
%timeit func_OP(f_arr, nech, n)
# 1 loop, best of 5: 1.5 s per loop
%timeit func(f_arr, nech, n)
# 1 loop, best of 5: 567 ms per loop
%timeit func2(f_arr, nech, n)
# 1 loop, best of 5: 352 ms per loop
%timeit func_OP_nb(f_arr, nech, n)
# 1000 loops, best of 5: 1.87 ms per loop
%timeit func_nb(f_arr, nech, n)
# 1000 loops, best of 5: 1.7 ms per loop
%timeit func2_nb(f_arr, nech, n)
# 1000 loops, best of 5: 768 μs per loop
%timeit func_pdisty(f_arr, nech, n)
# 10 loops, best of 5: 44.5 ms per loop
%timeit func2_cy(f_arr, nech, n)
# 1000 loops, best of 5: 1 ms per loop
nech = n = 2000
f_arr = np.arange(n)
%timeit func_OP(f_arr, nech, n)
# 1 loop, best of 5: 6.01 s per loop
%timeit func(f_arr, nech, n)
# 1 loop, best of 5: 2.3 s per loop
%timeit func2(f_arr, nech, n)
# 1 loop, best of 5: 1.42 s per loop
%timeit func_OP_nb(f_arr, nech, n)
# 100 loops, best of 5: 7.31 ms per loop
%timeit func_nb(f_arr, nech, n)
# 100 loops, best of 5: 6.82 ms per loop
%timeit func2_nb(f_arr, nech, n)
# 100 loops, best of 5: 3.05 ms per loop
%timeit func_pdisty(f_arr, nech, n)
# 1 loop, best of 5: 344 ms per loop
%timeit func2_cy(f_arr, nech, n)
# 100 loops, best of 5: 3.95 ms per loop
nech = n = 4000
f_arr = np.arange(n)
%timeit func_OP(f_arr, nech, n)
# 1 loop, best of 5: 24.3 s per loop
%timeit func(f_arr, nech, n)
# 1 loop, best of 5: 9.27 s per loop
%timeit func2(f_arr, nech, n)
# 1 loop, best of 5: 5.71 s per loop
%timeit func_OP_nb(f_arr, nech, n)
# 10 loops, best of 5: 29 ms per loop
%timeit func_nb(f_arr, nech, n)
# 10 loops, best of 5: 27.3 ms per loop
%timeit func2_nb(f_arr, nech, n)
# 100 loops, best of 5: 12.2 ms per loop
%timeit func_pdisty(f_arr, nech, n)
# 1 loop, best of 5: 706 ms per loop
%timeit func2_cy(f_arr, nech, n)
# 100 loops, best of 5: 15.9 ms per loop
...and with the input sizes you provided:
nech = 24 * 30 * 24
n = 70299
f_arr = np.random.random(n)
%timeit -n1 -r1 func_OP(f_arr, nech, n)
# 1 loop, best of 1: 1h 4min 50s per loop
%timeit -n1 -r1 func(f_arr, nech, n)
# 1 loop, best of 1: 21min 14s per loop
%timeit -n1 -r1 func2(f_arr, nech, n) # only one run / loop
# 1 loop, best of 1: 13min 9s per loop
%timeit func_OP_nb(f_arr, nech, n)
# 1 loop, best of 5: 4.74 s per loop
%timeit func_nb(f_arr, nech, n)
# 1 loop, best of 5: 4 s per loop
%timeit func2_nb(f_arr, nech, n)
# 1 loop, best of 5: 1.62 s per loop
# %timeit func_pdisty(f_arr, nech, n)
# -- MEMORY ERROR --
%timeit func2_cy(f_arr, nech, n)
# 1 loop, best of 5: 2.2 s per loop
uj5u.com熱心網友回復：
您已經得到了一些很好的答案，它們提供了快速的 numpy 和 numba 實作。我想介紹一個體面的 Cython 實作以進行公平比較。
首先，讓我們在我的機器上計時 norok2 最快的 numba 實作：
In [3]: %timeit func2_nb(f_arr, n)
3.4 s ± 129 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
為了在 Cython 中獲得快速的記憶體訪問，向其提供有關輸入陣列的所有必要資訊至關重要。在您的情況下，
f_arr是一個 C-contiguous np.ndarray，所以我們使用 C-contiguous memoryviewsdouble[::1]而不是普通的 memoryviewsdouble[:]。不同之處在于，索引 C 連續記憶體視圖減少為純 C 代碼f_arr[i]，而后者減少為f_arr[i f_arr.strides[0]].接下來值得一提的是，Python 的冪運算子
a**2將被 C 代碼替代pow(a,2)，即我們稱之為pow-function。在緊密回圈中呼叫函式會產生不必要的開銷，即使在 C 中也是如此。在 C 中，我們只需撰寫a*a. 所以讓我們在 Cython 中做同樣的事情。
在代碼中：
%%cython -c=-O3 -c=-march=native
# for MSVC use: -c=/Ox -c=/arch:AVX2
# or: -c=/Ox -c=/arch:AVX512 if your CPU supports AVX512
import numpy as np
cimport numpy as np
cimport cython
@cython.boundscheck(False)
@cython.wraparound(False)
def func2_cy(double[::1] f_arr, int nech):
cdef int n = f_arr.size
cdef double[::1] sf2 = np.zeros(nech)
cdef int i, j
for i in range(nech):
for j in range(n-i-2):
sf2[i] = (f_arr[i j 1]-f_arr[i])*(f_arr[i j 1]-f_arr[i])
sf2[i] /= (n - i - 2)
return np.asarray(sf2)
在我的機器上使用 macOS 上的 Apple Clang 13.1.6 編譯，這比前面提到的 numba 實作要慢：
In [5]: %timeit func2_cy(f_arr, n)
4.2 s ± 32 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
但是，應該知道 Cython 基本上只是一個Python 到 C的編譯器。這意味著，在我們洗掉所有 python 互動之后，我們可以嘗試使用 C 代碼時應用的相同性能技巧：傳遞優化標志-O3并啟用當前平臺上可用的所有 CPU 指令-march=native。請注意，這也意味著良好的 Cython 代碼（即在緊密回圈中沒有 python 互動的代碼）的性能在很大程度上取決于您的 C 編譯器。根據我的經驗，由于 MSVC 糟糕的自動矢量化，numba 通常比 Windows 上的 Cython 快。在 macOS/Linux 和 gcc 或 clang 上，通常情況正好相反。
One of these well-known performance tricks is loop-unrolling in order to give the compiler hints to SIMD-vectorize the specific loop. Unrolling the innermost loop, the function looks like this:
%%cython -c=-O3 -c=-march=native
# for MSVC use: -c=/Ox -c=/arch:AVX2
# or: -c=/Ox -c=/arch:AVX512 if your CPU supports AVX512
#cython: cdivision=True
import numpy as np
cimport numpy as np
cimport cython
@cython.boundscheck(False)
@cython.wraparound(False)
def func3_cy(double[::1] f_arr, int nech):
cdef int n = f_arr.size
cdef double[::1] sf2 = np.zeros(nech)
cdef double fi, fij, fij1, fij2, fij3
cdef int i, j, ntmp
for i in range(nech):
# unroll the inner loop for better SIMD vectorization
j = 1
ntmp = (n-i-2) - ((n-i-2) % 4)
# we use a while-loop since Cython doesn't support range(j, ntmp, 4)
while j < ntmp:
# unpack the values and hope the CPU keeps them in its register
fi = f_arr[i]
fij = f_arr[i j 1]
fij1 = f_arr[i j 2]
fij2 = f_arr[i j 3]
fij3 = f_arr[i j 4]
sf2[i] = ((fij-fi)*(fij-fi) (fij1-fi)*(fij1-fi) (fij2-fi)*(fij2-fi) (fij3-fi)*(fij3-fi))
j = 4
for j in range(ntmp, n - i - 2):
sf2[i] = (f_arr[i j 1] - f_arr[i]) * (f_arr[i j 1] - f_arr[i])
sf2[i] /= (n-i-2)
return np.asarray(sf2)
On my machine, this is nearly two times faster than numba:
In [7]: %timeit func3_cy(f_arr, n)
1.7 s ± 54.4 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
To take it one step further, we can parallelize the outer loop with the help of prange which is just a thin wrapper around OpenMP:
%%cython -c=-O3 -c=-march=native -c=-fopenmp --link-args=-fopenmp
# for MSVC use: -c=/Ox -c=/openmp -c=/arch:AVX2
# or: -c=/Ox -c=/openmp -c=/arch:AVX512 if your CPU supports AVX512
#cython: cdivision=True
import numpy as np
cimport numpy as np
cimport cython
from cython.parallel cimport prange
@cython.boundscheck(False)
@cython.wraparound(False)
def func4_cy(double[::1] f_arr, int nech):
cdef int n = f_arr.size
cdef double[::1] sf2 = np.zeros(nech)
cdef double fi, fij, fij1, fij2, fij3
cdef int i, j, ntmp
for i in prange(nech, nogil=True, schedule="static", chunksize=1):
# unroll the inner loop for better SIMD vectorization
j = 1
ntmp = (n-i-2) - ((n-i-2) % 4)
# we use a while-loop since Cython doesn't support range(j, ntmp, 4)
while j < ntmp:
# unpack the values and hope the CPU keeps them in its register
fi = f_arr[i]
fij = f_arr[i j 1]
fij1 = f_arr[i j 2]
fij2 = f_arr[i j 3]
fij3 = f_arr[i j 4]
sf2[i] = ((fij-fi)*(fij-fi) (fij1-fi)*(fij1-fi) (fij2-fi)*(fij2-fi) (fij3-fi)*(fij3-fi))
j = 4
for j in range(ntmp, n - i - 2):
sf2[i] = (f_arr[i j 1] - f_arr[i]) * (f_arr[i j 1] - f_arr[i])
sf2[i] /= (n-i-2)
return np.asarray(sf2)
On my machine with 8 CPU Cores, this is by far the fastest implementation:
In [9]: %timeit func4_cy(f_arr, n)
329 ms ± 15.5 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
However, it's worth mentioning that Numba also supports thread parallelism, so I'd expect a similar performance with Numba.
uj5u.com熱心網友回復：
通過一些基本的 python/numpy 重寫，我可以加快你的代碼速度
def structFunPython4(f,N,nech):
sf2 = np.zeros(N)
#count = np.zeros(N)
for i in range(nech):
# indN = np.arange(1,N-i-1)
sf2[i] = np.sum(np.power(f[i 1:N-1]-f[i],2))
#for j in range(1,N-i-1):
# sf2[i] = np.power((f[i j]-f[i]),2)
# #count[i] = 1
sf2[i] = sf2[i]/(N-i-2)
return sf2
對于適度的樣本量：
In [53]: f = np.arange(100); N=f.shape[0]; nech=98
他們匹配：
In [54]: np.allclose(structFunPython(f,N,nech),structFunPython4(f,N,nech))
Out[54]: True
和時間：
In [55]: timeit structFunPython(f,N,nech)
34.4 ms ± 109 μs per loop (mean ± std. dev. of 7 runs, 10 loops each)
In [56]: timeit structFunPython4(f,N,nech)
2.22 ms ± 35.9 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)
我首先替換for i in indN:為for i in range(1,N-1-i):. 對陣列的迭代比對串列或 的迭代慢range。
正如其他人指出的那樣，我們不需要迭代count.
但最大的變化是j用陣列切片代替迭代，對整個切片和陣列求和進行電源。
我還沒有看到i足夠的迭代來消除它。 f[i 1:N-1]slice 的長度不同，從nech低到 0。
MATLAB 做了很多 jit 編譯，所以你可以通過迭代來解決。早在 1990 年代，當我使用 MATLAB 時，它的形式會很糟糕——那些版本需要全矩陣計算才能獲得任何合理的速度。Python 級別的迭代很慢（解釋），并且在陣列上比在串列上慢。盡可能嘗試使用整體陣列方法。或用于numba編譯計算。
====
轉載請註明出處，本文鏈接：https://www.uj5u.com/ruanti/456633.html