线性代数之矩阵特征值与特征向量的数值求解方法

文章目录

前言1. 幂迭代法（Power Iteration）幂法与反幂法求解矩阵特征值幂法求最大特征值编程实现补充说明

2. 逆幂迭代法（Inverse Iteration）移位反幂法

3. QR 算法（QR Algorithm）——稠密矩阵理论推导编程实现

4. 雅可比方法（Jacobi Method）——对称矩阵编程实现

5. Lanczos 算法（稀疏矩阵）6. 分治法7. 方法选择指南8. 关键公式与说明参考文献

前言

定n×n维矩阵A，满足下式的数λ称作矩阵A的一个特征值：

A

u

=

λ

u

Au = \lambda u

Au=λu 推广形式

特征值问题可推广到更一般的形式。假设

u

=

f

(

x

)

u = f(x)

u=f(x) 是一个连续函数，

A

=

d

d

x

A = \frac{d}{dx}

A=dxd​ 表示微分运算，则二阶微分方程：

d

2

u

d

x

2

=

k

2

u

\frac{d^2u}{dx^2} = k^2u

dx2d2u​=k2u 可表示为：

A

2

u

=

k

2

u

A^2u = k^2u

A2u=k2u 这是特征值问题在微分算子中的表现形式。

特征方程与求解方法

根据定义

(

A

−

λ

I

)

u

=

0

(A - \lambda I)\mathbf{u} = 0

(A−λI)u=0 若

A

−

λ

I

A - \lambda I

A−λI 非奇异，则方程只有零解。因此，特征值需满足：

det

⁡

(

A

−

λ

I

)

=

0

\det(A - \lambda I) = 0

det(A−λI)=0 此方程称为特征方程，其根即为矩阵

A

A

A 的特征值。

示例

给定矩阵：

A

=

[

1

2

3

2

]

A = \begin{bmatrix} 1 & 2 \\ 3 & 2 \end{bmatrix}

A=[13​22​] 其特征方程为：

det

⁡

[

1

−

λ

2

3

2

−

λ

]

=

(

1

−

λ

)

(

2

−

λ

)

−

6

=

0

\det\begin{bmatrix} 1 - \lambda & 2 \\ 3 & 2 - \lambda \end{bmatrix} = (1 - \lambda)(2 - \lambda) - 6 = 0

det[1−λ3​22−λ​]=(1−λ)(2−λ)−6=0 展开并化简：

λ

2

−

3

λ

−

4

=

0

⟹

λ

1

=

4

,

λ

2

=

−

1

\lambda^2 - 3\lambda - 4 = 0 \implies \lambda_1 = 4,\ \lambda_2 = -1

λ2−3λ−4=0⟹λ1​=4, λ2​=−1

然而，对于高阶矩阵，特征值的解析解通常难以直接计算，需借助数值方法（如QR算法、幂迭代法等）进行求解。

1. 幂迭代法（Power Iteration）

目标：求解矩阵的模最大特征值及其对应特征向量。

幂法与反幂法求解矩阵特征值

本节介绍如何使用幂法和反幂法分别求解矩阵的模最大和模最小特征值。给定矩阵

A

A

A，假设其有

n

n

n 个实特征值：

∣

λ

1

∣

>

∣

λ

2

∣

>

⋯

>

∣

λ

n

∣

|λ_1| > |λ_2| > \cdots > |λ_n|

∣λ1​∣>∣λ2​∣>⋯>∣λn​∣ 对应的特征向量为

u

1

,

u

2

,

…

,

u

n

u_1, u_2, \ldots, u_n

u1​,u2​,…,un​。

幂法求最大特征值

步骤说明：

初始向量选取： 随机选取初始向量

x

1

x_1

x1​，可表示为特征向量的线性组合：

x

1

=

c

1

u

1

+

c

2

u

2

+

⋯

+

c

n

u

n

x_1 = c_1u_1 + c_2u_2 + \cdots + c_nu_n

x1​=c1​u1​+c2​u2​+⋯+cn​un​

迭代计算：

第一次迭代：

A

x

1

=

c

1

A

u

1

+

c

2

A

u

2

+

⋯

+

c

n

A

u

n

=

λ

1

c

1

x

2

Ax_1 = c_1Au_1 + c_2Au_2 + \cdots + c_nAu_n = λ_1c_1x_2

