CS205A Lecture 7 Eigenproblems; Algorithms

课程主页：https://graphics.stanford.edu/courses/cs205a-13-fall/schedule.html

这次回顾第七讲，这一讲主要介绍了计算特征值的方法。

首先回顾我们的问题条件：

$$\begin{aligned} A\in \mathbb R^{n\times n}是对称矩阵\\ \vec x_1,...,\vec x_n \in \mathbb R^n 为特征向量\\ |\lambda_1|\ge |\lambda_2| \ge ...\ge |\lambda_n |是特征值 \end{aligned}$$

幂迭代(Power Iteration)

注意到所有的特征向量张成$\mathbb R^n $，所以对于任意向量$\vec v$，我们有

$$\vec v = c_1 \vec x_1 +...+c_n \vec x_n$$

因此

$$\begin{aligned} A\vec v &= c_1 A\vec x_1 +....+c_n A\vec x_n\\ &=c_1 \lambda_1 \vec x _1 +...+ c_n \lambda_n \vec x_n \\ & =\lambda_1 \Big( c_1 \vec x_1 +\frac{\lambda_2} {\lambda_1} c_2 \vec x_2 +...+\frac{\lambda_n} {\lambda_1} c_n \vec x_n \Big) \\ &\ldots \\ A^k \vec v &=\lambda_1^k \Big( c_1 \vec x_1 + \Big (\frac{\lambda_2} {\lambda_1}\Big )^k c_2 \vec x_2 +...+\Big (\frac{\lambda_n} {\lambda_1}\Big )^k c_n \vec x_n \Big) \end{aligned}$$

这里假设特征值互不相同，那么

$$\lim_{k\to \infty} \Big (\frac{\lambda_i} {\lambda_1}\Big )^k=0, 2\le i \le n$$

所以当$k\to \infty $时，

$$A^k \vec v \approx \lambda_1^k \vec v$$

于是产生如下算法：

幂迭代(power iteration)

1. 选择$\vec v_1 \in \mathbb R^n$为任意非零向量。

2. 迭代直至收敛：

   $$\vec v_k =A\vec v_{k-1}$$

显然，如果$|\lambda_1 |>1$，那么

$$\lim_{k\to \infty} \Arrowvert \vec v_k \Arrowvert \to \infty$$

这样就会产生计算问题。注意到对于特征向量，我们只在意其方向，所以产生了如下算法：

正规化幂迭代(normalized power iteration)

1. 选择$\vec v_1 \in \mathbb R^n$为任意非零向量。

2. 迭代直至收敛：

   $$\begin{aligned} \vec w_k &=A\vec v_{k-1} \\ \vec v_k &=\frac{\vec w_k } {\Arrowvert \vec w_k \Arrowvert} \end{aligned}$$

注意，该算法返回模最大的特征值。

反向迭代(Inverse Iteration)

这里假设$A$可逆，那么如果

$$A\vec x =\lambda \vec x$$

我们有

$$\begin{aligned} \vec x &= \lambda A^{-1} \vec x \\ A^{-1} \vec x& =\frac 1 \lambda \vec x \end{aligned}$$

因此$A^{-1}$的特征值满足

$$\frac 1{|\lambda_n |} \ge \frac 1{|\lambda_{n-1} |} \ge \ldots \ge \frac 1{|\lambda_1 |}$$

所以利用上述事实找到$\frac 1 {\lambda_n }$，即模最小的特征值，算法如下：

反向幂迭代(inverse power iteration)

1. 选择$\vec v_1 \in \mathbb R^n$为任意非零向量。
2. 迭代直至收敛：
   • 求解$\vec w_k$：$A\vec w_k =\vec v_{k-1}$
   • 正规化：$\vec v_k =\frac{\vec w_k }
     {\Arrowvert \vec w_k \Arrowvert}$

利用$LU$分解，可以将提升上述算法的效率：

1. 分解$A=LU$
2. 选择$\vec v_1 \in \mathbb R^n$为任意非零向量。
3. 迭代直至收敛：
   • 用forward substitution求解$\vec y_k$：$L\vec y_k =\vec v_{k-1}$
   • 用back substitution求解$\vec w_k$：$U \vec w_k =\vec y_k$
   • 正规化：$\vec v_k =\frac{\vec w_k }{\Arrowvert \vec w_k \Arrowvert}$

注意，该算法返回模最小的特征值。

平移(Shifting)

首先回顾恒等式：

$$A^k \vec v =\lambda_1^k \Big( c_1 \vec x_1 + \Big (\frac{\lambda_2} {\lambda_1}\Big )^k c_2 \vec x_2 +...+\Big (\frac{\lambda_n} {\lambda_1}\Big )^k c_n \vec x_n \Big)$$

