下一节:基本算法
上一级:Jacobi-Davidson方法
上一节:Jacobi-Davidson方法
基本理论
Jacobi-Davidson方法基于两个基本原则的结合。第一个原则是针对特征问题A{x}=\lambda {x}应用Galerkin方法,相对于由正交基{v}_1,\ldots,{v}_m张成的给定子空间。使用不同于Krylov子空间的方法是由Davidson[99]提出的,他还提出了特定的选择,与我们即将采取的方法不同,用于构造正交基向量{v}_j。Galerkin条件是
AV_m{s}-\theta V_m{s} \perp \{{v}_1, \ldots, {v}_m\} ,
这导致简化的系统
V_m^\ast AV_m {s} - \theta {s} = 0, \tag{4.48}
其中V_m表示列从{v}_1到{v}_m的矩阵。这个方程有m个解(\theta_j^{(m)}, {s}_j^{(m)})。m对(\theta_j^{(m)}, {u}_j^{(m)}\equiv V_ms_j^{(m)})被称为{A}相对于由{V}_m列张成的子空间的Ritz值和Ritz向量。这些Ritz对是{A}的特征对的近似,我们的目标是通过精心选择的子空间扩展来获得更好的近似。
在此基础上,Jacobi-Davidson方法的另一个原则开始发挥作用。这一思想可以追溯到Jacobi[241]。假设我们有一个对应于给定特征值\lambda的特征向量{x}的近似{u}_j^{(m)}。那么Jacobi建议(在原始论文中针对强对角占优矩阵)计算{u}_j^{(m)}的正交校正{t},使得
A({u}_j^{(m)} + {t}) = \lambda ({u}_j^{(m)} +{t}) .
由于{t}\perp {u}_j^{(m)},我们可以将自己限制在与{u}_j^{(m)}正交的子空间中。在该子空间中,由{u}_j^{(m)}张成的子空间上的算子A由
\left(I-{u}_j^{(m)}{{u}_j^{(m)}}^\ast\right)A\left(I-{u}_j^{(m)}{{u}_j^{(m)}}^\ast\right) ,
给出,我们发现{t}满足方程
\left(I-{u}_j^{(m)}{{u}_j^{(m)}}^\ast\right)(A-\lambda I)\left(I-{u}_j^{(m)}{{u}_j^{(m)}}^\ast\right){t} =- (A-\theta_j^{(m)} I){u}_j^{(m)} .
在实际情况下,我们不知道\lambda,显然的解决方案是用其近似\theta_j^{(m)}替换它,这导致了针对更新t_j^{({m})}\perp {u}_j^{(m)}的Jacobi-Davidson校正方程:
\left(I-{u}_j^{(m)}{{u}_j^{(m)}}^\ast\right)\left(A-\theta_j^{(m)}I\right)\left(I-{u}_j^{(m)}{{u}_j^{(m)}}^\ast\right){t}_j^{({m})} = -r^{(m)}_j \equiv -(A-\theta_j^{(m)}I){u}_j^{(m)}. \tag{4.49}
这个校正方程只是近似求解,其近似解\tilde{t}^{(m)}用于子空间的扩展。这是与Krylov子空间方法的一个根本区别;我们不是选择一个作为算子作用于给定起始向量的幂的子空间,而是选择一些没有Krylov结构的子空间,并将给定矩阵投影到这个子空间上。
从(4.48)我们得出r_j^{({m})}\perp \{{v}_1, \ldots,{v}_m\},特别是r_j^{({m})}\perp {u}_j^{(m)},因此Jacobi-Davidson校正方程表示一个一致的线性系统。
可以证明,(4.49)的精确解导致最大的\theta_j^{(m)}对于增加的m向\lambda_{\max}(A)立方收敛(对于A的其他特征值的收敛也可以做出类似陈述,前提是每一步中适当选择Ritz值)。
在[411]中建议仅近似求解方程(4.49),例如,通过一些最小残差(MINRES)[350]步骤,如果有可用的适当预处理器{K}用于A-\theta_j^{(m)}I,但实际上任何近似技术对于t_j^{(m)}都是允许的,前提是考虑到投影(I-{u}_j^{(m)}{{u}_j^{(m)}}^{\ast})。在我们的模板中,我们将介绍使用Krylov子空间方法近似t_j^{(m)}的方法。
我们现在将讨论如何使用预处理为像广义最小残差(GMRES)或共轭梯度平方(CGS)这样的迭代求解器,应用于方程(4.49)。假设我们有一个可用的左预处理器{K}用于算子A-\theta_j^{(m)}I,在某种意义上{K}^{-1}(A-\theta_j^{(m)} I) \approx I。由于\theta_j^{(m)}随迭代计数{m}变化,预处理器也可能变化,尽管通常对于不同的\theta值使用相同的{K}是实用的。当然,预处理器{K}也必须限制在与{u}_j^{(m)}正交的子空间中,这意味着我们必须有效地使用
\widetilde{K}\equiv\left(I-{u}_j^{(m)}{{u}_j^{(m)}}^\ast\right){K}\left(I-{u}_j^{(m)}{{u}_j^{(m)}}^\ast\right) .
这可能看起来过于复杂,而且乍一看似乎涉及大量由于投影(I-{u}_j^{(m)}{{u}_j^{(m)}}^{\ast})包围矩阵向量操作的开销。但实际上,这归结为少数相当简单的操作,正如我们现在将展示的。
我们假设使用带有初始猜测{t}_0=0的Krylov求解器和左预处理来近似求解校正方程(4.49)。由于起始向量在{u}_j^{(m)}正交的子空间中,Krylov求解器的所有迭代向量都将在这个空间中。在这个子空间中,我们必须为Krylov求解器提供的向量{v}计算向量{z} \equiv\widetilde{K}^{-1}\widetilde{A}{v},并且
\widetilde{A}\equiv\left(I-{u}_j^{(m)}{{u}_j^{(m)}}^\ast\right)\left(A-\theta_j^{(m)} I\right)\left(I-{u}_j^{(m)}{{u}_j^{(m)}}^\ast\right).
我们将分两步进行。首先我们计算
\widetilde{A}v =\left(I-{u}_j^{(m)}{{u}_j^{(m)}}^\ast\right)\left(A-\theta_j^{(m)} I\right)\left(I-{u}_j^{(m)}{{u}_j^{(m)}}^\ast\right){v}=\left(I-{u}_j^{(m)}{{u}_j^{(m)}}^\ast\right){y}
其中y \equiv (A-\theta_j^{(m)} I){v}因为{{u}_j^{(m)}}^\ast {v}=0。然后,使用左预处理,我们解{z}\perp {u}_j^{(m)}从
\widetilde{K} {z}=\left(I-{u}_j^{(m)}{{u}_j^{(m)}}^\ast\right){y} .
由于{z}\perp {u}_j^{(m)},因此{z}满足K{z}={y}-\alpha {u}_j^{(m)}或{z}=K^{-1}{y} - \alpha K^{-1}{u}_j^{(m)}。条件{z}\perp {u}_j^{(m)}导致
\alpha = \frac{{{u}_j^{(m)}}^\ast K^{-1}y}{{{u}_j^{(m)}}^\ast K^{-1}{u}_j^{(m)}} .
向量\widehat{y}\equiv {K}^{-1}{y}从{K}\widehat{y}={y}解出,同样地,\widehat{u}\equiv{K}^{-1}{u}_j^{(m)}从{K}\widehat{u}={u}_j^{(m)}解出。最后一个方程在方程(4.49)的迭代过程中只需解一次,因此对于{i}_S次迭代的线性求解器,实际上需要{i}_S+1次预处理器的操作。同样地,Krylov求解器中左预处理算子的矩阵向量乘法在一次迭代中只需要一次内积和一次向量更新,而不是投影算子(I-{u}_j^{(m)}{{u}_j^{(m)}}^{\ast})的四次操作。这在算法4.15中给出的解决方案模板中得到了详细说明。沿着类似的思路,可以获得一个虽然比左预处理情况稍微昂贵但仍然高效的模板,用于右预处理算子。这个模板在算法4.16中描述。有关校正方程预处理的更多细节,请参见[412]。
如果我们对方程(4.49)中的t_j^{(m)}形成近似\widetilde{t}^{(m)}\equiv {K}^{-1}r -\alpha {K}^{-1} {u}_j^{(m)},其中\alpha使得\widetilde{t}^{(m)}\perp {u}_j^{(m)},并且没有通过迭代求解器加速,我们得到一个由Olsen、Jørgensen和Simons[344]提出的过程。
下一节:基本算法
上一级:Jacobi-Davidson方法
上一节:Jacobi-Davidson方法
Susan Blackford
2000-11-20