线性代数 --- Gauss消元的部分主元法和完全主元法

Gauss消元的部分主元法和完全主元法

心怀二意的人，在他一切所行的路上都没有定见。----雅各书1章8节

笔者的一些话：刚开始写这篇文章的时候，我觉得高斯消元很简单。因为，这时的我已经完成了我一直想写的一篇关于高斯消元的文章。

线性代数 --- 什么是高斯消元法，什么又是高斯-若尔当消元？_松下J27的博客-CSDN博客_高斯若尔当消元法Gauss Jordan Elimination高斯若尔当消元法https://blog.csdn.net/daduzimama/article/details/120486666 不仅如此，我还完成了我自认为比较满意的另一篇巨作 ---> 矩阵的LU分解，顺便说一下，完成这篇文章一直是我的一个心愿。

线性代数 --- LU分解（Gauss消元法的矩阵表示）_松下J27的博客-CSDN博客LU分解是一个伟大的分解。https://blog.csdn.net/daduzimama/article/details/124684810

因为，我觉得LU分解的重要性，在国内被低估了。不管是大学的教材还是考试，大家好像只是侧重于LU分解的计算过程，或者说是求解步骤，但是对于LU分解的真正要义，则不太关注。比如说，很少有人会问：

“为什么要有LU分解？”

“学那么多分解干啥，我就用高斯消元法打天下，可以吗？”

尤其是，很少有人会真正去关心LU分解每一步的由来，每一个矩阵的形成过程，等等。当然，没有这些想法和问题我也非常理解，毕竟让我现在去参加大学线性代数的考试，肯定也是不及格。此外，在我们平时的工作中，我们也会遇到一些问题，比如说，我在用matlab的时候，我就发现

“为什么有时候用matlab自带的LU分解函数算出来的L矩阵和U矩阵，和我自己手算的不一样？”

“为什么matlab算出来的LU分解中的L（很多时候）并不是标准的下三角矩阵？”

“为什么matlab的LU分解函数会提供那么多的用法？”

比如说，下面这个matlab官方的说明文件。

“为什么LU分解有这么多变种？”

“P矩阵是什么，D矩阵是什么？Q矩阵又是什么？”

等等这些问题，都应该引起或激发我们的兴趣。大家要相信一点，兴趣是最好的老师，但是，很多时候，兴趣并不是生来就有的，就拿我来说，我自己一开始对线性代数并不感兴趣，而是因为工作需求，确切的说，是因为，我在工作中，并不喜欢浅尝辄止。比如说，在公司里面做项目，正确的做法是满足公司和客户的需求，在预定的时间节点保质保量的完成。但我不会这样，我想的是，做这个项目我能从中间学到什么。所以，我不会在公司预定的时间节点保质保量的完成。而是按照我自己预定的时间节点完成任务(这一点也不自私)，不搞懂这个项目背后的问题，我是不会下手的。而且，我还不喜欢调别人的库，所以，我喜欢自己动手，自己实现(这会让你更有成就感，因为，你以后每次再调用函数的时候，都在调用你自己的库)。

在我其中的一个项目中，我就没用现成的库。为此，我自己用C语言实现高斯消元法和LU分解，我参考了一本叫《numerical recipes in C》的书。书里面提到了几个之前没有听过的，有趣的概念--->"不选主元高斯消元法"和“选主元高斯消元法”，而选主元的消元法里面又有“部分主元法”，“完全主元法”和“隐式主元法”。正如下面的几页截图里看到的。

乍一眼看上去，觉得豁然开朗，给人一种so easy的感觉。觉得原来这就是所谓的“部分主元法”和“完全主元法”啊！。我写这篇文章的初衷，就是想把这两个比较陌生的概念分享给大家(就是本文的标题)，如果今后大家在别的地方遇到了这两个陌生的词汇，马上就能明白。哪知道真正写起来才发现，并没有自己想象的那么简单。。。（以下才是正文）

不选主元法

在说部分主元法和完全主元法之前，我们要先说说高斯消元法的另一种，不选主元法。

不选主元法的做法就是：对于方阵而言，令主对角线上的元素a11,a22,....ann为主元，从第一列开始，依次消除该列主元下的所有元素，直到最后一列。例如，对于第一列而言，令a11处为主元，然后消除a11下面的所有元素。然后，再令a22处的元素为主元，消除a22下面的所有元素。以此类推。（下图中我用红笔圈出来的元素就是主元所在的位置。）他与部分主元法和完全主元法的最大区别就是：他既不进行行交换，也不进行列交换。