Ax1​=c1​Au1​+c2​Au2​+⋯+cn​Aun​=λ1​c1​x2​ 规范化后得到：

x

2

=

u

1

+

c

2

c

1

λ

2

λ

1

u

2

+

⋯

+

c

n

c

1

λ

n

λ

1

u

n

x_2 = u_1 + \frac{c_2}{c_1} \frac{λ_2}{λ_1}u_2 + \cdots + \frac{c_n}{c_1} \frac{λ_n}{λ_1}u_n

x2​=u1​+c1​c2​​λ1​λ2​​u2​+⋯+c1​cn​​λ1​λn​​un​

第二次迭代：

A

x

2

=

λ

1

u

1

+

c

2

c

1

λ

2

2

λ

1

u

2

+

⋯

+

c

n

c

1

λ

n

2

λ

1

u

n

=

λ

1

x

3

Ax_2 = λ_1u_1 + \frac{c_2}{c_1} \frac{λ_2^2}{λ_1}u_2 + \cdots + \frac{c_n}{c_1} \frac{λ_n^2}{λ_1}u_n = λ_1x_3

Ax2​=λ1​u1​+c1​c2​​λ1​λ22​​u2​+⋯+c1​cn​​λ1​λn2​​un​=λ1​x3​ 规范化后得到：

x

3

=

u

1

+

c

2

c

1

λ

2

2

λ

1

2

u

2

+

⋯

+

c

n

c

1

λ

n

2

λ

1

2

u

n

x_3 = u_1 + \frac{c_2}{c_1} \frac{λ_2^2}{λ_1^2}u_2 + \cdots + \frac{c_n}{c_1} \frac{λ_n^2}{λ_1^2}u_n

x3​=u1​+c1​c2​​λ12​λ22​​u2​+⋯+c1​cn​​λ12​λn2​​un​

通用迭代公式： 第

k

k

k 次迭代的通式为：

x

k

+

1

=

u

1

+

c

2

c

1

λ

2

k

λ

1

k

u

2

+

⋯

+

c

n

c

1

λ

n

k

λ

1

k

u

n

x_{k+1} = u_1 + \frac{c_2}{c_1} \frac{λ_2^k}{λ_1^k}u_2 + \cdots + \frac{c_n}{c_1} \frac{λ_n^k}{λ_1^k}u_n

xk+1​=u1​+c1​c2​​λ1k​λ2k​​u2​+⋯+c1​cn​​λ1k​λnk​​un​

收敛性分析： 由于

∣

λ

1

∣

>

∣

λ

i

∣

(

i

≥

2

)

|λ_1| > |λ_i| \, (i \geq 2)

∣λ1​∣>∣λi​∣(i≥2)，当

k

k

k 充分大时，高阶小项趋于零，可得：

A

x

k

+

1

≈

λ

1

u

1

,

x

k

+

1

≈

u

1

Ax_{k+1} \approx λ_1u_1, \quad x_{k+1} \approx u_1

Axk+1​≈λ1​u1​,xk+1​≈u1​

具体实现步骤：

随机初始化非零向量

v

0

\boldsymbol{v}_0

v0​。迭代计算：

v

k

+

1

=

A

v

k

∥

A

v

k

∥

\boldsymbol{v}_{k+1} = \frac{A\boldsymbol{v}_k}{\|A\boldsymbol{v}_k\|}

vk+1​=∥Avk​∥Avk​​估计特征值：

λ

≈

v

k

⊤

A

v

k

\lambda \approx \boldsymbol{v}_k^\top A \boldsymbol{v}_k

λ≈vk⊤​Avk​

编程实现

具体实现时，并没有λ1和u1的值，因此，迭代计算

x

k

+

1

=

A

x

k

x_{k+1}=Ax_k

xk+1​=Axk​后，规范化

x

k

+

1

x_{k+1}

xk+1​即可。注意：最大特征值是指模最大的那个特征值

% 幂法求最大特征值

clc;

clear;

close all;

% 第一种写法

A=[4 2 -2; -2 8 1 ; 2 4 -4];

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)

eig(A)

% 第二种写法

v0 = [1;1;1];

u0 = [1;1;1];

% A = [2,-1,0;-1,2,-1;0,-1,2];

v = A * u0;

u = v / norm(v, inf);

