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的条款加以总结。