数值计算day5-特征值与特征向量

上节课主要介绍了线性方程组的两种迭代求解算法，一个是Jacobi迭代（同步更新），一个是高斯塞德尔迭代（异步更新）。对于特殊的三对角系统，一种更简单快捷的Thomas算法也可以用来求解。之后介绍了向量范数与矩阵范数的概念，线性系统数值解的相对误差可以通过条件数来判定。本节课主要介绍矩阵的特征值，特征向量，以及其中涉及到的几种数值算法。

1. 特征值与特征向量

给定 $n \times n$ 维矩阵 $A$ ，满足下式的数 $\lambda$ 称作矩阵 $A$ 的一个特征值： $Au=\lambda u$ 而对应的向量 $u$ 称作对应特征值 $\lambda$ 的一个特征向量。上述运算可以推广到更一般的形式，而不仅仅是针对矩阵运算。假设 $u=f(x)$ 是一个关于 $x$ 的连续函数， $A=\frac{d}{dx}$ 表示微分运算，则 $\frac{d^2u}{dx^2}=k^2u$ 就可以表示为： $A^2u = k^2u$ 这是特征值问题的一种更广泛的表现形式。
特征值与特征向量在工程与科学中有着非常重要的作用。例如，在振动研究中，特征值表示系统或组件的固有频率，而特征向量表示这些振动的模式。
为了求解一个矩阵的特征值与特征向量，根据定义，可以得到： $(A-\lambda I)u=0$ 假设 $A-\lambda I$ 是非奇异矩阵，则上式只有零解，不满足求解要求，因此，想要求解特征值，需要 $det(A-\lambda I) = 0$ ，该式称作特征方程，特征方程的根即为矩阵 $A$ 的特征值。
例： $A=\begin{bmatrix}1 & 2\\3 & 2\end{bmatrix}$
其特征方程为 $det(A-\lambda I) = det\begin{bmatrix}1-\lambda & 2\\ 3 & 2-\lambda\end{bmatrix}=(1-\lambda)(2-\lambda)-6=0$ 求解可得特征值为 $\lambda_1=4,\lambda_2=-1$ 然而，对于阶数比较高的矩阵，特征值的求解不会这么简单，需要使用数值求解方法。

2. 幂法与反幂法

2.1 幂法

这节主要介绍如何使用幂法和反幂法分别求解一个矩阵所有特征值中模最大和模最小的值。给定矩阵 $A$ ，假设其有 $n$ 个实特征值 $|\lambda_1|>|\lambda_2|>\cdots>|\lambda_n|$ 其对应的特征向量为 $u_1,u_2,...,u_n$ 。幂法是通过迭代法来计算最大特征值 $\lambda_1$ ，首先随机选取初始向量 $x_1$ ， $x_1=c_1u_1+c_2u_2+...+c_nu_n$ 迭代计算 $x_2$ ， $Ax_1=c_1Au_1+c_2Au_2+...+c_nAu_n=\lambda_1c_1x_2$ $x_2=u_1+\frac{c_2}{c_1}\frac{\lambda_2}{\lambda_1}u_2+...+\frac{c_n}{c_1}\frac{\lambda_n}{\lambda_1}u_n$ 进一步计算 $x_3$ , $Ax_2=\lambda_1u_1+\frac{c_2}{c_1}\frac{\lambda^2_2}{\lambda_1}u_2+...+\frac{c_n}{c_1}\frac{\lambda^2_n}{\lambda_1}u_n=\lambda_1x_3$ $x_3=u_1+\frac{c_2}{c_1}\frac{\lambda^2_2}{\lambda^2_1}u_2+...+\frac{c_n}{c_1}\frac{\lambda^2_n}{\lambda^2_1}u_n$ 可以看到，其迭代公式为: $x_{k+1}=u_1+\frac{c_2}{c_1}\frac{\lambda^k_2}{\lambda^k_1}u_2+...+\frac{c_n}{c_1}\frac{\lambda^k_n}{\lambda^k_1}u_n$ 注意到 $\lambda_1$ 是最大的特征值，因此 $\frac{\lambda_i}{\lambda_i}<1$ ，当 $k$ 充分大时，有 $Ax_{k+1}=\lambda_1 u_1,x_{k+1}=u_1$ 。具体实现时，并没有 $\lambda_1$ 和 $u_1$ 的值，因此，迭代计算 $x_{k+1}=Ax_k$ 后，规范化 $x_{k+1}$ 即可（有待进一步验证，采用向量范数并不合理）。注意：最大特征值是指模最大的那个特征值
MATLAB实现：

