lanczos算法及C++实现（三）实对称三对角阵特征值分解的分治算法

1. 基本思想

在第一篇中，我们讨论了lanczos算法的基本框架。当我们用lanczos算法将一个实对称阵转化成三对角阵之后，我们可以用第二篇中的QR算法计算三对角阵的特征值特征向量。

本篇我们将讨论计算该三对角阵更加快速的算法——分治法（Divide and Conquer），该算法最早由Cuppen于1981年提出。

给定实对称三对角阵

$T=\left(\begin{matrix}\alpha_1 & \beta_1 & & & \dots\\ \beta_1 & \alpha_2 & \beta_2 & & \dots \\ & \beta_2 &\alpha_3 & \beta_3 & \dots \\ & & & & \dots \end{matrix}\right)$

分治法的基本思想是将T分成两块三对角阵，求得这两块三对角阵的特征值特征向量之后，再合并，求合并之后的特征值特征向量:

分（Divide）：

Block diagonal plus correction.png

求 $\hat{T_1},\hat{T_2}$ 的特征值分解：

$\hat{T_1}=Q_1 D_1 Q^T_1$

$\hat{T_2}=Q_2 D_2 Q^T_2$

合并（Conquer）

T可以表示为：

$T=\left(\begin{matrix}Q_1 & \\ & Q_2\end{matrix}\right)(D+\beta zz^T) \left(\begin{matrix}Q^T_1 & \\ & Q^T_2\end{matrix}\right)$

其中， $D=\left(\begin{matrix}D_1 & \\ & D_2\end{matrix}\right), z^{{T}}=(q_{{1}}^{{T}},q_{{2}}^{{T}})$ , $q_{{1}}^{{T}}$ 是 $Q_{{1}}$ 的最后一行， $q_{{2}}^{{T}}$ 是 $Q_{{2}}$ 的第一行

现在要求T的特征值分解， $T=Q \Lambda Q^T$ ，可以先求 $D+\beta zz^T$ 的特征值分解 $D+\beta zz^T=P \Lambda P^T$ ，再令 $Q=\left(\begin{matrix}Q_1 & \\ & Q_2\end{matrix}\right)P$

如何求 $D+\beta zz^T$ 的特征值分解呢？如果 $\beta=0$ 则显然P=I。下面假设 $\beta\ne 0$

假设 $\lambda\ne 0$ 是特征值， $p\ne 0$ 是对应的特征向量，根据特征值定义我们有

$(D+\beta zz^T)p=\lambda p$

$(D-\lambda I)p + \beta zz^Tp=0$

$p + (D-\lambda I)^{-1}\beta zz^Tp=0$

$z^Tp + z^T(D-\lambda I)^{-1}\beta zz^Tp=0$

假如D所有对角线元素各不相同，并且z的所有元素不为0 。可以证明 $z^Tp\ne 0$ 。

反证法，否则 $(D+\beta zz^T)p=Dp = \lambda p$ ，p是D的特征向量。因为D对角线元素各不相同，因此p只能有一个非0分量，而由于z所有分量非0，所以得到 $z^Tp\ne 0$

由此我们有

$1 + \beta z^T(D-\lambda I)^{-1} z=0$

或者写成如下形式：

$1 + \beta \sum_i \frac{z_i^2}{d_i-\lambda}=0$

其中 $d_i = D_{ii}$

因此求特征值问题最后归结为求解一个方程的根。这个方程被称为"Secular Equation"。该方程的解法下面会详细讨论，这里略过。得到特征值之后，根据上面第二行 $(D-\lambda I)p + \beta zz^Tp=0$ ，可以得出 $(D-\lambda I)^{-1}z$ 是特征值 $\lambda$ 对应的特征向量。

至此我们得到了分治法的框架，总结如下：

[Q,D] = DivideAndConquer(T)
//输入T实对称三对角阵，输出T的特征值分解T=Q D Q'，Q为正交阵，每列为T的特征向量，D为对角阵，对角线元素为T的特征值
if T is 1x1 matrix
         Q=1, D=T
