幂法在线视频

各位同学

大家好

我是北京理工大学朱国庆老师

接下来我们给大家讲解

关于矩阵特征值和特征向量的计算问题

首先 回顾一下

什么叫做矩阵的特征值与特征向量

我们假设矩阵A是一个n阶的方阵

如果存在这么一个数λ

和n维的非零列向量x

使得Ax=λx这个关系式成立

我们就说

这个数λ是这个矩阵A的特征值

然后这个非零向量x称为矩阵A

对应于特征值λ的特征向量

前面我们给大家讲

矩阵的计算问题的时候会说一句

这有一个叫做谱半径的概念

其实本质上就是矩阵的特征值

按模最大的那一个

包括矩阵的范数 二范数

也是需要计算这个特征值

还有那个二条件数也需要计算特征值

然后在工程当中有很多的振动问题

比如桥梁或者建筑物的振动

机械零件和飞机机翼的振动

以及一些稳定性分析和相关分析

在数学当中都可以转换成

这个矩阵的特征值和特征向量的计算的问题

这个线性代数当中告诉你了

就是说如果矩阵A是一个2阶的

包括3阶的方阵的话

那么按照行列式展开的方法求特征多项式

就是λE-A的行列式

求它的根的形式

我们可以把矩阵A的特征值给它计算出来

但是这里面牵扯一个问题就是这个n

如果足够大的时候

就是按照特征展开多项式的方法的话

这个特征多项式的根

它的这个表达式计算起来比较麻烦一点

然后求它的根

工作量就更加的大了

用这种方法求解矩阵特征值其实是不切实际的

因此我们需要研究

当这个矩阵A的阶数比较大的时候

这个矩阵A的特征值与特征向量的数值解法的问题

我们知道

现在的计算机得到了广泛的应用

那么在计算机当中

适用于求解高阶矩阵的特征值的问题

我们讲一个要叫做幂法

另外一个是反幂法

首先 来说一下什么叫做幂法

这解幂法的时候

我们首先假设矩阵A是n阶的实矩阵

也就是矩阵A的每一个元素都是实数

然后λi ui是矩阵A的n个特征值

这里面λi可能是有重根

也就是λ1跟λ2有可能相等

或者是某一个λ跟其它的特征值是一样的

ui是相应于λi的特征向量

我们假设这个所有的λi当中满足一个条件

也就是λ1的模要大于λ2的模

然后λ2的模要大于等于λ3的模

大于等于λn的模

这里面λ1 λ2 λ3

我们是按照按模递减的顺序做了一个排列

这个ui分别对应λi的特征向量

这里面我们假设u1 u2 un

它们是线性无关的

我们看一下计算机当中

适用的幂法是怎么推导出来的

首先任取一个非零的向量

这里面是初始向量V0

我们矩阵A构造一个向量序列Vk

V1是拿矩阵A乘上V0

这是一个V1

然后有V1了之后

我们是A乘上V1

构造一个量

然后实际上这个V1是AV0

所以是A的平方V0

那么Vk+1

实际上是拿矩阵A乘上Vk

也就是A的k+1次幂

然后乘上V0

这样构造一个迭代的向量的序列

根据前面的这个假设

我们说初始向量V0

其实可以写成

特征向量u1 u2 un的一个线性组合

我们可以写成这样的形式

把V0写成a1u1

加上a2u2加上anun的形式

这里面我们要假设A1是不等于零的

按照前面的假设

V0写成了

a1u1 a2u2加上anun的形式

那么Vk它是A乘上Vk

那么Vk它是A乘上Vk
也就是A的k次幂乘上V0的形式

也就是A的k次幂乘上V0的形式

V0写成了上面的形式之后

那么我们把那个AV0的形式

直接写成了a1λ1的k次幂u1

加上a2λ2的k次幂u2

再加上anλn的k次幂un的这个形式

这里面我们根据这个假设

λ1的k次幂 λ2的k次幂

λn的k次幂有一个大小关系

我们先把那个λ1的k次幂提出来

然后中括号内的这部分第一项a1u1

后面那里面每一项都是ai

然后λi