上图中，用蓝色圆珠笔标注的2.1，3.1，2.3(笔记中的2.3有误，应改为3.2，表示第3行第2列的元素)，表示这一步消除的元素在矩阵A中的位置。蓝色钢笔所圈出的数值表示为了消除该元素，主元所乘的倍数。此外，这个倍数的负数，就是用于记录每一步消元过程的消元矩阵E中对应位置的值。这些消元矩阵E在我的笔记中没有写出，我在这里补上。

$\small E_{21}=\begin{bmatrix} 1 & 0 & 0\\ -2& 1 &0 \\ 0&0 &1 \end{bmatrix}$ , $\small E_{31}=\begin{bmatrix} 1 &0 &0 \\ 0& 1 &0 \\ 1&0 &1 \end{bmatrix}$ , $\small E_{32}=\begin{bmatrix} 1&0 &0 \\ 0& 1& 0\\ 0& 3&1 \end{bmatrix}$

利用这三个消元矩阵，就可以用矩阵的计算来描述整个消元过程，即：

$\small E_{32}E_{31}E_{21}A=U$

或者，我们令A左边的所有消元矩阵的乘积为矩阵Z。即：

$\small Z=E_{32}E_{31}E_{21}$

$\small ZA=U$

这样一来，只需通过一次矩阵的操作ZA，就能完成之前三次的消元。

此外，我们也同时得到了U的还原矩阵L，所谓还原就是消元的逆过程：

$\large L={E_{21}}^{-1}{E_{31}}^{-1}{E_{32}}^{-1}$

幸运的是，由于你所使用的不选主元法，你只需把消元矩阵中的对应元素改变符号后放在单位矩阵I中，就能直接写出L矩阵：

$\large L=\begin{bmatrix} 1 & 0 &0 \\ 2 & 1 &0 \\ 0& 0& 1 \end{bmatrix}\begin{bmatrix} 1 & 0 &0 \\ 0 & 1 &0 \\ -1& 0& 1 \end{bmatrix}\begin{bmatrix} 1 & 0 &0 \\ 0 & 1 &0 \\ 0& -3& 1 \end{bmatrix}=\begin{bmatrix} 1 & 0 &0 \\ 2 & 1 &0 \\ -1& -3& 1 \end{bmatrix}$

下面我用Matlab来验证一下消元过程，这是第一步消元的代码：

%% 不选主元法---松下J27
A=[2 1 1;4 1 0;-2 2 1]
E21=[1 0 0;-2 1 0;0 0 1]
X=E21*A

注意，这一结果和我笔记中第一步消元的结果相同。

继续，第二部消元。

E31=[1 0 0;0 1 0;1 0 1]
X=E31*X

同样，得到相同的结果。

第三步，也是最后一步消元。

E32=[1 0 0;0 1 0;0 3 1]
X=E32*X

可以看的，用矩阵的乘法来表示消元过程所得到的结果和我在笔记本中自己手算得到的结果一致。

同样，用Matlab来验证消元矩阵Z，通过让Z左乘A，即，U=ZxA。看看他是不是能够一次性完成整个消元过程。实际上，这里的U就是LU分解中的U矩阵。此外，我们还可以看看他的逆矩阵，也就是LU分解中的L矩阵。

%消元矩阵的积Z矩阵
Z=E32*E31*E21
U=Z*A
%复原矩阵L就是Z的逆
L=inv(Z)
%Matlab自带的lu分解函数
[l u]=lu(A)

注意，如果用matlab自带的lu分解函数去计算这个矩阵，会得到不一样的结果。如下，这是因为，matlab的函数所使用的是分部主元法，而不是我这里所使用的不选主元法。

以下是方程组Ax=b剩余的求解过程：

事实上，我上面所写的每一步，都是对矩阵的LU分解的一个回顾。

部分主元法(Partial pivoting)

下图很好的说明了部分主元法的计算流程：

1，先在第一列中搜索，找到绝对值最大的元素r1。

2，然后通过行交换，把r1所在的行，放在a11的位置。

3，消除r1（a11）下面的所有元素，第一列消元结束。

4，在第二列中a22以下的元素中搜索，找到绝对值最大的元素r2。

5，通过行交换，把r2放在a22的位置。

6，消除a22下面的所有元素，至此，对于第二列的消元结束。

7，依此类推，直到最后一个元素ann。

注意：我这里说的a11,a22，。。。ann实际上就是主对角线上的元素。（见下图）

下面是我的手写笔记：

