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

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

1. 特征值与特征向量

给定(n imes n)维矩阵(A)，满足下式的数(lambda)称作矩阵(A)的一个特征值：$$Au=lambda u$$ 而对应的向量(u)称作对应特征值(lambda)的一个特征向量。上述运算可以推广到更一般的形式，而不仅仅是针对矩阵运算。假设(u=f(x))是一个关于(x)的连续函数，(A=frac{d}{dx})表示微分运算，则$$frac{d^2u}{dx2}=k^2u$$ 就可以表示为：$$A^2u = k^2u$$ 这是特征值问题的一种更广泛的表现形式。
特征值与特征向量在工程与科学中有着非常重要的作用。例如，在振动研究中，特征值表示系统或组件的固有频率，而特征向量表示这些振动的模式。
为了求解一个矩阵的特征值与特征向量，根据定义，可以得到：$$(A-lambda I)u=0$$ 假设(A-lambda I)是非奇异矩阵，则上式只有零解，不满足求解要求，因此，想要求解特征值，需要(det(A-lambda I) = 0)， 该式称作特征方程，特征方程的根即为矩阵(A)的特征值。
例： $$A=egin{bmatrix}1 & 23 & 2end{bmatrix}$$
其特征方程为$$det(A-lambda I) = detegin{bmatrix}1-lambda & 2 3 & 2-lambdaend{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}{lambda2_1}u_2+...+frac{c_n}{c_1}frac{lambda^2_n}{lambda2_1}u_n$$ 可以看到，其迭代公式为: $$x_{k+1}=u_1+frac{c_2}{c_1}frac{lambda^k_2}{lambdak_1}u_2+...+frac{c_n}{c_1}frac{lambda^k_n}{lambdak_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 imes 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-2frac{vv^T}{vTv}$$ 其中(v=c+ |c|_2e)，(c)与(e)为列向量，(|c|_2)为向量的2范数。(H)矩阵是对称的，是正交的。 这就意味着(HAH)与(A)是相似的。下面详细说明，如何通过(H)矩阵将(A)分解为一个正交阵和一个上三角阵的乘积。

• 取(c_1)为矩阵(A)的第一列元素：$$c_1=egin{bmatrix}a_{11}a_{21}\vdotsa_{n1}end{bmatrix},e_1=egin{bmatrix}pm1\vdotsend{bmatrix}$$ (e_1)的第一位元素的符号与(c_1)的第一个元素的符号保持一致, (v_1=c_1+|c_1|_2e)
• 令 (H_1=I-2frac{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=egin{bmatrix}H_1c_1 & H_1 a_2 & cdots & H_1a_nend{bmatrix}$$ $$H_1c_1=c_1-2frac{v_1v_1^Tc_1}{v_1Tv_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=egin{bmatrix}a'{11} & a'{12} & cdots & a'{1n} 0 & a'{22} & cdots & a'{2n} vdots & vdots & ddots & vdots & a'{n2} & cdots & a'_{nn}end{bmatrix}$$
• 取(c_2)为矩阵第二列：$$c_2=egin{bmatrix}0a'{21}\vdotsa'{n1}end{bmatrix},e_2=egin{bmatrix}0\pm1\vdotsend{bmatrix}$$ (e_2)的第二位元素的符号与(c_2)的第二个元素的符号保持一致, (v_2=c_2+|c_2|_2e_2)
• 令(H_2=I-2frac{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=egin{bmatrix}a''{11} & a''{12} & cdots & a''{1n} 0 & a''{22} & cdots & a''{2n} vdots & vdots & ddots & vdots & 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)的特征值。