求解非奇异线性方程组的正交降阶法

李胜利，王川龙，温瑞萍

（1.太原理工大学 数学学院，山西 太原 030024；2.工程科学计算山西省重点实验室，山西 太原 030012）

1 预备知识

线性方程组

是非奇异矩阵，其解法已有了广泛的研究，主要包括直接法［1-3］和迭代法［4-5］.其中直接法包括高斯消去法，Gram-Schmidt正交化QR方法，双对角化方法等.单从浮点数（Flop Numbers）角度看，高斯消去法是解线性方程组（1）最经济的方法［1-2，6-7］.但是高斯消去法数值稳定性较差［7］，我们不得不考虑“增长因子”的影响.Gram-Schmidt正交化QR［3，8-11］与双对角化方法［6］对于病态问题更可靠，但是计算量又过大.本文提出的正交降阶方法既能保证较好的数值稳定性，又比Gram-Schmidt正交化QR以及双对角化方法所需计算量要少.

为方便，用Rn×n表示n阶实矩阵全体，Rn表示n维实向量全体，AT和xT分别表示矩阵A和向量x的转置，（x，y）表示向量x和向量y的内积.A＝［v1，v2，…，vn］是A∈Rn×n的一个列划分，即vi＝（a1，i，a2，i，…，an，i）T，i＝1，2，…，n.

下面是本文将要用到的一个重要结果.

定理1 设A∈Rn×n为非奇异矩阵，将A的第1列中第1个不为零的元素记为ai，1，将A的后n-1 列分别与第1列正交后去掉第1列和第i行剩下的（n-1）×（n-1）矩阵仍然非奇异.

证明 设A的各列向量分别记为v1，v2，…，vn-1，vn，现在将后n-1列分别和第1列正交化，用Gram-Schmidt正交化得

用矩阵表示为

将矩阵分块，令p＝（a1，1，a2，1，…，ai-1，1）T，q＝（ai+1，1，ai+2，1，…，an，1）T，r＝（ai，2，ai，3，…，ai，n）T，l＝（l1，2，l1，3，…，l1，n）T，I表示（n-1）×（n-1）单位矩阵，则式（2）变成

经第1步变换后得到第1列，与后n-1 列分别正交，于是有

这里0表示n-1 维零向量，所以可得到

所以向量-ai，1lT+rT可由-plT+B和-qlT+C线性表示.又因为将A的后n-1列分别与第1列正交后的矩阵的秩仍然为n，去掉第1列后剩余的矩阵的秩为n-1.所以剩余矩阵的行秩也为n-1.由式（3）知第i行可由其余行线性表示，则去掉第i行后剩余的n-1行线性无关，所以剩余的n-1 阶矩阵非奇异.

2 正交降阶分解算法

2.1 算法的构造

为求解线性方程组（1），记A的第1列中第1个不为零的元素为ai，1，先将矩阵A的第1列v1与后n-1列v2，…，vn分别正交，这等价于将矩阵A右乘矩阵P，于是线性方程组变成下列的等价形式：

对矩阵A＝［v1，v2，…，vn-1，vn］∈Rn×n，现在将后n-1列v2，…，vn分别和第1列v1正交化，将Gram-Schmidt正交化后A的各列记为则原方程变为

将矩阵A第1次正交化后的n阶矩阵去掉第1列和第i行，剩下的n-1阶非奇异方阵记作A2.由定理1知，线性方程组（2）的系数矩阵中第i行与其它行线性相关，去掉第i行的n-1阶非奇异方阵即为A2.将b-v1y1去掉第i个元素后的向量记作，于是n阶线性方程组（1）转化为n-1阶线性方程组，以上过程称为正交降阶的过程.对得到的新的低阶线性方程组重复以上步骤，直到降为一阶线性方程.在重复正交降阶的过程中分别设yk＝xk+lk，k+1xk+1+…+lk，nxn，k＝1，…，n.最终可求解得到y＝（y1，y2，…，yn）T.因而有

最终将一个n阶线性方程组转化为一个n阶上三角线性方程组，利用回代法可求得其解.

算法1 正交降阶分解算法

第1步：令k＝1，＝x.

第2步：记A的第1列中第1个不为零的元素为ai1，将矩阵A的第1列与其它列分别正交（也就是对A作列变换），正交后的各列记作v1，v2，…，vn-k+1，得到即v1y（k）+v2xk+1+…+vn-k+1xn＝b，其中y（k）是xk，xk+1，…，xn的线性组合，且各项系数构成R的第k行元素.

