MATLAB学习笔记——2. 矩阵

概述

MATLAB 对于矩阵的支持非常好，以矩阵运算为代表的基本运算功能一直是 MATLAB 引以为自豪的核心与基础。我们可以把向量和矩阵都视作矩阵进行统一的操作。在下文中我们默认讨论二维矩阵，但 MATLAB 支持多维矩阵。行向量即行数为 1 的矩阵，列向量即列数为 1 的矩阵。

在内存中，MATLAB 使用列主序进行连续存储，与 Fortran 相同，与 C 语言是反的。在下标的使用中，MATLAB 默认下标从 1 开始，与 Fortran 相同，与 C 语言等绝大部分编程语言都不同。

矩阵的尺寸信息可以通过下面的语句获取：

• size(A) 行数和列数，返回一个1x2的行向量；
• size(A,1) 行数；
• size(A,2) 列数。

下面我们默认矩阵是由浮点数组成的二维数组，但是MATLAB也支持更高维度的数组，而且元素类型除了浮点数，还可以是字符、布尔值、结构体数组、元胞数组等。

矩阵的创建

矩阵的字面值输入：

• 空格或者逗号视作单行元素的分隔；
• 分号视作不同行元素的分隔，也可以输入回车进行换行。（输入分号+回车也只是一次换行效果，并不会叠加）

例如

>> [1 2 3 4]
ans =
     1     2     3     4
>> [1 2 3;4 5 6]
ans =
     1     2     3
     4     5     6
>> [1 2
3 4]
ans =
     1     2
     3     4
>> [5 6;
7 8]
ans =
     5     6
     7     8

(如果尺寸不合法，则会报错，不允许长短不一致的矩阵构造)

可以使用特殊矩阵的构造函数：

• zeros(m,n) 创造 m 行 n 列的全零矩阵；zeros(m) 创造 m 行 m 列的全零方阵；
• ones(m,n) 创造 m 行 n 列的全 1 矩阵；ones(m) 创造 m 行 m 列的全 1 方阵；
• rand(m,n) 创造 m 行 n 列的随机矩阵；rand(m) 创造 m 行 m 列的随机方阵；
• eye(m) 创造 m 行 m 列的单位矩阵。

可以从文本文件直接导入一个矩阵，例如

1 2 3
4 5 6
7 8 9

使用load a.txt就会创建矩阵 a，并且把 txt 文件中的数据赋值给 a。(MATLAB 导入数据文件并不要求文件的后缀名，建议导入的文本文件名符合 MATLAB 的变量命名规范)

冒号表达式用于产生等差数列：(使用冒号表达式可以避免循环，计算效率更高)

• a:b:c 首项为a，公差为b，尾项由c确定但不一定包含。

>> 1:2:10
ans =
     1     3     5     7     9

% for(i=a;i<=c;i++){
%   ...
% }

• a:b 相当于把中间的默认步长b=1省略。

与之类似的是linspace函数：

• linspace(a,b,n) 生成行向量，第一个值是 a，最后一个值是 b，中间为等差数列，一共 n 个数；
• 缺省 n 时，默认取 n=100。

矩阵的运算

矩阵可以直接和数相加减，默认是逐个元素进行的。

>> zeros(2,3)+2
ans =
     2     2     2
     2     2     2

可以对矩阵执行内置函数，默认会逐个元素执行。

>> a=[1 2 3;4 5 6];sin(a)
ans =
    0.8415    0.9093    0.1411
   -0.7568   -0.9589   -0.2794

矩阵的两种转置：共轭转置A'，转置A.'。(对于实数矩阵两者一样)

矩阵乘法：A*B执行标准的矩阵乘法。

>> A=[1 1;2 2];B=[1 2;3 4];A*B
ans =
     4     6
     8    12

两个矩阵进行逐个元素的乘法：A.*B

>> A=[1 1;2 2];B=[1 2;3 4];A.*B
ans =
     1     2
     6     8

同理，还有逐个元素的乘方.^。 尤其注意，这里的运算例如+或.*对于尺寸兼容的 A 和 B 也是可以的，

矩阵运算时的隐式扩展

两个矩阵称为兼容的，如果它们的尺寸分别为(m1,m2,m3)和(n1,n2,n3)，总是可以加 1 使得描述它们尺寸的整数个数一致。 要求要么 m1=n1，要么其中一个为 1，这样为 1 的那个矩阵可以在这个维度上简单复制达到相同的尺寸，对于 m2 与 n2，m3 与 n3 也同理。

