NumPy 矩阵乘法

矩阵乘法是线性代数和数据科学中的核心运算。NumPy提供了多种方式进行矩阵乘法，包括np.dot、np.matmul和@运算符。理解它们之间的差异以及如何在不同维度下正确使用，对于高效计算至关重要。

1. 矩阵乘法的数学定义

给定矩阵 A (形状 m×n) 和 B (形状 n×p)，它们的乘积 C = A·B 的形状为 m×p，其中每个元素 C[i,j] = Σ_{k} A[i,k] * B[k,j]。

2. 基本矩阵乘法函数

2.1 np.dot - 点积/矩阵乘法

np.dot 是最早的函数，行为取决于输入维度：

• 如果两个都是一维数组，则计算内积（返回标量）。
• 如果一个是二维，另一个是一维，则计算矩阵与向量的乘积（返回一维数组）。
• 如果两个都是二维，则执行矩阵乘法。
• 对于更高维度，它执行的是最后一个轴和倒数第二个轴的乘积（不推荐用于高维，建议使用 np.matmul 或 @）。

import numpy as np

# 一维点积
a = np.array([1, 2, 3])
b = np.array([4, 5, 6])
print("np.dot(a, b) =", np.dot(a, b))   # 1*4+2*5+3*6 = 32

# 二维矩阵乘法
A = np.array([[1, 2], [3, 4]])
B = np.array([[5, 6], [7, 8]])
print("np.dot(A, B):\n", np.dot(A, B))

# 矩阵与向量
x = np.array([1, 2])
print("np.dot(A, x) =", np.dot(A, x))   # [1*1+2*2, 3*1+4*2] = [5, 11]

输出：

np.dot(a, b) = 32
np.dot(A, B):
 [[19 22]
 [43 50]]
np.dot(A, x) = [5 11]

2.2 np.matmul 和 @ 运算符

np.matmul 和 Python 3.5+ 引入的 @ 运算符是推荐的矩阵乘法方式，它们的行为更加一致，且支持广播。

• 如果两个都是二维，执行矩阵乘法。
• 如果一个是二维，另一个是一维，则将一维数组视为行向量或列向量（根据位置），返回一维数组。
• 对于高维数组，它沿最后两个维度进行矩阵乘法，其他维度广播（类似堆叠的矩阵）。

# 使用 @ 运算符
print("A @ B:\n", A @ B)

# 矩阵与向量
print("A @ x =", A @ x)

# 高维数组示例
X = np.random.rand(3, 2, 5)
Y = np.random.rand(3, 5, 4)
Z = X @ Y   # 形状 (3, 2, 4) - 对每个批次进行矩阵乘法
print("X @ Y shape:", Z.shape)

输出：

A @ B:
 [[19 22]
 [43 50]]
A @ x = [5 11]
X @ Y shape: (3, 2, 4)

3. 一维数组与多维数组的乘法细节

当一维数组参与矩阵乘法时，它被解释为行向量还是列向量取决于位置。

A = np.array([[1, 2], [3, 4]])
v = np.array([1, 2])

# 左乘：v @ A  -> v 视为行向量 (1x2) 乘以 (2x2) 得到 (1x2) 但返回一维数组
print("v @ A =", v @ A)   # [1*1+2*3, 1*2+2*4] = [7, 10]

# 右乘：A @ v  -> v 视为列向量 (2x1) 得到 (2x1) 返回一维数组
print("A @ v =", A @ v)   # [1*1+2*2, 3*1+4*2] = [5, 11]

输出：

v @ A = [7 10]
A @ v = [5 11]

4. 批量矩阵乘法

在处理一批矩阵时（例如形状 (batch, m, n) 和 (batch, n, p)），@ 运算符会沿第一个维度广播，对每个批次进行矩阵乘法。

batch_A = np.random.rand(10, 3, 4)   # 10个3x4矩阵
batch_B = np.random.rand(10, 4, 5)   # 10个4x5矩阵
batch_C = batch_A @ batch_B           # 10个3x5矩阵
print("batch_C shape:", batch_C.shape)

# 如果两个批次维度不同但可广播，也会自动广播
A = np.random.rand(1, 3, 4)          # 可广播到 (10,3,4)
B = np.random.rand(10, 4, 5)
C = A @ B   # A 广播到 (10,3,4)，然后乘法
print("广播后 shape:", C.shape)

输出：

batch_C shape: (10, 3, 5)
广播后 shape: (10, 3, 5)

5. 外积与内积

除了标准矩阵乘法，NumPy还提供了专门的外积和内积函数。

• np.outer：两个一维数组的外积，返回二维数组。
• np.inner：两个一维数组的内积（返回标量），对于高维数组，沿最后一个轴求和（类似于 np.tensordot 的特例）。

u = np.array([1, 2, 3])
v = np.array([4, 5])

# 外积
outer = np.outer(u, v)
print("外积:\n", outer)   # shape (3,2)

# 内积（一维）
inner = np.inner(u, u)
print("内积:", inner)     # 1+4+9=14

# 二维内积：对应位置行向量点积
A = np.array([[1, 2], [3, 4]])
B = np.array([[5, 6], [7, 8]])
inner_2d = np.inner(A, B)   # A的每一行与B的每一行做点积，结果形状 (2,2)
print("二维内积:\n", inner_2d)  # [[1*5+2*6, 1*7+2*8], [3*5+4*6, 3*7+4*8]]

输出：

外积:
 [[ 4  5]
 [ 8 10]
 [12 15]]
内积: 14
二维内积:
 [[17 23]
 [39 53]]

6. 张量积与 np.einsum

对于更复杂的乘积（如张量收缩），np.einsum 提供了强大的下标表示法，可以表达任意维度的乘法求和。

A = np.array([[1, 2], [3, 4]])
B = np.array([[5, 6], [7, 8]])

# 矩阵乘法用 einsum 表示为 'ij,jk->ik'
C = np.einsum('ij,jk->ik', A, B)
print("einsum 矩阵乘法:\n", C)

# 矩阵的迹（对角线求和）
trace = np.einsum('ii->', A)
print("einsum 迹:", trace)

# 外积
u = np.array([1, 2])
v = np.array([3, 4])
outer = np.einsum('i,j->ij', u, v)
print("einsum 外积:\n", outer)

输出：

einsum 矩阵乘法:
 [[19 22]
 [43 50]]
einsum 迹: 5
einsum 外积:
 [[3 4]
 [6 8]]

einsum 非常灵活，但需要熟悉下标语法。对于简单矩阵乘法，推荐使用 @ 运算符。

7. 性能考虑

• NumPy 的矩阵乘法底层调用 BLAS 库（如 OpenBLAS、MKL），经过高度优化，对大型矩阵效率极高。
• 对于小矩阵，函数调用开销可能占主导，但通常仍比纯 Python 循环快很多。
• 如果需要重复计算大量小矩阵乘法，考虑使用批处理形式（利用广播）以减少循环开销。
• 对于稀疏矩阵，应使用 scipy.sparse 模块，而非密集矩阵乘法。