gaussjordan消元法

（Gauss-Jordan）消元法选主元的高斯-约当（Gauss-Jordan）消元法在很多地方都会用到，例如求一个矩阵的逆矩阵、解线性方程组（插一句：LM算法求解的一个步骤），等等。

它的速度不是最快的，但是它非常稳定（来自网上的定义：一个计算方法，如果在使用此方法的计算过程中，舍入误差得到控制，对计算结果影响较小，称此方法为数值稳定的），同时它的求解过程也比较清晰明了，因而人们使用较多。

下面我就用一个例子来告诉你Gauss-Jordan法的求解过程吧。

顺便再提及一些注意事项以及扩展话题。

对本文中所提到的“主元”等概念的解释，可以参考此链接。

假设有如下的方程组：写成矩阵形式就是：AX=B，其中：且X=(X1, X2, X3)T。

文章来源：/现对矩阵A作初等变换，同时矩阵B也作同样的初等变换，则当A化为单位矩阵的时候，有：显而易见，我们得到了方程组的解X=(1, 2, 4)T。

所以，我们要以一定的策略，对A和B施以一系列的初等变换，当A化为单位矩阵的时候，B就为方程组的解。

选主元的G-J消元法通过这样的方法来进行初等变换：在每一个循环过程中，先寻找到主元，并将主元通过行变换（无需列变换）移动到矩阵的主对角线上，然后将主元所在的行内的所有元素除以主元，使得主元化为1；然后观察主元所在的列上的其他元素，将它们所在的行减去主元所在的行乘以一定的倍数，使得主元所在的列内、除主元外的其他元素化为0，这样就使得主元所在的列化为了单位矩阵的形式。

这就是一个循环内做的工作。

然后，在第二轮循环的过程中，不考虑上一轮计算过程中主元所在的行和列内的元素，在剩下的矩阵范围内寻找主元，然后（如果其不在主对角线上的话）将其移动到主对角线上，并再次进行列的处理，将列化为单位矩阵的形式。

余下的步骤依此类推。

具体的计算过程的一个例子，请看下面我举的求逆矩阵的过程。

如果要解系数矩阵相同、右端向量不同的N个方程组，在设计程序的时候，没有必要”解N 次方程组“，我们完全可以在程序中，将所有的右端向量以矩阵的数据结构（类似于二维数组）来表示，在系数矩阵作行变换的时候，矩阵里的每一个右端向量也做同样的变换，这样，我们在一次求解运算的过程中，实际上就是同时在解N个方程组了，这是要注意的地方。

列主元素消去法列主元素消去法（Gauss-Jordan 消元法）是一种线性代数中常用的消元方法，用于求解线性方程组的解。

这种方法的基本思想是，将线性方程组的增广矩阵通过一系列的初等变换，化为一个阶梯矩阵或行简化阶梯矩阵，从而得到线性方程组的解。

具体步骤如下：构造增广矩阵，即将系数矩阵和常数矩阵组合成一个矩阵。

将增广矩阵转化为一个上三角矩阵（也叫阶梯矩阵）。

反向消元，将阶梯矩阵转化为一个行简化阶梯矩阵。

根据简化矩阵求解方程组。

这种方法的优点是计算简单、容易理解，且可避免误差的积累。

但是，如果矩阵的规模较大，运算量会很大，计算时间较长。

此时可以使用更高效的算法，如LU分解、QR分解等。

假设有一个 $n$ 个未知量和 $n$ 个方程的线性方程组，可以写成矩阵形式如下：$Ax = b$其中，$A$ 是一个 $n \times n$ 的系数矩阵，$x$ 是一个 $n \times 1$ 的未知量向量，$b$ 是一个 $n \times 1$ 的常数向量。

为了求解 $x$，可以将方程组的增广矩阵表示如下：$\begin{bmatrix} a_{11} & a_{12} & a_{13} & \cdots & a_{1n} & b_{1} \ a_{21} & a_{22} & a_{23} & \cdots & a_{2n} & b_{2} \ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \ a_{n1} & a_{n2} & a_{n3} & \cdots & a_{nn} & b_{n} \end{bmatrix}$ 其中，$a_{ij}$ 表示矩阵的第 $i$ 行第 $j$ 列的元素。