A=[4 2 -2; -2 8 1 ; 2 4 -4]
x=[1 1 1]'
for i=1:16
x=A*x; x=x/norm(x,Inf);
end
A*x./x
%%%%%%%%%%%%%%%%
A =

     4     2    -2
    -2     8     1
     2     4    -4

ans =

    7.7502
    7.7503
    7.7502

>> eig(A)

ans =

   -3.6102
    3.8599
    7.7504
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function e = MaxEig(A)
    x = ones(size(A));
    for i=1:40
        x=A*x; 
        [mx,id] = max(abs(x));
        x=x/x(id);
    end
    e = A*x./x; 
    [mx,id] = max(abs(e));
    e = e(id);
end

2.2 反幂法

给定矩阵 $A$ ，假设其有 $n$ 个实特征值 $|\lambda_1|>|\lambda_2|>\cdots>|\lambda_n|$ 其对应的特征向量为 $u_1,u_2,...,u_n$ 。 $\lambda_n$ 是最小特征值，首先注意到如果 $Au=\lambda u$ ，则 $A^{-1}Au=A^{-1}\lambda u$ , 即 $u=A^{-1}\lambda u$ ，因此有 $A^{-1}u=\frac{1}{\lambda}u$ 。可以看到，当 $\lambda_n$ 为矩阵 $A$ 的最小特征值时， $\frac{1}{\lambda_n}$ 将是 $A^{-1}$ 的最大特征值，此时运用幂法求解 $A^{-1}$ 的最大特征值，取倒数，即为 $A$ 的最小特征值。反幂法中，需要注意的是，当最小特征值为0时，其倒数是没有定义的，反幂法求解的是第二小的特征值，且需要采用移位反幂法。
MATLAB实现

function e = MinEig(A)
    invA = inv(A);
    x = ones(size(A));
    for i=1:40
        x=invA*x; 
        [mx,id] = max(abs(x));
        x=x/x(id);
    end
    e = invA*x./x; 
    [mx,id] = max(abs(e));
    e = 1/e(id);
end
%%%%%%%%%%%%%%%%%%%%
>> A = [4 0 0;0 -1 0;0 0 -9]

A =

     4     0     0
     0    -1     0
     0     0    -9

>> MinEig(A)

ans =

    -1

>> eig(A)

ans =

    -9
    -1
     4
>> MaxEig(A)

ans =

    -9

3. QR分解法

幂法与反幂法是用来求解矩阵的最大特征值与最小特征值，想要求解矩阵所有特征值，可以使用QR分解法。假设 $A$ 是一个 $n\times n$ 的矩阵，有 $n$ 个互不相同的实特征值。QR分解的理论保证是，如果对矩阵 $A$ 做如下相似变换： $B=C^{-1}AC$ ，则特征值保持不变。这是因为如果 $Au=\lambda u$ ，则取 $v=C^{-1}u$ ，就有 $ACv = Au = \lambda Cv$ ，进一步有 $C^{-1}ACv=\lambda v$ ，因此 $\lambda$ 也是 $C^{-1}AC$ 的特征值。对 $A$ 进行QR分解的算法流程：

$A_1=A$ ，对 $A_1$ 做QR分解 $A_1=Q_1R_1$ 其中 $Q_1$ 是一个正交矩阵( $Q_1Q^T_1=I$ )， $R_1$ 是一个上三角矩阵
令 $A_2=R_1Q_1$ ，即 $A_2=Q^{-1}_1A_1Q_1$ ，特征值与 $A$ 是一致的。进一步对 $A_2$ 做QR分解： $A_2=Q_2R_2$
$A_3=R_2Q_2=Q^{-1}_2A_2Q_2=Q^{-1}_2Q^{-1}_1A_1Q_1Q_2$
重复上述过程，直到 $A_n$ 成为一个上三角矩阵，其特征值即为对角线上的元素。

MATLAB实现

A=[ 6 -7 2 ; 4 -5 2; 1 -1 1]
A0=A;
for i=1:40
    [Q R]=qr(A);
    A=R*Q;
