Skip to content

Approximating Eigenvalues

本页面还在施工中

The Power Methods

Principle

幂法的目标是求得矩阵 \(A\in \mathbb{R}^{n\times n}\) 模最大的特征值 \(\lambda_1\) 和其对应的特征向量。作为一个简化，我们只考虑

\[ |\lambda_1| > |\lambda_2| \geqslant \cdots \geqslant |\lambda_n|\geqslant 0 \]

并且 \(A\) 可对角化的情况，这样每个特征值对应的特征向量分别是

\[ v_1, v_2, \cdots, v_n \]

事实上，对于其它一些情况也可以修改幂法后使用，但是这种情况的处理和推导是最方便的，更多的讨论我将留坑到数值代数。

这样，对于 \(\forall x\in \mathbb{R}^n\)，有

\[ x=\sum_{i=1}^n \beta_i v_i \]

令 \(x^{(0)}=x\)，然后不断进行 \(x^{(n+1)}=Ax^{(n)}\) 的迭代，可以得到

\[ \begin{gathered} x^{(1)}=\sum_{i=1}^n \beta_i A v_i=\sum_{i=1}^n \beta_i \lambda_i v_i\\ \cdots\\ x^{(n)}=\sum_{i=1}^n \beta_i \lambda_i^n v_i=\lambda_1^n \left[\beta_1v_1 + \sum_{i=2}^n\beta_i \left(\frac{\lambda_i}{\lambda_1}\right)^n v_i\right] \end{gathered} \]

这样，如果 \(n\) 足够大，就会有

\[ \left.\begin{matrix} x^{(n)}\approx\lambda_1^n \beta_1v_1\\ x^{(n+1)}\approx\lambda_1^{n+1} \beta_1v_1 \end{matrix}\right\} \Rightarrow \lambda_1\approx \frac{x^{(n+1)}_i}{x^{(n)}_i}, \quad i=1,\cdots,n \]

对应的特征向量就选用 \(v_1\approx x^{(n)}\)。这里的下标 \(i\) 是取第 \(i\) 个分量，似乎在理论上取 \(1\) 到 \(n\) 都可以。

Implementation: Normalization with Infinity Norm

实际算法中，会选用下标 \(i\) 使得 \(\left|x_i^{(n)}\right|\) 取到 \(\left\|x^{(n)}\right\|_{\infty}\)，即

\[ i=\mathop{\arg\max}_{1\leqslant i\leqslant n} \left|x_i^{(n)}\right| \]

这里为了让除法更加稳定。并且为了在迭代中尽量稳定，每一步都会使用无穷范数进行规范化 (normalization)。举一个例子，已有 \(x^{(k)}\)，欲求 \(x^{(k+1)}\)，则取

\[ p_k=\mathop{\arg\max}_{1\leqslant p_k\leqslant n} \left|x_{p_k}^{(k)}\right| \]

然后计算规范化向量

\[ u^{(k)} = \frac{x^{(k)}}{\left|x^{(k)}_{p_k}\right|} \]

迭代计算 \(x^{(k+1)}\) 的公式改为

\[ x^{(k+1)}=Au^{(k)} \]

如果此时停止，那么就取 \(v_1=u^{(k)}\)（或者取 \(x^{(k)}\) 也行，然后取

\[ \lambda_1=\frac{x^{(k+1)}_{p_k}}{u^{(k)}_{p_k}}=\left|x^{(k)}_{p_k}\right| \]

Remark

  • 随机初始化的 \(x^{(0)}\) 都可以，即使 \(\beta_1=0\)，也会奏效

    为什么呢？读者可以自己思考一下

  • 对于最大模特征值是重根的情况也适用，但是要求代数重数等于几何重数
  • 如果 \(\lambda_1=-\lambda_2\)，直接套用这个方法将会失效，但是在数值代数中可以扩展修改该方法后继续奏效

Accelerate Power Method

Aitken's \(\Delta^2\) Method 可以应用于加速幂法。但是存在一种对原始的 \(A\) 进行简单变换得到 \(B\)，再对 \(B\) 应用幂法的预处理方法可以简单地加速幂法，其原理是进行特征值的整体平移。

首先关注幂法的收敛速率，通过对幂法的原理分析，很容易可以得到最关键的常数是

\[ \left|\frac{\lambda_1}{\lambda_2}\right| \]

为了加速幂法，一个简单的方法就是尽量减小 \(\left|\frac{\lambda_1}{\lambda_2}\right|\)。我们考虑一种简单情况，\(\lambda_1=1.8\)，\(\lambda_2=1.6\)，\(\lambda_n=-0.4\)，这样就有

\[ \left|\frac{\lambda_1}{\lambda_2}\right|=\frac{9}{8}=1.125 \]

现在，我们取

\[ p=\frac{\lambda_2+\lambda_n}{2}=0.6 \]

这样取 \(B=A-pI\)，这样 \(B\) 的特征值 \(\lambda_1', \lambda_2',\cdots, \lambda_n'\) 有

\[ \lambda_i'=\lambda_i-p,\quad i=1,\cdots,n \]

我们会发现有

\[ \lambda_2'=\frac{\lambda_2-\lambda_n}{2},\quad \lambda_n'=\frac{\lambda_n-\lambda_2}{2} \]

两者的模长相同，则保证了 \(\lambda_1'\) 仍然是模长最大的特征值。

这就是为什么 \(p\) 取 \(\frac{\lambda_2+\lambda_n}{2}\) 而不是 \(\frac{\lambda_1+\lambda_n}{2}\) 的原因

我们会发现，此时有

\[ \left|\frac{\lambda_1'}{\lambda_2'}\right|= \left|\frac{1.2}{1.0}\right|=1.2>1.125 \]

在特征值全是实数的特定条件下观察，这是因为 \(\lambda_1-\lambda_2\) 不变，但是 \(\lambda_2\) 变小了，所以 \(\left|\frac{\lambda_1}{\lambda_2}\right|\) 就减小了。

因此对 \(B\) 应用幂法求得 \(\lambda_1'\) 和特征向量 \(v_1'\)，则可以得到 \(\lambda_1=\lambda_1'+p\)，\(v_1=v_1'\)，而且整个过程可以很好地加速。

Variants of Original Power Method

原始的幂法可以衍生出反幂法 (inverse power method) 和移位幂法等变种。

  • 反幂法：求模长最小的特征值和其对应的特征向量
  • 移位幂法：求和所给的值最接近的特征值和其对应的特征向量

反幂法的原理很简单，\(A\) 的特征值是 \(|\lambda_1|\geqslant |\lambda_2|\geqslant \cdots > |\lambda_n|\)，那么可得 \(A^{-1}\) 的特征值是

\[ \left|\frac{1}{\lambda_n}\right|>\left|\frac{1}{\lambda_{n-1}}\right|\geqslant \cdots \geqslant \left|\frac{1}{\lambda_1}\right| \]

所以对 \(A^{-1}\) 应用幂法，就可以求得 \(A\) 的模长最小的特征值和其对应的特征向量。实际实现时，并不用显式地求出 \(A^{-1}\)，而是将 \(x^{(n+1)}=A^{-1}x^{(n)}\) 改成求解线性方程组

\[ Ax^{(n+1)}=x^{(n)} \]

移位幂法则是基于反幂法再进行一点修改，设要求 \(A\) 和 \(p\) 最接近的特征值和其特征向量，那么可以对 \(A-pI\) 应用反幂法。通过理论分析，我们可以得知 \(A-pI\) 模长最小的特征值再加上 \(p\) 还原之后就是和 \(p\) 最接近的特征值。

Comments