幂法和反幂法(C++)

幂法

设矩阵$A\in {\mathbb{R}}^{n\times n}$满足
$$|{\lambda }_1|>|{\lambda }_2|???|{\lambda }_n|$$
对应的$n$个特征向量$x_1,x_2,?,x_n$线性无关. 则称模最大的特征值${\lambda }_1$为矩阵$A$的主特征值, 称对应的特征向量$x_1$为主特征向量.

幂法就是用于求主特征值和主特征向量的一种数值方法.

算法描述

首先任取一个非零的初始向量$v_0$, 由矩阵$A$构造一个向量序列如下:
$$v_k=A{v}_{k?1}=?=A^k{v}_0,k=1,2,?$$
由假设, $v_0$可以表示为
$$v_0={\alpha }_1{x}_1+{\alpha }_2{x}_2+?+{\alpha }_n{x}_n$$
若记$v_k^{(i)}$为$v_k$的第$i$个分量, 则有
$$\begin{array}{rl}{v}_k& =A^k{v}_0\\ & =\sum _{i=1}^n{\alpha }_i{\lambda }_i^k{x}_i\\ & ={\lambda }_1^k\left({\alpha }_1{x}_1+\sum _{i=2}^n{\alpha }_i{\left(\frac{{\lambda }_i}{{\lambda }_1}\right)}^k{x}_i\right)\\ & ={\lambda }_1^k({\alpha }_1{x}_1+{\epsilon }_k)\end{array}$$
其中, ${\epsilon }_k={\sum }_{i=2}^n{\alpha }_i{\left(\frac{{\lambda }_i}{{\lambda }_1}\right)}^k{x}_i$. 若${\alpha }_1\ne 0$, $x_1^{(i)}\ne 0$, 由$\underset{k\to \infty }{\lim }{\epsilon }_k=0$知,
$$?\begin{array}{l}\underset{k\to \infty }{\lim }\frac{v_k}{{\lambda }_1^k}={\alpha }_1{x}_1\\ \underset{k\to \infty }{\lim }\frac{v_{k+1}^{(i)}}{v_k^{(i)}}={\lambda }_1\end{array}$$
因此, 当$k$充分大的时候, $v_k$近似于主特征向量, $v_{k+1}$与$v_k$对应的非零分量的比值近似于主特征值.

在实际计算中, 我们还需要对计算结果进行规范化. 因为当$|{\lambda }_1|>1$时, $v_k$的非零分量趋于无穷. 从而造成了计算的不稳定, 会出现上溢或者下溢. 为此, 我们每次计算后都将结果进行规范化:
$$?\begin{array}{l}{v}_0=u_0\ne 0\\ v_k=A{u}_{k?1}\\ u_k=\frac{v_k}{\Vert v_k{\Vert }_{\infty }}\end{array}$$
此时, ${\lambda }_k=\Vert v_k{\Vert }_{\infty }$.

代码实现

根据上述公式, 我们可以编写如下程序:

#include <armadillo>
#include <list>
using std::list;
using namespace arma;
/*
 * 幂法
 * A :待求矩阵
 * x0:初始向量
 * X :迭代过程保存
 * e :求解精度
 */
double power_method(const mat &A, const vec &x0, list<vec> &X, const double &e)
{
    if (A.n_cols != A.n_rows)
        throw "非方阵不能求特征值！";
    if (A.n_cols != x0.n_elem)
        throw "初始向量与矩阵维数不匹配！";
    if (A.empty())
        return 0.;
    X.clear();
    vec x(A * x0);
    double lambda1(x.max());
    X.push_back(x);
    double lambda2((x = A * (x /= lambda1)).max());
    X.push_back(x);
    while ((lambda1 > lambda2 ? lambda1 - lambda2 : lambda2 - lambda1) > e)
    {
        lambda2 = (x = A * (x /= lambda1 = lambda2)).max();
        X.push_back(x);
    }
    return lambda2;
}

例 用幂法求矩阵
$$A=\left(\begin{array}{ccc}1& 1& 0.5\\ 1& 1& 0.25\\ 0.5& 0.25& 2\end{array}\right)$$
的主特征值和主特征向量.

首先编写main函数如下:

#include <stdio.h>
int main()
{
    mat A{{1, 1, 0.5}, {1, 1, 0.25}, {0.5, 0.25, 2}};
    vec x0{1., 1., 1.};
    list<vec> X;
    char str[6];
    printf("%.10f\n", power_method(A, x0, X, 1e-5));
    auto p(X.begin());
    int i(0);
    while (p != X.end())
    {
        sprintf(str, "x_%d=", ++i);
        p->print(str);
        ++p;
    }
    return 0;
}

最后得到结果如下:

2.5365374322
x_1=        
   2.5000   
   2.2500   
   2.7500
x_2=
   2.2273
   1.9773
   2.6591
x_3=
   2.0812
   1.8312
   2.6047
x_4=
   2.0021
   1.7521
   2.5753
x_5=
   1.9578
   1.7078
   2.5588
x_6=
   1.9325
   1.6825
   2.5494
x_7=
   1.9180
   1.6680
   2.5440
x_8=
   1.9096
   1.6596
   2.5409
x_9=
   1.9047
   1.6547
   2.5391
x_10=
   1.9019
   1.6519
   2.5380