请注意，部分主元法与不选主元法的不同之处是，部分主元法进行了行交换，并引入了置换矩阵P。这里我再插播一条广告，不太了解置换矩阵的同学可以看看。线性代数 --- 置换矩阵 - permutation matrix（个人笔记扫描版）_松下J27的博客-CSDN博客_置换矩阵置换矩阵 - permutation matrixhttps://blog.csdn.net/daduzimama/article/details/120494274

计算到这里，已经得到了我们想要的U,以及我们前面消元过程中的全部消元矩阵E。我们还能根据每一步消元的乘数，直接写出还原矩阵L。这种做法，在我的另外一篇LU分解的文章中提过(也就是文中开篇处提到的文章)，但在这里不灵了。以至于我虽然写出了L矩阵，最终的解是错误的。

这里，我对前面的错误进行了总结，由于置换矩阵的引入，无法直接写出L矩阵，最终还是根据还原矩阵L的定义算出了L。即，还原是消元的逆过程，还原矩阵L是消元矩阵Z的逆矩阵。

$\small L=Z^{-1}$

这里我再说一下，由于网上现有的一些资源和资料十分有限，我暂时没找到直接写出L矩阵的方法。但是我总觉得先求Z，再根据Z的逆求得L的过程是十分繁琐的，我想可能存在直接写出L的方法，希望知道的朋友能够在下面的评论区告诉我和大家，万分感谢。

注：上面得笔记中，通过求逆所得到的L，是我用计算机通过matlab算出来的。

用Matlab来验证我们之前的计算结果：

%% 部分主元法PA=LU
A=[2 1 1;4 1 0;-2 2 1]

%先根据消元过程中的各个矩阵求Z（L-1）
E31=[1 0 0;0 1 0;0.5 0 1]
E21=[1 0 0;-0.5 1 0;0 0 1]
X=E31*E21
P32=[1 0 0;0 0 1;0 1 0]
X=P32*X
P21=[0 1 0;1 0 0;0 0 1]
X=X*P21
E32=[1 0 0;0 1 0;0 -0.2 1]
Z=E32*X
%用消元矩阵求得U
U=Z*A
%求出非标准下三角矩阵，记为LL
LL=inv(Z)

%用matlab自带的函数来验证结果
[Matlab_L Matlab_U]=lu(A)

（在matlab中，我用LL来表示非标准下三角矩阵L）

从非标准下三角矩阵到标准下三角矩阵：

分解进行到这里，又要注意了。按照我之前对于LU分解的固有理解，“LU分解是矩阵A的一种三角分解”。也就是说，矩阵A可以用一个标准的下三角矩阵L和一个标准的上三角矩阵U的乘积来表示。既然这里的L和U都是特指的标准三角矩阵，为什么这里我算出来的L不是一个标准三角矩阵，是不是哪里出问题了？

通过仔细观察我们上面得到的非标准三角矩阵L，我们可以看到，只要经几次行交换，就能得到我们期望的标准三角矩阵L。为了表示这两者的区别，并且记录行交换的过程。我这里做了两个操作：

1，把之前得到的矩阵L重命名为L’，因为之前的L是不标准的三角矩阵，所以，我这里用L'表示，并且用L表示经过换行后得到的，也就是我们一直想要的标准下三角矩阵。

2，同时，我这里延续高斯消元的传统，用矩阵P来表示和记录行交换的过程。

这里，我们看到，经过一系列的推导，最终得到了一个重要结论，那就是：

可以把非标准三角矩阵L'变成标准下三角矩阵L的置换矩阵P，正是我们苦苦寻求的PA=LU分解中的矩阵P。

用矩阵来表示就是：

$\small PL'=L, PA=LU$

最终，按照部分主元法，求得了举世闻名的LU分解的另一种表达形式 ---> PA=LU。这里我们还要提一下，用PA=LU的来求解方程组的过程和用A=LU的求解过程，有一些细微差异。 A=LU的求解过程可以概括为，先用正向回代求解Ly=b，得到y。然后再用反向回代求解Ux=y，得到x。而，PA=LU的求解过程变成，先用正向回代求解Ly=Pb,然后再用反向回代求解Ux=y。

结合上面已有的code，最终我们可以用matlab自带的LU分解函数来验证我们的计算结果。这里我用P4L表示PA=LU中的P。

%用一系列的置换矩阵P把LL变成标准型的L
P31=[0 0 1;0 1 0;1 0 0]
P32=[1 0 0;0 0 1;0 1 0]
P4L=P31*P32
L=P4L*LL
%至此，已经求出了PA=LU分解中的全部矩阵P,L,U
%但是要证明PA中的P和把LL变成L的P是同一个置换矩阵

