我正在用 Python 撰寫代碼來計算離散拉普拉斯算子作為二維稀疏矩陣。代碼如下:
def discretizeLaplacian(N, step):
NN = N**2
dxx = (float) (1/((step)**2))
i_loc, j_loc, vals = np.array([]), np.array([]), np.array([])
for i in range(1, N-1):
for j in range(1, N-1):
iv0 = j*N i
iv1 = (j - 1)*N i
iv2 = (j 1)*N i
iv3 = j*N i - 1
iv4 = j*N i 1
i_loc = np.concatenate((i_loc,[iv0, iv0, iv0, iv0, iv0]))
j_loc = np.concatenate((j_loc,[iv0, iv1, iv2, iv3, iv4]))
vals = np.concatenate((vals,dxx*np.array([-4, 1, 1, 1, 1])))
sparseMatrix = csc_matrix((vals, (i_loc, j_loc)), shape = (NN, NN))
return sparseMatrix
與 Matlab 相比,for 回圈運行的時間要長得多。我想知道是否可以將其撰寫為嵌套串列理解以使其更快。
uj5u.com熱心網友回復:
性能緩慢來自演算法的糟糕復雜性。確實,原始代碼的復雜度是O(N**4)! 這是因為np.concatenate它通過復制舊陣列并在新陣列的末尾添加一些專案來創建一個新陣列。這意味著O(N**2)執行 3 個增長陣列的副本。通常,您應該避免np.concatenate在回圈中生成不斷增長的陣列。在這種情況下,您應該使用 Python 串列。
請注意,您可以使用np.tile重復陣列的值并預先計算常量dxx * np.array([-4, 1, 1, 1, 1])。
這是更正后的代碼:
def discretizeLaplacian(N, step):
NN = N**2
dxx = (float) (1/((step)**2))
tmp = dxx * np.array([-4, 1, 1, 1, 1])
i_loc, j_loc = [], []
for i in range(1, N-1):
for j in range(1, N-1):
iv0 = j*N I
iv1 = (j - 1)*N i
iv2 = (j 1)*N i
iv3 = j*N i - 1
iv4 = j*N i 1
i_loc.extend([iv0, iv0, iv0, iv0, iv0])
j_loc.extend([iv0, iv1, iv2, iv3, iv4])
i_loc = np.array(i_loc, dtype=np.float64)
j_loc = np.array(j_loc, dtype=np.float64)
vals = np.tile(tmp, (N-2)**2)
sparseMatrix = csc_matrix((vals, (i_loc, j_loc)), shape = (NN, NN))
return sparseMatrix
以下是 N=200 在我的機器上的性能結果:
Original code: 22.510 s
Corrected code: 0.068 s
雖然您可以使用串列理解來設定i_locwith [j*N I for i in range(1, N-1) for j in range(1, N-1)],但j_loc由于要附加多個專案,這并不容易。您可以構建一個串列/元組串列,然后將其展平,但這會降低代碼的可讀性并且不會更快。
如果這還不夠快,你可以使用numpy的矢量化功能(例如np.meshgrid,np.arange和基本的數學運算),以構建i_loc和j_loc。請考慮閱讀有關此的Numpy 教程。Python 回圈通常很慢,因為 Python 代碼通常被解釋。Matlab 在內部使用即時編譯器(JIT) 來快速執行這種基于回圈的代碼。您可以使用Numba或Cython在 Python 中執行相同的操作。
轉載請註明出處,本文鏈接:https://www.uj5u.com/shujuku/399153.html
