2.9 使用其他矩阵代数方法求解
前述已经介绍了逆矩阵的求解方法,以及利用LU分解、Cholesky分解和QR分解求解线性方程组。量化投资分析师应充分了解上述内容,以解决矩阵中包含大量金融数据所带来的问题。
某些情况下,我们求的解不收敛,可以使用迭代法解决此类问题,如Jacobi迭代、Gauss-Seidel迭代和SOR迭代法。
2.9.1 Jacobi迭代法
Jacobi迭代法通过对矩阵的对角元迭代求解线性方程组,计算结果收敛时终止迭代。方程Ax = B中,矩阵A = D + R,矩阵D为对角矩阵。以一个4×4的矩阵A为例:
通过迭代得出答案:
Ax = B
(D + R)x = B
Dx = B – Rx
xn+1 = D–1(B – Rxn)
与Gauss-Seidel迭代法不同,Jacobi方法欲求xn+1必须先求出xn,这将占用两倍的内存。然而,矩阵中每个元素的计算都以并行方式完成会显著提高运算速度。
如果矩阵A是一个不可约严格对角占优(strictly irreducibly diagonally dominant)矩阵,通过Jacobi迭代法得到的解一定是收敛的。不可约严格对角占优矩阵是每个对角元的绝对值都大于所在行非对角元绝对值之和的矩阵。
某些情况下,即使矩阵不满足上述条件也能通过Jacobi迭代得到收敛解:
本题矩阵A使用与Cholesky分解中相同的值,利用jacobi
函数进行25次迭代:
调用jacobi
函数求解x
:
最终x的值与通过Cholesky分解所求相同。
2.9.2 Gauss-Seidel迭代法
Gauss-Seidel迭代法与Jacobi迭代法很相似。方程Ax = B中,矩阵A = L + U,矩阵L为下三角矩阵,矩阵U为上三角矩阵。以一个4×4的矩阵A为例:
迭代得:
Ax = B
(L + U)x = B
Lx = B – Ux
xn+1 = L–1(B – Uxn)
利用下三角矩阵L计算xn+1,不必先求出xn,这相比Jacobi迭代法将节省一半的存储空间。
利用Gauss-Seidel迭代法求解的收敛速度很大程度取决于矩阵性质,需要严格对角占优或正定矩阵。即使这些条件没有满足,Gauss-Seidel迭代的结果仍可能收敛。
用Python实现Gauss-Seidel迭代法:
利用NumPy模块的tril()
函数通过下三角矩阵U返回下三角矩阵A,利用tol
定义的公差迭代求得x。
使用与Jacobi迭代和Cholesky分解中相同的矩阵,我们将gauss()
函数n的最大值设为100来计算x:
验证求得x值是否与前述相等:
结果显示求得的x与通过Jacobi迭代和Cholesky分解的结果相同。