高斯消元法，是线性代数中的一个算法，可用来求解线性方程组，并可以求出矩阵的秩，以及求出可逆方阵的逆矩阵。

高斯消元法的原理是：若用初等行变换将增广矩阵化为，则AX = B与CX = D是同解方程组。

所以我们可以用初等行变换把增广矩阵转换为行阶梯阵，然后回代求出方程的解。

以上是线性代数课的回顾，下面来说说高斯消元法在编程中的应用。

首先，先介绍程序中高斯消元法的步骤：(我们设方程组中方程的个数为equ，变元的个数为var，注意：一般情况下是n个方程，n个变元，但是有些题目就故意让方程数与变元数不同)1. 把方程组转换成增广矩阵。

2. 利用初等行变换来把增广矩阵转换成行阶梯阵。

枚举k从0到equ – 1，当前处理的列为col(初始为0) ，每次找第k行以下(包括第k行)，col列中元素绝对值最大的列与第k行交换。

如果col列中的元素全为0，那么则处理col + 1列，k不变。

3. 转换为行阶梯阵，判断解的情况。

①无解当方程中出现(0, 0, …, 0, a)的形式，且a != 0时，说明是无解的。

②唯一解条件是k = equ，即行阶梯阵形成了严格的上三角阵。

利用回代逐一求出解集。

③无穷解。

条件是k < equ，即不能形成严格的上三角形，自由变元的个数即为equ – k，但有些题目要求判断哪些变元是不缺定的。

这里单独介绍下这种解法：首先，自由变元有var - k个，即不确定的变元至少有var - k个。

我们先把所有的变元视为不确定的。

在每个方程中判断不确定变元的个数，如果大于1个，则该方程无法求解。

如果只有1个变元，那么该变元即可求出，即为确定变元。

以上介绍的是求解整数线性方程组的求法，复杂度是O(n3)。

浮点数线性方程组的求法类似，但是要在判断是否为0时，加入EPS，以消除精度问题。

高斯消元法简介在信息学竞赛中，很多问题都可以转化成线性方程组或者与之相关的问题。

因此，我们需要了解线性方程组的各种解法。

数值分析1顺序消去法、列主元、列主元Gauss-Jordan消去法function x = Gauss (A, b)% 求解方程组的Gauss消去法，调用方法为% x = Gauss (A, b)% 其中% A 为方程组的系数矩阵，b为方程组的右端项% x 为方程组的解[n,m] = size (A); nb = length (b);if n~=merror ('% 系数矩阵必须为方的！');endif m~=nberror ('% b 的维数与方程组的行数不匹配！'); endfor k = 1:n-1% 消元过程for i = k+1:nm = A (i,k)/A(k,k);for j = k+1:nA (i,j) = A (i,j)-m*A (k,j);endb (i) = b (i)-m*b (k);endendx=zeros (size (b));for k = n:-1:1for j = k+1:nb (k) = b (k)-A (k,j)*x (j);endx (k) = b (k)/A(k,k);endendfunction x = Gauss_Elim (A, b)% 求解方程组的列主元Gauss消去法，调用方法为% x = Gauss_Elim (A, b)% 其中% A 为方程组的系数矩阵，b为方程组的右端项% x 为方程组的解[n,m] = size (A); nb = length (b);error ('% 系数矩阵必须为方的！');endif m~=nberror ('% b 的维数与方程组的行数不匹配！');endfor k = 1:n-1% 选主元a_max = 0;for i = k:nif abs (A (i,k))>a_maxa_max = A (i,k); r=i;endendif abs(a_max)<1e-15error ('% 系数矩阵奇异，无法求解方程组！');end% 交换两行if r>kfor j = k:nz=A (k,j); A (k,j)=A (r,j);A (r,j)=z;endz=b (k);b (k)=b (r);b (r)=z;end% 消元过程for i = k+1:nm = A (i,k)/A(k,k);for j = k+1:nA (i,j) = A (i,j)-m*A (k,j);endb (i) = b (i)-m*b (k);endend% 回代过程if abs (A (n,n))<1e-15error ('% 系数矩阵奇异，无法求解方程组！'); endx=zeros (size (b));for k = n:-1:1for j = k+1:nb (k) = b (k)-A (k,j)*x (j);endx (k) = b (k)/A(k,k);endendfunction x = Gauss_Jordan (A, b)% 求解方程组的列主元Gauss-Jordan消去法，调用方法为% x = Gauss_Jordan (A, b)% 其中% A 为方程组的系数矩阵，b为方程组的右端项% x 为方程组的解[n,m] = size (A); nb = length (b);error ('% 系数矩阵必须为方的！');endif m~=nberror ('% b 的维数与方程组的行数不匹配！'); endfor k = 1:n% 选主元a_max = 0;for i = k:nif abs (A (i,k))>a_maxa_max = A (i,k); r=i;endendif abs(a_max)<1e-15error ('% 系数矩阵奇异，无法求解方程组！'); end% 交换两行if r>kfor j = k:nz=A (k,j); A (k,j)=A (r,j);A (r,j)=z;endz=b (k);b (k)=b (r);b (r)=z;end% 消元计算b (k) = b (k)/A(k,k);for j = k+1:nA (k,j) = A (k,j)/A(k,k);endfor i=1:nfor j=k+1:nA (i,j) = A (i,j)-A (i,k)*A (k,j); endb (i)=b (i)-A (i,k)*b (k); endendendx = b; % 输出bend。

