Jacobi-Davidson方法

线性化方法的可能缺点是问题的维度加倍；即，一个包含n维矩阵M、C和K的问题被转化为一个包含2n维矩阵A和B的广义问题。这一点在接下来要讨论的雅可比-戴维森方法中得以避免。在该方法中，二次特征值问题(QEP)首先被投影到一个低维子空间，从而得到一个维度适中的QEP。这个低维投影QEP可以使用任意选择的方法求解。子空间的扩展通过雅可比-戴维森校正方程实现。对于多项式特征值问题，这一技术首次在[408, sec. 8]中提出并讨论。

正如我们将在下文看到的，该方法也可以直接应用于多项式特征值问题

(\lambda^\ell C_\ell+\cdots +\lambda C_1 +C_0 ) x=0, \tag{9.19}

并且无需转化为广义“线性”特征值问题，其中C_i（i = 0, 1, \ldots, \ell）是给定的n \times n矩阵。为简单起见，我们仅介绍\ell = 2的情况。

在雅可比-戴维森迭代步骤的第一部分，用于求解多项式特征值问题

\Psi(\lambda)x =0 \tag{9.20}

其中

(\theta^2 V^\ast C_2 V +\theta V^\ast C_1 V + V^\ast C_0V){s}=0 \tag{9.21}

矩阵V的列{v}_i构成了搜索子空间的基。为了稳定性考虑，这些列被构造为正交归一的。投影问题与原问题具有相同的阶数，但维度要小得多，通常{m}\lln。我们假设解向量{s}是归一化的，\Vert{s}\Vert _2 = 1。首先，选择具有所需特性的Ritz值\theta，例如实部最大或最接近某个目标\tau，并计算相关的特征向量{s}。然后计算Ritz向量{u}\equiv V {s}和残差r\equiv \Psi(\theta){u}。为了扩展搜索空间，计算向量{p}，

{p} \equiv \Psi'(\theta){u},

其中

\Psi'(\theta)= 2 \theta C_2+C_1

在雅可比-戴维森迭代步骤的第二部分，搜索子空间\mathrm{span}({V})通过一个与{u}正交的向量{t}进行扩展，该向量近似地满足校正方程

\left(I-\frac{{p}\, {u}^\ast}{{u}^\ast {p}}\right)\Psi(\theta) \left(I-{u}\, {u}^\ast\right){t}=-r. \tag{9.22}

下一列{V}通过对近似解与先前计算的列{v}_1, \ldots, {v}_{m}进行正交归一化获得。

此过程重复进行，直到检测到特征对(\lambda,x)，即直到残差向量r足够小。算法的基本形式在算法9.1中呈现。我们参考[413]中的例子，一个来自声学问题的二次特征值问题，通过这种简化技术求解，n高达250,000。

算法 9.1：用于二次特征值问题的雅可比-戴维森方法

(1)  选择一个 n \times m 正交归一矩阵 V
(2)  对于 i=0,1,2
(3)      计算 W_{i} = C_{i} V 和 M_{i} = V^{} W_{i}
(4)  结束循环
(5)  迭代直到收敛
(6)      计算特征值问题 (\theta^{2} M_{2} + \theta M_{1} + M_{0}) s = 0 的特征对 (\theta, s)
(7)      选择满足 |s|{2} = 1 的所需特征对 (\theta, s)
(8)      计算 u = V s，p = \Psi^{\prime}(\theta) u，r = \Psi(\theta) u
(9)      如果 |r|{2} < \epsilon，则 \lambda = \theta，x = u，停止
(10)     在 u 的正交补空间中近似求解 \left(I - \frac{p u^*}{u^* p}\right) \Psi(\theta) \left(I - u u^*\right) t = -r
(11)     将 t 对 V 进行正交化，v = t / |t|{2}
(12)     对于 i = 0,1,2
(13)         计算 w{i} = C_{i} v
(12)         M_{i} = V^{} w_{i}，M_{i} = \begin{bmatrix} M_{i} & V^{} w_{i} \\ v^{} W_{i} & v^{*} w_{i} \end{bmatrix}，W_{i} = \left[ W_{i}, w_{i} \right]
(15)     结束循环
(16)     扩展 V = [V, v]

如果搜索子空间的维度变得过大，例如，当V的列数达到{m}={m}_{\max}时，可以继续使用维度较小的搜索子空间。选择由上一步中最佳的{m}_{\min}个Ritz对所张成的空间作为新的搜索子空间；即，取V=[u_1,\ldots,u_{{m}_{\min}}]，其中u_1,\ldots, u_{{m}_{\min}}是与最佳的{m}_{\min}个Ritz值\theta_1,\ldots,\theta_{{m}_{\min}}相关的Ritz向量。然后对V应用修正的Gram-Schmidt正交归一化，并以此矩阵重新开始。注意，新的搜索矩阵V=V_{{m}_{\min}}可以用旧的搜索矩阵V=V_{{m}_{\max}}表示为V_{{m}_{\min}}=V_{{m}_{\max}} T，其中T是一个{m}_{\max}\times {m}_{\min}的矩阵。变换矩阵T可以显式计算，并用于更新辅助矩阵W_i（=W_i T）和M_i（=T^\ast M_i T）。

已经检测到的特征向量可以保留在搜索子空间中（显式降阶），如果需要更多特征向量的话。保留特征向量在搜索子空间中可以防止过程重建已知的特征向量。

在§4.7.3和§8.4.2中，我们建议使用“隐式”降阶，因为该方法基于Schur形式，具有更稳定的校正方程。然而，在这种一般多项式设置中，如何结合隐式降阶尚不明确：Schur形式和广义Schur形式难以轻易推广。