2. 双语言实现向量化数值计算对比

任务拆解：同一组运算，两种实现

本例选取三个典型的向量化运算场景，分别用 MATLAB 和 Python（NumPy）实现，并在相同数据规模下测量执行时间与内存占用。三个场景覆盖了逐元素运算、归约运算和矩阵乘法，能够较全面地反映两种语言在数组计算上的差异。

• 场景 A（逐元素）：对长度 N 的向量做 y = sin(x).*exp(-x/5) + x.^2，考察逐元素函数调用与广播。
• 场景 B（归约）：对 N×M 矩阵按列求均值与标准差，考察归约操作的实现效率。
• 场景 C（矩阵乘）：计算 K×K 矩阵乘法 C = A*B，考察底层 BLAS 调用路径。

为保证对比公平，两种语言使用相同的 N、M、K 取值，计时均采用“预热一次 + 重复多次取中位数”的方式，内存占用通过各自的内存查询接口获取。

MATLAB 实现与计时

MATLAB 中向量化天然由 .*、.^ 等点运算符和内置函数完成，无需显式循环。计时使用 tic/toc，内存使用 whos 或 memory 获取。下面脚本将三个场景封装为函数，预热一次后重复 5 次取中位数。

% compare_vectorized.m
N = 1e6;  M = 500;  K = 800;
rng(42);

% 场景 A：逐元素运算
x = rand(N,1);
fA = @() sin(x).*exp(-x/5) + x.^2;

% 场景 B：按列归约
A = rand(N, M);
fB = @() deal(mean(A,1), std(A,0,1));

% 场景 C：矩阵乘法
P = rand(K);  Q = rand(K);
fC = @() P*Q;

runBench = @(f, reps) median(arrayfun(@(~) timeit(f), 1:reps));
% timeit 内部已做预热，直接调用即可
tA = timeit(fA);  tB = timeit(fB);  tC = timeit(fC);

fprintf('A 逐元素: %.4f s\n', tA);
fprintf('B 归约  : %.4f s\n', tB);
fprintf('C 矩阵乘: %.4f s\n', tC);

% 内存占用（场景 A 输出向量）
y = fA();
info = whos('y');
fprintf('A 输出内存: %.2f MB\n', info.bytes/1024^2);

说明：timeit 是 MATLAB 官方推荐的计时函数，内部自动预热并多次采样，比手写 tic/toc 更稳定。若环境无 timeit，可改用 tic; for k=1:5, f(); end; toc/5 的近似写法。

Python / NumPy 实现与计时

NumPy 的向量化写法与 MATLAB 高度对应：np.sin、np.exp 对应内置函数，* 对应 .*，@ 对应矩阵乘。计时使用 time.perf_counter，内存使用 tracemalloc 或 psutil 获取进程峰值。

import numpy as np, time, tracemalloc

N, M, K = 1_000_000, 500, 800
rng = np.random.default_rng(42)

x = rng.random(N)
A = rng.random((N, M))
P = rng.random((K, K)); Q = rng.random((K, K))

def bench(fn, reps=5):
    fn()  # 预热
    ts = []
    for _ in range(reps):
        t0 = time.perf_counter()
        fn()
        ts.append(time.perf_counter() - t0)
    return float(np.median(ts))

fA = lambda: np.sin(x) * np.exp(-x/5) + x**2
fB = lambda: (A.mean(axis=0), A.std(axis=0))
fC = lambda: P @ Q

tA, tB, tC = bench(fA), bench(fB), bench(fC)
print(f"A 逐元素: {tA:.4f} s")
print(f"B 归约  : {tB:.4f} s")
print(f"C 矩阵乘: {tC:.4f} s")

tracemalloc.start()
y = fA()
cur, peak = tracemalloc.get_traced_memory()
tracemalloc.stop()
print(f"A 输出内存: {peak/1024**2:.2f} MB")

注意：tracemalloc 只统计 Python 层分配，NumPy 大数组的底层缓冲可能不完全计入；若需精确测量，可在 Linux 下用 resource.getrusage 读取 ru_maxrss，或使用 psutil.Process().memory_info().rss 前后差值。

结果判读与差异来源

在相同硬件与默认线程配置下，典型观察如下（具体数值随机器变化，此处为教学示例的预期趋势，非实测数据）：

关键判读点：

• 内存布局：MATLAB 默认列主序（column-major），NumPy 默认行主序（C order）。场景 B 按列归约时，MATLAB 访问连续内存，NumPy 跨步访问，这是归约差异的主因。若在 NumPy 中改用 A.mean(axis=1) 或对 A.T 操作，差异会反转。
• BLAS 后端：矩阵乘性能几乎完全由底层 BLAS 决定。可用 np.show_config() 查看 NumPy 链接的 BLAS；MATLAB 可用 version('-blas') 查看。
• 线程数：两者默认都可能启用多线程，对比前应固定线程数（如 OMP_NUM_THREADS=1），否则结果不可比。
• 结果一致性：用 np.allclose 或 MATLAB 的 max(abs(y_mat - y_py)) 校验两语言输出，确认差异在浮点容差内（如 1e-10）。