算法

带状Lanczos算法的完整描述如下。

算法 7.16：NHEP的带状Lanczos方法

(1)  对于 k = 1, 2, \ldots, m，设置 \widehat{v}_k = b_k
      对于 k = 1, 2, \ldots, p，设置 \widehat{w}_k = c_k
      设置 m_c = m, p_c = p，并且 \mathcal{I}_v = \mathcal{I}_w = \emptyset
(2)  对于 j = 1, 2, \ldots，直到收敛或 m_c = 0 或 p_c = 0 执行以下步骤：
(3)      计算 \left\|\widehat{v}_j\right\|_2
(4)      判断是否应该对 \widehat{v}_j 进行消去，如果是，执行以下操作：
             (a) 如果 j - m_c > 0，设置 \mathcal{I}_w = \mathcal{I}_w \cup \{j - m_c\}
             (b) 设置 m_c = m_c - 1。如果 m_c = 0，设置 j = j - 1 并停止
             (c) 对于 k = j, j + 1, \ldots, j + m_c - 1，设置 \widehat{v}_k = \widehat{v}_{k + 1}
             (d) 返回步骤(3)
(5)      计算 \left\|\widehat{w}_j\right\|_2
(6)      判断是否应该对 \widehat{w}_j 进行消去，如果是，执行以下操作：
             (a) 如果 j - p_c > 0，设置 \mathcal{I}_v = \mathcal{I}_v \cup \{j - p_c\}
             (b) 设置 p_c = p_c - 1。如果 p_c = 0，设置 j = j - 1 并停止
             (c) 对于 k = j, j + 1, \ldots, j + p_c - 1，设置 \widehat{w}_k = \widehat{w}_{k + 1}
(7)      设置 t_{j, j - m_c} = \left\|\widehat{v}_j\right\|_2 和 \tilde{t}_{j, j - p_c} = \left\|\widehat{w}_j\right\|_2
          对 v_j = \widehat{v}_j / t_{j, j - m_c} 和 w_j = \widehat{w}_j / \tilde{t}_{j, j - p_c} 进行归一化
(8)      计算 \delta_j = w_j^T v_j。如果 \delta_j = 0，停止
(9)      计算 \widehat{v}_{j + m_c} = A v_j
(10)     设置 k_v = \max\{1, j - p_c\}
           设置 \mathcal{I}_v^{(e)} = \mathcal{I}_v \cup \{k_v, k_v + 1, \ldots, j - 1\}
           对于 \mathcal{I}_v^{(e)} 中的 k（按升序），设置
               t_{k j} = w_k^T \widehat{v}_{j + m_c} / \delta_k 和 \widehat{v}_{j + m_c} = \widehat{v}_{j + m_c} - v_k t_{k j}
(11)     计算 \widehat{w}_{j + p_c} = A^T w_j
(12)     设置 k_w = \max\{1, j - m_c\}
           设置 \mathcal{I}_w^{(e)} = \mathcal{I}_w \cup \{k_w, k_w + 1, \ldots, j - 1\}
           对于 \mathcal{I}_w^{(e)} 中的 k（按升序），设置
               \widetilde{t}_{k j} = \widehat{w}_{j + p_c}^T v_k / \delta_k 和 \widehat{w}_{j + p_c} = \widehat{w}_{j + p_c} - w_k \widetilde{t}_{k j}
(13)    对于 k = j + 1, j + 2, \ldots, j + m_c，设置
               t_{j, k - m_c} = w_j^T \widehat{v}_k / \delta_j 和 \widehat{v}_k = \widehat{v}_k - v_j t_{j, k - m_c}
(14)    对于 k = j + 1, j + 2, \ldots, j + p_c，设置
               \widetilde{t}_{j, k - p_c} = \widehat{w}_k^T v_j / \delta_j 和 \widehat{w}_k = \widehat{w}_k - w_j \widetilde{t}_{j, k - p_c}
