3. 矩阵算术与广播机制

在矩阵计算中，最容易出错的地方往往不是算法本身，而是运算符的语义。Python（NumPy）与MATLAB都支持矩阵运算，但同一个符号在两边可能代表完全不同的操作。本节围绕三类核心问题展开：逐元素运算、矩阵乘法（*、@、dot）以及广播机制，并通过一个实际的数据变换任务把差异与等价写法串起来。

一、逐元素运算：语义一致，写法略有差别

逐元素运算指两个形状相同的数组对应位置元素直接相加减乘除。NumPy中 a * b 默认就是逐元素乘法（Hadamard积），MATLAB中对应的是 a .* b。也就是说，NumPy的 * 等价于MATLAB的 .*，而MATLAB的 * 才是矩阵乘法。这是跨语言迁移时第一个必须建立的映射。

import numpy as np
A = np.array([[1, 2], [3, 4]])
B = np.array([[5, 6], [7, 8]])
print(A * B)        # 逐元素乘：[[5,12],[21,32]]
print(A + B)        # 逐元素加
print(A ** 2)       # 逐元素平方

% MATLAB 等价写法
A = [1 2; 3 4];
B = [5 6; 7 8];
disp(A .* B)        % 逐元素乘，结果同上
disp(A + B)
disp(A .^ 2)

注意MATLAB中 A^2 表示矩阵幂（A*A），而 A.^2 才是逐元素平方；NumPy中 A**2 是逐元素平方，矩阵幂需用 np.linalg.matrix_power(A, 2)。这类“同符号不同义”的点必须逐个厘清。

二、矩阵乘法：*、@ 与 dot 的分工

NumPy中矩阵乘法有三种常见写法：np.dot(A, B)、A @ B 以及（不推荐的）np.matmul(A, B)。对于二维数组，三者结果一致；对于高维数组，@ 与 matmul 按“最后两维做矩阵乘、前面维度按广播处理”，而 dot 的规则不同，容易产生意外形状。因此工程上推荐统一使用 @。

C = A @ B           # 推荐：矩阵乘法
C2 = np.dot(A, B)   # 二维时等价
# 验证：[[1*5+2*7, 1*6+2*8],[3*5+4*7, 3*6+4*8]]
# = [[19,22],[43,50]]

% MATLAB 中矩阵乘法直接用 *
C = A * B;          % 结果 [[19 22];[43 50]]

由此得到清晰的等价关系：NumPy的 @ 对应 MATLAB的 *；NumPy的 * 对应 MATLAB的 .*。把这条映射记牢，绝大多数迁移错误都能避免。此外，矩阵乘法要求左矩阵列数等于右矩阵行数，形状不匹配时NumPy抛出 ValueError，MATLAB报“内部矩阵维度必须一致”，两者都会在运算前拦截，属于可预期的错误。

三、广播机制：NumPy自动扩展，MATLAB需显式扩展

广播（broadcasting）是NumPy在形状不同的数组间做逐元素运算时自动扩展维度的规则。其判定从尾部维度向前对齐：两维相等，或其中一维为1，则可广播；否则报错。MATLAB从R2016b起也引入了隐式扩展（implicit expansion），规则与NumPy基本一致，但更早版本必须用 bsxfun 或 repmat 显式扩展。

# 场景：给每列减去该列均值（按列中心化）
X = np.array([[1., 2., 3.],
              [4., 5., 6.],
              [7., 8., 9.]])
col_mean = X.mean(axis=0)      # 形状 (3,)
Xc = X - col_mean              # (3,3) - (3,) 广播成功
print(Xc)

% MATLAB 等价实现
X = [1 2 3; 4 5 6; 7 8 9];
colMean = mean(X, 1);          % 1x3
Xc = X - colMean;              % 隐式扩展（R2016b+）
% 旧版本：Xc = bsxfun(@minus, X, colMean);

这里有一个方向约定需要先明确：NumPy的 axis=0 表示沿行方向压缩，得到“每列一个值”，形状为 (3,)；MATLAB的 mean(X, 1) 表示沿第1维（行）压缩，同样得到每列均值，形状为 1×3。两者语义一致，只是维度编号习惯不同。若要做行中心化，NumPy用 axis=1 得到形状 (3,)，需写成 X - row_mean[:, None] 才能正确对齐到列方向；MATLAB用 mean(X, 2) 得到 3×1，直接相减即可。

四、用一张图看清三条运算路径

下面的信号链式示意图把同一组输入 A、B 在两种语言下的三条运算路径并列呈现：逐元素乘、矩阵乘、以及广播减法，并标注了形状变化与结果，便于对照记忆。

五、实际数据变换任务：标准化与加权求和

假设有一份3个样本、3个特征的矩阵 X，任务分两步：先对每个特征做标准化（减均值、除标准差），再用权重向量 w 对特征加权求和得到每个样本的得分。这个任务同时用到广播和矩阵乘法，是检验理解的好例子。

X = np.array([[2., 4., 6.],
              [4., 6., 8.],
              [6., 8., 10.]])
mu = X.mean(axis=0)          # [4, 6, 8]
sd = X.std(axis=0)           # [1.633, 1.633, 1.633]
Z = (X - mu) / sd            # 广播：两次逐元素运算
w = np.array([0.5, 0.3, 0.2])
scores = Z @ w               # 矩阵乘：(3,3)@(3,) → (3,)
print(scores)

% MATLAB 等价实现
X = [2 4 6; 4 6 8; 6 8 10];
mu = mean(X, 1);
sd = std(X, 0, 1);           % 0 表示按 N 归一化，与 NumPy 默认一致
Z = (X - mu) ./ sd;          % 注意用 ./ 逐元素除
w = [0.5; 0.3; 0.2];
scores = Z * w;              % 矩阵乘，结果 3x1
disp(scores)

这里有两个细节值得强调。第一，NumPy的 std 默认按 N 归一化（总体标准差），MATLAB的 std(X,0,1) 同样按 N 归一化，二者一致；若MATLAB用默认的 std(X,1,1) 则按 N-1 归一化，结果会有细微差异。第二，MATLAB中除以标准差必须用 ./，若误写 / 会触发矩阵右除，语义完全不同。NumPy中 / 本身就是逐元素除，无需额外符号。

运行后 scores 的形状在NumPy中为 (3,)，在MATLAB中为 3×1。若后续要与形状为 (3,1) 的列向量拼接，NumPy需先 scores[:, None] 升维，否则广播会把 (3,) 与 (3,1) 扩展成 (3,3)，得到意料之外的矩阵。这类形状陷阱在数据流水线中非常常见，建议在每个变换步骤后打印 .shape 或 size() 做快速核对。