Solving linear equation systems

现代数值计算方法一般都会转化为求解线性方程组. 线性方程组的矩阵形式如下所示:

[A][x]=[b]

where A is m by n matrix,m≥n, x∈Rn, and b∈Rm.

由于求解线性方程组在数值分析中基础性核心地位, 线性方程组求解算法和技术一直层出不穷: 直接法,迭代法,并行算法,分块技术,超线程技术……

如果矩阵A是稀疏矩阵,或者是对角阵, 利用矩阵A的结构特点, 波前法,共轭法等可以高效求解大规模线性方程组.

对于普通的稠密矩阵A, The most directly solution of the equation system is

[x]=[A]−1[b]

Calculating the inverse of matrix [A] is easy to be implemented, however it requires approximately n! operations. 对于大规模线性方程组, 未知变量较多时,n! 是一个天文数字. 但当A是低阶矩阵时,这种直接求解算法简易,可靠. 计算量较小.例如下列4X4矩阵求逆算法[3]:

bool gluInvertMatrix(const double m[16], double invOut[16]){    double inv[16], det;    int i;    inv[0] = m[5]  * m[10] * m[15] -              m[5]  * m[11] * m[14] -              m[9]  * m[6]  * m[15] +              m[9]  * m[7]  * m[14] +             m[13] * m[6]  * m[11] -              m[13] * m[7]  * m[10];    inv[4] = -m[4]  * m[10] * m[15] +               m[4]  * m[11] * m[14] +               m[8]  * m[6]  * m[15] -               m[8]  * m[7]  * m[14] -               m[12] * m[6]  * m[11] +               m[12] * m[7]  * m[10];    inv[8] = m[4]  * m[9] * m[15] -              m[4]  * m[11] * m[13] -              m[8]  * m[5] * m[15] +              m[8]  * m[7] * m[13] +              m[12] * m[5] * m[11] -              m[12] * m[7] * m[9];    inv[12] = -m[4]  * m[9] * m[14] +                m[4]  * m[10] * m[13] +               m[8]  * m[5] * m[14] -                m[8]  * m[6] * m[13] -                m[12] * m[5] * m[10] +                m[12] * m[6] * m[9];    inv[1] = -m[1]  * m[10] * m[15] +               m[1]  * m[11] * m[14] +               m[9]  * m[2] * m[15] -               m[9]  * m[3] * m[14] -               m[13] * m[2] * m[11] +               m[13] * m[3] * m[10];    inv[5] = m[0]  * m[10] * m[15] -              m[0]  * m[11] * m[14] -              m[8]  * m[2] * m[15] +              m[8]  * m[3] * m[14] +              m[12] * m[2] * m[11] -              m[12] * m[3] * m[10];    inv[9] = -m[0]  * m[9] * m[15] +               m[0]  * m[11] * m[13] +               m[8]  * m[1] * m[15] -               m[8]  * m[3] * m[13] -               m[12] * m[1] * m[11] +               m[12] * m[3] * m[9];    inv[13] = m[0]  * m[9] * m[14] -               m[0]  * m[10] * m[13] -               m[8]  * m[1] * m[14] +               m[8]  * m[2] * m[13] +               m[12] * m[1] * m[10] -               m[12] * m[2] * m[9];    inv[2] = m[1]  * m[6] * m[15] -              m[1]  * m[7] * m[14] -              m[5]  * m[2] * m[15] +              m[5]  * m[3] * m[14] +              m[13] * m[2] * m[7] -              m[13] * m[3] * m[6];    inv[6] = -m[0]  * m[6] * m[15] +               m[0]  * m[7] * m[14] +               m[4]  * m[2] * m[15] -               m[4]  * m[3] * m[14] -               m[12] * m[2] * m[7] +               m[12] * m[3] * m[6];    inv[10] = m[0]  * m[5] * m[15] -               m[0]  * m[7] * m[13] -               m[4]  * m[1] * m[15] +               m[4]  * m[3] * m[13] +               m[12] * m[1] * m[7] -               m[12] * m[3] * m[5];    inv[14] = -m[0]  * m[5] * m[14] +                m[0]  * m[6] * m[13] +                m[4]  * m[1] * m[14] -                m[4]  * m[2] * m[13] -                m[12] * m[1] * m[6] +                m[12] * m[2] * m[5];    inv[3] = -m[1] * m[6] * m[11] +               m[1] * m[7] * m[10] +               m[5] * m[2] * m[11] -               m[5] * m[3] * m[10] -               m[9] * m[2] * m[7] +               m[9] * m[3] * m[6];    inv[7] = m[0] * m[6] * m[11] -              m[0] * m[7] * m[10] -              m[4] * m[2] * m[11] +              m[4] * m[3] * m[10] +              m[8] * m[2] * m[7] -              m[8] * m[3] * m[6];    inv[11] = -m[0] * m[5] * m[11] +                m[0] * m[7] * m[9] +                m[4] * m[1] * m[11] -                m[4] * m[3] * m[9] -                m[8] * m[1] * m[7] +                m[8] * m[3] * m[5];    inv[15] = m[0] * m[5] * m[10] -               m[0] * m[6] * m[9] -               m[4] * m[1] * m[10] +               m[4] * m[2] * m[9] +               m[8] * m[1] * m[6] -               m[8] * m[2] * m[5];    det = m[0] * inv[0] + m[1] * inv[4] + m[2] * inv[8] + m[3] * inv[12];    if (det == 0)        return false;    det = 1.0 / det;    for (i = 0; i < 16; i++)        invOut[i] = inv[i] * det;    return true;}

为了提高求解效率, 利用变换等多种方法简化求解线性方程组, 是降低计算量的重要途径. 常用的算法主要有: 高斯消元法(Gaussian elimination), Gauss-Jordan 法, LU分解法(LU decomposition), QR分解法(QR decomposition)等.

[8]In Gaussian elimination is unstable when meeting zero pivots. As small pivots operate similar to zeros during Gaussian elimination, some LU decompositions become numerically unstable if relatively small pivots are used. In order to prevent this instability, realtively large entries are selected as pivot entries. This prevents large factors from appearing in the computed L and U, which reduces roundoff errors during compuation. [9] provided some case studies of pivoting in LU decomposition. There are three kinds of pivoting: partial pivoting, complete pivoting and rook pivoting.

Using QR factorization is much better to solve Ax=b than using LU factorization even with partial pivoting[10].

[1] 列出了几种主要算法的简要对比：

Method Computing Complexity Notes Laplace expansion n! numerical errors Gaussian elimination n3/3+n2/2 numerical errors LU decomposition n3/3+n2−n/3 numerical errors QR decomposition 2n3+3n3 Good Precision

[1] http://qucs.sourceforge.net/tech/node99.html
[2] https://en.wikipedia.org/wiki/System_of_linear_equations
[3] https://stackoverflow.com/questions/1148309/inverting-a-4x4-matrix
[4] http://fourier.eng.hmc.edu/e176/lectures/NM/node10.html
[5] http://www.math.iit.edu/~fass/477577_Chapter_4.pdf
[6] http://www.aaronschlegel.com/qr-decomposition-householder-reflections/
[7] http://www.math.usm.edu/lambers/mat610/sum10/lecture9.pdf
[8] http://buzzard.ups.edu/courses/2014spring/420projects/math420-UPS-spring-2014-reid-LU-pivoting.pdf
[9] http://www.cnblogs.com/lacozhang/p/3693203.html
[10] https://www-old.math.gatech.edu/academic/courses/core/math2601/Web-notes/3num.pdf