《用Python学习数值分析--解方程组》

2017-08-25

  解一元方程，我们通常面临的问题是非线性方程。对于解方程组，我们一般解决线性方程组，如

2 * x1 + 3 * x2 = 4

5 * x1 + 6 * x 2 = 7 

  很少会涉及到非线性方程组（ Systems of Nonlinear Equations），例如简单的一种

2 * x12  +  3 * x2  = 4

5 * x1 + 6 * x2      = 7

  它所对应的模型会复杂的多，较少见，暂时不谈。光是线性方程组的解题，就足够我们学习研究的了。深入的课程还有《矩阵分析》、《矩阵计算》，我觉得它们的难度就大了不少。关于这部分知识对于做算法开发人员的重要性，就不需要再多说了。矩阵在解方程中极为重要。

  同样，解方程组也有”直接法“和”迭代法“两种方式。直接法，就是我们高中时学到的消元法，也称高斯消元法，解析的方式。无论有多少个变量，理论上都是可以通过高斯消元法来求解的，只是复杂度为O(n3)，更为严重的问题是内存需求，若方阵维数10w，计算所需内存为76GB。而且，手动计算时，我们可以使用分数来存储中间结果，但是，用计算机来计算时，我们一般也只能使用float或者double来存储中间值，若再使用自定义的Scalar量来存分数，这代价又不知增加了多少，由于多次计算，前面计算的一点点误差将导致在后续过程中被无限的累积放大，算出来的结果也许对我们没什么意义了。如果通过数值迭代方式来求解，永远只能取得近似解。

一，解析方式。线性代数中LU分解可以看成就是高斯消元法。前面提到了算法复杂度为O(n3)。

  关于LU分解，不得不提到Alan Turning，就是那个图灵奖所纪念的人，wiki上也说是 Tadeusz Banachiewicz在1938年提出的。不管是谁的贡献吧。问题的关键是把matrix和Gaussian Elimination联系起来的时间才不过百年，matrix概念提出时间也不过两百年。这里，术语”高斯消元法“虽然是几年高斯1810年提出的算法，更早的在公元三世纪就被中国的数学家刘徽所使用。

  LU分解的概念是相当简单的。参考wiki即可。其中，具体算法实现时，需要一个P的置换矩阵，A = PA。防止换行时出现除0的情况。算法复杂度为2 * n3/3。其中 Ux = L-1 * b， Ux部分计算非常简单，n2 复杂度， 但是L-1 * b 计算复杂度也是O(n3)。2w×2w的矩阵，如果使用scipy.linalg.lu() 来求解，需要15G左右的内存，E5 单核需要计算70s。实在可怕。

  因为LU分解所对应的A矩阵必须是square matrix，所以便有了推广的LDU分解，LDU算法中，A矩阵可以是方形矩阵，中间的D为对角方阵。

二，迭代法

  迭代法和单变量方程定点迭代法类似。迭代法主要分为两类：

定常迭代法
雅可比法
高斯-赛德尔迭代
逐次超松弛迭代法（SOR）

Krylov子空间法（原型是共轭梯度法）
广义最小残量法（GMRES）
双共轭梯度方法（BiCG）

1.1，Jacobi method
  算法复杂度为 O(n2)，

x0 = initial value vector

xk+1 = D-1(b -  (L + U ) xk )，其中 k = 0,1,2,...

1.2，Gauss-Seidel Method

  可以说这个方法也是Jacobi方法的变体，这里，x= [u, v, w, a, b, c...]，计算v时使用到了u的值，

x0 = initial value vector

xk+1 = D-1(b -  Uxk  - Lxk+1 )，其中 k = 0,1,2,...

   对于这两种方法，如果A是严格对角占优矩阵，那么算法是收敛的。

1.3，SOR

  SOR 是 Gauss-Seidel 方法的一种变体。

x0 = initial value vector

xk+1 = (wL + D-1)[(1-w)Dxk - wUxk] + w(D + wL)-1b，其中 k = 0,1,2,...

  当w=1时，就是Gauss-Seidel方法了。当w>1时，就是over-relaxtion；当w<1时，这个方法就是Successive Under-Relaxtion。

  在线性系统分析中，对称矩阵格外受到喜爱。相较于一般矩阵，对称矩阵只有一半的无关量。当解决Ax = b时，不需要采用矩阵分解，新的解决方法是“共轭梯度法”，这对于大型的稀疏矩阵特别有效。而共轭梯度法是Krylov空间法的原型。只要满足特定条件，就能够解决其他方法不能攻克的难题。

  不同于解方程时scipy提供能不同的算法和接口，但是，解线性系统的由 numpy.linalg.solve() 提供，对于稀疏矩阵，由scipy.sparse.linalg 模块提供。

三，非线性系统求解

  一些情况下，我们还是需要求解非线性系统的。在这里对于多变量的非线性系统，我们依然可以使用强大的牛顿法，都依赖于泰勒多项式展开。一般而言，我们也不关注二次项。函数在x0 附近的值为 F(x) = F(x0) + DF(x0)*(x-x0) + O(x - x0)2。

x0 = initial value

xk+1 = xk - (DF(xk))-1 F(xk)， 其中 k = 0,1,2,...

  DF(x) 一般写作J(x)，就是雅可比矩阵。

  另外一种方法是Broyden's Method，如解方程一文中所讲过的，当我们不知道一阶导数时，就可以采用割线法来近似代替，在这里不能简单的采用割线法替代，单依然可以用近似项来替代delta项。

  对于非线性系统的求解，由scipy.optimize 模块提供。以此可知非线性系统求解等价于最优化问题。

四，总结。

  在写本文时，我参考了中文wiki的词条“迭代法”，但是，这里的“最常见的迭代法是牛顿法。其他还包括最速下降法、共轭迭代法、变尺度迭代法、最小二乘法、线性规划、非线性规划、单纯型法、惩罚函数法、斜率投影法、遗传算法、模拟退火等等”并不那么准确。最优化问题等价于方程组求解，所以，最优解也可以被认为是方程解。但是，需要说明的是这里的误差可能相当大。

ref：

https://en.wikipedia.org/wiki/Nonlinear_system
https://en.wikipedia.org/wiki/Matrix_decomposition
https://en.wikipedia.org/wiki/LU_decomposition
https://en.wikipedia.org/wiki/Gaussian_elimination
https://en.wikipedia.org/wiki/QR_decomposition
https://en.wikipedia.org/wiki/Computational_complexity_of_mathematical_operations
https://en.wikipedia.org/wiki/Iterative_method
https://en.wikipedia.org/wiki/Sparse_matrix

如果有任何意见，欢迎留言讨论。 

[ 主页 ]