MATLAB 踩坑记录

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

矩阵元素个数

除了最常见的size(A)可以获取数组的尺寸，还有如下函数：

numel(A) 可以获取（任意维度数组）所有元素个数；
length(A) 可以获取向量的长度（也就是元素个数），但是扩展到高维数组的行为是非常反直觉的——返回最大维数长度，相当于max(size(A))。

矩阵偏移

函数 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

矩阵连乘的效率问题

对于多个矩阵的乘法，MATLAB默认从左到右计算，也可以使用小括号改变运算顺序，虽然矩阵乘法满足结合律，但是运算顺序仍可能极大地影响运算效率！易知\(m\times n\)矩阵乘以\(n \times p\)矩阵的运算量为\(O(mnp)\)。

例如A*B*C，A 为 500×2，B 为 2×500，C 为 500×2，则下面两组运算量有巨大差异（但是MATLAB不会自动进行这种优化）

A*B*C的运算量为\(O(500\times 2 \times 500) + O(500 \times 500 \times 2)\)；
A*(B*C)的运算量为\(O(2\times 500 \times 2) + O(500 \times 2 \times 2)\)。

这里\(O\)记号的比例常数与矩阵规模无关。因此对于连续的矩阵乘法，考虑在适当位置加上小括号，使得小括号内的运算结果是一个小矩阵。

二维网格索引问题

记录一下MATLAB在生成二维网格和绘图时索引顺序问题，索引和直观的想法很不同，如果取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);

linspace得到的是一个等差数列，我们同时限制了两端的值（闭区间）和总的元素个数。

接下来基于它们生成矩阵网格

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

这里得到的X和Y都是二维矩阵，它们的尺寸是Ny\(\times\)Nx！而不是直观上的Nx\(\times\)Ny。

实际上，X的每一行都是x的副本，总行数为Ny；Y的每一列都是y的副本，总列数为Nx。

1
2
3

X(*,j) = x(j);

Y(i,*) = y(i);

通常我们使用矩阵X和Y来对网格函数赋值，得到初值矩阵，例如

1 2	U = sin(X) + cos(Y) V = X.*Y;

这实质上是 \[ \begin{aligned} U(i,j) &= \sin(X(i,j)) + \cos(Y(i,j)) = \sin(x_j) + \cos(y_i)\\ V(i,j) &= X(i,j) \times Y(i,j) = x_j\,y_i \end{aligned} \] 因此得到的矩阵U也是一个Ny\(\times\)Nx的矩阵，它的第一个分量代表y坐标，第二个分量代表x坐标，矩阵与网格的对应关系如下图所示

但是我们希望的是\(U(i,j)=u(x_i,y_j)\)，对应为Nx\(\times\)Ny的矩阵，矩阵与网格的对应关系如下图所示

下一个环节通常是绘图，例如surf或mesh

1
2
3

surf(X, Y, U);
% or
mesh(X, Y, U);

它们的逻辑也是类似的，绘制的每一个数据点坐标为 \[ P_{ij}=(X(i,j),Y(i,j),Z(i,j)) = (x_j,y_i,u_{ij}) \]

使用surf默认会显示网格，如果网格过密则会全黑，此时可以改成mesh。

注意到这里的索引顺序之后，有若干解决办法：

第一种：接受MATLAB的设定，采用第一个分量为y，第二个分量为x的表示
第二种：在meshgrid和mesh/surf函数的处理时加上转置，即算法内部仍然采用第一个分量为x，第二个分量为y的形式
第三种：拒绝使用meshgrid，自行配置网格的坐标矩阵，保证X(i,j)=x(i)和Y(i,j)=y(j)，参考代码如下

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);

矩阵可视化

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

1	imagesc(B)