NumPy的矩陣運算

Question

我想計算以下值對所有i和j：

M_ki = Sum[A_ij - A_ik - A_kj + A_kk, 1 <= j <= n] 

我如何使用NumPy的（Python）的它沒有明確的循環呢？

謝謝！

Answer 1

下面是解決這類問題的一般策略。

首先，寫一個小的腳本，在兩個不同的功能明確寫入循環，並在最後的測試確保這兩個功能是完全一樣的：

import numpy as np 
from numpy import newaxis 

def explicit(a): 
    n = a.shape[0] 
    m = np.zeros_like(a) 
    for k in range(n): 
     for i in range(n): 
      for j in range(n): 
       m[k,i] += a[i,j] - a[i,k] - a[k,j] + a[k,k] 
    return m 

def implicit(a): 
    n = a.shape[0] 
    m = np.zeros_like(a) 
    for k in range(n): 
     for i in range(n): 
      for j in range(n): 
       m[k,i] += a[i,j] - a[i,k] - a[k,j] + a[k,k] 
    return m 

a = np.random.randn(10,10) 
assert np.allclose(explicit(a), implicit(a), atol=1e-10, rtol=0.) 

然後，矢量化功能一步一步通過編輯implicit，在每一步運行腳本，以確保他們繼續保持相同：

步驟1

def implicit(a): 
    n = a.shape[0] 
    m = np.zeros_like(a) 
    for k in range(n): 
     for i in range(n): 
      m[k,i] = (a[i,:] - a[k,:]).sum() - n*a[i,k] + n*a[k,k] 
    return m 

步驟2

def implicit(a): 
    n = a.shape[0] 
    m = np.zeros_like(a) 
    m = - n*a.T + n*np.diag(a)[:,newaxis] 
    for k in range(n): 
     for i in range(n): 
      m[k,i] += (a[i,:] - a[k,:]).sum() 
    return m 

步驟3

def implicit(a): 
    n = a.shape[0] 
    m = np.zeros_like(a) 
    m = - n*a.T + n*np.diag(a)[:,newaxis] 
    m += (a.T[newaxis,...] - a[...,newaxis]).sum(1) 
    return m 

的Et瞧'！最後一個沒有循環。爲了矢量化這種方程，broadcasting是要走的路！

警告：確保explicit是您想要矢量化的公式。我不確定是否也應該將不依賴於j的條款加以總結。