% P4L*A
% L*U
%用matlab自带的函数来验证结果
L
U
P4L
[l u p]=lu(A)

补充: 手动计算非标准三角矩阵L'的方法

前面说过，在使用部分主元法对矩阵A进行消元的过程中。由于出现了行交换，导致无法根据每次消元所乘的倍数（或根据消元矩阵E）直接写出L矩阵，后面有高人在评论区回复了一种方法。我发现，他所提供的方法虽然不能直接写出标准的下三角矩阵L，但是，依然可以在“进行了行交换”的基础上直接写出非标准三角矩阵L'。

首先，因为基于我们之前得出的结论：

$\small L=Z^{-1}$

即，还原矩阵L，是多个消元矩阵的乘积Z矩阵的逆(或者说是逆过程更好)。根据逆矩阵的定义，我们有：

$\large ZZ^{-1}=ZL=I$

因此，我们可以得到如下公式：

$\large ZL=E_{32}P_{32}E_{31}E_{21}P_{21}\cdot P_{21}^{-1}E_{21}^{-1}E_{31}^{-1}P_{32}^{-1}E_{32}^{-1}=I$

上方公式中，从点乘符号“ $\LARGE \mathbf{\cdot}$ ”开始，每一次相乘都会得到单位矩阵I，最终结果就是n个单位矩阵相乘，最终得到一个单位矩阵I。

其中：

$\large Z=E_{32}P_{32}E_{31}E_{21}P_{21}$ ， $\large L={P_{21}}^{-1}{E_{21}}^{-1}{E_{31}}^{-1}{P_{32}}^{-1}{E_{32}}^{-1}$

又因为，置换矩阵P的逆等于他自己，我们把L改写成如下形式：

$\large L=P_{21}{E_{21}}^{-1}{E_{31}}^{-1}P_{32}{E_{32}}^{-1}$

我们基于上面的公式逐步求出L(非标准下三角矩阵)：

(注意: 这里是不需要用到计算器的，只需要知道如何写出消元矩阵E的逆和如何进行行交换就行)

1，先写 ${E_{32}}^{-1}$ ，把 $E_{32}$ 中的元素{3,2}(第三行，第二列)改变符号后放在相同尺寸单位矩阵中对应的位置即可。

$\large {E_{32}^{-1}}=\begin{bmatrix} 1 &0 &0 \\ 0& 1 &0 \\ 0& 0.2 &1 \end{bmatrix}$

2, 交换 ${E_{32}}^{-1}$ 的2，3行

$\large P_{32}E{_{32}}^{-1}=\begin{bmatrix} 1 &0 &0 \\ 0& 0.2 &1 \\0& 1 &0 \end{bmatrix}$

3, 把 $E_{31}$ 中的元素{3，1}改变符号后，放在上一步计算结果对应的位置上。

$\large E{_{31}}^{-1}P_{32}E{_{32}}^{-1}=\begin{bmatrix} 1 &0 &0 \\ 0& 0.2 &1 \\-0.5& 1 &0 \end{bmatrix}$

4，把 $E_{21}$ 中的元素{2，1}改变符号后，放在上一步计算结果对应的位置上。

$\large E{_{21}}^{-1}E{_{31}}^{-1}P_{32}E{_{32}}^{-1}=\begin{bmatrix} 1 &0 &0 \\ 0.5& 0.2 &1 \\-0.5& 1 &0 \end{bmatrix}$

5，交换上一步计算结果的12行

$\large P_{21}E{_{21}}^{-1}E{_{31}}^{-1}P_{32}E{_{32}}^{-1}=\begin{bmatrix} 0.5& 0.2 &1\\ 1 &0 &0 \\-0.5& 1 &0 \end{bmatrix}=L$

这与上文中，通过Matlab求解消元矩阵Z的逆矩阵，所得到的非标准下三角矩阵L'(即，上文中的L'矩阵或matlab中的LL矩阵)的结果一模一样。

我们依然要通过行交换把L'变成标准的下三角矩阵L。在这个变换的过程中，我们所使用的置换矩阵P就是，PA=LU分解中的矩阵P。

这样一来就能极大的缩减L矩阵的运算量和运算时间，即使是编程实现，也可以使用这种“取巧”的办法。而这里所说的“取巧”实际上只不过是充分利用了P矩阵和E矩阵的一些性质。

在这里我要特别感谢网友voodoocjl的回复和帮助，万分感谢，解答了我长久以来的一个疑惑。

完全主元法(Full pivoting)