除上λ1的k次幂

然后是ui

这里边我们做一个简单的记号

写成λ1的k次幂a1u1

加上εk

这个εk其实是一个附加项

这里根据前面的假设

λi除上λ1的时候

它是严格小于1的

这里面那个i是从2一直到n的

然后如果这个k你确定有无穷大的时候

我们这个附加项εk

它的极限就等于零的

所以如果是这个k足够大的话

我们假设Vk除上λ1的k次幂

k趋近无穷大的时候

这里面其实ε趋近于零

就剩下一个a1u1了

当k充分大的时候

我们说这个Vk是λ1的k次幂

这个a1u1

我们就说它其实就是

矩阵对应于特征值λ1的一个近似特征向量

根据前面的讨论

Vk+1是等于AVk的

所以这个λ1的k+1次幂a1u1

就等于λ1倍的Vk

我们把这个V用Vk第i个下标

下面加了个下标 i

表示Vk这个向量的第i个分量

那么k充分大的时候

Vk+1的第二个分量

除上Vk第二个分量

它就近似的等于λ1了

这里面的λ1就是我们所要求的

矩阵A的主特征值λ1的近似值了

这里边总结一下

前面的第四个式子

其实Vk实际上是等于λ1的k次幂a1u1的

那么第六个式子λ1相当于是

近似地通过Vk+1的迭代向量

Vk+1的第二个分量

除上第k个迭代向量

Vk的第二个分量得到的

这种通过由已知向量V0以及矩阵A的乘幂

A的k次幂构造出来的迭代序列Vk

来计算矩阵A的特征值λ1

最大特征值λ1

以及相应的特征向量u1的这个过程

我们一般都称之为幂法

书上把前面的这个推导过程做了一个总结

这就是下面这个定理

假设矩阵A是n乘n阶的矩阵

特征值λi满足

λ1大于λ2大于λ3大于等于λn

这里面都是按模递减的顺序排列

而且与λi对应的特征向量u1 u2 un

它是线性无关的

那么我们对任意的非零的数值向量V0而言

这里面Vk构造出来这个序列A的k次幂V0

它是趋近于u1的

也就是我们要求的那个对应λ1特征向量

然后Vk+1的第二个分量

除上Vk的第二个分量

它就是收敛到了λ1的

这个里面需要给大家注意的就是

Vk+1形式上写成了λ1的k次幂 a1u1

实际上计算的时候会遇到这样一个问题

如果λ1我们大家知道

λ1如果是大于1的一个数

那么k+1次幂

这个k如果是趋近无穷的话

那么λ1的k+1次幂

这将是一个非常大的一个数

也就是有可能在计算机讨论当中

它是一个溢出的一个数

我们没办法对它进行实际的运算

那么我们讨论的时候

实际上是要对Vk+1做一个标准化的处理

也叫做归一化的处理

具体是这样的

假设与λ1对应的特征向量是u1

我们说如果一个向量x写成x1 x2 xn的形式

那么我们把xi

就是x的分量当中最大的那个拿出来

按模最大的那个拿出来记成xr

然后xr是取成了这个向量当中最大的那个分量

我们取初始向x0

把x(0)标准化为y(0)

实际上是拿x(0)

除上max(x0)

也就是最大的那个分量

使得标准化之后那个向量

我们记为y(0)

这个y(0)其实就是

我们要计算的x(0)这个表达式当中

按这个分量最大

归一化之后得到一个结果

这个里面

这个y(k)记为u1

除上这个max(u1)

然后max(xk)其实最后就收敛到

这个λ1的情况

MatLab在计算的时候

它的语言实际上是这么来算的

x是一个n维的列向量

x是一个n维的列向量

m r中括号

[m,r]

然后等于max

然后是(abs)(x)

这里面x的第r个分量xr

是需要取最大

y就等于x除上xr

我们具体来看一下我们的幂法的算法

首先需要输入矩阵A

然后一个初始向量x

然后这里面有一个误差限ε

还有一个最大的迭代次数N

这里面最大迭代次数N

这里面最大迭代次数N

实际上是防止我们

进行迭代过程当中无限循环的

