zoukankan      html  css  js  c++  java
  • [矩阵计算]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=egin{bmatrix} alpha_1 & eta_1 & & & 0 \ eta_1 & alpha_2 &ddots & & \ & ddots & ddots & ddots &  \ & & ddots & alpha_{n-1} & eta_{n-1} \ 0 & & & eta_{n-1} & alpha_n  end{bmatrix} $$

    比较

    $$ AQ=QT $$

    两边矩阵的每一列,得到

    $$ Aq_i = eta_{i-1}q_{i-1}+a_iq_i+eta_iq_{i+1}, quad i=1,2,cdots, n$$

    由于$ eta_0, eta_n, q_0, q_{n+1} $未定义,补充定义$ eta_0q_0 = eta_nq_{n+1} = 0 $,这样上式两边乘$q_i^T$得到

    $$ alpha_i = q_i^TAq_i, qquad eta_i = q_{i+1}^TAq_i=||Aq_i-eta_{i-1}q_{i-1}-alpha_iq_i||_2$$

    可以看出只要给定任意的$q_1 in mathbb{R}^n$且$||q_1||_2=1$,就能够利用递推关系得到全部$ q_i, alpha_i, eta_i $,迭代格式为

    $ alpha_1 = q_1^TAq_1, $       //初始值

    $ r_i=Aq_i-alpha_iq_i-eta_{i-1}q_{i-1}, $

    $ eta_i=||r_i||_2, $         //由上一轮$alpha$,$q$产生新的值

    $ q_{i+1}=r_i/eta_1 (eta_i eq 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子空间。

    注意其中若某一步$eta_i=0$,则此时得到的$T_j$的特征值将都是$A$的特征值。

    二、具体算法

    Lanczos算法求大型稀疏矩阵$A$特征值(三对角矩阵)

    1. 输入$A in Smathbb{R}^{n imes 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, eta_j:=||r_j||_2$

    4. 若$eta_j=0$则结束,否则

        $q_{j+1}:=r_j/eta_j, u_{j+1}:=Aq_{j+1}-eta_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$,在$kll n$时$T_k$的特征值内就含有$lambda$的很好的近似值;位于区间内部而又与其他特征值分离得不好的特征值$lambda$,需$kgg n$,$T_k$的特征值中才会有$lambda$的较好地近似值。

    因此,当我们只需计算大型稀疏对称矩阵$A$的少数几个两端特征值时,通常只需迭代很少几步($kll n$),$T_k$的两端特征值就是$A$两端特征值的很好的近似值。

    一个直观的例子是Parllet的数据,当$n=10^4$时,取$k=300$就可以求出10个$A$的两端特征值和对应的特征向量的很好的近似值。

    若求$A$全部的特征值时,一般$k$要远大于$n$。Cullum和Willoughby的经验是对于绝大多数矩阵来讲,只需$kleqslant 3n$就可以求出其几乎全部的不同特征值达到机器精度的近似值。

    显然,对于$ extbf{HC}=E extbf{C}$问题的精确基态和几个低能态来说,Lanczos方法是有用而且可以进一步发展的。在这个问题上最高能级一般不要求求出,那么得到低能级的精确值的方法还可以进一步加速。具体的方法就是计算化学中常用的Davidson对角化。

  • 相关阅读:
    Django源码解析(1):启动程序
    python之importlib模块
    Django中间件:CsrfViewMiddleware
    Django的admin组件
    Linux学习之CentOS--CentOS6.4下Mysql数据库的安装与配置【转】
    C#读取Xml【转】
    在eclipse导入项目的步骤【转】
    Spring学习(一)——Spring中的依赖注入简介【转】
    Spring学习(二)——Spring中的AOP的初步理解[转]
    Spring之AOP
  • 原文地址:https://www.cnblogs.com/fnight/p/5720072.html
Copyright © 2011-2022 走看看