MATLAB 踩坑记录

记录一下 MATLAB 踩过的小坑，MATLAB 的各种内置函数的用法实在是太奇怪了，一个小细节就可以搞出各种问题。

矩阵偏移

函数 circshift 可以用于向量或矩阵的循环偏移，不会修改原始数据。

向量偏移

A = [1, 2, 3, 4, 5];
circshift(A, 1)
% [5, 1, 2, 3, 4]

circshift(A, -1)
% [2, 3, 4, 5, 1]

矩阵的偏移就比较复杂了，按照前面的做法会以行为整体进行偏移！！！

A = [1 2 3;
    4 5 6;
    7 8 9];
circshift(A, 1);
%     7     8     9
%     1     2     3
%     4     5     6

如果希望以列进行偏移，有两种方式，第一种方式是额外指定维数（传递两个标量）

A = [1 2 3;
    4 5 6;
    7 8 9];
circshift(A, 1, 2)
%     3     1     2
%     6     4     5
%     9     7     8

第二种方式是指定每一个维度的偏移量，传递一个向量

A = [1 2 3;
    4 5 6;
    7 8 9];
circshift(A, [0, 1])
%     3     1     2
%     6     4     5
%     9     7     8

max/min 函数

MATLAB 对于内置函数的重载真是五花八门，连最基础的最大值最小值函数也有 n 种用法，会根据输入参数的不同进行不同的运算，判断逻辑不够直观，不看官方文档很容易造成误用。

下面只以最大值函数 sum 为例，最小值函数 min 的用法也是一样的。
需要预先说明的是，这里讨论的最大值的含义：对于实数直接比较大小，对于复数则是比较模的大小。

第一种用法是传入一个变量，接收一个变量的情况，此时会沿着输入数组的大小不等于 1 的第一个维度进行比较。 返回值的尺寸与输入数组基本相同，除了将第一个不为1的维度变成1。

具体来说：

如果传入的是一个向量（无论是行向量还是列向量），返回最大值，即输入尺寸 [1,n] 或 [n,1]，输出结果尺寸[1,1]，也就是标量；
如果传入的是一个矩阵，返回一个行向量，每一个元素是矩阵对应列的最大值，即输入尺寸 [m, n]（$m > 1$），输出结果尺寸[1,n]；
如果传入的是一个多维数组，在第一个不为1的维度取最大值，例如输入尺寸[3,4,5]，输出结果尺寸[1,4,5]。

例如

>> max([1,2,3])
ans =
     3
>> max([1;2;3])
ans =
     3
>> max([1,2;3,4])
ans =
     3     4

第二种用法是传入两个变量，那么会对这两个变量的对应元素进行比较，返回最大值所组成的数组。
理想情况是两个变量的尺寸相同，那么返回值的尺寸也与之相同，例如

1
2
3

>> max([1,2,4],[3,0,5])
ans =
     3     2     5

如果两个变量的尺寸不同，那么就会对它们进行隐式扩展（就像对它们直接加减所进行的操作一样），然后对扩展后的两个数组进行逐个元素比较，例如

>> max([1,2,4],2)
ans =
     2     2     4
>> max([1,2,4],[3;0;5])
ans =
     3     3     4
     1     2     4
     5     5     5

第三种用法则比较奇怪，为了与其它情况区分，必须使用 [] 作为第二个实参，可以传入一个整数以指定比较的维度。
例如 1 代表按列比较，2 代表按行比较。

>> max([1,2;3,4],[],1)
ans =
     3     4
>> max([1,2;3,4],[],2)
ans =
     2
     4

这里指定 1 就和默认的 max([1,2;3,4]) 结果一致，但是指定 2 就会返回一个列向量，由对应行的最大值组成。

其实还可以指定维度向量（行向量列向量均可），对多个维度求最大值，例如对于矩阵来说，[1,2] 就代表对所有元素求最大值。

1
2
3

>> max([1,2;3,4],[],[1,2])
ans =
     4

第四种用法则是最常见的需求：返回最大值，无论它是向量、矩阵还是多维数组，只要所有元素中的最大值，结果一定是一个标量。例如

