有沒有辦法用numpy有效地反轉矩陣數組？

通常我會在for循環中反轉3x3矩陣的數組，如下例所示。不幸的是for循環很慢。有沒有更快，更有效的方法來做到這一點？有沒有辦法用numpy有效地反轉矩陣數組？

import numpy as np 
A = np.random.rand(3,3,100) 
Ainv = np.zeros_like(A) 
for i in range(100): 
    Ainv[:,:,i] = np.linalg.inv(A[:,:,i])

來源

2012-08-15 katrasnikj

請參閱這裏：http://stackoverflow.com/questions/211160/python-inverse-of-a-matrix – 2012-08-15 15:18:41

另外，你有看看嗎？ http://docs.scipy.org/doc/numpy/reference/routines.array-manipulation.html – 2012-08-15 15:20:25

我不認爲你正確理解我的問題。我在問如何倒置不是一個，而是多個矩陣 - 一組矩陣。 – katrasnikj 2012-08-15 15:21:49

事實證明，你在numpy.linalg代碼中被燒燬了兩層。如果你看看numpy.linalg.inv，你可以看到它只是對numpy.linalg.solve（A，inv（A.shape [0]）的調用。這會在每次迭代中重新創建單位矩陣因爲所有的數組都是相同的大小，這是浪費時間。通過預先分配單位矩陣來削減這個步驟可以減少約20％的時間（fast_inverse）。我的測試表明，預先分配數組或者分配數組它從一個結果列表並沒有太大的區別

看更深的一層，你會發現拉包例程的調用，但它包裹在幾個健全性檢查，如果你去掉所有這些，只是調用lapack在你的for循環（因爲你已經知道你的矩陣的尺寸，也許知道它是真實的，並不複雜），事情運行得更快（請注意，我已經使我的陣列更大）：

import numpy as np 
A = np.random.rand(1000,3,3) 
def slow_inverse(A): 
    Ainv = np.zeros_like(A) 

    for i in range(A.shape[0]): 
     Ainv[i] = np.linalg.inv(A[i]) 
    return Ainv 

def fast_inverse(A): 
    identity = np.identity(A.shape[2], dtype=A.dtype) 
    Ainv = np.zeros_like(A) 

    for i in range(A.shape[0]): 
     Ainv[i] = np.linalg.solve(A[i], identity) 
    return Ainv 

def fast_inverse2(A): 
    identity = np.identity(A.shape[2], dtype=A.dtype) 

    return array([np.linalg.solve(x, identity) for x in A]) 

from numpy.linalg import lapack_lite 
lapack_routine = lapack_lite.dgesv 
# Looking one step deeper, we see that solve performs many sanity checks. 
# Stripping these, we have: 
def faster_inverse(A): 
    b = np.identity(A.shape[2], dtype=A.dtype) 

    n_eq = A.shape[1] 
    n_rhs = A.shape[2] 
    pivots = zeros(n_eq, np.intc) 
    identity = np.eye(n_eq) 
    def lapack_inverse(a): 
     b = np.copy(identity) 
     pivots = zeros(n_eq, np.intc) 
     results = lapack_lite.dgesv(n_eq, n_rhs, a, n_eq, pivots, b, n_eq, 0) 
     if results['info'] > 0: 
      raise LinAlgError('Singular matrix') 
     return b 

    return array([lapack_inverse(a) for a in A]) 


%timeit -n 20 aI11 = slow_inverse(A) 
%timeit -n 20 aI12 = fast_inverse(A) 
%timeit -n 20 aI13 = fast_inverse2(A) 
%timeit -n 20 aI14 = faster_inverse(A)

結果令人印象深刻：

20 loops, best of 3: 45.1 ms per loop 
20 loops, best of 3: 38.1 ms per loop 
20 loops, best of 3: 38.9 ms per loop 
20 loops, best of 3: 13.8 ms per loop

編輯：我沒有什麼獲取返回的解決看起來不夠緊密。事實證明，'b'矩陣被覆蓋並且最終包含結果。此代碼現在提供一致的結果。

來源

2012-08-17 02:34:07

numpy數組是否必須連續並按特定順序（'C'或'F'）？ – 2012-08-17 12:44:30

非常好。你能否像eig一樣:-) – 2012-08-17 12:44:49

非常感謝你的回答。 – katrasnikj 2012-08-17 17:51:35

For循環確實不一定比替代方法慢得多，在這種情況下，它也不會幫你太多。但這裏是一個建議：

import numpy as np 
A = np.random.rand(100,3,3) #this is to makes it 
          #possible to index 
          #the matrices as A[i] 
Ainv = np.array(map(np.linalg.inv, A))

定時這個解決方案與您的解決方案產生了一個小而明顯的區別：

# The for loop: 
100 loops, best of 3: 6.38 ms per loop 
# The map: 
100 loops, best of 3: 5.81 ms per loop

我試圖用numpy的常規「矢量化」與創建的希望甚至更清潔的解決方案，但我不得不再看一下。數組A中排序的變化可能是最重要的變化，因爲它利用了numpy數組按列順序排列的事實，因此數據的線性讀出以這種方式稍微快一點。

來源

2012-08-17 01:10:49 Dikkemik

有沒有辦法用numpy有效地反轉矩陣數組？

回答

相關問題