i = 0;

while norm(u - u0, inf) >= 1e-6

u0 = u;

v = A * u0;

u = v / norm(v, inf);

i = i+1;

end

norm(v, inf)

i

u

补充说明

最大特征值： 幂法求得的是模最大的特征值

λ

1

λ_1

λ1​。

2. 逆幂迭代法（Inverse Iteration）

目标：求解靠近

μ

\mu

μ 的最小模特征值。

给定矩阵 ( A )，假设其有 ( n ) 个实特征值：

∣

λ

1

∣

>

∣

λ

2

∣

>

⋯

>

∣

λ

n

∣

|\lambda_1| > |\lambda_2| > \cdots > |\lambda_n|

∣λ1​∣>∣λ2​∣>⋯>∣λn​∣

其对应的特征向量为(

u

1

,

u

2

u_1, u_2

u1​,u2​,

…

,

u

n

\ldots, u_n

…,un​)。(

λ

n

\lambda_n

λn​ ) 是最小特征值。首先注意到如果 (

A

u

=

λ

u

Au= \lambda u

Au=λu )，则：

A

−

1

A

u

=

A

−

1

λ

u

⟹

u

=

A

−

1

λ

u

A^{-1}Au = A^{-1}\lambda u \implies u = A^{-1}\lambda u

A−1Au=A−1λu⟹u=A−1λu

因此有：

A

−

1

u

=

1

λ

u

A^{-1}u = \frac{1}{\lambda}u

A−1u=λ1​u

可以看到，当

λ

n

\lambda_n

λn​为矩阵 A 的最小特征值时，(

1

λ

n

\frac{1}{\lambda_n}

λn​1​ ) 将是

A

−

1

A^{-1}

A−1的最大特征值。此时运用幂法求解

A

−

1

A^{-1}

A−1 的最大特征值，取倒数，即为 A 的最小特征值。反幂算法中需要注意的是，当最小特征值为 0 时，其倒数是没有定义的，此时反幂法求解的是第二小的特征值，且需要采用移位反幂法。

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

−

μ

I

)

(A - \mu I)

(A−μI) 进行 LU 分解。随机初始化向量

v

0

\boldsymbol{v}_0

v0​。迭代求解：

(

A

−

μ

I

)

v

k

+

1

=

v

k

⇒

v

k

+

1

=

v

k

+

1

∥

v

k

+

1

∥

(A - \mu I)\boldsymbol{v}_{k+1} = \boldsymbol{v}_k \quad \Rightarrow \quad \boldsymbol{v}_{k+1} = \frac{\boldsymbol{v}_{k+1}}{\|\boldsymbol{v}_{k+1}\|}

(A−μI)vk+1​=vk​⇒vk+1​=∥vk+1​∥vk+1​​

A = [3,0,-10;-1,3,4;0,1,-2];

I = eye(3,3);

p = 4.3;

u0 = [1;1;1];

v = inv(A - p * I) * u0;

u = v / norm(v, inf);

i = 0;

while norm(u - u0, inf) > 1e-5

u0 = u;

v = inv(A - p * I) * u0;

u = v / norm(v, inf);

i ++;

end;

i

u

x = p + 1 / norm(v, inf)

综述所述，可以总结反幂法求解特征向量的特点如下：

位移技术： 对每个已求得的特征值

λ

i

\lambda_i

λi​，构造矩阵

A

0

−

λ

i

I

A_0-\lambda_iI

A0​−λi​I，使其接近奇异；

加速收敛： 反幂法迭代公式为

x

k

+

1

=

(

A

−

σ

I

)

−

1

x

k

x_{k+1}=(A-\sigma I)^{-1}x_k

xk+1​=(A−σI)−1xk​，其中

σ

\sigma

σ接近特征值。此时

(

A

−

σ

I

)

−

1

(A-\sigma I)^{-1}

(A−σI)−1的模最大特征值对应的特征向量即为A的

σ

\sigma

σ附近特征值的特征向量。

高精度优势： 当

λ

i

\lambda_i

λi​精度较高时，反幂法可以在少量迭代内快速收敛到对应特征向量。

3. QR 算法（QR Algorithm）——稠密矩阵

目标：求解所有特征值（稠密矩阵）。 步骤：

将

A

A

A 转化为上 Hessenberg 矩阵。迭代 QR 分解：