第3步：利用各列与第1 列的正交性，求得y（k），得到v2xk+1+…+vn-k+1xn＝b-v1y（k），去掉其系数矩阵的第i行，去掉右端b-v1y（k）的第i个元素，得到新的n-k阶方程.系数矩阵、未知向量和右端项仍记作和b.k＝k+1，若k＜n返回第2步，若k＝n则进入第4步.

第4步：得到等价上三角线性方程组Rx＝y，求解上三角线性方程组.

2.2 算法的复杂性分析

下面考虑正交降阶分解算法的运算量.运用正交降阶法，运算量为个flops；而利用Gram-Schmidt正交化，将n阶系数矩阵QR 分解需要个flops.

由表1可知，从运算量来看，高斯消去法是解线性方程组最经济的方法.尽管如此，有下面理由说明为什么仍然要考虑正交化方法.

对于病态问题，正交化方法更加可靠，也就是说正交化方法能更好地保持数值稳定性，不必像高斯消去法那样担心“增长因子”.

表1 各种直接法运算量比较Tab.1 Computational comparison of various direct methods

3 数值实验

例1 首先在MATLAB 上随机生成n阶非奇异的矩阵，然后将运用Gram-Schmidt正交化进行QR 分解，将所占用的CPU 的时间与正交降阶分解算法所占用的CPU 的时间进行对比，结果如表2 所示.

表2 Gram-Schmidt QR 分解法与算法1CPU 占用时间对比Tab.2 Comparison of CPU time between Gram-Schmidt QR decomposition method and algorithm 1 （s）

从CPU 的占用时间来看，当n＝500，1 000，2 000，5 000时，Gram-Schmidt QR 算法大约是算法1所占用时间的1.23，2.31，2.82，3.61倍，且阶数越大，优势越明显.

例2 设矩阵

a＝ones（1 000，1），A＊a＝b，Ax＝b的精确解为x＝ones（1 000，1）.

LU 分解方法和算法1求解求线性方程组（1）的过程对比如图1 所示.

在图1中，C表示解的分量，S表示解.

例2 表明与运算量最少的高斯消去法相比，运算量次少的算法1具有较强的数值稳定性.

例3A取病态矩阵的典型Hilbert矩阵.A＝hilb（8），x＝ones（8，1），b＝A＊x，此方程的精确解为x＝ones（8，1）.用不同方法求解该方程组，得到的结果如表3 所示.

表3 算法1与QR 方法的求解数值结果Tab.3 Comparison result for algorithm 1and QR decomposition method

例3 表明对病态Hilbert 矩阵，算法1 较Gram-Schmidt正交化QR 方法数值更稳定.

从上述数值实验的结果可以看出，本文提出的正交降阶法不仅运算速度得到了很好的改善，运算精度也得到了提高，而且具有良好的数值稳定性.

图1 算法1和LU 分解法的对比Fig.1 Comparison between LU decomposition method and algorithm 1

［1］Golub G H，Van Loan C F.Matrix Computations［M］.Baltimore：The Johns Hopkins University Press，1996.

［2］徐树方.数值线性代数［M］.北京：北京大学出版社，2011.

［3］徐树方.矩阵计算的理论与方法［M］.北京：北京大学出版社，1995.

［4］Strang G.Linear Algebra and Its Applications［M］.2nd ed.New York：Academic Press，1980.

［5］Ruhe A.Numerical aspects of gram-schmidt orthogonal-ization of vectors［J］.Linear Algebra and its Applica-tions，1983，52-53：591-601.

［6］Giraud L，Langou J.Rounding error analysis of the classical gram-schmidt orthogonalization process［J］.Numer.Math.，2005，101（1）：87-100.

［7］Rice J R.Experiments on gram-schmidt orthogonalization［J］.Math.Comput.，1966，20（94）：325-328.

［8］Daniel J W，Gragg W B，Kaufman L，et al.Reorthogonalizationand stable algorithms for updating the gramschmidt QR factorization［J］.Math.Comput.，1976，136（30）：772-795.

［9］Giraud L，Langou J.When modified gram-schmidt generates a well-conditioned set of vectors［J］.IMA Journal of Numerical Analysis，2002，22（4）：521-528.

［10］Wang C L，Yan X H.On convergence of splitting iteration methods for non-Hermitian positive definite linear systems［J］.International Journal of Computer Mathematics，2013，90（2）：292-305.

［11］Meng G Y，Wang C L，Yan X H.Self-Adaptive Non-Stationary Parallel Multisplitting Two-Stage Iterative Methods for Linear Systems［M］.Berlin Heidelberg：Springer-Verlag，2012：38-47.