end
A
ev=diag(A)
eig(A0)
%%%%%%%%%%%%%%%%
>> QRiteration

A =

     6    -7     2
     4    -5     2
     1    -1     1

A =

    2.0000   -8.6603   -7.4066
    0.0000   -1.0000   -1.0690
    0.0000   -0.0000    1.0000

ev =

    2.0000
   -1.0000
    1.0000

ans =

   -1.0000
    2.0000
    1.0000

在每次迭代中，对矩阵进行QR分解的操作是通过使用一种特殊的矩阵Householder矩阵来实现的，其形式为 $H=I-2\frac{vv^T}{v^Tv}$ 其中 $v=c+ \|c\|_2e$ ， $c$ 与 $e$ 为列向量， $\|c\|_2$ 为向量的2范数。 $H$ 矩阵是对称的，是正交的。 这就意味着 $HAH$ 与 $A$ 是相似的。下面详细说明，如何通过 $H$ 矩阵将 $A$ 分解为一个正交阵和一个上三角阵的乘积。

取 $c_1$ 为矩阵 $A$ 的第一列元素： $c_1=\begin{bmatrix}a_{11}\\a_{21}\\\vdots\\a_{n1}\end{bmatrix},e_1=\begin{bmatrix}\pm1\\0\\\vdots\\0\end{bmatrix}$ $e_1$ 的第一位元素的符号与 $c_1$ 的第一个元素的符号保持一致, $v_1=c_1+\|c_1\|_2e$
令 $H_1=I-2\frac{v_1v_1^T}{v_1^Tv_1}$ , $Q_1=H_1,R_1=H_1A$ ，则 $Q_1$ 是对称正交矩阵（HouseHolder, $A=Q_1R_1$ ）,注意 $v_1^T v_1=(c_1+\|c_1\|_2e_1)^T (c_1+\|c_1\|_2 e_1)=2\|c_1\|^2_2+2a_{11}\|c_1\|_2$ $R_1=H_1A=\begin{bmatrix}H_1c_1 & H_1 a_2 & \cdots & H_1a_n\end{bmatrix}$ $H_1c_1=c_1-2\frac{v_1v_1^Tc_1}{v_1^Tv_1}=c_1-\frac{(\|c_1\|_2^2+a_{11}\|c_1\|_2)v_1}{\|c_1\|^2_2+a_{11}\|c_1\|_2}=-\|c_1\|_2e_1$ 因此 $R_1=H_1A$ 的第一列除了第一个元素，其他均为 $0$ $R_1=\begin{bmatrix}a'_{11} & a'_{12} & \cdots & a'_{1n}\\ 0 & a'_{22} & \cdots & a'_{2n} \\ \vdots & \vdots & \ddots & \vdots\\0 & a'_{n2} & \cdots & a'_{nn}\end{bmatrix}$
取 $c_2$ 为矩阵第二列： $c_2=\begin{bmatrix}0\\a'_{21}\\\vdots\\a'_{n1}\end{bmatrix},e_2=\begin{bmatrix}0\\\pm1\\\vdots\\0\end{bmatrix}$ $e_2$ 的第二位元素的符号与 $c_2$ 的第二个元素的符号保持一致, $v_2=c_2+\|c_2\|_2e_2$
令 $H_2=I-2\frac{v_2v_2^T}{v_2^Tv_2}$ ， $Q_2=Q_1H_2=H_1H_2$ , $R_2=H_2R_1=H_2H_1A$ , $H_2$ 也为对称正交矩阵（ $A=H^T_1H^T_2R_2=Q_2R_2$ ），类似地， $R_2=H_2R_1$ 的第二列第二个元素一下均为 $0$ $R_2=\begin{bmatrix}a''_{11} & a''_{12} & \cdots & a''_{1n}\\ 0 & a''_{22} & \cdots & a''_{2n} \\ \vdots & \vdots & \ddots & \vdots\\0 & 0 & \cdots & a''_{nn}\end{bmatrix}$
重复上述过程，直到 $R_{n-1}$ 成为一个上三角矩阵，矩阵 $A$ 分解为正交矩阵 $Q_{n-1}$ 和上三角矩阵 $R_{n-1}$ 的乘积。

MATLAB实现：

