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）：
求$\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 （这两个假设来自wiki，但是感觉没什么道理，anyway先pass）。可以证明$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
[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加速，此外分治算法还有稳定性等问题。这里就不详细探讨了。
3. 参考文献
https://en.wikipedia.org/wiki/Divide-and-conquer_eigenvalue_algorithm
A. Melman Numerical solution of a secular equation
J Rutter A Serial Implementation of Cuppen's Divide and Conquer Algorithm for the Symmetric Eigenvalue Problem Chapter 3。本文的一些插图来自这篇文章。
本文属作者原创，转载请注明出处
http://www.cnblogs.com/qxred/p/dcalgorithm.html
本系列目录：