1
2
3

>> max([1,2;3,4],[],'all')
ans =
     4

除此之外，还有很多用法，例如：

使用两个返回值接收，那么还会同时返回最大值的索引；
可以传入 missingflag，这个参数可以控制如何处理缺失值的比较。

但是目前用不到，暂时不管这些细节了。

sum 函数

求和函数 sum 有非常多的用法，和前面的最值比较相似，但是也需要整理一下。

第一种用法是传入一个变量，此时会沿着输入数组的大小不等于 1 的第一个维度进行求和。 返回值的尺寸与输入数组基本相同，除了将第一个不为1的维度变成1。

具体来说：

如果传入的是一个向量（无论是行向量还是列向量），返回所有元素之和，即输入尺寸 [1,n] 或 [n,1]，输出结果尺寸 [1,1]，也就是标量；
如果传入的是一个矩阵，返回一个行向量，每一个元素是矩阵对应列的求和，即输入尺寸 [m, n]（$m > 1$），输出结果尺寸 [1,n]；
如果传入的是一个多维数组，在第一个不为 1 的维度求和，例如输入尺寸 [3,4,5]，输出结果尺寸 [1,4,5]。

例如

>> sum([1,2,3])
ans =
     6
>> sum([1;2;3])
ans =
     6
>> sum([1,2;3,4])
ans =
     4     6

第二种用法可以传入一个整数以指定求和的维度，例如 1 代表按列求和，2 代表按行求和。（与 max 函数不同，不需要 [] 占位）

例如

>> sum([1,2;3,4],1)
ans =
     4     6
>> sum([1,2;3,4],2)
ans =
     3
     7

可以传入维数向量（行向量列向量均可），对多个维度求和，例如对于矩阵来说，[1,2] 就代表对所有元素求和。

1
2
3

>> sum([1,2;3,4],[1,2])
ans =
    10

第三种用法则是最常见的需求：返回所有元素之和，无论它是向量、矩阵还是多维数组，都对所有元素求和，结果一定是一个标量。

1
2
3

>> sum([1,2;3,4],'all')
ans =
    10

diag 函数

这个函数创建对角矩阵或获取矩阵的对角元素。（MATLAB 就那么喜欢把几个功能合到一个内置函数中才舒服吗，一点儿也不考虑可读性）

包括四种具体用法：

D = diag(v) 返回包含主对角线上向量 v 的元素的对角方阵。
D = diag(v,k) 将向量 v 的元素放置在方阵的第 k 条对角线上。（生成方阵不会丢失 v 的数据）
x = diag(A) 返回矩阵 A 的主对角线元素组成的列向量。
x = diag(A,k) 返回 A 的第 k 条对角线上元素组成的列向量。

第 k 条对角线的含义为：k=0 表示主对角线，k>0 位于主对角线上方，k<0 位于主对角线下方。

生成对角阵例如

>> v = [2 1 -1 -2 -5];
>> diag(v)
ans =
     2     0     0     0     0
     0     1     0     0     0
     0     0    -1     0     0
     0     0     0    -2     0
     0     0     0     0    -5
>> diag(v,1)
ans =
     0     2     0     0     0     0
     0     0     1     0     0     0
     0     0     0    -1     0     0
     0     0     0     0    -2     0
     0     0     0     0     0    -5
     0     0     0     0     0     0

获取对角线对应的列向量例如

>> A = reshape(1:30, 5, 6)
A =
     1     6    11    16    21    26
     2     7    12    17    22    27
     3     8    13    18    23    28
     4     9    14    19    24    29
     5    10    15    20    25    30
>> diag(A)
ans =
     1
     7
    13
    19
    25
>> diag(A,1)
ans =
     6
    12
    18
    24
    30

二维网格索引问题

MATLAB 在二维网格生成时存在两套索引约定。如果不加区分，在数值计算（NPDE）中容易产生混淆，尤其在采用 (Nx = Ny) 时，问题会更加隐蔽。

