《用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为对角方阵。

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

  1. 定常迭代法
    1. 雅可比法
    2. 高斯-赛德尔迭代
    3. 逐次超松弛迭代法(SOR)
  2. Krylov子空间法(原型是共轭梯度法
    1. 广义最小残量法GMRES
    2. 双共轭梯度方法(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 模块提供。

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

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

  DF(x) 一般写作J(x),就是雅可比矩阵。
  另外一种方法是Broyden's Method,如解方程一文中所讲过的,当我们不知道一阶导数时,就可以采用割线法来近似代替,在这里不能简单的采用割线法替代,单依然可以用近似项来替代delta项。
  对于非线性系统的求解,由scipy.optimize 模块提供。以此可知非线性系统求解等价于最优化问题。

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

ref:

  1. https://en.wikipedia.org/wiki/Nonlinear_system
  2. https://en.wikipedia.org/wiki/Matrix_decomposition
  3. https://en.wikipedia.org/wiki/LU_decomposition
  4. https://en.wikipedia.org/wiki/Gaussian_elimination
  5. https://en.wikipedia.org/wiki/QR_decomposition
  6. https://en.wikipedia.org/wiki/Computational_complexity_of_mathematical_operations
  7. https://en.wikipedia.org/wiki/Iterative_method
  8. https://en.wikipedia.org/wiki/Sparse_matrix
如果有任何意见,欢迎留言讨论。


[ 主页 ]
COMMENTS
POST A COMMENT

(optional)



(optional)