文章目录
- 幂法
- 算法描述
- 代码实现
- 算法改进
- Aitken外推法
- Rayleigh商加速
- 反幂法与原点位移
- 反幂法
- 原点位移
- 实例分析
幂法
设矩阵
A
∈
R
n
×
n
A\in\mathbb R^{n\times n}
A∈Rn×n满足
∣
λ
1
∣
>
∣
λ
2
∣
⩾
⋯
⩾
∣
λ
n
∣
|\lambda_1|>|\lambda_2|\geqslant\cdots\geqslant|\lambda_n|
∣λ1∣>∣λ2∣⩾⋯⩾∣λn∣
对应的
n
n
n个特征向量
x
1
,
x
2
,
⋯
,
x
n
x_1,x_2,\cdots,x_n
x1,x2,⋯,xn线性无关. 则称模最大的特征值
λ
1
\lambda_1
λ1为矩阵
A
A
A的主特征值, 称对应的特征向量
x
1
x_1
x1为主特征向量.
幂法就是用于求主特征值和主特征向量的一种数值方法.
算法描述
首先任取一个非零的初始向量
v
0
v_0
v0, 由矩阵
A
A
A构造一个向量序列如下:
v
k
=
A
v
k
−
1
=
⋯
=
A
k
v
0
,
k
=
1
,
2
,
⋯
v_k=Av_{k-1}=\cdots=A^kv_0\quad,k=1,2,\cdots
vk=Avk−1=⋯=Akv0,k=1,2,⋯
由假设,
v
0
v_0
v0可以表示为
v
0
=
α
1
x
1
+
α
2
x
2
+
⋯
+
α
n
x
n
v_0=\alpha_1x_1+\alpha_2x_2+\cdots+\alpha_nx_n
v0=α1x1+α2x2+⋯+αnxn
若记
v
k
(
i
)
v_k^{(i)}
vk(i)为
v
k
v_k
vk的第
i
i
i个分量, 则有
v
k
=
A
k
v
0
=
∑
i
=
1
n
α
i
λ
i
k
x
i
=
λ
1
k
(
α
1
x
1
+
∑
i
=
2
n
α
i
(
λ
i
λ
1
)
k
x
i
)
=
λ
1
k
(
α
1
x
1
+
ε
k
)
\begin{aligned} v_k&=A^kv_0\\ &=\sum_{i=1}^n\alpha_i\lambda_i^kx_i\\ &=\lambda_1^k\left(\alpha_1x_1+\sum_{i=2}^n\alpha_i\left(\frac{\lambda_i}{\lambda_1}\right)^kx_i\right)\\ &=\lambda_1^k(\alpha_1x_1+\varepsilon_k) \end{aligned}
vk=Akv0=i=1∑nαiλikxi=λ1k(α1x1+i=2∑nαi(λ1λi)kxi)=λ1k(α1x1+εk)
其中,
ε
k
=
∑
i
=
2
n
α
i
(
λ
i
λ
1
)
k
x
i
\varepsilon_k=\sum\limits_{i=2}^n\alpha_i\left(\dfrac{\lambda_i}{\lambda_1}\right)^kx_i
εk=i=2∑nαi(λ1λi)kxi. 若
α
1
≠
0
\alpha_1\ne0
α1=0,
x
1
(
i
)
≠
0
x_1^{(i)}\ne0
x1(i)=0, 由
lim
k
→
∞
ε
k
=
0
\lim\limits_{k\to\infin}\varepsilon_k=0
k→∞limεk=0知,
{
lim
k
→
∞
v
k
λ
1
k
=
α
1
x
1
lim
k
→
∞
v
k
+
1
(
i
)
v
k
(
i
)
=
λ
1
\begin{cases} \lim\limits_{k\to\infin}\dfrac{v_k}{\lambda_1^k}=\alpha_1x_1\\ \lim\limits_{k\to\infin}\dfrac{v_{k+1}^{(i)}}{v_k^{(i)}}=\lambda_1 \end{cases}
⎩⎪⎪⎨⎪⎪⎧k→∞limλ1kvk=α1x1k→∞limvk(i)vk+1(i)=λ1
因此, 当
k
k
k充分大的时候,
v
k
v_k
vk近似于主特征向量,
v
k
+
1
v_{k+1}
vk+1与
v
k
v_k
vk对应的非零分量的比值近似于主特征值.
在实际计算中, 我们还需要对计算结果进行规范化. 因为当
∣
λ
1
∣
>
1
|\lambda_1|>1
∣λ1∣>1时,
v
k
v_k
vk的非零分量趋于无穷. 从而造成了计算的不稳定, 会出现上溢或者下溢. 为此, 我们每次计算后都将结果进行规范化:
{
v
0
=
u
0
≠
0
v
k
=
A
u
k
−
1
u
k
=
v
k
∥
v
k
∥
∞
\begin{cases} v_0=u_0\ne0\\ v_k=Au_{k-1}\\ u_k=\dfrac{v_k}{\|v_k\|_\infin} \end{cases}
⎩⎪⎪⎨⎪⎪⎧v0=u0=0vk=Auk−1uk=∥vk∥∞vk
此时,
λ
k
=
∥
v
k
∥
∞
\lambda_k=\|v_k\|_\infin
λk=∥vk∥∞.
代码实现
根据上述公式, 我们可以编写如下程序:
#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
=
(
1
1
0.5
1
1
0.25
0.5
0.25
2
)
A=\begin{pmatrix} 1&1&0.5\\ 1&1&0.25\\ 0.5&0.25&2 \end{pmatrix}
A=⎝⎛110.5110.250.50.252⎠⎞
的主特征值和主特征向量.
首先编写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
算法改进
由于 ∣ ∥ v k ∥ ∞ − λ 1 ∣ ≈ c ∣ λ 2 λ 1 ∣ k \big|\|v_k\|_\infin-\lambda_1\big|\approx c\left|\dfrac{\lambda_2}{\lambda_1}\right|^k ∣∣∥vk∥∞−λ1∣∣≈c∣∣∣∣λ1λ2∣∣∣∣k, 因此幂法是线性收敛的算法. 但当 λ 2 \lambda_2 λ2很接近 λ 1 \lambda_1 λ1时, 收敛很慢, 这时就需要一些加速收敛的方法:
Aitken外推法
记
m
k
:
=
∥
v
k
∥
∞
m_k:=\|v_k\|_\infin
mk:=∥vk∥∞, 我们有如下外推加速:
m
~
k
=
m
k
−
(
m
k
−
m
k
−
1
)
2
m
k
−
2
m
k
−
1
+
m
k
−
2
,
k
⩾
3
\tilde m_k=m_k-\frac{(m_k-m_{k-1})^2}{m_k-2m_{k-1}+m_{k-2}}\quad,k\geqslant3
m~k=mk−mk−2mk−1+mk−2(mk−mk−1)2,k⩾3
u ~ k ( j ) = u k ( j ) − ( u k ( j ) − u k − 1 ( j ) ) 2 u k ( j ) − 2 u k − 1 ( j ) + u k − 2 ( j ) , u k ( j ) ≠ 1 \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)}}\quad,u_k^{(j)}\ne1 u~k(j)=uk(j)−uk(j)−2uk−1(j)+uk−2(j)(uk(j)−uk−1(j))2,uk(j)=1
Rayleigh商加速
若
A
∈
R
n
×
n
A\in\mathbb R^{n\times n}
A∈Rn×n为对称矩阵, 则幂法所得的规范化向量
u
k
u_k
uk的Rayleigh商加速给出特征值
λ
1
\lambda_1
λ1的较好的近似值, 即:
⟨
A
u
k
,
u
k
⟩
⟨
u
k
,
u
k
⟩
=
λ
1
+
O
(
∣
λ
2
λ
2
∣
2
k
)
\frac{\lang Au_k,u_k\rang}{\lang u_k,u_k\rang}=\lambda_1+O\left(\left|\frac{\lambda_2}{\lambda_2}\right|^{2k}\right)
⟨uk,uk⟩⟨Auk,uk⟩=λ1+O(∣∣∣∣λ2λ2∣∣∣∣2k)
其中,
⟨
⋅
,
⋅
⟩
\lang\cdot,\cdot\rang
⟨⋅,⋅⟩表示
R
n
\mathbb R^n
Rn上的内积.
反幂法与原点位移
反幂法
反幂法, 顾名思义, 就是与幂法相反的算法. 它用来用来计算矩阵按模最小的特征值及其特征向量. 设
A
∈
R
n
×
n
A\in\mathbb R^{n\times n}
A∈Rn×n为非奇异矩阵, 即它的特征值满足
∣
λ
1
∣
⩾
∣
λ
2
∣
⩾
⋯
⩾
∣
λ
n
−
1
∣
>
∣
λ
n
∣
>
0
|\lambda_1|\geqslant|\lambda_2|\geqslant\cdots\geqslant|\lambda_{n-1}|>|\lambda_n|>0
∣λ1∣⩾∣λ2∣⩾⋯⩾∣λn−1∣>∣λn∣>0
则
A
−
1
A^{-1}
A−1的特征值
λ
1
−
1
,
λ
2
−
1
,
⋯
,
λ
n
−
1
\lambda_1^{-1},\lambda_2^{-1},\cdots,\lambda_n^{-1}
λ1−1,λ2−1,⋯,λn−1满足
∣
λ
n
−
1
∣
>
∣
λ
n
−
1
−
1
∣
⩾
⋯
⩾
∣
λ
1
−
1
∣
|\lambda_n^{-1}|>|\lambda_{n-1}^{-1}|\geqslant\cdots\geqslant|\lambda_1^{-1}|
∣λn−1∣>∣λn−1−1∣⩾⋯⩾∣λ1−1∣
因此,
λ
1
−
1
\lambda_1^{-1}
λ1−1为
A
−
1
A^{-1}
A−1的主特征值. 因此, 我们只需要对
A
−
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 − p I A-pI A−pI考虑反幂法1, 我们可以得到 A − p I A-pI A−pI按模最小的特征值及其特征向量. 也就是说, 我们可以求得矩阵 A A A的最接近 p p 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
=
(
2
1
0
1
3
1
0
1
4
)
A=\begin{pmatrix} 2&1&0\\1&3&1\\0&1&4 \end{pmatrix}
A=⎝⎛210131014⎠⎞
的特征值.
首先修改主函数如下并运行:
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
假设这是可行的. ↩︎