function [Q R] = QRFactorization(R)
% The function factors a matrix [A] into an orthogonal matrix [Q]
% and an upper-triangular matrix [R].
% Input variables:
% A  The (square) matrix to be factored.
% Output variables:
% Q  Orthogonal matrix.
% R  Upper-triangular matrix.

nmatrix = size(R);
n = nmatrix(1);
I = eye(n);
Q = I;
for j = 1:n-1
    c = R(:,j);
    c(1:j-1) = 0;
    e(1:n,1)=0;
    if c(j) > 0
        e(j) = 1;
    else
        e(j) = -1;
    end
    clength = sqrt(c'*c);
    v = c + clength*e;
    H = I - 2/(v'*v)*v*v';
    Q = Q*H;
    R = H*R;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
>> A=[ 6 -7 2 ; 4 -5 2; 1 -1 1]

A =

     6    -7     2
     4    -5     2
     1    -1     1

>> [Q R] = qr(A)

Q =

   -0.8242    0.3925   -0.4082
   -0.5494   -0.7290    0.4082
   -0.1374    0.5608    0.8165

R =

   -7.2801    8.6537   -2.8846
         0    0.3365   -0.1122
         0         0    0.8165

>> [Q1 R1] = QRFactorization(A)

Q1 =

   -0.8242    0.3925   -0.4082
   -0.5494   -0.7290    0.4082
   -0.1374    0.5608    0.8165

R1 =

   -7.2801    8.6537   -2.8846
    0.0000    0.3365   -0.1122
   -0.0000    0.0000    0.8165

4. 总结

这节课主要介绍了矩阵特征值与特征向量的概念，对于低阶矩阵，可以通过特征方程求解出矩阵特征值；对于高阶矩阵，可以使用幂法与反幂法求解出矩阵的最大特征值与最小特征值，要求出矩阵的所有特征值，则可以对矩阵进行QR分解，将矩阵 $A$ 相似化为一个上三角矩阵，上三角矩阵的对角线元素即为矩阵 $A$ 的特征值。

人面猴
序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
沈念sama阅读 216,496评论 6赞 501
死咒
序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
沈念sama阅读 92,407评论 3赞 392
救了他两次的神仙让他今天三更去死
文/潘晓璐我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
开封第一讲书人阅读 162,632评论 0赞 353
道士缉凶录：失踪的卖姜人
文/不坏的土叔我叫张陵，是天一观的道长。经常有香客问我，道长，这世上最难降的妖魔是什么？我笑而不...
开封第一讲书人阅读 58,180评论 1赞 292
港岛之恋（遗憾婚礼）
正文为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
茶点故事阅读 67,198评论 6赞 388
恶毒庶女顶嫁案：这布局不是一般人想出来的
文/花漫我一把揭开白布。她就那样静静地躺着，像睡着了一般。火红的嫁衣衬着肌肤如雪。梳的纹丝不乱的头发上，一...
开封第一讲书人阅读 51,165评论 1赞 299
城市分裂传说
那天，我揣着相机与录音，去河边找鬼。笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。我是一名探鬼主播，决...
沈念sama阅读 40,052评论 3赞 418
双鸳鸯连环套：你想象不到人心有多黑
文/苍兰香墨我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
开封第一讲书人阅读 38,910评论 0赞 274
万荣杀人案实录
序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
沈念sama阅读 45,324评论 1赞 310
护林员之死
正文独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋以下内容为张勋视角年9月15日...
茶点故事阅读 37,542评论 2赞 332
白月光启示录
正文我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
茶点故事阅读 39,711评论 1赞 348
活死人
序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
沈念sama阅读 35,424评论 5赞 343
日本核电站爆炸内幕
正文年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
茶点故事阅读 41,017评论 3赞 326
男人毒药：我在死后第九天来索命
文/蒙蒙一、第九天我趴在偏房一处隐蔽的房顶上张望。院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
开封第一讲书人阅读 31,668评论 0赞 22
一桩弑父案，背后竟有这般阴谋
文/苍兰香墨我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。一阵脚步声响...
开封第一讲书人阅读 32,823评论 1赞 269
情欲美人皮
我被黑心中介骗来泰国打工，没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。一个月前我还...
沈念sama阅读 47,722评论 2赞 368
代替公主和亲
正文我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
茶点故事阅读 44,611评论 2赞 353