Matlab快速入门之线性代数：线性方程组

本文属于Matlab快速入门之线性代数的第二篇，即线性方程组，内容主要包括：计算注意事项、通解、方阵方程组、超定方程组、欠定方程组、多右端线性方程组的求解、迭代法、多线程计算等。

计算注意事项

进行科学计算时，最重要的一个问题是对联立线性方程组求解。

在矩阵表示法中，常见问题采用以下形式：给定两个矩阵 A 和 b，是否存在一个唯一矩阵x使Ax = b或xA = b？

考虑1×1示例具有指导意义。例如，方程7x = 21是否具有唯一解？

答案当然是肯定的。方程有唯一解x = 3。通过除法很容易求得该解：x = 21/7 = 3。

该解通常不是通过计算 7 的倒数求得的，即先计算 7^–1 = 0.142857...，然后将 7^–1 乘以21。这将需要更多的工作，而且如果7^–1以有限位数表示时，准确性会较低。类似注意事项也适用于多个未知数的线性方程组；MATLAB 在解此类方程时不会计算矩阵的逆。

尽管这不是标准的数学表示法，但MATLAB使用标量示例中常见的除法术语来描述常规联立方程组的解。斜杠 / 和反斜杠 \ 这两个除号分别对应 MATLAB 函数 mrdivide 和 mldivide。两种运算符分别用于未知矩阵出现在系数矩阵左侧或右侧的情况：

`x = b/A`	表示使用 `mrdivide` 获得的矩阵方程 xA = b 的解。
`x = A\b`	表示使用 `mldivide` 获得的矩阵方程 Ax = b 的解。

考虑将方程Ax = b或xA = b的两端“除以”A。系数矩阵 A 始终位于“分母”中。

x = A\b的维度兼容性条件要求两个矩阵A和b的行数相同。这样，解 x 的列数便与 b 的列数相同，并且其行维度等于A的列维度。对于x = b/A，行和列的角色将会互换。

实际上，Ax=b 形式的线性方程组比 xA=b 形式的线性方程组更常见。因此，反斜杠的使用频率要远高于斜杠的使用频率。本节其余部分将重点介绍反斜杠运算符；斜杠运算符的对应属性可以从以下恒等式推知：(b/A)' = (A'\b')。

系数矩阵A不需要是方阵。如果A的大小为 m×n，则有三种情况：

m = n	方阵方程组。求精确解。
m > n	超定方程组，即方程个数多于未知数个数。求最小二乘解。
m < n	欠定方程组，即方程个数少于未知数个数。使用最多 m 个非零分量求基本解。

mldivide 算法. mldivide运算符使用不同的求解器来处理不同类型的系数矩阵。通过检查系数矩阵自动诊断各种情况。有关详细信息，请参阅mldivide参考页的“算法”小节。

>> help mldivide
 \   Backslash or left matrix divide.
    A\B is the matrix division of A into B, which is roughly the
    same as INV(A)*B , except it is computed in a different way.
    If A is an N-by-N matrix and B is a column vector with N
    components, or a matrix with several such columns, then
    X = A\B is the solution to the equation A*X = B. A warning 
    message is printed if A is badly scaled or nearly singular.
    A\EYE(SIZE(A)) produces the inverse of A.
 
    If A is an M-by-N matrix with M < or > N and B is a column
    vector with M components, or a matrix with several such columns,
    then X = A\B is the solution in the least squares sense to the
    under- or overdetermined system of equations A*X = B. The
    effective rank, K, of A is determined from the QR decomposition
    with pivoting. A solution X is computed which has at most K
    nonzero components per column. If K < N this will usually not
    be the same solution as PINV(A)*B.  A\EYE(SIZE(A)) produces a
    generalized inverse of A.
 
    C = mldivide(A,B) is called for the syntax 'A \ B' when A or B is an
    object.

通解

线性方程组Ax = b的通解描述了所有可能的解。可以通过以下方法求通解：

求对应的齐次方程组Ax = 0的解。使用null命令通过键入null(A)来执行此操作。这会将解空间的基向量恢复为Ax = 0。任何解都是基向量的线性组合。
求非齐次方程组Ax = b的特定解。

然后，可将Ax = b的任何解写成第 2 步中的Ax = b的特定解加上第 1 步中的基向量的线性组合之和。

本节其余部分将介绍如何使用 MATLAB 求 Ax = b 的特定解，如第 2 步中所述。

方阵方程组

最常见的情况涉及到一个方阵系数矩阵 A 和一个右侧单列向量 b。

非奇异系数矩阵. 如果矩阵 A 是非奇异矩阵，则解 x = A\b 的大小与 b 的大小相同。例如：

>> A = pascal(3);
>> u = [3; 1; 4];
>> x = A\u

x =

    10
   -12
     5

可以确认 A*x 恰好等于 u。