第二步 我们取k=1 λ=0

然后y作为x的一个归化

或者叫做标准化的一个处理

第三步 计算x等于Ay

然后β是等于max(x)

然后y就等于x除上β

还是y作为x的一个标准化处理

这里面如果λ减掉β是小于ε的

我们就得到了

这个β就是我们的特征值

y就是我们的特征向量

这个时候我们就停止计算

否则的话 我们转入下一步

我们这个k如果还小于迭代次数的话

我们把k+1赋给新的k

然后把那个β赋给新的λ

然后重新进行第三步的一个循环运算

否则的话 我们就输出所谓的失败信息

这个算法就是这样设计的

这个里面我给大家做一个简单的推广

前面的假设当中

λ1 λ2 λn

它们之间有一个按模降幂排列

最大特征值只有一个

其实这里面最大特征值

按模最大特征值

如果有m个的话

这个算法也是可行的

首先我们来看一下Vk

如果我们写成a的k次幂V0的形式

这里面V0写成了a1u1 a2u2

amum anun形式的话

我们把前面m个a1u1一直加到amum的话

我们把它们整理合在一起

去乘λ1的k次幂

然后是λ1的k次幂a1u1

加上a2u2 amum的形式

然后其它的

后面λm+1到λn的时候

λm+1的k次幂除上λ1的k次幂

λn的k次幂除上λ1的k次幂

当这个k足够大的时候

就是这个所谓的附加项

它是要趋近于零的

第二种情况

如果我们讨论的这个按模最大的特征值

是一个正负的

也就是说λ1等于-λ2的时候

当然这个按模的话

这两个还是一样的

那么我们可以取Vk序列的偶数列

也可以取Vk序列的奇数列

也就是说对于Vk这个序列的偶数列而言

偶数列当中的后位项的第二个分量

除上偶数列的前一个分量当中的那个偶数项

得到的就是λ1的平方

这个也是可以进行计算的

第三个

如果我们的V0是任意的

但是这个时候A0取成零了

那么x0的写法当中

前面的a1u1那部分可能是会出现问题的

按照舍入误差的讨论

其实这里面这个迭代了k次之后

a1实际上是可以不等于零的

照样也可以求出来λk和uk的运算

这个地方 舍入误差其实是在

正向地帮助我们迭代运算的

下面看一个具体的例子

用幂法来计算矩阵这个A

它的按模最大的特征值

以及对应的特征向量

首先我们取这个V0的时候

取的是0 0 1

这个也是为了我们计算的简化

我们把这个V0直接写成了0 0 1

那么它的归化也好

标准化也好

都是u0的形式

它们两个是一样的

我们把Vk是等于Auk-1的

然后μk实际上取的是

Vk当中的那个最大元素

Vk的标准化就是拿Vk除上μk

就是我们的那个uk了

通过计算V1是等于2 4 1

当然是个列向量

然后μ1取的是

这个2 4 1当中最大的那个

就是那个4

然后u1当然就是μ1分之V1了

这里面2除上4是0.5

4除上4是1

然后1除上4是0.25

这个V1就是0.5 1 0.25

我们把这个所有的这个运算程序

当然是通过计算机来算的

我们把这个第k次的运算的结果

Vk的结果 μk的结果

以及这个uk的结果给大家罗列下来

其实通过这个可以看得出来

随着k的次数的增加

然后这个μk实际上是一个趋近于 10.9999

当然你就可以理解为接近11就行了

然后那个Vk的结果也是逐渐的趋于固定

不再变动

所以我们按照这个运算的话

我们可以知道

这个λ1其实就约等于是11的

10.9999当然就近似等于11了

然后x1的取值

这里面当k取7的时候

和k取8的时候都是0.5 1.0 0.75

我们就可以近似地把它理解成

特征向量就是0.5 1.0和0.75了

关于幂法的讨论其实还有一个问题

就是如果是这个初始向量选的不好的话

它可能会牵扯到这个运算速度过慢

我们有一个幂法的加速的讨论

这里面介绍两个

第一个就是原点移位法

我们假设λi是A的特征值

那么λi-λ0