gauss_jordan法求矩阵的逆先来点预备知识。

矩阵的3种运算我们称之为“⾏初等变换”：1. 交换任意2⾏2. 某⼀⾏的元素全部乘以⼀个⾮0数3. 某⼀⾏的元素加上另⼀⾏对应元素的N倍，N不为0以矩阵实施⾏初等变换等同于在矩阵左边乘以⼀个矩阵。

当要求矩阵A的逆时，在A的右边放⼀个单位矩阵，我们称[A|I]为增⼴矩阵。

对增⼴矩阵实施⾏初等变换，即左乘⼀个矩阵P，如果使得P[A|I]= [PA|P]=[I|P]，则P就是A−1。

通过⼀系列的⾏初等变换把[A|I]变成[I|P]的形式，有很多种途径，⽽数值计算就是要找到⼀种确定性的便于计算机执⾏的⽅法，gauss jordan消元法就是这样⼀种⽅法，第i次迭代时，它让增⼴矩阵的第i⾏乘以⼀个系数，使得增⼴矩阵的第i⾏第i列上的元素变为1，然后让第i⾏以外的其他⾏加上第i⾏上对应元素的N倍，使得其他⾏的第i列上的元素变为0。

下⾯举例说明gauss_jordan消元法的计算过程。

上例中A是⼀个3阶矩阵，所以经过3次迭代得到了它的逆矩阵，每次迭代增⼴矩阵中的每个元素（共3⾏6列18个元素）都要变⼀次，所以算法的时间复杂度为O(2∗n3)观察上⾯增⼴矩阵的变化过程，我们发现在每⼀步迭代的结果中，增⼴矩阵左侧有x列已化为单位矩阵时，右侧就有n-x列保持着单矩阵的样⼦，即总能从增⼴矩阵中抽出n列组成⼀个单矩阵。

同时左侧已化为单位矩阵的那⼏列在以后的⾏初等变换为始终保持不变。

所以，可以把右侧不再是单位矩阵的列存储到左侧已变为单矩阵的列上，这样就不需要额外的内存来存储整个增⼴矩阵了，内存开销减少了⼀半，同时算法的时间复杂度也降为O(n3)（虽然量级上没有变化）。

数值计算的迭代过程往往都伴随着舍⼊误差的累积，所以最终的结果也会有误差，如果这个最终的误差在⼀个可控的范围内，则称该算法为数值稳定的算法，否则为数值不稳定的算法。