完全主元法和部分主元法的最大不同之处在于，部分主元法是按列选择绝对值最大的元素，而完全主元法是在“整个矩阵”中去寻找。只是这里所说的“整个矩阵”的搜索范围，会随着消元的进行，越来越小。

下面我详细说一下：

1，在整个矩阵中找到绝对值最大的元素r1。

2，通过行交换和列交换，把r1放在a11的位置。（注意：这里增加了列交换）

3，对于第一例，令r1为主元，消除a11下面的所有元素。

4，缩小搜索范围，在除了第一行和第一列的矩阵中，找到绝对值最大的元素r2。

5，通过行交换和列交换，把r2放在a22的位置。

6，令r2为第二列的主元，消除a22下面的所有元素。

7，依此类推，直到ann。

下面是我手写的个人笔记：

为了要把最大的元素9换到a11的位置上去，除了使用行置换矩阵P，我们首次引入了列置换矩阵Q。这里，我要重复之前的一个概念《向量与矩阵的乘法》，那就是，如果要对一个矩阵进行行操作，也就是行与行之间的线性组合，需要用一个行向量v左乘这个矩阵A，即，vA。而如果要对一个矩阵进行列操作，也就是列与列之间的线性组合，需要用一个列向量v右乘这个矩阵A，即Av。当时我写了一个口诀，叫：“前（左）乘行，行操作。后（右）乘列，列操作。”

现在看来如果要用行置换矩阵P和列置换矩阵Q对矩阵进行操作，依然可以用我们前面的口诀，只是这里变成了矩阵和矩阵相乘，就不存在什么行向量和列向量了。同时，还要牢记，我们是怎么构建行置换矩阵P的，也用同样的方式构建置换矩阵Q。即，置换矩阵就是改变顺序的单位矩阵。

笔记中有误：上面笔记中最后一个等式的左边中A的左边并不等不之前部分消元法的Z。在之前的部分消元法中Z=E32 P32 E31 E21 P21，而这里A的左边是E32 P32 E31 E21 P31。

注意：这里我们得到的LL矩阵(这是用matlab对Z求逆算出来的)，依然是我们之前提到的非标准下三角矩阵。

我们依然要通过行交换把LL变成标准的下三角矩阵。在这个变换过程中，我们记录下来的置换矩阵P就是，PAQ=LU分解中的矩阵P。（在部分主元法的笔记中，我把非标准三角矩阵L，重命名为L'，并进行了推导得到PL'=L。这里，我把非标准三角矩阵L，直接命名为LL，得到PLL=L。）

最后，同样，我们可以用matlab来验证我们最终的计算结果。可惜的是matlab自带函数的PAQ=LU分解，只针对稀疏矩阵。但我们可以用他的PA=LU分解函数验证，只是这里，我们要把A矩阵变成Q*A。

%% 完全主元法，PAQ=LU
A=[2 3 4;4 7 5;4 9 5]
%逐步消元直到U
Q21=[0 1 0;1 0 0;0 0 1]
P31=[0 0 1;0 1 0;1 0 0]
X=A*Q21
X=P31*X
E21=[1 0 0;-7/9 1 0;0 0 1]
X=E21*X
E31=[1 0 0;0 1 0;-1/3 0 1]
X=E31*X
Q32=[1 0 0;0 0 1;0 1 0]
P32=[1 0 0;0 0 1;0 1 0]
X=X*Q32
X=P32*X
E32=[1 0 0;0 1 0;0 -10/21 1]
X=E32*X

%用ZAQ表示消元过程
%E32P32E31E21P31AQ21Q32=U
Q=Q21*Q32
Z=E32*P32*E31*E21*P31
%ZAQ=U
Z*A*Q
%求出非标准下三角矩阵，记为LL
LL=inv(Z)

%用一系列的置换矩阵P把LL变成标准型的L
P21=[0 1 0;1 0 0;0 0 1]
P21*LL
P31=[0 0 1;0 1 0;1 0 0]
P4L=P31*P21
L=P4L*LL

[l u p]=lu(A*Q)

补充：

使用完全主元法对稀疏矩阵（sparse matrix）进行LU分解时，能极大的减少计算内存资源的利用。

也就是说，如果我们把用完全主元法对稀疏矩阵做lu分解得到的LU矩阵和部分主元法所求得的LU矩阵相比，前者的结果中会有很多0元素（较少占用计算机资源），而后者则不然。

例1：这是MATLAB自带的例子

已知一个479x479的稀疏矩阵，并且我们用spy函数画出这个矩阵。其中，nz表示的是非零元素的个数。