numpy Matrix-Arithmetik mit seinem Diagonalelement

Question

Ich liebe numpy weil es vektorisiert Betrieb ermöglicht, wie:

mat1 = np.array([[1,2],[3,4]]) 
mat2 = np.array([[10,20],[30,40]]) 
mat3 = (mat1 + mat2)*2.0 # vectorization way. nice. 

Aber ich kann nicht finden, wie man diese Art von Operation zu tun, mit Diagonalelementen. Was Ich mag würde zu tun ist,

Ist es möglich, über eine Vektorisierung Weise mit numpy zu betreiben?

Answer 1

Zum ersten exemple:

mit:

In [3]: A 
""" 
array([[1, 3, 4, 0, 4], 
     [2, 3, 3, 3, 0], 
     [1, 0, 4, 1, 0], 
     [0, 3, 3, 2, 0], 
     [2, 1, 0, 3, 2]]) 
""" 

In [4]: Aii=vstack((diag(A),)*A.shape[0]) 
""" 
array([[1, 3, 4, 2, 2], 
     [1, 3, 4, 2, 2], 
     [1, 3, 4, 2, 2], 
     [1, 3, 4, 2, 2], 
     [1, 3, 4, 2, 2]]) 
""" 

In [5]: Ajj=Aii.T # transpose 
In [6]: B= 1/ (Aii+Ajj-2*A)  

Oder mit abstrakteren Tools:

B1 = 1/(np.add.outer(diag(A),diag(A))-2*A) 
B2 = A/np.sqrt(np.multiply.outer(diag(A),diag(A))) 

Answer 2

Für Fälle wie diejenigen, die keine Abhängigkeit zwischen Iterationen, scheint broadcasting wie die offensichtlicher Ansatz nach Erweiterung der Dimensionen der diagonalen 1D-Array zu einer 2D-Version mit None/np.newaxis und Durchführung der Berechnung s mit seiner 1D-Version, um A[i,i] bzw. A[j,j] zu simulieren. So diese beiden Fälle zu lösen, könnte man wie so vorgehen -

Ad = np.diag(A) 
case1_out_vectorized = 1/(Ad[:,None] + Ad - 2*A) 
case2_out_vectorized = A/np.sqrt(Ad[:,None]*Ad) 

Probelauf -

In [33]: # Random input array 
    ...: A = np.random.rand(4,4) 
    ...: 
    ...: # Naive loopy implmentation (used here for verification) 
    ...: m = A.shape[0] 
    ...: case1_out_loopy = np.zeros((m,m)) 
    ...: case2_out_loopy = np.zeros((m,m)) 
    ...: for i in range(m): 
    ...:  for j in range(m): 
    ...:   case1_out_loopy[i,j] = 1/(A[i,i] + A[j,j] - 2*A[i,j]) 
    ...:   case2_out_loopy[i,j] = A[i,j]/np.sqrt(A[i,i]*A[j,j]) 
    ...:   

In [34]: # Proposed approach    
    ...: Ad = np.diag(A) 
    ...: case1_out_vectorized = 1/(Ad[:,None] + Ad - 2*A) 
    ...: case2_out_vectorized = A/np.sqrt(Ad[:,None]*Ad) 
    ...: 

In [35]: np.allclose(case1_out_loopy,case1_out_vectorized) 
Out[35]: True 

In [36]: np.allclose(case2_out_loopy,case2_out_vectorized) 
Out[36]: True