A

k

=

Q

k

R

k

,

A

k

+

1

=

R

k

Q

k

A_k = Q_k R_k, \quad A_{k+1} = R_k Q_k

Ak​=Qk​Rk​,Ak+1​=Rk​Qk​当

A

k

A_k

Ak​ 收敛为上三角矩阵时，对角线元素即为特征值。

理论推导

幂法与反幂法用于求解矩阵的最大特征值与最小特征值。若想求解矩阵的所有特征值，可以使用QR分解法。假设矩阵 A 是

n

×

n

n \times n

n×n 的方阵，且其 n 个特征值均为互不相同的实数。QR分解法的理论保证如下：

若对矩阵 A 进行相似变换

B

=

C

−

1

A

C

B = C^{-1}AC

B=C−1AC ，则变换后的矩阵 B 的特征值与 A 一致。这是因为： 若

A

u

=

λ

u

Au = \lambda u

Au=λu ，令

v

=

C

−

1

u

v = C^{-1}u

v=C−1u，则有

A

C

v

=

A

u

=

λ

C

v

ACv = Au = \lambda Cv

ACv=Au=λCv 进一步可得

C

−

1

A

C

v

=

λ

v

C^{-1}ACv = \lambda v

C−1ACv=λv 因此，

λ

\lambda

λ 也是

C

−

1

A

C

C^{-1}AC

C−1AC 的特征值。

据此，可以通过以下步骤实现特征值和特征向量的求解：

初始化

令 ( A_1 = A )，并对 ( A_1 ) 进行QR分解：

A

1

=

Q

1

R

1

A_1 = Q_1R_1

A1​=Q1​R1​ 其中，( Q_1 ) 是正交矩阵（满足 ( Q_1Q_1^T = I )），( R_1 ) 是上三角矩阵。 迭代生成新矩阵

计算 ( A_2 = R_1Q_1 )，即：

A

2

=

Q

1

−

1

A

1

Q

1

A_2 = Q_1^{-1}A_1Q_1

A2​=Q1−1​A1​Q1​ 此时 ( A_2 ) 的特征值与 ( A ) 一致。继续对 ( A_2 ) 进行QR分解：

A

2

=

Q

2

R

2

A_2 = Q_2R_2

A2​=Q2​R2​ 重复迭代

计算 ( A_3 = R_2Q_2 )，即：

A

3

=

Q

2

−

1

A

2

Q

2

=

Q

2

−

1

Q

1

−

1

A

1

Q

1

Q

2

A_3 = Q_2^{-1}A_2Q_2 = Q_2^{-1}Q_1^{-1}A_1Q_1Q_2

A3​=Q2−1​A2​Q2​=Q2−1​Q1−1​A1​Q1​Q2​ 终止条件

重复上述步骤，直至 ( A_n ) 收敛为一个上三角矩阵。此时，矩阵对角线上的元素即为 ( A ) 的所有特征值。

注：QR分解通过不断迭代将原矩阵相似变换为上三角矩阵，从而直接读取对角线元素作为特征值。此方法适用于实对称矩阵或具有实特征值的方阵。

编程实现

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)

特点：

复杂度

O

(

n

3

)

O(n^3)

O(n3)，LAPACK 的核心算法。结合位移（如 Wilkinson 位移）优化收敛。

‌特征值修正‌ QR方法得到的特征值可能存在微小误差，反幂法可进一步修正

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)

Q

eig(A0)

% 使用反幂法求特征向量，并对特征值进行修正

a = ev;

n = size(a,1);

x = zeros(n);

for i = 1:n

x0 = ones(n,1);

[b,x0] = MinEig(A0-a(i,1)*eye(n));

x(:,i) = x0;

a(i,:) = a(i,:) + b;

end

a

x

4. 雅可比方法（Jacobi Method）——对称矩阵

目标：求解对称矩阵的所有特征值和特征向量。

Jacobi方法的基本思想是通过一次正交变换,将A中的一对非零的非对角元素化成零并且使得非对角元素的平方和减小。反复进行上述过程,使变换后的矩阵的非对角元素的平方和趋于零，从而使该矩阵近似为对角矩阵，得到全部特征值和特征向量。

步骤：

通过 Givens 旋转矩阵

G

k

G_k

Gk​ 逐步对角化：

A