else
         Divide
                  Block diagonal plus correction.png
         [Q1,D1]=DivideAndConquer( $\hat{T_1}$ )
         [Q2,D2]=DivideAndConquer( $\hat{T_2}$ )
         let $D=\left(\begin{matrix}D_1 & \\ & D_2\end{matrix}\right), z^{{T}}=(q_{{1}}^{{T}},q_{{2}}^{{T}}), Q_{{1}}$ 的最后一行和 $Q_{{2}}$ 的第一行拼接成的向量
         Find eigen values $\lambda=[\lambda_1, \lambda_2 \dots ]^T$ of $D+\beta zz^T$ by solving Secular Equation
                  $1 + \beta \sum_i \frac{z_i^2}{d_i-\lambda}=0$
         Calculate corresponding eigen vectors
                  $p_i=(D-\lambda_i I)^{-1}z$
         Update eigen vectors of T:
                 $Q=\left(\begin{matrix}Q_1 & \\ & Q_2\end{matrix}\right)P$
        where $P=[p_1,p_2,\dots ]$
        return D=diag( $\lambda$ ), Q
end

2. Secular Equation的解法

这节讨论secular方程 $1 + \beta \sum_{i=1}^n \frac{z_i^2}{d_i-\lambda}=0$ 的解法，这里假设 $\forall i, z_i \ne 0,\qquad \forall i\ne j, d_i\ne d_j$

任选一个i,作变量代换：

$\lambda = d_i + \beta \mu$

$\delta_j=\frac{d_j-d_i}{\beta}$

这样secular方程就变为

$h_i(\mu)=1 + \sum_{i=j}^n \frac{z_j^2}{\delta_j-\mu}=0$

这个方程有n个根，且根各不相同，按照从小到大的顺序，记这n个根为 $\mu_1< \mu_2 < \dots <\mu_n$ 。

不失一般性，假设 $\delta_1< \delta_2 < \dots < \delta_{i-1} < 0=\delta_{i} < \delta_{i+1} < \dots <\delta_n$ ，则有 $\delta_j < \mu_j < \delta_{j+1}$ ，对于j=n的情况， $\delta_n < \mu_n < +\infty =\delta_{n+1}$

举个例子，函数 $h_3(\mu) = 1+\frac{0.1}{-2-\mu}+\frac{0.5}{-0.5-\mu}+\frac{0.2}{-\mu}+\frac{0.2}{1-\mu}+\frac{0.4}{2-\mu}$ 的图像如下

可以看到 $0 < \mu_i < \delta_{i+1}$ ，我们可以在区间 $(0,\delta_{i+1})$ 内求解第i个根。我们把这个区间放大，得到如下图像。

可以看到这段区间内，函数h单调增，与x轴有1个交点。如何求这个点呢，最简单的是2分法，可以在log时间内收敛。

更快的是牛顿迭代法。但是这里注意一个问题，就是如果初始点选的不好，牛顿法会发散，找不到根。如何保证牛顿法收敛呢？这里我们需要对变量作一个替换，

令 $\gamma=\frac1{\mu}$

这样，原来的定义域 $0<\mu_i<\delta_{i+1}$ 就变成了 $\gamma_i>\frac1{\delta_{i+1}}$ ，对于i=n的情况， $\gamma_n>0$

现在我们看一下新的函数 $h_i(\gamma)$ 在区间 $(\frac1{\delta_{i+1}},\infty)$ 上的图像.

可见，在区间 $(1,\infty)$ 上，函数 $h_i(\gamma)$ 单调递减。假设 $h_i(\gamma)$ 的根为 $\gamma^*$ ，论文Numerical solution of a secular equation证明了，如果初始点选在 $\gamma^*$ 的左侧，则牛顿法必定收敛。换句话说，如果初始点 $\gamma$ 满足 $h_i(\gamma)>0$ ，则用牛顿法得到的 $\gamma$ 序列收敛到 $\gamma^*$

