在科学研究和工程技术中，许多问题都可以归结为求解如下的线性方程组:
$$
\left\{\begin{aligned}
a_{11} x_{1}+a_{12} x_{2}+\cdots+a_{1 n} x_{n} &=b_{1} \\
a_{21} x_{1}+a_{22} x_{2}+\cdots+a_{2 n} x_{n} &=b_{2} \\
\cdots \cdots & \\
a_{n 1} x_{1}+a_{n 2} x_{2}+\cdots+a_{n n} x_{n} &=b_{n}
\end{aligned}\right.
$$
方程组常写为矩阵形式，即
$$
Ax=b
$$
其中
$$
\begin{gathered}
A=\left(\begin{array}{cccc}
a_{11} & a_{12} & \cdots & a_{1 n} \\
a_{21} & a_{22} & \cdots & a_{2 n} \\
\vdots & \vdots & & \vdots \\
a_{n 1} & a_{n 2} & \cdots & a_{n n}
\end{array}\right) \\
x=\left(x_{1}, x_{2}, \cdots, x_{n}\right)^{\mathrm{T}}, \quad b=\left(b_{1}, b_{2}, \cdots, b_{n}\right)^{\mathrm{T}} .
\end{gathered}
$$
从求解线性方程组问题提出以来，人们便不断地探索研究高效实用的求解算法。本篇博文，
从一个计算需要30万年的方法说起
克莱姆法则，又译克拉默法则（Cramer's Rule）是线性代数中一个关于求解线性方程组的定理。它适用于变量和方程数目相等的线性方程组，是瑞士数学家克莱姆（1704-1752）于1750年，在他的《线性代数分析导言》中发表的。其实莱布尼兹〔1693〕，以及马克劳林〔1748〕亦知道这个法则，但他们的记法不如克莱姆。（摘自百度百科）
即对于系数矩阵的行列式$|A|\neq0$ 的情况下，方程组的解存在且唯一，由克拉默法则可知方程组$Ax=b$ 的解为
$$
x_{i}=\frac{\left|A_{i}\right|}{|A|}, \quad i=1,2, \cdots, n
$$
其中$A_i(i=1,2,cdots,n)$ 为系数矩阵$A$ 的第$i$ 列用方程组的右端项$b$ 替换后所得到的矩阵。用克拉默方法求解$n$ 阶线性方程组，乘除法的运算量为$(n+1)!(n-1)+n$。当$n=20$时，此数约为$9.7\times10^{20}$. 用每秒钟作$10^8$ 次乘除法运算的计算机求解这个方程组，大约需要30.8万年。这说明即使在求解阶数较小的方程组时，该方法也是不能使用的。鉴于这个原因，人们提出了许多求解线性方程组的数值解法，而且至今仍在不断地研究高效实用的新算法。
高斯消去法（高斯：为什么要消去我！）
算法原理
高斯消去法，本质上就是初等数学中的消元法，只不过是将求解步骤规范化，便于编写计算机程序 。高斯消去法首先是将方程组进行消元运算，将其化为一个等价的同解的上三角方程组（这一过程称为消元过程），然后通过求解上三角方程组得到原方程组的解（这一过程称之为回代）。(挖坑多线程的高斯消元法)
高斯消去法可以算得上计算方法这门课里面接触的第一个深入讨论的算法，算法符合尚且学过初等数学的我们的常规认识，然后又在该认识的基础上，介绍可以被计算机理解并执行的算法，可以说高斯消去法作为众多及算法方法书介绍的第一个算法是非常合适的。让我们再说回到算法上。当系数举证的各阶顺序主子式都不为0时，高斯消去法可以顺利进行，下面用一般式介绍高斯消去法消元的一般过程。
第1步 假定$|A|\neq0$ （$|A|\neq0$ 的情况下可通过换行保证新的$a_{11}^{(0)}\neq0$），将第1个方程保留不动，将$2\sim n$个方程中$x_1$的系数全部消为零，即第$i$个方程减去第1个方程乘以$\frac{a_{i1}^{(0)}}{a_{11}^{(0)}}(i=2,3,\cdots,n)$得
$$
\left\{\begin{aligned}
a_{11}^{(0)} x_{1}+& a_{12}^{(0)} x_{2}+a_{13}^{(0)} x_{3}+\cdots+a_{1 n}^{(0)} x_{n}=b_{1}^{(0)} \\
& a_{22}^{(1)} x_{2}+a_{23}^{(1)} x_{3}+\cdots+a_{2 n}^{(1)} x_{n}=b_{2}^{(1)} \\
& a_{32}^{(1)} x_{2}+a_{33}^{(1)} x_{3}+\cdots+a_{3 n}^{(1)} x_{n}=b_{3}^{(1)} \\
&\qquad\qquad\cdots \cdots & \\
& a_{n 2}^{(1)} x_{2}+a_{n 3}^{(1)} x_{3}+\cdots+a_{n n}^{(1)} x_{n}=b_{n}^{(1)}
\end{aligned}\right.
$$
其中
$$
\begin{aligned}
&a_{i j}^{(1)}=a_{i j}^{(0)}-\frac{a_{i 1}^{(0)}}{a_{11}^{(0)}} a_{1 j}^{(0)}=a_{i j}^{(0)}-l_{i 1} a_{1 j}^{(0)}, \quad i=2,3, \cdots, n ; j=2,3, \cdots, n, \\
&b_{i}^{(1)}=b_{i}^{(0)}-\frac{a_{i 1}^{(0)}}{a_{11}^{(0)}} b_{1}^{(0)}=b_{i}^{(0)}-l_{i 1} b_{1}^{(0)}, \quad i=2,3, \cdots, n, \\
&l_{i 1} \triangleq \frac{a_{i 1}^{(0)}}{a_{11}^{(0)}}, \quad i=2,3, \cdots, n .
\end{aligned}
$$
将方程组（5）改写为矩阵形式，即
$$
A^{(1)}x=b^{(1)}
$$
其中
$$
A^{(1)}=\left(\begin{array}{ccccc}
a_{11}^{(0)} & a_{12}^{(0)} & a_{13}^{(0)} & \cdots & a_{1 n}^{(0)} \\
0 & a_{22}^{(1)} & a_{23}^{(1)} & \cdots & a_{2 n}^{(1)} \\
0 & a_{32}^{(1)} & a_{33}^{(1)} & \cdots & a_{3 n}^{(1)} \\
\vdots & \vdots & \vdots & & \vdots \\
0 & a_{n 2}^{(1)} & a_{n 3}^{(1)} & \cdots & a_{n n}^{(1)}
\end{array}\right), \quad b^{(1)}=\left(\begin{array}{c}
b_{1}^{(0)} \\
b_{2}^{(1)} \\
b_{3}^{(1)} \\
\vdots \\
b_{n}^{(1)}
\end{array}\right)
$$
第2步 与第一步类似，保留方程组1，2个方程不动，将$3\sim n$个方程中$x_2$的系数消为零，即第$i$个方程减去第2个方程乘以$\frac{a_{i2}^{(0)}}{a_{22}^{(1)}}(i=3,4,\cdots,n)$得
$$
\left\{\begin{aligned}
a_{11}^{(0)} x_{1}+& a_{12}^{(0)} x_{2}+a_{13}^{(0)} x_{3}+\cdots+a_{1 n}^{(0)} x_{n}=b_{1}^{(0)} \\
& a_{22}^{(1)} x_{2}+a_{23}^{(1)} x_{3}+\cdots+a_{2 n}^{(1)} x_{n}=b_{2}^{(1)} \\
& a_{32}^{(1)} x_{2}+a_{33}^{(1)} x_{3}+\cdots+a_{3 n}^{(1)} x_{n}=b_{3}^{(1)} \\
&\qquad\qquad\cdots \cdots & \\
& a_{n 2}^{(1)} x_{2}+a_{n 3}^{(1)} x_{3}+\cdots+a_{n n}^{(1)} x_{n}=b_{n}^{(1)}
\end{aligned}\right.
$$
其中
$$
\begin{aligned}
&a_{i j}^{(2)}=a_{i j}^{(1)}-\frac{a_{i 2}^{(1)}}{a_{22}^{(1)}} a_{2 j}^{(0)}=a_{i j}^{(1)}-l_{i 2} a_{2 j}^{(1)}, \quad i=3,4, \cdots, n ; j=3,4, \cdots, n, \\
&b_{i}^{(2)}=b_{i}^{(1)}-\frac{a_{i 1}^{(1)}}{a_{22}^{(1)}} b_{2}^{(1)}=b_{i}^{(1)}-l_{i 2} b_{1}^{(0)}, \quad i=3,4, \cdots, n, \\
&l_{i 2} \triangleq \frac{a_{i 2}^{(1)}}{a_{22}^{(1)}}, \quad i=3,4, \cdots, n .
\end{aligned}
$$
将方程组（9）改写为矩阵形式，即
$$
A^{(2)}x=b^{(2)}
$$
其中
$$
A^{(1)}=\left(\begin{array}{ccccc}
a_{11}^{(0)} & a_{12}^{(0)} & a_{13}^{(0)} & \cdots & a_{1 n}^{(0)} \\
0 & a_{22}^{(1)} & a_{23}^{(1)} & \cdots & a_{2 n}^{(1)} \\
0 & 0 & a_{33}^{(2)} & \cdots & a_{3 n}^{(2)} \\
\vdots & \vdots & \vdots & & \vdots \\
0 & 0 & a_{n 3}^{(2)} & \cdots & a_{n n}^{(2)}
\end{array}\right), \quad b^{(2)}=\left(\begin{array}{c}
b_{1}^{(0)} \\
b_{2}^{(1)} \\
b_{3}^{(2)} \\
\vdots \\
b_{n}^{(2)}
\end{array}\right)
$$
重复以上操作，进行$n-1$次消元后，原方程组可以化为如下等价的同解上三角方程组：
$$
A^{(n-1)}x=b^{(n-1)}
$$
其中
$$
A^{(n-1)}=\left(\begin{array}{ccccccc}
a_{11}^{(0)} & a_{12}^{(0)} & \cdots & a_{1 k}^{(0)} & a_{1, k+1}^{(0)} & \cdots & a_{1 n}^{(0)} \\
& a_{22}^{(1)} & \cdots & a_{2 k}^{(1)} & a_{2, k+1}^{(1)} & \cdots & a_{2 n}^{(1)} \\
& & \ddots & \vdots & \vdots & & \vdots \\
& & & a_{k k}^{(k-1)} & a_{k, k+1}^{(k-1)} & \cdots & a_{k n}^{(k-1)} \\
& & & & a_{k+1, k+1}^{(k)} & \cdots & a_{k+1, n}^{(k)} \\
& & & & &\ddots & \vdots \\
& & & & & &a_{n n}^{(n-1)}
\end{array}\right), b^{(n-1)}=\left(\begin{array}{l}
b_{1}^{(0)} \\
b_{k}^{(1)} \\
b_{k+1}^{(k-1)} \\
\vdots \\
b_{n}^{(n-1)}
\end{array}\right)
$$
下面介绍回代过程，使用初等数学的方法，从最后一个方程开始依次求得$x_n,x_{n-1},\cdots,x_1$，具体公式如下
$$
\left\{\begin{aligned}
x_{n}=& \frac{b_{n}^{(n-1)}}{a_{n n}^{(n-1)}} \\
x_{k}=& \frac{b_{k}^{(k-1)}-\sum_{j=k+1}^{n} a_{k j}^{(k-1)} x_{j}}{a^{(k-1)}}, \quad k=n-1, n-2, \cdots, 2,1
\end{aligned}\right.
$$
算法实现
为了便于高斯消去法在计算机上的实现以及减少计算机的储存量，在开始计算时将方程组的系数矩阵存放在二维数组$A$中，方程的右端项存放在一维数组$b$中，消元过程中各步被消除的元用消元相乘的系数代替，以节省存储空间。消元过程结束后，二维数组$A$和一维数组$b$的内容如下：
$$
A=\left(\begin{array}{cccccc}
a_{11}^{(0)} & a_{12}^{(0)} & a_{13}^{(0)} & \cdots & a_{1, n-1}^{(0)} & a_{1 n}^{(0)} \\
l_{21} & a_{22}^{(1)} & a_{23}^{(1)} & \cdots & a_{2, n-1}^{(1)} & a_{2 n}^{(1)} \\
l_{31} & l_{32} & a_{33}^{(2)} & \cdots & a_{3, n-1}^{(2)} & a_{3 n}^{(2)} \\
\vdots & \vdots & \vdots & & \vdots & \vdots \\
l_{n 1} & l_{n 2} & l_{n 3} & \cdots & l_{n, n-1} & a_{n n}^{(n-1)}
\end{array}\right), \quad b=\left(\begin{array}{c}
b_{1}^{(0)} \\
b_{2}^{(1)} \\
b_{3}^{(2)} \\
\vdots \\
b_{n}^{(n-1)}
\end{array}\right)
$$
高斯消元法对应matlab代码如下：
function x=ge(A,b,options)
% solve a linear system using Gaussian Elimination
% specify options to be one of the following:
% ge: naive Gaussian Elimination
% pp: Gaussian Elimination with partial pivoting
% spp: Gaussian Elimination with scaled partial pivoting
% Step 1: Gaussian Elimination with scaled partial pivoting.
% Overwrite the matrix.
n=length(b);
p = (1:n)’; % initialize the pivoting vector
if strcmp(options,’pp’)
s=ones(1,n);
elseif strcmp(options,’spp’)
s = max(abs(A’));% compute the scale of each row
end
for k = 1:(n-1)
if strcmp(options,’pp’) || strcmp(options,’spp’)
r = abs(A(p(k),k)/s(p(k)));
kp = k;
for i = (k+1):n
t = abs(A(p(i),k)/s(p(i)));
if t > r, r = t; kp = i; end
end
l = p(kp); p(kp) = p(k); p(k) = l; % interchange p(kp) and p(k)
end
for i = (k+1):n
A(p(i),k) = A(p(i),k)/A(p(k),k);
for j = (k+1):n
A(p(i),j) = A(p(i),j)-A(p(i),k)*A(p(k),j);
end
end
end
%p % output the pivoting vector
%A % output the overwritten matrix
% Step 2: Forward substitution to solve L*y = b, where
% L(i,j) = 0 for j>i; L(i,i) = 1; L(i,j) = A(p(i),j) for i>j.
y = zeros(n,1); % initialize y to be a column vector
y(1) = b(p(1));
for i = 2:n
y(i) = b(p(i));
for j = 1:(i-1)
y(i) = y(i)-A(p(i),j)*y(j);
end
end
%y % output y
% Step 3: Back substitution to solve U*x = y
% U(i,j) = A(p(i),j) for j>=i; U(i,j) = 0 for i>j.
x = zeros(n,1); % initialize x to be a column vector
x(n) = y(n)/A(p(n),n);
for i = (n-1): -1 : 1
x(i) = y(i);
for j = (i+1):n
x(i) = x(i)-A(p(i),j)*x(j);
end
x(i) = x(i)/A(p(i),i);
end
end列主元高斯消去法
在高斯消去法的消元过程中主元$a_{kk}^{(k-1)}=0$将会导致无法消元或者$\left|a_{kk}^{(k-1)}\right|$很小导致其他元素数量级剧增和舍入误差的增加，在消去的每一步进行按列选主元的操作，故被称为列主元高斯消去法。列主元消去法代码如下：
%求解线性方程组AX=b;
function x=Column_Gauss(A,b)
%准备工作：获取方程组的信息；
[m,n]=size(A); %获取矩阵的行和列；
x=zeros(n,1);
%按列选主元；
for k=1:n-1 %因为在主元（主对角线上的元素），所以k=行号=列号；--理解为列号；
max=abs(A(k,k)); %先假设对角元的元素均为最大值；
r=k; %记录此时的行号；
for i=k+1:n %比较同一列其他元素值与主对角元元素的大小；i是行号；
if max<abs(A(i,k))
max=abs(A(i,k));
r=i; %更新同一列最大值元素所在的行号；
end
end
%换行；要将主对角元的行的元素全部与列主元最大元素所在行的元素交换；
if r~=k
for c=k:n %只需交换同一行中对角元列号之后的元素；--理解为列号；
a(c)=A(k,c); %将原主对角元元素所在位置之后的行元素存到a中；--把位置空出来；
A(k,c)=A(r,c); %将列主元最大元素所在行的其他元素交换到主对角元元素所在行的想对应的位置上；--交换
A(r,c)=a(c); %将存在a中的主对角元的元素放到原来列主元最大元素所在行对应的位置上；--把空出来的元素放到移出元素的行的位置；---终于实现交换了！！！
end %对应的b也要按照行的交换而交换；
d=b(k); %把原主对角元所在行的b移出来；
b(k)=b(r); %再把最大元素所在行对应的b分量放到主对角元的位置；
b(r)=d; %把原来移出来的元素放到列主元最大元素所在行对应的b的位置；---实现b的交换；
end
%消元；就是把主对角元下面的元素全部划为0---目的是要系数矩阵A化为一个上三角的系数矩阵；
for i=k+1:n
for j=i:n
A(i,j)=A(i,j)-A(i,k)*A(k,j)/A(k,k);
end
b(i)=b(i)-A(i,k)*b(k)/A(k,k);
A(i,k)=0;
end
%回代求解；此时已经化为上三角型的系数矩阵了；主要从下往上代入求解即可；
x(n)=b(n)/A(n,n);
for i=n-1:(-1):1
sum=0;
for j=i+1:n
sum=sum+A(i,j)*x(j);
end
x(i)=(b(i)-sum)/A(i,i);
end
end矩阵的三角分解法
在实际问题中，经常遇到系数矩阵$A$相同，但右端项不同的多个线性方程组$Ax=b^{(i)}(i=1,2,\cdots,m)$. 对高斯消去法作进一步的输入分析，由矩阵运算到处矩阵的三角分解法。利用矩阵的三角分解求解上述的一类线性方程组，以达到减少运算量的目的。其次，针对具有特殊系数矩阵的线性方程组，由矩阵的三角分解建立一些特殊的求解方法。
从矩阵运算的观点来看，高斯消去法的消元过程实质上是将增广矩阵$(A^{(0)},b^{(0)})$通过逐步左乘一些列的初等下三角矩阵$L_k(k=1,2,\cdots,n-1)$最终变为矩阵$(A^{(n-1)},b^{(n-1)})$.
一般地，
$$
\begin{aligned}
&A^{(k-1)} \triangleq\left(\begin{array}{ccccc}
a_{11}^{(0)} & a_{12}^{(0)} & & \ldots & \cdots & a_{1 n}^{(0)} \\
& a_{22}^{(1)} & & \cdots & \cdots & a_{2 n}^{(1)} \\
& & \ddots & & & \vdots \\
& & & a_{k k}^{(k-1)} & \cdots & a_{k n}^{(k-1)} \\
& & & & \vdots \\
& & & & a_{n k}^{(k-1)} & \cdots & a_{n n}^{(k-1)}
\end{array}\right)\\
&\qquad L_{k} \triangleq\left(\begin{array}{cccccc}
1 & & & & & \\
& \ddots & & & & \\
& & 1 & & & \\
& & -l_{k+1, k} & 1 & & \\
& & -l_{k+2, k} & & 1 & \\
& & \vdots & & \ddots \\
& & -l_{n k} & & & 1
\end{array}\right)\\
&\quad A^{(k)} \triangleq\left(\begin{array}{cccccc}
a_{11}^{(0)} & a_{12}^{(0)} & & \ldots & a_{1, k+1}^{(0)} & \cdots & a_{1 n}^{(0)} \\
& a_{22}^{(1)} & & \cdots & a_{2, k+1}^{(1)} & \cdots & a_{2 n}^{(1)} \\
& & \ddots & & \vdots & & \vdots \\
& & & a_{k k}^{(k-1)} & a_{k, k+1}^{(k-1)} & \cdots & a_{k n}^{(k-1)} \\
& & & & a_{k+1, k+1}^{(k)} & \cdots & a_{k+1, n}^{(k)} \\
& & & & \vdots & & \vdots \\
& & & & a_{n, k+1}^{(k)} & \cdots & a_{n n}^{(k)}
\end{array}\right)
\end{aligned}
$$
Comments NOTHING