什么时候会造成数值不稳定？⽐如算法某⼀步要除以⼀个很⼩的数，⼩到绝对值趋近于0，商趋于⽆穷⼤，此时舍⼊误差⼤到不可控。

2012-2013（1）专业课程实践论文列主元Gauss-Jordan消去法求矩阵的逆夏文春，0818180226，R数学08-2班董蒙蒙，0818180104，R数学08-1班Gauss 消元法始终是消去对角线下方的元素，而Gauss-Jordan 列主元消去法是消去对角线上方和下方的元素，是对Gauss 消元法的一种修正。

Gauss-Jordan 列主元消去法在算法复杂度上虽然没有Gauss 消元法小，但是它在矩阵求逆理论有重要的应用。

Gauss-Jordan 列主元消去法求方阵的逆，步骤如下：1 detA ←1；对于k=1,2,…,n 做到步82 按列选主元素|k i k a |=max ni k ≤≤|ik a |；0c ←k i k a ，)(k Ip ←k i .3 如果0c =0，则计算停止（此时A 为奇异矩阵）.4 如果k i =k ，则转不5，否则换行：kj a ↔j i k a （j=1,2,…,n ）,detA ←—detA.5 detA ←0c *detA ；6 计算h ←kk a ←1/0c ;ik a ←ik m =-ik a ⋅h(i=1,2,…,n;i ≠k ).7 消元计算ij ij a a ←+kj ik a m1,2,3......,;1,2,3......,;i n i k j n j k =≠⎧⎫⎨⎬=≠⎩⎭. 8 计算主行);,,2,1(k j n j h a a ij kj ≠⋅⋅⋅=⋅←.9 交换列对于k=n-1,n-2,.......,2,1;1) t=Ip(k);2) 如果t ≤k ，则转至3），否则换列：ik a ↔it a (i=1,2,…,n ).3) 继续循环否是否是开始 输入原始方阵所在文件夹名 将初始矩阵的信息数据导入程序 按列选主元素|k i k a |=max n i k ≤≤|ik a |；k i k a c ←0 0c >=0.0000001 如果k i =k ，则转不5，否则换行：kj a ↔j i k a （j=1,2,…,n ）；detA ←—detA 按列逐次消元，并且记录每次选取列元所在 i>n 输出#include <iostream.h>#include <fstream.h>#include <stdlib.h>#include <math.h>#define a(u,v) a[(u)*n+(v)-w]void main(){//-----------------------文件流---------------------------- char filename1[256],filename2[256];cout<<"Enter the path to input file:";cin>>filename1;cout<<"Enter the path to output file:";cin>>filename2;ifstream infile(filename1,ios::in|ios::nocreate); ofstream outfile(filename2);if(!infile){cout<<"Can not open input file:"<<filename1<<' '; exit(1);}//----------------------变量定义---------------------------- int n,w,k,i,j;int *p;double *a;double t;//------------------------输入------------------------------ infile>>n;w=n+1;p=new int[w];a=new double[n*n];for(i=0;i<n*n;i++)infile>>a[i];//------------------------算法------------------------------for(k=1;k<w;k++) //step 1{ p[k]=k;for(i=k;i<w;i++) //step 2{if(fabs(a(i,k))>fabs(a(p[k],k))) p[k]=i;}if(a(p[k],k)<1e-10) //step 3{cout<<"A is singular!"<<endl;exit(1);}if(p[k]!=k) //step 4for(i=1;i<w;i++){t=a(k,i);a(k,i)=a(p[k],i);a(p[k],i)=t;}a(k,k)=1/a(k,k); //step 5for(i=1;i<w;i++)if(i!=k)a(i,k)*=(-a(k,k)); //step 6for(i=1;i<w;i++) //step 7if(i!=k)for(j=1;j<w;j++)if(j!=k)a(i,j)+=a(i,k)*a(k,j);for(j=1;j<w;j++) //step 8if(j!=k)a(k,j)*=a(k,k);}for(k=n;k>0;k--) //step 9if(k!=p[k])for(i=1;i<w;i++){t=a(i,k);a(i,k)=a(i,p[k]);a(i,p[k])=t;}outfile<<"各步主行行号依次为::"<<endl; for(i=1;i<w;i++)outfile<<p[i]<<endl;' ';outfile<<"A的逆矩阵为:"<<endl; //step 10 for(i=1;i<w;i++){for(j=1;j<w;j++)outfile<<a(i,j)<<' ';outfile<<endl;}delete a,p;infile.close();outfile.close();}四、算法实现一）利用高斯约当列主元消去法求方阵A=⎪⎪⎪⎭⎫ ⎝⎛121061084642 的逆矩阵解：运行程序(1) 显示出 Enter the path to input file:输入1.txt ，回车。

