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
本系列目錄：