如果 A 和 b 为方阵并且大小相同，则 x= A\b 也具有相同大小：

>> b = magic(3);
>> X = A\b

X =

    19    -3    -1
   -17     4    13
     6     0    -6

可以确认 A*x 恰好等于 b。

以上两个示例具有确切的整数解。这是因为系数矩阵选为 pascal(3)，这是满秩矩阵（非奇异的）。

奇异系数矩阵. 如果方阵 A 不包含线性无关的列，则该矩阵为奇异矩阵。如果 A 为奇异矩阵，则 Ax = b 的解将不存在或不唯一。如果 A 接近奇异或检测到完全奇异性，则反斜杠运算符 A\b 会发出警告。

如果 A 为奇异矩阵并且 Ax = b 具有解，可以通过键入以下内容求不是唯一的特定解：

>> P = pinv(A)*b

pinv(A) 是 A 的伪逆。如果 Ax = b 没有精确解，则 pinv(A) 将返回最小二乘解。例如：

>> A = [ 1     3     7
     -1     4     4
      1    10    18 ]

A =

     1     3     7
    -1     4     4
     1    10    18

为奇异矩阵，可以通过键入以下内容进行验证：

>> rank(A)

ans =

     2

由于 A 不是满秩，它有一些等于零的奇异值。

精确解。对于 b =[5;2;12]，方程 Ax = b 具有精确解，给定

>> b =[5;2;12]

b =

     5
     2
    12

>> pinv(A)*b

ans =

    0.3850
   -0.1103
    0.7066

通过键入以下内容验证 pinv(A)*b 是否为精确解：

>> A*pinv(A)*b

ans =

    5.0000
    2.0000
   12.0000

上述结果，我们发现返回值恰好为原始向量b

最小二乘解。但是，如果 b = [3;6;0]，则 Ax = b 没有精确解。在这种情况下，pinv(A)*b 会返回最小二乘解。键入：

>> b = [3;6;0]

b =

     3
     6
     0

>> A*pinv(A)*b

ans =

   -1.0000
    4.0000
    2.0000

上述结果不会返回原始向量 b。

通过得到增广矩阵 [A b] 的简化行阶梯形式，可以确定 Ax = b 是否具有精确解。为此，对于此示例请输入

>> rref([A b])

ans =

    1.0000         0    2.2857         0
         0    1.0000    1.5714         0
         0         0         0    1.0000

由于最下面一行全部为零（最后一项除外），因此该方程无解。在这种情况下，pinv(A) 会返回最小二乘解。

超定方程组

此示例说明在对试验数据的各种曲线拟合中通常会如何遇到超定方程组。

在多个不同的时间值 t 对数量 y 进行测量以生成以下观测值。可以使用以下语句输入数据并在表中查看该数据。

>> t = [0 .3 .8 1.1 1.6 2.3]';
>> y = [.82 .72 .63 .60 .55 .50]';
>> B = table(t,y)

B = 

     t      y  
    ___    ____

      0    0.82
    0.3    0.72
    0.8    0.63
    1.1     0.6
    1.6    0.55
    2.3     0.5

尝试使用指数衰减函数对数据进行建模:

y(t)=c₁+c₂e^−t。

上一方程表明，向量 y 应由两个其他向量的线性组合来逼近。一个是元素全部为 1 的常向量，另一个是带有分量 exp(-t) 的向量。未知系数 c1 和 c2 可以通过执行最小二乘拟合来计算，这样会最大限度地减小数据与模型偏差的平方和。在两个未知系数的情况下有六个方程，用 6×2 矩阵表示。

>> E = [ones(size(t)) exp(-t)]

E =

    1.0000    1.0000
    1.0000    0.7408
    1.0000    0.4493
    1.0000    0.3329
    1.0000    0.2019
    1.0000    0.1003

使用反斜杠运算符获取最小二乘解：

>> c = E\y

c =

    0.4760
    0.3413

也就是说，对数据的最小二乘拟合为：y(t)=0.4760+0.3413e^−t.

以下语句按固定间隔的 t 增量为模型求值，然后与原始数据一同绘制结果：

>> T = (0:0.1:2.5)';
>> Y = [ones(size(T)) exp(-T)]*c;
>> plot(T,Y,'-',t,y,'o')

E*c 与 y 不完全相等，但差值可能远小于原始数据中的测量误差。

如果矩形矩阵 A 没有线性无关的列，则该矩阵秩亏。如果 A 秩亏，则 AX = B 的最小二乘解不唯一。如果 A 秩亏，则 A\B 会发出警告，并生成一个最小二乘解。您可以使用 lsqminnorm 求在所有解中具有最小范数的解 X。

欠定方程组

本例演示了欠定方程组的解不唯一的情况。欠定线性方程组包含的未知数比方程多。MATLAB 矩阵左除运算求基本最小二乘解，对于 m×n 系数矩阵，它最多有 m 个非零分量。