Gauss-Jordan消去法考虑初二时所学的消元法。

所谓消元法，就是把若干个方程分别乘上常数然后加起来，达到消去某些元的目的。

我们把它看做这些方程的一个线性叠加，在这个过程中，最重要的是怎么选择叠加的系数。

比如在Gauss消去法中，通常是拿两个方程(子)和(丑)叠加，消去一个元(甲)，叠加系数由甲的系数决定，设(子)中甲的系数为甲子，(丑)中甲的系数为甲丑，那么可以用(子)乘上甲丑减去(丑)乘上甲子来消去甲。

两个方程叠加一般可以消去一个元。

类似的，对于一个N*N线性方程组，N个方程叠加，理应存在N个叠加系数，使得叠加后消去N-1个元，那么这些系数是怎么样的呢？如果能得到这些系数，就可以直接叠加N个方程，得到只含一个元的一次方程（因已消去N-1个元），这样，对每个未知元，都可以通过这个叠加变成仅含该元的一次方程，从而获得解。

下面讨论这些叠加系数应满足的条件。

设N*N的方程如下：a(1,1)*x(1)+a(1,2)*x(2)+…+a(1,n)*x(n)=b(1)(1)a(2,1)*x(1)+a(2,2)*x(2)+…+a(2,n)*x(n)=b(2)(2) …………………………………………………………a(n,1)*x(1)+a(n,2)*x(2)+…+a(n,n)*x(n)=b(n)(n)我们假设消去N-1个元后剩下的未知数为x(i)，叠加系数分别为c(i,1)，c(i,2)，…，c(i,n)。

那么叠加后的方程为：c(i,1)*(1) +c(i,2)*(2) +…+c(i,n)*(n)。

即：[c(i,1)*a(1,1)+c(i,2)*a(2,1)+…+c(i,n)*a(n,1)]x(1)+[c(i,1)*a(1,2)+c(i,2)*a(2,2)+…+c(i,n)*a(n,2)]x(2)+……………………………………………………+[c(i,1)*a(1,i)+c(i,2)*a(2,i)+…+c(i,n)*a(n,i)]x(i)+…………………………………………………+[c(i,1)*a(1,n)+c(i,2)*a(2,n)+…+c(i,n)*a(n,n)]x(n)=c(i,1)*b(1) +c(i,2)*b(2) +…+c(i,n)*b(n)按要求，除x(i)的系数不为零之外，其他未知数的系数应该为零，另外，考虑到叠加系数全部乘上一个非零常数之后的效果一样，于是可假设x(i)的系数为1，于是：c(i,1)*a(1,1)+c(i,2)*a(2,1)+…+c(i,n)*a(n,1)=0c(i,1)*a(1,2)+c(i,2)*a(2,2)+…+c(i,n)*a(n,2)=0…………………………………………………c(i,1)*a(1,i)+c(i,2)*a(2,i)+…+c(i,n)*a(n,i)=1…………………………………………………c(i,1)*a(1,n)+c(i,2)*a(2,n)+…+c(i,n)*a(n,n)=0这样相当于一个N*N的线性方程组，它仅跟原方程的系数矩阵A={a(i,j)}有关，对每个i，都对应于这样一个向量{c(i,j)，j从1到n}，这些向量（共N个）可以组成一个新的矩阵C={c(i,j)}，不难知道，C的第i行乘以A的第j列刚好是delta符号：δ(i,j)(i等于j时为1，i不等于j时为0)，其实就是说：CA=E为单位矩阵。