如何找这个初始点呢？论文Numerical solution of a secular equation构造了 $h_i(\gamma)$ 的一个下界函数 $g_i(\gamma)\le h_i(\gamma)$ ，通过求解 $g_i(\gamma)=0$ 得到根 $\gamma$ ，则有 $0=g_i(\gamma)\le h_i(\gamma)$ ，从而求得初始点。

现在我们推出这个下界函数。先讨论i<n的情况，然后讨论i=n的情况。

我们直接通过函数 $h_i(\mu)$ 来推出下界函数。根据前面所述，

$h_i(\mu)=1 + \sum_{i=j}^n \frac{z_j^2}{\delta_j-\mu}$

$=1 + \sum_{j\ne i}^n \left(\frac{z_j^2}{\delta_j} + \frac{z_j^2/\delta_j^2}{1/\mu - 1/\delta_j} \right) -\frac{z_i^2}{\mu}$

$=1 + \sum_{j<i}^n \left(\frac{z_j^2}{\delta_j} + \frac{z_j^2/\delta_j^2}{1/\mu - 1/\delta_j} \right) + \sum_{j>i+1}^n \left(\frac{z_j^2}{\delta_j} + \frac{z_j^2/\delta_j^2}{1/\mu - 1/\delta_j} \right) -\frac{z_i^2}{\mu} + \frac{z_{i+1}^2}{\delta_{i+1}} + \frac{z_{i+1}^2/\delta_{i+1}^2}{1/\mu - 1/\delta_{i+1}}$

当j<i时， $\delta_j < 0$ ，因此 $1/\mu-1/\delta_j>0$

当j>i时， $\mu<\delta_j$ ，同样有 $1/\mu-1/\delta_j>0$

所以

$h_i(\mu)\ge 1+\sum_{j\ne i}\frac{z_j^2}{\delta_j} + \frac{z_{i+1}^2/\delta_{i+1}^2}{1/\mu - 1/\delta_{i+1}} - \frac{z_i^2}{\mu}$

这时，用 $\gamma=\frac1{\mu}$ 带入，得到

$h_i(\gamma)\ge 1+\sum_{j\ne i}\frac{z_j^2}{\delta_j} + \frac{z_{i+1}^2/\delta_{i+1}^2}{\gamma - 1/\delta_{i+1}} - z_i^2\gamma$

令右边等于0，我们得到一个方程，通过通分，去分母，最后转成一个二次方程：

$z_i^2\gamma^2 - \left(1+\frac{z_i^2}{\delta_{i+1}} +\sum_{j\ne i}\frac{z_j^2}{\delta_j}\right)\gamma + \left(1+\sum_{j\ne i}\frac{z_j^2}{\delta_j}\right)/\delta_{i+1} - \frac{z_{i+1}^2}{\delta_{i+1}^2}=0$

用求根公式可以得到两个解，注意定义域 $\gamma_i>\frac1{\delta_{i+1}}$ ，因此取大根即可。

当i=n时，此时的定义域为 $\gamma_n>0$ ,其实 $h_n(\gamma=0)=h_n(\mu=\infty)=1$ 也是有定义的，而且此时 $\gamma=0$ 是该函数的最左边的点，因此直接选取 $\gamma=0$ 作为初始点即可。

之后牛顿迭代：

$\gamma(k+1)=-\frac{h_i(\gamma(k))}{h'_i(\gamma(k))}+\gamma(k)$

一般两到三次迭代即可收敛，可以认为是O(1) 时间。

这里只讨论了分治算法最基本的实现，论文A Serial Implementation of Cuppen's Divide and Conquer Algorithm for the Symmetric Eigenvalue Problem还提到了用deflation加速，此外分治算法还有稳定性等问题。这里就不详细探讨了。