以下是一个简单的随机示例：

>> R = [6 8 7 3; 3 5 4 1]
>> rng(0);
>> b = randi(8,2,1)

R =

     6     8     7     3
     3     5     4     1


b =

     7
     8

线性方程组 Rp = b 有两个方程，四个未知数。由于系数矩阵包含较小的整数，因此适合使用 format 命令以有理格式显示解。通过以下命令可获取特定解

>> format rat
>> p = R\b

p =

       0       
      17/7     
       0       
     -29/7

其中一个非零分量为 p(2)，因为 R(:,2) 是具有最大范数的 R 的列。另一个非零分量为 p(4)，因为 R(:,4) 在消除 R(:,2) 后起控制作用。

欠定方程组的完全通解可以通过 p 加上任意零空间向量线性组合来表示，可以使用 null 函数（使用请求有理基的选项）计算该空间向量。

>> Z = null(R,'r')

Z =

      -1/2           -7/6     
      -1/2            1/2     
       1              0       
       0              1

可以确认 R*Z 为零，并且残差 R*x - b 远远小于任一向量 x（其中

x = p + Z*q

由于 Z 的列是零空间向量，因此 Z*q 是以下向量的线性组合：

为了说明这一点，选择任意 q 并构造 x。

q = [-2; 1];
x = p + Z*q;

计算残差的范数。

>> format short
norm(R*x - b)

ans =

     0

如果有无限多个解，则最小范数解具有特别意义。您可以使用 lsqminnorm 计算最小范数最小二乘解。该解具有 norm(p) 的最小可能值。我这里用的Matlab R2016a，没有lsqminnorm函数，下方运行结果来自官方文档。

p = lsqminnorm(R,b)

p =

    -207/137   
     365/137   
      79/137   
    -424/137

多右端线性方程组的求解

某些问题涉及求解具有相同系数矩阵 A 但具有不同右端 b 的线性方程组。如果可以同时使用不同的 b 值，则可以将 b 构造为多列矩阵，并使用单个反斜杠命令求解所有方程组：X = A\[b1 b2 b3 …]。

但是，有时不同的 b 值并非全部同时可用，也就是说，您需要连续求解若干方程组。如果使用斜杠 (/) 或反斜杠 (\) 求解其中一个方程组，则该运算符会对系数矩阵 A 进行分解，并使用此矩阵分解来求解。然而，随后每次使用不同的 b 求解类似方程组时，运算符都会对 A 进行同样的分解，而这是一次冗余计算。

此问题的求解是预先计算 A 的分解，然后重新使用因子对 b 的不同值求解。但是，实际上，以这种方式预先计算分解可能很困难，因为需要知道要计算的分解（LU、LDL、Cholesky 等）以及如何乘以因子才能对问题求解。例如，使用 LU 分解，您需要求解两个线性方程组才能求解原始方程组 Ax = b：

[L,U] = lu(A);
x = U \ (L \ b);

对于具有若干连续右端的线性方程组，建议使用 decomposition 对象求解。借助这些对象，您可利用预先计算矩阵分解带来的性能优势，而不必了解如何使用矩阵因子。您可以将先前的 LU 分解替换为：

dA = decomposition(A,'lu');
x = dA\b;

如果您不确定要使用哪种分解，decomposition(A) 会根据 A 的属性选择正确的类型，类似于反斜杠的功能。

以下简单测试验证了此方法可能带来的性能优势。该测试分别使用反斜杠 (\) 和 decomposition 对同一稀疏线性方程组求解 100 次。

n = 1e3;
A = sprand(n,n,0.2) + speye(n);
b = ones(n,1);

% Backslash solution
tic
for k = 1:100
    x = A\b;
end
toc

Elapsed time is 9.006156 seconds.

% decomposition solution
tic
dA = decomposition(A);
for k = 1:100
    x = dA\b;
end
toc

Elapsed time is 0.374347 seconds.

对于这个问题，decomposition 求解比单独使用反斜杠要快得多，而语法仍然很简单。

迭代法

如果系数矩阵 A 很大并且是稀疏矩阵，分解方法一般情况下将不会有效。迭代方法可生成一系列近似解。MATLAB 提供了多个迭代方法来处理大型的稀疏输入矩阵。

函数	说明
`pcg`	预条件共轭梯度法。此方法适用于 Hermitian 正定系数矩阵 A。
`bicg`	双共轭梯度法
`bicgstab`	双共轭梯度稳定法
`bicgstabl`	双共轭梯度稳定法(l)
`cgs`	共轭梯度二乘法
`gmres`	广义最小残差法
`lsqr`	LSQR 方法
`minres`	最小残差法。此方法适用于 Hermitian 系数矩阵 A。
`qmr`	拟最小残差法
`symmlq`	对称的 LQ 方法
`tfqmr`	无转置 QMR 方法