(15)    对于 \mathcal{I}_w 中的 k，设置 t_{j k} = \tilde{t}_{k j} \delta_k / \delta_j
(16)    设置 T_j^{(pr)} = [t_{i k}]_{i, k = 1, 2, \ldots, j}
(17)    测试收敛性
(18) 结束循环

我们强调，非厄米带状Lanczos方法可以直接按照算法7.16的描述来实现，并结合下面给出的进一步细节。为了将描述长度控制在一页内，在算法7.16中，T_j和\tilde{T}_j的带状部分的所有可能非零项都计算为内积。

然而，实际上大约只有这些项的一半需要明确计算为内积。其余的项可以通过以下关系获得：

t_{ik}=\tilde{t}_{ki}\, \delta_k / \delta_i, \tag{7.79}

这一关系源自T_j^{\rm (pr)}和\tilde{T}_j^{\rm (pr)}的连接7.71。特别是，在步骤(10)、(12)、(13)和(14)的讨论中，我们给出了使用7.79尽可能多的T_j和\tilde{T}_j项的公式。采用这些公式可以最小化内积的数量，但会牺牲一些数值稳定性。

接下来，我们将更详细地讨论算法7.16的一些步骤。

(4)

通常，是否需要对\hat{v}_j进行收缩的决策应基于检查是否

\Vert\hat{v}_j\Vert _2 \leq {\tt dtol}, \tag{7.80}

其中{\tt dtol}是一个适当选择的小收缩容差。如果7.80满足，则对\hat{v}_j进行收缩。在这种情况下，如果正值j-m_c，则添加到包含\tilde{T}_j的\tilde{T}_j^{\rm {(d)}}部分中非零行的索引集{\mathcal I}_w中，并且当前右块大小更新为m_c=m_c-1。如果m_c=0，则右块Krylov序列7.60耗尽，算法自然终止。如果m_c>0，删除向量\hat{v_j}，剩余右候选向量\hat{v}_{k+1}的索引k+1重置为k，最后，算法返回到步骤(3)。如果不满足7.80，则不需要收缩，算法继续步骤(5)。

(6)

类似于7.80，是否需要对\hat{w}_j进行收缩的决策基于检查是否

\Vert\hat{w}_j\Vert _2 \leq {\tt dtol}. \tag{7.81}

如果7.81满足，则对\hat{w}_j进行收缩。在这种情况下，如果正值j-p_c，则添加到包含T_j的T_j^{\rm {(d)}}部分中非零行的索引集{\mathcal I}_v中，并且当前左块大小更新为p_c=p_c-1。如果p_c=0，则左块Krylov序列7.61耗尽，算法自然终止。如果p_c>0，删除向量\hat{w}_j，剩余左候选向量\hat{w}_{k+1}的索引k+1重置为k，最后，算法返回到步骤(5)。如果不满足7.81，则不需要收缩，算法继续步骤(7)。

(7)

向量\hat{v}_j和\hat{w}_j都通过了收缩检查，现在被归一化以成为下一个右和左Lanczos向量v_j和w_j。归一化使得

\left\Vert v_j \right\Vert _2 = \left\Vert w_j \right\Vert _2 = 1\quad \mathrm{for all}\quad j.

(8)

在这里，我们计算\delta_j并检查是否发生中断。如果\delta_j=0，则需要前瞻以继续算法。

(9)

在这一步中，我们通过计算新的候选向量\hat{v}_{j+m_c}，作为最新右Lanczos向量v_j的A倍数，来推进右块Krylov序列。

(10)

向量\hat{v}_{j+m_c}与左Lanczos向量w_k，k\in {\mathcal I}_v^{(\rm e)}进行双正交化。注意，这一双正交化是通过TSMGS过程完成的。通过使用7.79，可以在计算t_{kj}时节省一个内积。更确切地说，应使用以下t_{kj}的公式：

t_{kj}= \begin{cases} \tilde{t}_{jk}\, \delta_j/\delta_k & \text{如果 } k = j- p_c; \\ {w_k^T \hat{v}_{j+m_c}}/{\delta_k} & \text{其他情况}.\end{cases} \tag{7.82}