x_11=
   1.9002
   1.6502
   2.5374
x_12=
   1.8992
   1.6492
   2.5370
x_13=
   1.8987
   1.6487
   2.5368
x_14=
   1.8983
   1.6483
   2.5367
x_15=
   1.8982
   1.6482
   2.5366
x_16=
   1.8980
   1.6480
   2.5366
x_17=
   1.8980
   1.6480
   2.5366
x_18=
   1.8979
   1.6479
   2.5365
x_19=
   1.8979
   1.6479
   2.5365

算法改进

由于$|\Vert v_k{\Vert }_{\infty }?{\lambda }_1|\approx c{\left|\frac{{\lambda }_2}{{\lambda }_1}\right|}^k$, 因此幂法是线性收敛的算法. 但当${\lambda }_2$很接近${\lambda }_1$时, 收敛很慢, 这时就需要一些加速收敛的方法:

Aitken外推法

记$m_k:=\Vert v_k{\Vert }_{\infty }$, 我们有如下外推加速:
$${\tilde{m}}_k=m_k?\frac{(m_k?m_{k?1}{)}^2}{m_k?2m_{k?1}+m_{k?2}},k?3$$

$${\tilde{u}}_k^{(j)}=u_k^{(j)}?\frac{(u_k^{(j)}?u_{k?1}^{(j)}{)}^2}{u_k^{(j)}?2u_{k?1}^{(j)}+u_{k?2}^{(j)}},u_k^{(j)}\ne 1$$

Rayleigh商加速

若$A\in {\mathbb{R}}^{n\times n}$为对称矩阵, 则幂法所得的规范化向量$u_k$的Rayleigh商加速给出特征值${\lambda }_1$的较好的近似值, 即:
$$\frac{?A{u}_k,u_k?}{?u_k,u_k?}={\lambda }_1+O\left({\left|\frac{{\lambda }_2}{{\lambda }_2}\right|}^{2k}\right)$$
其中, $??,??$表示${\mathbb{R}}^n$上的内积.

反幂法与原点位移

反幂法

反幂法, 顾名思义, 就是与幂法相反的算法. 它用来用来计算矩阵按模最小的特征值及其特征向量. 设$A\in {\mathbb{R}}^{n\times n}$为非奇异矩阵, 即它的特征值满足
$$|{\lambda }_1|?|{\lambda }_2|???|{\lambda }_{n?1}|>|{\lambda }_n|>0$$
则$A^{?1}$的特征值${\lambda }_1^{?1},{\lambda }_2^{?1},?,{\lambda }_n^{?1}$满足
$$|{\lambda }_n^{?1}|>|{\lambda }_{n?1}^{?1}|???|{\lambda }_1^{?1}|$$
因此, ${\lambda }_1^{?1}$为$A^{?1}$的主特征值. 因此, 我们只需要对$A^{?1}$应用幂法即可得到矩阵按模最小的特征值及其特征向量, 这种方法被称为反幂法.

在幂法的基础上, 我们可以编写反幂法的程序如下:

/*
 * 反幂法
 * A :待求矩阵
 * x0:初始向量
 * X :迭代过程保存
 * e :求解精度
 */
double reverse_power_method(const mat &A, const vec &x0, list<vec> &X, const double &e)
{
    mat A_inv(inv(A));
    return 1 / power_method(A_inv, x0, X, e);
}

原点位移

在反幂法的基础上, 如果我们对矩阵$A?pI$考虑反幂法¹, 我们可以得到$A?pI$按模最小的特征值及其特征向量. 也就是说, 我们可以求得矩阵$A$的最接近$p$的特征值及其特征向量. 我们同样可以实现原点位移法如下:

/*
 * 原点位移法
 * A :待求矩阵
 * x0:初始向量
 * X :迭代过程保存
 * e :求解精度
 * p :待求的中心点
 */
double p_power_method(const mat &A, const vec &x0, list<vec> &X, const double &e, const double &p)
{
    mat A_p(inv(A - p * eye(A.n_rows, A.n_cols)));
    return 1 / power_method(A_p, x0, X, e) + p;
}

实例分析

用原点位移法求矩阵
$$A=\left(\begin{array}{ccc}2& 1& 0\\ 1& 3& 1\\ 0& 1& 4\end{array}\right)$$
的特征值.

首先修改主函数如下并运行:

int main()
{
    mat A{{2., 1., 0.}, {1., 3., 1.}, {0., 1., 4.}};
    vec x0{1., 1., 1.};
    list<vec> X;
    char str[6];
    printf("%.10f\n", p_power_method(A, x0, X, 1e-5, 1.2679));
    auto p(X.begin());
    int i(0);
    while (p != X.end())
    {
        sprintf(str, "x_%d=", ++i);
        p->print(str);
        ++p;
    }
    return 0;
}

得到如下运行结果:

1.2679491924 
x_1=
  6.7764e+003
 -4.9600e+003
  1.8158e+003
x_2=
  2.0327e+004
 -1.4881e+004
  5.4467e+003
x_3=
  2.0328e+004
 -1.4881e+004
  5.4470e+003
x_4=
  2.0328e+004
 -1.4881e+004
  5.4470e+003
x_5=
  2.0328e+004
 -1.4881e+004
  5.4470e+003

---

1. 假设这是可行的. ??