3. MATLAB内存管理与向量化

在MATLAB中，性能瓶颈往往不在算法本身，而在内存的分配、复制与释放。理解数组在内存中的布局方式，是写出高效并行与GPU代码的前提。MATLAB的数值数组默认采用列优先（column-major）存储：一个 m×n 的矩阵在内存中按列连续排列，第1列的所有元素先存放，接着是第2列，依此类推。这与C/C++的行优先习惯相反，直接决定了遍历顺序对缓存命中率的影响。

考虑一个 10000×10000 的 double 矩阵，占用约 800 MB 内存。若按行遍历：

A = rand(10000);
tic;
for i = 1:10000
    for j = 1:10000
        A(i,j) = A(i,j) + 1;
    end
end
toc;

由于内层循环沿行方向跳跃访问，每次跨越一整列（10000个元素），缓存行利用率极低。交换循环顺序，让内层沿列方向连续访问，速度通常可提升数倍。这是列优先布局带来的最直接优化：最内层循环应遍历第一个维度。

写时复制：MATLAB的惰性复制机制

MATLAB采用写时复制（Copy-on-Write, COW）策略管理数组内存。当你执行 B = A; 时，MATLAB并不会立即复制A的数据，而是让B与A共享同一块内存，仅增加一个引用计数。只有当B或A被修改时，才会真正分配新内存并复制数据。这一机制在函数传参时尤为关键：

function y = processData(x)
    x(1) = x(1) * 2;   % 此处触发复制
    y = x;
end

调用 processData(A) 时，若函数内部修改了输入参数x，MATLAB必须为x分配新内存，因为A仍在外部被引用。这会导致一次完整的数据复制。若函数不修改输入，则不会复制。因此，避免在函数内部修改输入参数是减少内存复制的有效手段。

写时复制也影响变量赋值。以下代码会触发复制：

A = rand(5000);
B = A;        % 共享内存，无复制
B(1) = 0;     % 触发复制，B获得独立内存

而如果先清除A再修改B，则可能避免复制：

A = rand(5000);
B = A;
clear A;      % 引用计数降为1
B(1) = 0;     % 就地修改，无复制

这种“先clear再修改”的模式在内存敏感的场景中非常实用，尤其是在处理大型数组时。

就地操作与临时变量

MATLAB的许多运算会隐式创建临时数组。例如 C = A + B; 会先分配C的内存，再逐元素相加。若A、B、C都是大矩阵，峰值内存可能达到三倍。更隐蔽的是链式表达式：

D = A .* B + C .* D - E ./ F;

该表达式会依次创建多个中间临时变量，每个都占用完整数组大小的内存。在GPU上，这些临时变量还会导致多次内核启动和数据传输，严重拖慢速度。

减少临时变量的策略包括：

• 使用就地操作函数：如 A = A + 1; 可改写为 A(:) = A(:) + 1;，但更推荐使用支持就地输出的函数，如 A = plus(A, 1); 在某些情况下可减少临时变量。
• 分块处理：将大数组切分为小块，逐块计算并写回，避免一次性分配巨大临时数组。
• 预分配输出数组：在循环前用 zeros 或 nan 预分配结果数组，避免循环中动态增长。
• 利用 bsxfun 或隐式扩展：MATLAB R2016b之后支持隐式扩展，但底层仍可能产生临时变量，需结合 pagefun 等GPU函数优化。

以下是一个分块处理的示例，用于计算两个大矩阵的逐元素乘积并求和，避免生成完整临时矩阵：

A = rand(20000, 20000);
B = rand(20000, 20000);
total = 0;
blockSize = 2000;
for i = 1:blockSize:20000
    for j = 1:blockSize:20000
        blockA = A(i:i+blockSize-1, j:j+blockSize-1);
        blockB = B(i:i+blockSize-1, j:j+blockSize-1);
        total = total + sum(blockA .* blockB, 'all');
    end
end

该代码每次只处理 2000×2000 的块，临时变量大小可控，峰值内存显著降低。在GPU上，类似的分块策略可配合 gpuArray 和 pagefun 进一步加速。

内存布局与向量化的协同

向量化不仅是“少写循环”，更是“按内存布局组织计算”。例如，计算矩阵每列的均值：

% 低效：逐列循环
for j = 1:n
    colMean(j) = mean(A(:,j));
end

% 高效：向量化
colMean = mean(A, 1);

mean(A, 1) 沿第一个维度（列方向）归约，数据连续访问，且MATLAB内部可调用优化过的BLAS/LAPACK例程。若改为 mean(A, 2)，则沿行方向归约，访问跳跃，速度下降。因此，向量化时优先沿列方向操作，与列优先布局保持一致。

在GPU上，这一原则同样适用。GPU内存访问以合并访问（coalesced access）为最优，即相邻线程访问相邻内存地址。由于MATLAB的gpuArray同样采用列优先布局，沿列方向的向量化操作能获得更高的内存带宽利用率。

面向GPU的实践建议

将上述原则迁移到GPU时，需额外注意三点。第一，减少主机与设备间的数据传输：一旦数据传入GPU，尽量在GPU上完成所有计算，避免频繁的 gather 和 gpuArray 转换。第二，利用 pagefun 处理批量矩阵运算：对于三维数组，pagefun 可对每个页（page）并行执行矩阵乘法、求逆等操作，避免循环。第三，预分配GPU数组：使用 gpuArray.zeros 预分配，避免在GPU上动态增长数组导致反复分配和复制。

以下代码展示了在GPU上分块计算并避免临时变量的模式：

A = gpuArray.rand(10000, 10000);
B = gpuArray.rand(10000, 10000);
C = gpuArray.zeros(10000, 10000);
blockSize = 1000;
for i = 1:blockSize:10000
    for j = 1:blockSize:10000
        idx = i:i+blockSize-1;
        jdx = j:j+blockSize-1;
        C(idx, jdx) = A(idx, jdx) .* B(idx, jdx);
    end
end
result = sum(C, 'all');

该模式将大矩阵切分为 1000×1000 的块，每块在GPU上独立计算并写回，避免了生成完整临时矩阵，同时保持了合并访问。实际运行时，可通过 gpuDevice 查看显存占用，确认峰值内存是否在预期范围内。