注意，7.82中所需的项\tilde{t}_{jk}可从步骤(7)获得。

(11)

在这一步中，我们通过计算新的候选向量\hat{w}_{j+p_c}，作为最新左Lanczos向量w_j的A^T倍数，来推进左块Krylov序列。

(12)

向量\hat{w}_{j+p_c}与右Lanczos向量v_k，k\in {\mathcal I}_w^{(\rm e)}进行双正交化。注意，这一双正交化是通过TSMGS过程完成的。同样，为了节省一个内积，应使用以下\tilde{t}_{kj}的公式：

\tilde{t}_{kj}= \begin{cases} t_{jk}\, \delta_j/\delta_k & \text{如果 } k = j - m_c; \\ {\hat{w}_{j+p_c}^T v_k}/{\delta_k} & \text{其他情况}.\end{cases}

(13)

右候选向量\hat{v}_{j+1},\hat{v}_{j+2},\ldots,\hat{v}_{j+m_c}与最新左Lanczos向量w_j进行双正交化。为了节省内积，可以使用以下t_{j,k-m_c}的公式：

t_{j,k-m_c} = \begin{cases} {w_j^T \hat{v}_{k}}/{\delta_j} & \text{如果 } k - m_c \le 0 \text{ 或 } k - m_c = j; \\ \tilde{t}_{k-m_c,j}\, \delta_{k-m_c}/{\delta_j} & \text{其他情况}.\end{cases}

然而，使用这一公式会牺牲一些数值稳定性。^[1]

(14)

左候选向量\hat{w}_{j+1},\hat{w}_{j+2},\ldots,\hat{w}_{j+p_c}与最新右Lanczos向量v_j进行双正交化。为了节省内积，可以使用以下\tilde{t}_{j,k-p_c}的公式：

\tilde{t}_{j,k-p_c} = \begin{cases} {\hat{w}_{k}^T v_j}/{\delta_j} & \text{如果 } k - p_c \le 0; \\ t_{k-p_c,j}\, \delta_{k-p_c}/{\delta_j} & \text{其他情况}.\end{cases}

然而，使用这一公式会牺牲一些数值稳定性。

(15)

在这一步中计算的t_{jk}项是由于\hat{v}_k向量收缩引起的垂直尖刺导致的T_j^{\rm {(pr)}}的第j行中的潜在非零项。注意，再次使用公式7.79。

(16)

现在，T_j^{\rm {(pr)}}的第j行和第j列中的所有潜在非零元素都已计算完毕，并将其添加到前一次迭代j-1的矩阵T_{j-1}^{(\rm pr)}中，以产生当前矩阵T_{j}^{(\rm pr)}。这里，我们约定在算法7.16中未明确定义的t_{ik}项设置为零。我们注意到，在初始迭代中，即只要j\leq m_c，分别j\leq p_c，算法7.16还产生负索引k\leq 0的t_{jk}项，分别\tilde{t}_{jk}项。这些项来自起始向量的双正交化，它们不是矩阵T_j^{\rm {(pr)}}的一部分。特别是，在进行特征值计算时，这些项不需要。然而，当算法7.16应用于降阶建模时，这些项至关重要，我们将在下面的§7.10.4中讨论。

(17)

对于特征值计算，通过计算j\times j矩阵T_j^{\rm {(pr)}}的特征值\theta_i^{(j)}，i=1,2,\ldots,j，来测试收敛性。如果某些\theta_i^{(j)}足够好地近似于A的期望特征值，则停止算法。对于降阶建模，如果算法生成的j阶模型足够好地近似于原始线性动力系统，则停止算法。

1. 在R. W. 弗罗因德(R. W. Freund)于2000年8月撰写的《使用不完全缩减的Krylov子空间方法进行降阶建模》一文中，描述了一种增强数值稳定性但需额外进行一些内积运算的算法版本。该文作为贝尔实验室数值分析手稿的一部分，详细探讨了如何在牺牲一定计算效率的前提下，提升算法在处理复杂问题时的稳定性能。