例如

A=[1 2 3]; %行向量
B=[10;20;30;40]; %列向量

它们做逐个元素的乘法时，会自动对尺寸进行扩展，在这个例子中就是都变成尺寸(3,4)的矩阵，在扩展的维度上只是单纯的复制

A->
1 2 3
1 2 3
1 2 3

B->
10 10 10
20 20 20
30 30 30
40 40 40

然后进行运算，计算A+B和A.*B得到的结果分别为

A+B =
    11    12    13
    21    22    23
    31    32    33
    41    42    43

A.*B =
    10    20    30
    20    40    60
    30    60    90
    40    80   120

矩阵的拼接

我们可以很方便地把两个矩阵拼起来，包括水平拼接和竖直拼接，它们分别要求两个矩阵有相同的行数或列数。

>> a=[1 2];b=[3 4]
b =
     3     4
>> [a b]
ans =
     1     2     3     4
>> [a;b]
ans =
     1     2
     3     4

矩阵的扩充赋值

MATLAB的矩阵尺寸非常自由，可以在赋值时直接扩充，例如扩充向量

a = [];
a(1) = 1;
a(2) = 2;

a % [1 2]

扩充矩阵

b = zeros(2,2);
b(2,:) = [2 2];
b(3,:) = [3 3];

b %
% 0 0
% 2 2
% 3 3

扩充后对于没有提供值的部分会自动补0。

矩阵的索引

注意索引全部使用小括号。

线性索引

首先，无论是矩阵还是退化的行向量列向量，或者高维矩阵，每一个元素都可以使用一个正整数进行索引，称为线性索引。(注意，在 C 语言等编程语言中，对二维数组用一个整数索引会得到一个一维数组，)

线性索引从 1 开始，完全对应于矩阵在内存中列主序的排列顺序。(在同一列中，相邻的元素就在内存中相邻)

>> a=[1 2 3;4 5 6]
a =
     1     2     3
     4     5     6
>> a(1)
ans =
     1
>> a(2)
ans =
     4
>> a(3)
ans =
     2
>> a(4)
ans =
     5

一般索引

更直观的索引是直接指定行数和列数(用逗号分隔)，对于更高维的数组同理。对于行向量，则和线性索引没有区别。

切片

矩阵支持丰富的切片操作，可以提取指定的部分元素，例如

• A(m:n,p) 第 m 到 n 行，第 p 列；
• A(m,:) 第 m 行，所有列；
• A(m:end,:) 第 m 到最后一行，所有列；

在对切片进行赋值时，切片可以是当前矩阵不存在的部分，赋值完成后会自动对矩阵进行扩充。

>> A=[1 2]
A =
     1     2
>> A(2,:)=1 % 这里是把切片的每一个元素都赋值1
A =
     1     2
     1     1

例如对行向量的操作：行向量x

• x(m) 是第 m 个元素，注意任何索引从 1 开始！
• x(m:n) 是第 m 到 n 个元素；
• x(m:end) 是第 m 个元素到最后；
• x(:) 取出向量的所有元素全部输出，并且作为列向量输出；

>> x=1:10
x =
     1     2     3     4     5     6     7     8     9    10
>> x(3)
ans =
     3
>> x(1:5)
ans =
     1     2     3     4     5
>> x(5:end)
ans =
     5     6     7     8     9    10
>> x(:)
ans =
     1
     2
     3
     4
     5
     6
     7
     8
     9
    10

矩阵的逆和左除右除

MATLAB提供了直接对矩阵求逆的语句：inv(A)，但是一般不会对大矩阵直接求逆，效率很低。 在科学计算编程中主要会求解线性方程组，包括下面的左除和右除语句：(两者一般不相等)

• 左除：A\B 求解 AX=B，实质结果为inv(A)*B；
• 右除：A/B 求解 XA=B，实质结果为B*inv(A)；

在使用A\B在求解线性方程组时会根据线性方程组的特点灵活使用各种数值算法，比直接求逆inv(A)快得多！

补充

矩阵连乘的效率问题

对于多个矩阵的乘法，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得到的是一个等差数列，我们同时限制了两端的值（闭区间）和总的元素个数。

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

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

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

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

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

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

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

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

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