[矩阵计算]Lanczos方法:求稀疏矩阵特征值
更新: 29 JUL 2016
由QR方法知,求矩阵$A$的特征值,大多需要先将其三对角化(详细方法见徐树方先生的教材。此处外链一个例子),即
$$ T=Q^TAQ $$
即找到正交矩阵$Q$使得$T$成为三对角矩阵。然而若$A$为大型稀疏矩阵,常用的方法如Householder和Givens变换都无法充分利用$A$的稀疏性,因此考虑直接计算$T$和$Q$的矩阵元以利用$A$的稀疏性加速运算。
一、Lanczos方法基本原理
将以上分解式中的$Q$写成
$$ Q=[q_1,q_2,\cdots,q_n] $$
其中$ q_i $为$Q$的列向量。$T$写成
$$ T=\begin{bmatrix} \alpha_1 & \beta_1 & & & 0 \\ \beta_1 & \alpha_2 &\ddots & & \\ & \ddots & \ddots & \ddots & \\ & & \ddots & \alpha_{n-1} & \beta_{n-1} \\ 0 & & & \beta_{n-1} & \alpha_n \end{bmatrix} $$
比较
$$ AQ=QT $$
两边矩阵的每一列,得到
$$ Aq_i = \beta_{i-1}q_{i-1}+a_iq_i+\beta_iq_{i+1}, \quad i=1,2,\cdots, n$$
由于$ \beta_0, \beta_n, q_0, q_{n+1} $未定义,补充定义$ \beta_0q_0 = \beta_nq_{n+1} = 0 $,这样上式两边乘$q_i^T$得到
$$ \alpha_i = q_i^TAq_i, \qquad \beta_i = q_{i+1}^TAq_i=||Aq_i-\beta_{i-1}q_{i-1}-\alpha_iq_i||_2$$
可以看出只要给定任意的$q_1 \in \mathbb{R}^n$且$||q_1||_2=1$,就能够利用递推关系得到全部$ q_i, \alpha_i, \beta_i $,迭代格式为
$ \alpha_1 = q_1^TAq_1, $ //初始值
$ r_i=Aq_i-\alpha_iq_i-\beta_{i-1}q_{i-1}, $
$ \beta_i=||r_i||_2, $ //由上一轮$\alpha$,$q$产生新的值
$ q_{i+1}=r_i/\beta_1 \ (\beta_i \neq 0)$
$ \alpha_{i+1} = q_{i+1}^TAq_{i+1}, \ i=1,2,\cdots, n-1$
此即Lanczos迭代。其中$q_i$称为Lanczos向量。在第$j$步产生的矩阵$T_j$称为j阶Lanczos矩阵,其特征值可能是$A$的部分特征值的很好的近似。详细内容参考Krylov子空间。
注意其中若某一步$\beta_i=0$,则此时得到的$T_j$的特征值将都是$A$的特征值。
二、具体算法
Lanczos算法求大型稀疏矩阵$A$特征值(三对角矩阵)
1. 输入$A \in S\mathbb{R}^{n\times n}, q_1 \in \mathbb{R}^n\ (||q_1||_2=1)$
2. $u_1:=Aq_1,\ j:=1$
3. $a_j:=q_j^Tu_j,\ r_j:=u_j-a_jq_j,\ \beta_j:=||r_j||_2$
4. 若$\beta_j=0$则结束,否则
$q_{j+1}:=r_j/\beta_j,\ u_{j+1}:=Aq_{j+1}-\beta_jq_j$
$j:=j+1$,转步3
这一算法的主要工作量集中在计算矩阵$A$与向量$v$的乘积$Av$上。在实际使用时,应根据$A$的具体特点,设计一个计算$Av$的子程序,使算法的运算量尽可能少。
三、Lanczos现象
如果Lanczos算法是在没有摄入误差的情况下执行的,则所得到的Lanczos向量$q_1,\cdots,q_j$是相互正交的,而且之多$n$步必然终止。但是,在误差出现的情况下,计算得到的Lanczos向量的正交性很快就会失去,有时甚至还是线性相关的。C.C.Paige指出失去正交性恰与近似特征值的精度提高有关。
在Lanczos矩阵$T_j$中,其特征值$\mu_j$只要与其他$\mu_i$不要太靠近,特征向量$||z_j||_2$就和1很接近。而迭代次数$k$越大(可以超过$n$),则$T_k$含有越多的$A$的较好地近似特征值。当$k$充分大时,$T_k$就含有$A$的所有相异的特征值,此即Lanczos现象。
粗略地讲,位于$A$谱区间两端且分离较好地特征值$\lambda$,在$k\ll n$时$T_k$的特征值内就含有$\lambda$的很好的近似值;位于区间内部而又与其他特征值分离得不好的特征值$\lambda$,需$k\gg n$,$T_k$的特征值中才会有$\lambda$的较好地近似值。
因此,当我们只需计算大型稀疏对称矩阵$A$的少数几个两端特征值时,通常只需迭代很少几步($k\ll n$),$T_k$的两端特征值就是$A$两端特征值的很好的近似值。
一个直观的例子是Parllet的数据,当$n=10^4$时,取$k=300$就可以求出10个$A$的两端特征值和对应的特征向量的很好的近似值。
若求$A$全部的特征值时,一般$k$要远大于$n$。Cullum和Willoughby的经验是对于绝大多数矩阵来讲,只需$k\leqslant 3n$就可以求出其几乎全部的不同特征值达到机器精度的近似值。
显然,对于$\textbf{HC}=E\textbf{C}$问题的精确基态和几个低能态来说,Lanczos方法是有用而且可以进一步发展的。在这个问题上最高能级一般不要求求出,那么得到低能级的精确值的方法还可以进一步加速。具体的方法就是计算化学中常用的Davidson对角化。
转载于:https://www.cnblogs.com/fnight/p/5720072.html
[矩阵计算]Lanczos方法:求稀疏矩阵特征值相关推荐
- 求矩阵特征值的方法和性质
求矩阵特征值的方法 Ax=mx,等价于求m,使得(mE-A)x=0,其中E是单位矩阵,0为零矩阵. |mE-A|=0,求得的m值即为A的特征值.|mE-A| 是一个n次多项式,它的全部根就是n阶方阵A ...
- 利用 Lanczos 方法实现张量的 HOSVD 分解
1. 特征值分解(EVD) 如果 AAA 是一个 m×mm\times mm×m 的 实对称矩阵(A=ATA=A^TA=AT) ,如果存在 mmm 维列向量 qqq 和实数 λ\lambdaλ 满足 ...
- 基于Lanczos方法的矩阵双对角化或三对角化
矩阵对角化通过简单线性运算,将复杂矩阵变为简单矩阵.例如传统的特征值和特征向量对角化可以将矩阵n阶方阵A变为A=T−1diag(a1,a2,...an)TA=T^{-1}diag(a_1,a_2,.. ...
- 求矩阵特征值和特征向量
求矩阵特征值和特征向量的一个小程序 代码较长,如果不能执行,就是要建立结构体,大家试试吧,希望能用. // // 实对称三对角阵的全部特征值与特征向量的计算 // // 参数: // 1. doubl ...
- c语言求矩阵特征值的程序,如何用C语言编写求对称矩阵的特征值和特征向量的程序编写对称矩阵的特征值和特征向量,其中矩阵用二维数组保存.特征向量要求有大到小放到数组里....
优质解答 //数值计算程序-特征值和特征向量 // //约化对称矩阵为三对角对称矩阵 //利用Householder变换将n阶实对称矩阵约化为对称三对角矩阵 //a-长度为n*n的数组,存放n阶实对称 ...
- 移位取逆迭代(shifted inverse iteration)求最近特征值和特征向量(python,数值积分)
第二十七篇 移位取逆迭代求最近特征值和特征向量 移位逆迭代 一种比"最大"特征值法更直接实现向量迭代收敛的的特征值方法是将移位向量迭代法改写为下面形式 其中p是一个标量" ...
- 【数学知识】三种方法求 [1,n] 中所有数欧拉函数(线性筛欧拉函数优化至 O(n) )
整理的算法模板合集: ACM模板 ①直接求小于或等于n,且与n互质的数个数(求[1,n]中所有数的欧拉函数时间复杂度:O(nn)O(n\sqrt{n})O(nn)) ②求[1,n]之间每个数的质因数 ...
- 用子函数的方法求一个3*4的数组的转置数组
<程序设计基础实训指导教程-c语言> ISBN 978-7-03-032846-5 p142 7.1.2 上级实训内容 [实训内容3]用子函数的方法求一个3*4的数组的转置数组 #incl ...
- 用子函数的方法求一维数组中所有元素之和
<程序设计基础实训指导教程-c语言> ISBN 978-7-03-032846-5 p142 7.1.2 上级实训内容 [实训内容2]用子函数的方法求一维数组中所有元素之和 #includ ...
最新文章
- python简单代码画皮卡丘-实现童年宝可梦,教你用Python画一只属于自己的皮卡丘...
- 【Python】监控视频中运动目标检测的代码实现及效果展示
- OS- -请求分页系统、请求分段系统和请求段页式系统(二)
- 启动linux_使用 UEFI 双启动 Windows 和 Linux | Linux 中国
- 【转】C# HttpWebRequest提交数据方式
- NameNode之文件系统目录树
- 乱码解决方案SecureCRT中文乱码解决方案
- 虚拟机中PXE-MOF:Exiting intel PXE ROM.Operating system not found解决方法
- 配置Kafka集群和zookeeper集群
- 我去,这么简单的条件表达式竟然也有这么多坑
- AutoLine源码之RobotFramework运行器
- 马里兰大计算机专业学phd博士,美国纽约州立大学石溪分校计算机专业博士CS PHD全奖OFFER...
- Java--实现简单的音频(mp3格式)播放
- 公有继承中 构造函数和析构函数的调用(包含内嵌子对象)
- 准备女儿的学前班毕业典礼
- 请求报错Required String parameter 'id' is not present
- 华为数通HCIA学习笔记之OSI参考模型TCP/IP模型
- 思考的梯子 | 黄金圈法则What-How-Why(超干货)
- 离线数仓12—— 数仓开发之DWD层
- Web 前端学习 之servlet技术(一)