三对角矩阵特征值数值求解 |
您所在的位置：网站首页 › 怎样求矩阵的对角矩阵 › 三对角矩阵特征值数值求解 |
|
背景
诸如以下的矩阵被称为三对角矩阵，也称为上下 Hessenberg 矩阵. 该矩阵除了三条对角上有数值以外，其余的矩阵元都为零. \[\begin{bmatrix} a_1 & b_1 & & & & &\\ c_1 & a_2 & b_2 & & & &\\ & c_2 & a_3 & b_3 & & &\\ & & \cdots & \cdots & \cdots & &\\ & & & c_{n-3} & a_{n-2} & b_{n-2} & \\ & & & & c_{n-2} & a_{n-1} & b_{n-1}\\ & & & & & c_{n-1} & a_n\\ \end{bmatrix} _{n\times n} \]这种矩阵不论是在工程应用，还是在数学物理等方面都十分常见，例如求解薛定谔方程的矩阵对角化等，因此掌握三对角矩阵的一些相关知识是十分必要的，本文主要是介绍三对角矩阵的一些解析解法以及相关数值计算。 三对角矩阵特征值求解 QR法在这个方法中，矩阵 \(A\) 被分解成以下形式：\(A=QR\)，其中 \(Q\) 是正交矩阵，\(R\) 是上三角矩阵。产生如下的矩阵序列： \[A_0 = Q_0R_0\\ A_1 = R_0Q_0 = Q_1R_1\\ A_2 = R_1Q_1 = Q_2R_2\\ \cdots\\ A_n = R_{n-1}Q_{n-1} = Q_{n}R_n \]因为 \(A_{i+1} = R_{i}Q_{i} = Q_{i}^{-1}Q_{i}R_{i}Q_{i} = Q_{i}^{-1} A_{i} Q_{i}\)，不难发现 \(A_{i+1}\) 与 \(A_{i}\) 相似，它们应该有相同的特征值。通过多次迭代，最终 \(A_{n}\) 趋于变为对角阵，其对角线上的元素给出原矩阵的特征值。 参考: https://blog.csdn.net/weixin_42304279/article/details/113314453 https://www.zhihu.com/question/23905796 https://zhuanlan.zhihu.com/p/102376857 特征多项式法特征多项式法可以像特征多项式 \(\phi(\lambda)\) 的根一样确定特征值。有一种有效的方法来构造三对角矩阵的特征多项式。使用符号法可以求特征值的归类，从而形成一个Sturmian序列。然后用对分法或试位法来求精确的特征值。由Householder变换得到的对称三对角矩阵的特征多项式为： \[\phi(\lambda)= \begin{bmatrix} a_1-\lambda & b_1 & & & & &\\ c_1 & a_2-\lambda & b_2 & & & &\\ & c_2 & a_3-\lambda & b_3 & & &\\ & & \cdots & \cdots & \cdots & &\\ & & & c_{n-3} & a_{n-2}-\lambda & b_{n-2} & \\ & & & & c_{n-2} & a_{n-1}-\lambda & b_{n-1}\\ & & & & & c_{n-1} & a_n-\lambda\\ \end{bmatrix} _{n\times n} \]得: \[\phi_0(\lambda) = 1\\ \phi_1(\lambda) = a_1-\lambda\\ \phi_2(\lambda) = (a_1-\lambda)(a_2-\lambda) - b_1c_1 = (a_2-\lambda)\phi_1(\lambda)-b_1c_1\phi_0(\lambda) \\ \phi_i(\lambda) = (a_i-\lambda)\phi_{i-1}(\lambda)-b_{i-1}c_{i-1}\phi_{i-2}(\lambda) \\ \]求解 \(\phi_{i}(\lambda)=0\) 的根即为矩阵的特征值。 综上，对于高次根 \(\phi(\lambda)\) 求解十分复杂。通常，对于很多程序里的大维度的矩阵特征值问题求解，QR 分解是常采用的方法，对于不同类型的矩阵 QR 分解，也存在相应的算法优化。 三对角矩阵方程组求解 Ax = d通常来说，求解下列矩阵方程组的解 \(x\) 的时间复杂度为 \(O(n^3)\)，但通过 Thomas Algorithm (又叫The tridiagonal matrix algorithm, TDMA)可以实现 \(O(n)\) 的时间复杂度。它是一种基于高斯消元法的算法，具体如下, \[\begin{bmatrix} a_1 & b_1 & & & & &\\ c_1 & a_2 & b_2 & & & &\\ & c_2 & a_3 & b_3 & & &\\ & & \cdots & \cdots & \cdots & &\\ & & & c_{n-3} & a_{n-2} & b_{n-2} & \\ & & & & c_{n-2} & a_{n-1} & b_{n-1}\\ & & & & & c_{n-1} & a_n\\ \end{bmatrix} _{n\times n} \begin{bmatrix} x_1\\ x_2\\ x_3\\ \cdots\\ x_{n-2}\\ x_{n-1}\\ x_n\\ \end{bmatrix} _{n\times n} = \begin{bmatrix} d_1\\ d_2\\ d_3\\ \cdots\\ d_{n-2}\\ d_{n-1}\\ d_n\\ \end{bmatrix} _{n\times n} \]不难发现, 核心公式为 \(c_{i-1}x_{i-1} + a_ix_i + b_ix_{i+1}=d_i\), 展开为 \[a_1x_1 + b_1x_2 = d_1 \\ c_1x_1 + a_2x_2 + b_2x_3 = d_2\\ c_2x_2 + a_3x_3 + b_3x_4 = d_3\\ c_{i-1}x_{i-1} + a_ix_i + b_ix_{i+1}=d_i\\ c_{n-2}x_{n-2} + a_{n-1}x_{n-1} + b_{n-1}x_{n}=d_{n-1}\\ c_{n-1}x_{n-1} + a_nx_n =d_n\\ \]迭代联立求解该式子,便可得到该方程组的解，详细参考： https://www.cnblogs.com/xpvincent/archive/2013/01/25/2877411.html 经典三对角矩阵本征值求解本部分以下面特殊的三对角矩阵为例 \[D_n = \begin{vmatrix} a & b & & & & &\\ c & a & b & & & &\\ & c & a & b & & &\\ & & \cdots & \cdots & \cdots & &\\ & & & c & a & b & \\ & & & & c & a & b\\ & & & & & c & a\\ \end{vmatrix} _{n\times n} \]把该矩阵行列式写成如下形式，其中 \(x+y=a, xy=bc\). \[D_n = \begin{vmatrix} x+y & xy & & & & &\\ 1 & x+y & xy & & & &\\ & 1 & x+y & xy & & &\\ & & \cdots & \cdots & \cdots & &\\ & & & 1 & x+y & xy & \\ & & & & 1 & x+y & xy\\ & & & & & 1 & x+y\\ \end{vmatrix} _{n\times n} \]我们先来计算 \(D_n\)，显然 \(D_1 = x+y, D_2 = (x+y)^2 -xy\). 由前面可知，\(D_n = (x+y)D_{n-1}-xyD_{n-2}\), \(D_n - xD_{n-1}= y(D_{n-1}-xD_{n-2})\)，反复利用这个递推关系可得 \(D_n - xD_{n-1}= y^{n-1}(D_{1}-xD_{0}) = y^n\). 因此 \(D_{n} = xD_{n-1} + y^n\), 再利用这个递推关系得: \[D_n = y^n + xD_{n-1} = y^n + x(y^{n-1} +xD_{n-2})\\ =y^n + xy^{n-1} + x^2y^{n-2} + \cdots + x^{n-1}y + x^n= \left\{ \begin{aligned} &\frac{x^{n+1}-y^{n+1}}{x-y}, & x\neq y \\ &(n+1)x^n, & x =y \end{aligned} \right. \]结合 \(xy = bc,x+y=a\)，可知: \(x=\frac{a+\sqrt{a^2-4bc}}{2},y=\frac{a-\sqrt{a^2-4bc}}{2}\)，因此, \[D_n = \left\{ \begin{aligned} &\frac{({\frac{a+\sqrt{a^2-4bc}}{2}})^{n+1}-({\frac{a-\sqrt{a^2-4bc}}{2}})^{n+1}}{\sqrt{a^2-4bc}}, & a^2-4bc\neq 0 \\ &(n+1)(a/2)^n, & a^2-4bc = 0 \end{aligned} \right. \]因为该矩阵的特征多项式为: \[\phi(\lambda) = \begin{vmatrix} a-\lambda & b & & & & &\\ c & a-\lambda & b & & & &\\ & c & a-\lambda & b & & &\\ & & \cdots & \cdots & \cdots & &\\ & & & c & a-\lambda & b & \\ & & & & c & a-\lambda & b\\ & & & & & c & a-\lambda\\ \end{vmatrix} _{n\times n} \]那么 \[\phi(\lambda)=\left\{ \begin{aligned} &\frac{({\frac{a-\lambda+\sqrt{(a-\lambda)^2-4bc}}{2}})^{n+1}-({\frac{(a-\lambda)-\sqrt{(a-\lambda)^2-4bc}}{2}})^{n+1}}{\sqrt{(a-\lambda)^2-4bc}}, & (a-\lambda)^2-4bc\neq 0 \\ &(n+1)((a-\lambda)/2)^n, & (a-\lambda)^2-4bc = 0 \end{aligned} \right. \]令 \(\phi(\lambda)=0\)，则 \(({\frac{a-\lambda+\sqrt{(a-\lambda)^2-4bc}}{2}})^{n+1}-({\frac{(a-\lambda)-\sqrt{(a-\lambda)^2-4bc}}{2}})^{n+1} = 0\)， 即 \(\frac{a-\lambda+\sqrt{(a-\lambda)^2-4bc}}{a-\lambda-\sqrt{(a-\lambda)^2-4bc}} = e^{\frac{2k\pi}{n+1}}\), 因此, 该解为 \(\lambda = a + 2\sqrt{bc}\cos\frac{k\pi}{n+1}\), 其中 \(k=1,2,\cdots,n\). 各类三对角矩阵的严格解 矩阵1 \[\begin{bmatrix} a & b & & & & &\\ b & a & b & & & &\\ & b & a & b & & &\\ & & \cdots & \cdots & \cdots & &\\ & & & b & a & b & \\ & & & & b & a & b\\ & & & & & b & a\\ \end{bmatrix} _{n\times n} \]该矩阵的本征值为: \[\lambda_k = a + 2b\cos{\frac{k\pi}{n+1}}, k = 1,2,3,\cdots, n. \]更一般地, 对于如下三对角矩阵， \[\begin{bmatrix} a & b & & & & &\\ c & a & b & & & &\\ & c & a & b & & &\\ & & \cdots & \cdots & \cdots & &\\ & & & c & a & b & \\ & & & & c & a & b\\ & & & & & c & a\\ \end{bmatrix} _{n\times n} \]该矩阵的本征值为: \[\lambda_k = a + 2\sqrt{bc}\cos{\frac{k\pi}{n+1}}, k = 1,2,3,\cdots, n. \]矩阵2 \[\begin{bmatrix} a-b & b & & & & &\\ b & a & b & & & &\\ & b & a & b & & &\\ & & \cdots & \cdots & \cdots & &\\ & & & b & a & b & \\ & & & & b & a & b\\ & & & & & b & a\\ \end{bmatrix} _{n\times n} \]该矩阵的本征值为: \[\lambda_k = a + 2b\cos{\frac{2k\pi}{2n+1}}, k = 1,2,3,\cdots, n. \]矩阵3 \[\begin{bmatrix} a-b & b & & & & &\\ b & a & b & & & &\\ & b & a & b & & &\\ & & \cdots & \cdots & \cdots & &\\ & & & b & a & b & \\ & & & & b & a & b\\ & & & & & b & a+b\\ \end{bmatrix} _{n\times n} \]该矩阵的本征值为: \[\lambda_k = a + 2b\cos{\frac{(2k-1)\pi}{2n}}, k = 1,2,3,\cdots, n. \]矩阵4 \[\begin{bmatrix} a+b & b & & & & &\\ b & a & b & & & &\\ & b & a & b & & &\\ & & \cdots & \cdots & \cdots & &\\ & & & b & a & b & \\ & & & & b & a & b\\ & & & & & b & a+b\\ \end{bmatrix} _{n\times n} \]该矩阵的本征值为: \[\lambda_k = a + 2b\cos{\frac{k\pi}{n}}, k = 1,2,3,\cdots, n. \]矩阵5 \[\begin{bmatrix} a & 1 & & & & &\\ 1 & b & 1 & & & &\\ & 1 & a & 1 & & &\\ & & \cdots & \cdots & \cdots & &\\ & & & 1 & b & 1 & \\ & & & & 1 & a & 1\\ & & & & & 1 & b\\ \end{bmatrix} _{2n\times 2n} \]该矩阵的本征值为: \[\lambda_k = \frac{a+b\pm[(a-b)^2+16\cos^2\frac{k\pi}{2n+1}]^{1/2}}{2}, k = 1,2,3,\cdots, n. \]矩阵6 \[\begin{bmatrix} a & 1 & & & & &\\ 1 & b & 1 & & & &\\ & 1 & a & 1 & & &\\ & & \cdots & \cdots & \cdots & &\\ & & & 1 & a & 1 & \\ & & & & 1 & b & 1\\ & & & & & 1 & a\\ \end{bmatrix} _{2n+1\times 2n+1} \]该矩阵的本征值为: \[\lambda_k = \frac{a+b\pm[(a-b)^2+16\cos^2\frac{k\pi}{2n+2}]^{1/2}}{2}, k = 1,2,3,\cdots, n.\\ \lambda = a. \]数值求解三对角矩阵对于求解高维度的三对角矩阵, 例如 10000*10000 维度，如果不是一些特殊的三对角矩阵，那么就不存在上面所述的严格解析解，此时则需要严格对角化求解矩阵的本征值。因为三对角矩阵是稀疏矩阵，为节省内存资源，通常不用常规的矩阵方法去存储矩阵元。另外，因为在很多实际问题中，我们往往只需要求解特征值的前几个值而已，并不要求全部得到所有的特征值，此时如果采用严格对角化往往会浪费很多时间，是否存在比较快速的方法求解三对角矩阵的特征值？ 矩阵幂法可以快速求解最大的特征值和相应的特征向量。 https://baike.baidu.com/item/幂法求矩阵特征值/3082486?fr=aladdin Lanczos 算法Lanczos 方法是求解大规模对称矩阵端部特征问题的一种正交投影方法. 此处先mark，有时间再看！ 参考: https://www.cnblogs.com/qxred/p/lanczos.html Python 三对角矩阵 scipy.sparse.linalg.eigs(A, k=6, M=None, sigma=None, which='LM', v0=None, ncv=None, maxiter=None, tol=0, return_eigenvectors=True, Minv=None, OPinv=None, OPpart=None)Python能否进行大规模数值计算？ - 幽州狼的回答 - 知乎 https://www.zhihu.com/question/68171807/answer/2186194693 Matlab 三对角矩阵利用 eigs 函数计算大型稀疏矩阵的特征值。 d = eigs(A) %返回一个向量，其中包含矩阵 A 的六个模最大的特征值。当使用 eig 计算所有特征值的计算量很大时（例如对于大型稀疏矩阵来说），这是非常有用的。 参考文献[1] Master thesis 1953: The Characteristic Roots of Certain Real Symmetric Matrices. https://trace.tennessee.edu/cgi/viewcontent.cgi?article=3834&context=utk_gradthes [2] Devadatta et al., Eigenvalues of tridiagonal pseudo-Toeplitz matrices. Linear Algebra and its Applications, 297, 63-80, 1999. https://doi.org/10.1016/S0024-3795(99)00114-7 [3] Moawwad EI-Mikkawy. A note on a three-term recurrence for a tridiagonal matrix. Applied Mathematics and Computation, 139, 503-511, 2003. https://doi.org/10.1016/S0096-3003(02)00212-6 [4] 杨胜良. 三对角行列式及其应用[J]. 工科数学, 2002, (02):102-104. [5] PhD thesis 1997: A New \(O(n^2)\) Algorithm for the Symmetric Tridiagonal Eigenvalue Eigenvector Problem.https://www.cs.utexas.edu/users/inderjit/public_papers/thesis.pdf 附录文件 |
今日新闻 |
推荐新闻 |
| CopyRight 2018-2019 办公设备维修网 版权所有 豫ICP备15022753号-3 |