所以收敛速度主要由$\Big (\frac{\lambda_2} {\lambda_1}\Big )$控制，如果这个比值接近$1$，那么收敛速度会大大减慢，但是注意到我们有

$$(A-\sigma) \vec x_i =(\lambda_i- \sigma) \vec x_i$$

所以我们可以计算$A-\sigma$的特征值

$$\lambda_1-\sigma,...,\lambda_n -\sigma$$

如果

$$\Big |\frac{\lambda_1 -\sigma}{\lambda_2 -\sigma} \Big | < \Big |\frac{\lambda_1 }{\lambda_2}\Big |$$

那么我们可以通过计算$A-\sigma$的特征值来更快得到$A$的特征值，但由于$\lambda_i $未知，所以$\sigma$的选择很困难。

接着用来另一种角度看待特征值问题，实际上我们是在找$\lambda $，使得

$$A\vec x \approx \lambda \vec x$$

所以对于$\vec v_0$，我们去猜测$\vec v_0$对应的特征值，即找到

$$\arg \min_{\lambda} \Arrowvert A\vec v_0- \lambda \vec v_0 \Arrowvert^2$$

记

$$\begin{aligned} f(\lambda) &\triangleq \Arrowvert A\vec v_0- \lambda \vec v_0 \Arrowvert^2 \\ &=( A\vec v_0- \lambda \vec v_0 )^T( A\vec v_0- \lambda \vec v_0 )\\ &=\vec v_0^T A^T A\vec v_0 -2\lambda \vec v_0^T A\vec v_0 +\lambda^2 \vec v_0^T \vec v_0 \end{aligned}$$

那么求梯度并令其为$0$可得

$$\begin{aligned} \nabla_{\lambda} f(\lambda) &= 2\lambda \vec v_0^T \vec v_0 -2 \vec v_0^T A\vec v_0=0\\ \lambda&= \frac{ \vec v_0^T A\vec v_0}{ \vec v_0^T \vec v_0} =\frac{ \vec v_0^T A\vec v_0}{ \Arrowvert \vec v_0 \Arrowvert^2_2} \end{aligned}$$

上述量被称为瑞利商(Rayleigh quotient)，有了这个量，我们可以认为$A-\lambda I_n$的最小特征值接近于$0$，于是用反幂法可以高效求解该特征值，于是产生了如下算法：

瑞利商迭代(Rayleigh quotient iteration)

1. 选择$\vec v_1 \in \mathbb R^n$为任意非零向量或者猜测一个特征向量。

2. 迭代直至收敛：

   • 计算对当且特征值的估计：

     $$\sigma_k = \frac{ \vec v_{k-1}^T A\vec v_{k-1}}{ \Arrowvert \vec v_{k-1} \Arrowvert^2_2}$$

   • 求解$\vec w_k$：$(A-\sigma_k I_n) \vec w_k =\vec v_{k-1}$

   • 正规化：$\vec v_k =\frac{\vec w_k}{\Arrowvert \vec w_k \Arrowvert}$

注意上述方法有一个缺点，因为$A-\sigma_k I_n$会改变，所以无法使用$LU$分解提升运行速度。

找到多个特征值

之前的算法都是介绍如何找到一个特征值，但是实际问题中往往要求解多个特征值，这一部分将介绍这点。

收缩法(Deflation)

考虑power iteration的过程，假设$\lambda_1 $为模最大的特征值，其对应的特征向量为$\vec x_1$，选择的初值为$\vec v$，注意到

$$\vec v = c_1 \vec x_1 +...+c_n \vec x_n$$

如果

$$\vec v ^T \vec x_1 =0$$

那么

$$c_1=0$$

所以

$$\begin{aligned} A^k \vec v &= c_2{\lambda_2^k} \vec x_2 +...+c_n{\lambda_n^k} \vec x_n \\ &=\lambda_2^k \Big( c_2 \vec x_2 +...+\Big (\frac{\lambda_n} {\lambda_2}\Big )^k c_n \vec x_n \Big) \end{aligned}$$

因此无论迭代多少次，都无法求出$\lambda_1 $，但是可以求出$\lambda_2$，这就提醒我们可以通过投影的方法逐步求出每个特征值，于是产生如下算法：

• 对$l=1,2,…$

  • 选择$\vec v_1 \in \mathbb R^n$为任意非零向量。

  • 迭代直至收敛：

    • 减去特征值上的投影：

      $$\vec u_k =\vec v_{k-1}-\text{proj}_{\text{span}\{\vec x_1,...,\vec x_{l-1}\}} \vec v_{k-1}$$

    • 计算$\vec w_k =A\vec u_k$

    • 正规化：$\vec v_k =\frac{\vec w_k}{\Arrowvert \vec w_k \Arrowvert}$

  • 将迭代的结果加入$\vec x_i$的集合中。