考虑对矩形区域 $[a,b]\times[c,d]$ 进行均匀剖分，首先是每个维度的节点坐标：

a = 0; b = 2;
c = 0; d = 1;
nx = 40; ny = 30;

x = linspace(a, b, nx);
y = linspace(c, d, ny);

然后需要生成二维网格对应的矩阵数据，MATLAB 提供两种二维网格对应矩阵的生成函数：

meshgrid（适合图形展示）
ndgrid（适合数值计算）

meshgrid 会返回两个 Ny × Nx 矩：

1	[X, Y] = meshgrid(x, y);

解释：X 的每一行都是 x 的副本，总行数为 Ny；Y 的每一列都是 y 的副本，总列数为 Nx。
索引关系：$X(i,j) = x(j), Y(i,j) = y(i)$
维度对应：第一维对应 y，第二维对应 x（反过来的！）
尺寸：[ny, nx]

矩阵与网格的对应关系如下图所示（第一象限沿着 $y=x$ 镜像对称）

等价实现：

X = zeros(ny,nx);
Y = zeros(ny,nx);

for i = 1:ny
    for j = 1:nx
        X(i,j) = x(j);
        Y(i,j) = y(i);
    end
end

例如

>> [X,Y] = meshgrid([1,2,3],[4,5,6,7])
X =
     1     2     3
     1     2     3
     1     2     3
     1     2     3
Y =
     4     4     4
     5     5     5
     6     6     6
     7     7     7

ndgrid 则会返回两个 Nx × Ny 矩阵

1	[X, Y] = ndgrid(x, y);

解释：X 的每一列都是 x 的副本，总列数为 Ny；Y 的每一行都是 y 的副本，总行数为 Nx。
索引关系：$X(i,j) = x(i), Y(i,j) = y(j)$
维度对应：第一维对应 x，第二维对应 y
尺寸：[nx, ny]

矩阵与网格的对应关系如下图所示（第一象限）

等价实现：

X = zeros(nx,ny);
Y = zeros(nx,ny);

for i = 1:nx
    for j = 1:ny
        X(i,j) = x(i);
        Y(i,j) = y(j);
    end
end

例如

>> [X,Y] = ndgrid([1,2,3],[4,5,6,7])
X =
     1     1     1     1
     2     2     2     2
     3     3     3     3
Y =
     4     5     6     7
     4     5     6     7
     4     5     6     7

显然 ndgrid 更符合 NPDE 的习惯，此时可以直接用 X 和 Y 进行函数赋值，得到的矩阵还是 [Nx,Ny] 尺寸的，例如

1
2
3

[X, Y] = ndgrid(x, y);
U = sin(X) + cos(Y)
V = X.*Y;

实质上是

1 2	U(i,j) = sin(X(i,j)) + cos(Y(i,j)) = sin(x(i)) + cos(y(j)) V(i,j) = X(i,j) * Y(i,j) = x(i) * y(j)

可以保证 $U(i,j)=u(x_i,y_j)$。如果用 meshgrid，这种写法的结果就是错的。

除了这两种内置函数，其实自行配置网格的坐标矩阵也是可以接受的，例如

hx = 2/(nx+1); hy = 2/(ny+1);

x = -1 + (1:nx)'*hx;
y = -1 + (1:ny)'*hy;

% Set up mesh
X = zeros(nx,ny);
Y = zeros(nx,ny);

for i = 1:nx
    for j = 1:ny
        X(i,j) = x(i);
        Y(i,j) = y(j);
    end
end

n = nx*ny;
X_vector = reshape(X,n,1);
Y_vector = reshape(Y,n,1);

虽然这段代码使用了两层 for 循环，看起来低效，但是好处是可读性非常高，保证符合 NPDE 的习惯。这段网格生成代码不是程序的性能瓶颈，因此完全没有问题，写成向量化的简洁语句反而会降低代码的可读性，尤其是这里还可能涉及到维度对应和转置的问题。

