如何用Numpy无显式循环计算满足双重求和的m×m矩阵B？

这题的关键是先从数学上简化计算逻辑，再利用Numpy的向量化操作来实现，完全不用写任何显式循环，而且效率拉满！

先搞懂数学推导 #

首先咱们把B的元素公式拆解开：

B[i, j] = ∑_{k=1,...,n} ∑_{l=1,...,n} A[i, k] * A[j, l]

注意到A[i,k]只和第i行有关，A[j,l]只和第j行有关，这两个求和项是独立的，所以可以拆成两个行求和的乘积：

B[i, j] = (∑_{k=1,...,n} A[i, k]) * (∑_{l=1,...,n} A[j, l])

换句话说，B其实就是A的行和数组与自身的外积。

Numpy实现代码 #

基于这个推导，我们只需要两步就能得到B：

1. 计算A每一行的和：

row_sums = A.sum(axis=1)

这里axis=1表示沿着列的方向求和，得到一个形状为(m,)的一维数组，每个元素对应A某一行的总和。

2. 计算行和数组的外积，有三种等价的写法，选你习惯的就行：

• 用Numpy自带的outer函数（最直观）：

B = np.outer(row_sums, row_sums)

• 用广播机制（手动扩展维度实现外积）：

B = row_sums[:, np.newaxis] * row_sums[np.newaxis, :]

• 用矩阵乘法（适合习惯线性代数写法的同学）：

B = row_sums.reshape(-1, 1) @ row_sums.reshape(1, -1)

验证正确性 #

咱们用小矩阵测试一下，确保结果对得上：

import numpy as np
m = 2
n = 3
A = np.array([[1,2,3], [4,5,6]])
row_sums = A.sum(axis=1)  # 输出: array([6, 15])
B = np.outer(row_sums, row_sums)
print(B)

输出结果：

[[ 36  90]
 [ 90 225]]

手动计算的话，B[0,0] = (1+2+3)(1+2+3)=66=36，B[0,1]=(1+2+3)(4+5+6)=615=90，完全一致，没问题！

为什么这个方案好？ #

• 完全没有显式循环，利用Numpy的底层C优化实现，速度比Python循环快几个数量级，尤其是当m、n很大的时候优势更明显；
• 代码极简，逻辑清晰，一眼就能看懂背后的数学意义；
• 内存效率高，不需要中间生成大的临时矩阵。

内容的提问来源于stack exchange，提问作者Simon Parker