k

+

1

=

G

k

⊤

A

k

G

k

A_{k+1} = G_k^\top A_k G_k

Ak+1​=Gk⊤​Ak​Gk​重复直到非对角元素接近零。

特点：

稳定但收敛慢，特征向量通过旋转矩阵累积。

编程实现

function [D,V,iter]=Jacobi_classical(A,maxIter,tol)

n = size(A, 1); % 矩阵的大小

V = eye(n); % 初始化特征向量矩阵为单位矩阵

iter = 0; % 初始化迭代次数

% 设置最大迭代次数和误差精度

if nargin < 3 || isempty(tol)

tol = 1e-9; % 默认误差精度

end

if nargin < 2 || isempty(maxIter)

maxIter = 1000; % 默认最大迭代次数

end

while(iter < maxIter)

iter=iter+1;

D=A;

n=size(D,1);

p=1;q=2;

for i=1:n

for j=i+1:n

if(abs(D(i,j))>abs(D(p,q)))%找到对称矩阵的上三角矩阵中最大的元素的下标

p=i;q=j;

end

end

end

if(abs(D(p,q))

break;

end

if(A(p,q)~=0)

d=(A(q,q)-A(p,p))/(2*A(p,q));

if(d>0)

t=1/(d+sqrt(d^2+1));

else

t=-1/(-d+sqrt(d^2+1));

end

c=1/sqrt(t^2+1);s=c*t;

else

c=1;s=0;

end

R=[c s;-s c];

A([p,q],:)=R'*A([p,q],:);

A(:,[p,q])=A(:,[p,q])*R;

V(:, [p, q]) = V(:,[p,q])*R;

end

D = diag(diag(D)); % 提取特征值

end

结果测试：

clc;

clear;

close all;

A=[ 6 -7 2 ; 4 -5 2; 1 -1 1]

[D, V, iter] = Jacobi_classical(A, 2000)

eig(A)

5. Lanczos 算法（稀疏矩阵）

目标：求解稀疏矩阵的部分极端特征值。 步骤：

生成 Krylov 子空间的正交基底。投影到三对角矩阵

T

k

T_k

Tk​：

T

k

=

V

k

⊤

A

V

k

T_k = V_k^\top A V_k

Tk​=Vk⊤​AVk​对

T

k

T_k

Tk​ 应用 QR 算法求特征值。 特点：

仅需矩阵-向量乘法，适合大规模稀疏矩阵。

6. 分治法

目标：高效求解对称三对角矩阵的所有特征值。 步骤：

将矩阵分解为子矩阵。递归求解子矩阵特征值。合并子问题解并修正。

特点：

复杂度

O

(

n

2

)

O(n^2)

O(n2)，适合大规模三对角矩阵。

7. 方法选择指南

场景推荐方法中小规模稠密矩阵QR 算法对称矩阵Jacobi 或 QR 算法稀疏矩阵的极端特征值Lanczos/Arnoldi 迭代最小/靠近

μ

\mu

μ 的特征值逆幂迭代法 + 位移工程问题中的部分特征值子空间迭代法

8. 关键公式与说明

特征方程 矩阵

A

A

A 的特征值满足：

det

⁡

(

A

−

λ

I

)

=

0

\det(A - \lambda I) = 0

det(A−λI)=0

2x2 矩阵：

λ

2

−

tr

(

A

)

λ

+

det

⁡

(

A

)

=

0

\lambda^2 - \text{tr}(A)\lambda + \det(A) = 0

λ2−tr(A)λ+det(A)=0n 阶矩阵：

P

(

λ

)

=

(

−

1

)

n

λ

n

+

⋯

+

det

⁡

(

A

)

P(\lambda) = (-1)^n \lambda^n + \dots + \det(A)

P(λ)=(−1)nλn+⋯+det(A)

提示：

高阶矩阵避免解析法，优先使用数值库（如 LAPACK、ARPACK）。对称矩阵的特征向量可正交化，提升计算稳定性。

参考文献

[1] *数值计算day5-特征值与特征向量 [2] 数值计算方法 Chapter7. 计算矩阵的特征值和特征向量 [3] 数值线性代数：Arnoldi求解特征值/特征向量 [4] 使用Matlab实现：幂法、反幂法(原点位移) [5] MATLAB求解矩阵特征值的六种方法