还有一个环节是二维的绘图，例如使用 surf 或 mesh，它们支持直接传入网格矩阵进行绘图，但是 MATLAB 绘图函数采用与 meshgrid 一致的索引约定（第一维对应 y，第二维对应 x），因此如果使用 ndgrid 生成的网格，绘图时通常需要进行转置以匹配该约定。

矩阵可视化

使用 imagesc 函数可以绘制二维矩阵

1	imagesc(B)

例如

1	imagesc([1,2,3;4,4,4])

imagesc 还是和前面的 meshgrid 一样，绘制的矩阵元素示意图的位置关系是符合 MATLAB 习惯的：

左上角：$A(1,1)$，右上角：$A(1,end)$；
左下角：$A(end,1)$，右下角：$A(end,end)$；
水平方向是矩阵的每一行，竖直方向是矩阵的每一列。

但是通常用这个函数只是为了检查矩阵本身，而不是展示网格上的值，因此没啥问题。

矩阵右除的歧义

MATLAB 对矩阵方程 $XA=B$ 提供了右除运算：$X = A / B$，但是这很可能与通常的除法或逐元素除法混淆，例如

>> 1/[1;2;3]
ans =
         0         0    0.3333

>> 1 /[1;2;3]
ans =
         0         0    0.3333

>> 1. /[1;2;3]
ans =
         0         0    0.3333

>> 1./[1;2;3]
ans =
    1.0000
    0.5000
    0.3333

>> 1 ./ [1;2;3]
ans =
    1.0000
    0.5000
    0.3333

这是毫无疑问的垃圾语法。

底层运算多线程导致的显著差异

即使是同一份 MATLAB 源文件，并且所有算法步骤在数学意义上都是确定性的（不含随机数），在相同 MATLAB 版本下分别运行于不同机器（例如大型 Linux 服务器和 Windows 本地机器）时，数值结果也可能不完全一致。在某些对舍入误差敏感的问题中，这种差异甚至可能表现为较明显的最终结果差异。这类现象不一定意味着源代码存在逻辑错误，而可能来自底层浮点计算路径的差异，尤其容易出现在 GMRES、稀疏矩阵求解、预条件器、SVD、QR、rank 截断等操作中。

核心原因在于，MATLAB 的许多线性代数操作并不是由 MATLAB 逐项解释执行，而是调用底层高性能数值库，例如 BLAS、LAPACK、MKL 以及稀疏矩阵相关的内部求解器。这些库通常会根据硬件环境自动选择不同的实现路径，并使用多线程并行计算来提高性能。并行计算本身并不会引入“随机误差”，但它会改变浮点运算的执行顺序，尤其是矩阵乘法、向量内积、范数计算、稀疏矩阵分解等操作中的分块方式、求和顺序和归约顺序。

由于浮点加法不满足严格的结合律，不同的求和顺序会导致不同的舍入误差。因此，即使 MATLAB 版本相同，不同操作系统、不同 CPU 指令集、不同线程数、不同线程调度方式，以及底层数值库针对硬件采用的不同优化策略，都可能使实际计算路径发生差异。对于条件数较大、迭代过程较长，或者包含阈值判断的算法，这种底层差异可能会被逐步放大。

这个问题在 GMRES 这类 Krylov 子空间迭代法中尤其明显。GMRES 的 Arnoldi 正交化过程包含大量内积和范数计算，例如

$$
h_{ij}=q_i^T A q_j, \qquad |v|_2 .
$$

这些量的最后若干位差异会影响后续 Krylov 基向量的构造，进而影响 Hessenberg 矩阵、残差历史以及停止准则。也就是说，早期很小的舍入误差并不一定只停留在最后几位；经过多步 Krylov 迭代后，它可能导致迭代次数、最终残差，甚至后续算法分支出现差异。

如果需要避免这种多线程并行导致的差异，可以在 MATLAB 中显式设置

1	maxNumCompThreads(1);

或者在启动时加上选项 -singleCompThread。

还可以通过环境变量直接禁止各种底层库使用多线程，例如

export OMP_NUM_THREADS=1
export MKL_NUM_THREADS=1
export OPENBLAS_NUM_THREADS=1
export BLIS_NUM_THREADS=1