其实那是A-λ0I的特征值

这个在高中代数当中

其实是通过简单的运算是可以得到的

我们把Vk+1写成（A-λ0I）Vk这种形式

这里面A-λ0I

k+1次幂V0

它可以写成λ1-λ0的k+1次幂

后面中括号这些部分

第一项a1u1其实保留了

后面的λ2-λ0

除上λ1-λ0的k+1次幂

我们通过原点移位

其实是希望λ2-λ0除上λ1-λ0

要比λ2除上λ1要更小才行

这样的话我们就可以保证这个幂法的加速了

一般来说是这样的

假设矩阵A有特征值λi

这里面满足λ1大于λ2

跟前面的假设是一样的

大于等于λn

取一个λ0

使得λ1-λ0要大于λi-λ0

而且后面λi-λ0

除上λ1-λ0这一部分

它的最大值

其实要小于λ2除上λ1

这样才能加速

用幂法求A-λ0I的按模最大特征值λ1*

那么这个幂法计算的时候

要比直接算可能速度要更快一点

那么我们这个λ1*算出来之后

真正的λ1就等于λ1*

加上λ0就可以了

这种方法我们都称之为叫原点移位法

本质上而言就是把我们通常的那个零点

移动到了λ0这个点

其实在实际应用的时候

那个矩阵A的特征值

可能我们是不知道的

也就说我们算λ的时候

可能没办法来直接确定

到底怎么去找它的一个初值

这个时候如果我们给大家建议

就是当收敛速度过慢的时候

可以适当的移动原点

通过这个式算

我们找到一个合适的运算的初值

这里面还有幂法的加速当中的第二个算法

就是所谓的埃特肯加速

它的原理是这样的

如果有一个数列ak它收敛到a

而且ak+1减掉a

除上ak-a

k趋近无穷的时候

它这个极限如果等于c

这个c是不等于零的

也就是我们通常说的一个序列

它收敛的时候

它是一个线性收敛的关系

那么当我们取极限的运算

当k趋近无穷大

我没办法取到无穷大

但是我可以假设它充分大的时候

这个表达式ak+1减a

除上ak-a这个表达式

我可以把它的后一项

也就是ak+2减a

除上ak加1减a

除上ak+1减a

这两个式子

我让它近似相等

记住 这里面约等号其实是一个形式的写法

我们实际上计算的时候

我们可以按照直等

也就是相等的情况来计算

这里面这个约等号当中所含的那个a

把a算出来那个结果作为a的一个近似值

我们一般是记作ak上面带了一个帽子

这样的一个形式的写法

可以证明的就是说

戴了帽子之后的ak这个序列

它要收敛的

要比那个a

原来那个ak收敛趋近a的时候

收敛的要快一点

这种用ak上面带帽的这种形式

收逼近a的这种方法

就叫做埃特肯的加速算法

具体的算法设计是这样的

首先输入矩阵A的元素

然后给定一个初始向量

给定这个误差限

给定最大的迭代次数

我们取k等于1 α0取成0 α1取成0

λ0取成1.0

λ0取成1.0

y还是作为x的一个标准化的处理

这里面x等于Ay的

然后把max

也就是x这个分量当中最大的一个记成α2

然后把x再做一个带入化处理

记成新的y

那么讨论的时候

我们要算那个λ

实际上等于α0减掉α1减掉α0的平方

除上α2减掉2倍α1

加上α0的这种形式

这里面这个λ0

这个一开始的时候等于1

如果我们算出来这个λ减掉λ0是小于这个ε

那么实际上就说

我们给定的算出来这个λ

就是我们要的那个真实的计算的λ的值了

这个时候那个y就是我们要的特征向量

这个运算就停止了

否则的话

我们把α1记成α0

把α2赋给α1

然后

把λ赋给λ0

把k+1再赋给k

也就是那个运算的迭代步数再增加1

接着再从第三步接着开始算起来

这样的话

我们可以一步一步地通过迭代

然后得到满足λ减λ0小于ε

的这样的一个特征值的一个运算的过程

好 这个加速度算法就讲到这个地方

这一部分的分享就到此结束