注意到如果$A$非对称，那么其特征向量不正交，上述算法会失效，此时就要使用别的算法。

现在假设

$$\begin{eqnarray*} A\vec x_1 &&=\lambda_1 \vec x_1 \tag 1\\ \Arrowvert \vec x_1\Arrowvert &&=1 \end{eqnarray*}$$

记$H$是Householder矩阵，使得

$$H\vec x_1 =\vec e_1 \tag 2$$

注意到相似变换不会影响特征值，而$H$是对称正交矩阵，所以我们考虑

$$HAH^T$$

右乘$\vec e_1$，可得

$$\begin{aligned} HAH^T \vec e_1&= HAH \vec e_1 &H对称\\ &=HAHH\vec x_1 &由(2)\\ &=HA\vec x_1&由 HH=HH^T=I_n \\ &=\lambda_1 H\vec x_1 &由(1)\\ &=\lambda_1 \vec e_1 &由(2) \end{aligned}$$

所以$HAH^T \vec e_1$的第一列是$\lambda_1 \vec e_1 $，即

$$HAH^T = \left( \begin{matrix} \lambda_1 & \vec b^T \\ \vec 0 & B \end{matrix} \right)$$

其中$B\in \mathbb R^{(n-1)\times (n-1)}$并且其特征值为$\lambda_2,…,\lambda_n$，所以可以对$B$使用power iteration，但算法依然是一次只计算一个特征值，后续将介绍如何一次计算多个特征值的算法。

$QR$迭代

依然利用相似矩阵不改变特征值，特别地，取相似矩阵为正交阵，那么

$$Q^{-1} AQ =Q^TAQ$$

那么该选择什么样的$Q$呢？考虑$QR$分解：

$$A=QR$$

那么

$$\begin{aligned} Q^T AQ &=Q^T QR Q\\ &=RQ \end{aligned}$$

这就产生如下算法：

1. 选择$A_1= A$
2. 对$k=1,2,…$
   • 分解$A_k =Q_kR_k$
   • 令$A_{k+1}=R_k Q_k$

由之前的推导可得$A_k$的特征值和$A$相同，如果上述算法收敛，那么

$$A_{\infty} =Q_{\infty}R_{\infty}=R_{\infty}Q_{\infty}$$

所以

$$A_{\infty}Q_{\infty}=Q_{\infty}R_{\infty}Q_{\infty}=Q_{\infty}A_{\infty}$$

为方便讨论，这里假设$A_{\infty }$和$Q_{\infty }$各自的特征值互不相同。

现在任取$Q_{\infty }$的特征值$\lambda $以及对应的特征向量$\vec x$，我们有

$$\begin{aligned} Q_{\infty}A_{\infty} \vec x &=A_{\infty}Q_{\infty}\vec x\\ &= \lambda A_{\infty}\vec x \end{aligned}$$

如果

$$A_{\infty}\vec x = \vec 0$$

那么$\vec x$是矩阵$A_{\infty}$属于特征值$0$的特征向量。

如果

$$A_{\infty}\vec x \neq \vec 0$$

$A_{\infty}\vec x$是矩阵$R_{\infty}$属于$\lambda$的特征向量，即

$$A_{\infty}\vec x =\lambda_1 \vec x$$

无论哪种情形，都有$\vec x $是$A_{\infty}$的特征向量，所以$Q_{\infty }$的特征向量都是$A_{\infty }$的特征向量，同理可得$A_{\infty }$的特征向量都是$Q_{\infty }$的特征向量，即两者的特征向量相同。

现在任取$Q_{\infty }$的特征向量$\vec x$，由正交性可得

$$Q_{\infty} \vec x =\pm \vec x$$

我们取$+1$对应的特征向量，所以

$$Q_{\infty} \vec x = \vec x \tag 3$$

假设该向量对应于$A$的特征值为$\lambda$，即

$$A_{\infty} \vec x =\lambda \vec x$$

因此

$$\begin{aligned} \lambda \vec x &=A_{\infty}\vec x \\ &=R_{\infty}Q_{\infty}\vec x \\ &=R_{\infty} \vec x&由(3) \end{aligned}$$

所以$R_{\infty}$的特征值和$A_{\infty }$的特征值相同，注意到$R_{\infty } $为上三角矩阵，所以其特征值为对角线元素，因此我们只要计算$R_{\infty } $的对角元即可。