前言

定n×n維矩陣A，滿足下式的數λ稱作矩陣A的一個特征值：
$$Au=\lambda u$$
推廣形式

特征值問題可推廣到更一般的形式。假設 $$u=f(x)$$ 是一個連續函數，$$A=\frac{d}{dx}$$ 表示微分運算，則二階微分方程：
$$\frac{d^{2}u}{d{x}^2}=k^{2}u$$
可表示為：
$$A^{2}u=k^{2}u$$
這是特征值問題在微分算子中的表現形式。

---

特征方程與求解方法

根據定義 $$(A?\lambda I)\mathbf{u}=0$$
若 $A?\lambda I$ 非奇異，則方程只有零解。因此，特征值需滿足：
$$\det (A?\lambda I)=0$$
此方程稱為特征方程，其根即為矩陣 $A$ 的特征值。

---

示例

給定矩陣：
$$A=\left[\begin{array}{cc}1& 2\\ 3& 2\end{array}\right]$$
其特征方程為：
$$\det \left[\begin{array}{cc}1?\lambda & 2\\ 3& 2?\lambda \end{array}\right]=(1?\lambda )(2?\lambda )?6=0$$
展開并化簡：
$${\lambda }^2?3\lambda ?4=0?{\lambda }_1=4,{\lambda }_2=?1$$

然而，對于高階矩陣，特征值的解析解通常難以直接計算，需借助數值方法（如QR算法、冪迭代法等）進行求解。

---

1. 冪迭代法（Power Iteration）

目標：求解矩陣的模最大特征值及其對應特征向量。

冪法與反冪法求解矩陣特征值

本節介紹如何使用冪法和反冪法分別求解矩陣的模最大和模最小特征值。給定矩陣 $A$，假設其有 $n$ 個實特征值：
$$|{\lambda }_1|>|{\lambda }_2|>?>|{\lambda }_n|$$
對應的特征向量為 $u_1,u_2,\dots ,u_n$。

---

冪法求最大特征值

步驟說明：

1. 初始向量選取：
   隨機選取初始向量 $x_1$，可表示為特征向量的線性組合：
   $$x_1=c_1{u}_1+c_2{u}_2+?+c_n{u}_n$$

2. 迭代計算：

   • 第一次迭代：
     $$A{x}_1=c_{1}A{u}_1+c_{2}A{u}_2+?+c_{n}A{u}_n={\lambda }_1{c}_1{x}_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$$

   • 第二次迭代：
     $$A{x}_2={\lambda }_1{u}_1+\frac{c_2}{c_1}\frac{{\lambda }_2^2}{{\lambda }_1}{u}_2+?+\frac{c_n}{c_1}\frac{{\lambda }_n^2}{{\lambda }_1}{u}_n={\lambda }_1{x}_3$$
     規范化后得到：
     $$x_3=u_1+\frac{c_2}{c_1}\frac{{\lambda }_2^2}{{\lambda }_1^2}{u}_2+?+\frac{c_n}{c_1}\frac{{\lambda }_n^2}{{\lambda }_1^2}{u}_n$$

3. 通用迭代公式：
   第 $k$ 次迭代的通式為：
   $$x_{k+1}=u_1+\frac{c_2}{c_1}\frac{{\lambda }_2^k}{{\lambda }_1^k}{u}_2+?+\frac{c_n}{c_1}\frac{{\lambda }_n^k}{{\lambda }_1^k}{u}_n$$

4. 收斂性分析：
   由于 $|{\lambda }_1|>|{\lambda }_i|(i\ge 2)$，當 $k$ 充分大時，高階小項趨于零，可得：
   $$A{x}_{k+1}\approx {\lambda }_1{u}_1,x_{k+1}\approx u_1$$

---

具體實現步驟：

1. 隨機初始化非零向量 ${\mathbf{v}}_0$。
2. 迭代計算：
   $${\mathbf{v}}_{k+1}=\frac{A{\mathbf{v}}_k}{\Vert A{\mathbf{v}}_k\Vert }$$
3. 估計特征值：
   $$\lambda \approx {\mathbf{v}}_k^{?}A{\mathbf{v}}_k$$

編程實現

具體實現時，并沒有λ1和u1的值，因此，迭代計算$x_{k+1}=A{x}_k$后，規范化$x_{k+1}$即可。注意：最大特征值是指模最大的那個特征值

% 冪法求最大特征值
clc;
clear;
close all;
% 第一種寫法
A=[4 2 -2; -2 8 1 ; 2 4 -4];
x = ones(size(A));
for i=1:40x=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-6u0 = u;v = A * u0;u = v / norm(v, inf);i = i+1;
end
norm(v, inf)
i
u

補充說明

• 最大特征值： 冪法求得的是模最大的特征值 ${\lambda }_1$。

---

2. 逆冪迭代法（Inverse Iteration）

目標：求解靠近 $\mu$ 的最小模特征值。

給定矩陣 ( A )，假設其有 ( n ) 個實特征值：

$$|{\lambda }_1|>|{\lambda }_2|>?>|{\lambda }_n|$$

其對應的特征向量為( $u_1,u_2$, $\dots ,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 時，其倒數是沒有定義的，此時反冪法求解的是第二小的特征值，且需要采用移位反冪法。

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

移位反冪法

步驟：

1. 對 $(A?\mu I)$ 進行 LU 分解。
2. 隨機初始化向量 ${\mathbf{v}}_0$。
3. 迭代求解：
   $$(A?\mu I){\mathbf{v}}_{k+1}={\mathbf{v}}_k?{\mathbf{v}}_{k+1}=\frac{{\mathbf{v}}_{k+1}}{\Vert {\mathbf{v}}_{k+1}\Vert }$$

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-5u0 = u;v = inv(A - p * I) * u0;u = v / norm(v, inf);i ++;
end;
i
u
x = p + 1 / norm(v, inf)

綜述所述，可以總結反冪法求解特征向量的特點如下：

位移技術： 對每個已求得的特征值 ${\lambda }_i$，構造矩陣$A_0?{\lambda }_{i}I$，使其接近奇異；

加速收斂： 反冪法迭代公式為$x_{k+1}=(A?\sigma I{)}^{?1}{x}_k$，其中$\sigma$接近特征值。此時$(A?\sigma I{)}^{?1}$的模最大特征值對應的特征向量即為A的$\sigma$附近特征值的特征向量。

高精度優勢： 當${\lambda }_i$精度較高時，反冪法可以在少量迭代內快速收斂到對應特征向量。

---

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

目標：求解所有特征值（稠密矩陣）。
步驟：

1. 將 $A$ 轉化為上 Hessenberg 矩陣。
2. 迭代 QR 分解：
   $$A_k=Q_k{R}_k,A_{k+1}=R_k{Q}_k$$
3. 當 $A_k$ 收斂為上三角矩陣時，對角線元素即為特征值。

理論推導

冪法與反冪法用于求解矩陣的最大特征值與最小特征值。若想求解矩陣的所有特征值，可以使用QR分解法。假設矩陣 A 是 $n\times n$ 的方陣，且其 n 個特征值均為互不相同的實數。QR分解法的理論保證如下：

若對矩陣 A 進行相似變換 $B=C^{?1}AC$ ，則變換后的矩陣 B 的特征值與 A 一致。這是因為：
若 $Au=\lambda u$ ，令 $v=C^{?1}u$，則有
$$ACv=Au=\lambda Cv$$
進一步可得
$$C^{?1}ACv=\lambda v$$
因此，$\lambda$ 也是 $C^{?1}AC$ 的特征值。

據此，可以通過以下步驟實現特征值和特征向量的求解：

1. 初始化

   • 令 ( A_1 = A )，并對 ( A_1 ) 進行QR分解：
     $$A_1=Q_1{R}_1$$
     其中，( Q_1 ) 是正交矩陣（滿足 ( Q_1Q_1^T = I )），( R_1 ) 是上三角矩陣。

2. 迭代生成新矩陣

   • 計算 ( A_2 = R_1Q_1 )，即：
     $$A_2=Q_1^{?1}{A}_1{Q}_1$$
     此時 ( A_2 ) 的特征值與 ( A ) 一致。繼續對 ( A_2 ) 進行QR分解：
     $$A_2=Q_2{R}_2$$

3. 重復迭代

   • 計算 ( 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$$

4. 終止條件

   • 重復上述步驟，直至 ( 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)$，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:nx0 = ones(n,1);[b,x0] = MinEig(A0-a(i,1)*eye(n));x(:,i) = x0;a(i,:) = a(i,:) + b;
enda
x

---

4. 雅可比方法（Jacobi Method）——對稱矩陣

目標：求解對稱矩陣的所有特征值和特征向量。

Jacobi方法的基本思想是通過一次正交變換,將A中的一對非零的非對角元素化成零并且使得非對角元素的平方和減小。反復進行上述過程,使變換后的矩陣的非對角元素的平方和趨于零，從而使該矩陣近似為對角矩陣，得到全部特征值和特征向量。

步驟：

1. 通過 Givens 旋轉矩陣 $G_k$ 逐步對角化：
   $$A_{k+1}=G_k^{?}{A}_k{G}_k$$
2. 重復直到非對角元素接近零。

特點：

• 穩定但收斂慢，特征向量通過旋轉矩陣累積。

編程實現

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:nfor j=i+1:nif(abs(D(i,j))>abs(D(p,q)))%找到對稱矩陣的上三角矩陣中最大的元素的下標p=i;q=j;endendendif(abs(D(p,q))<tol)break;endif(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));elset=-1/(-d+sqrt(d^2+1));endc=1/sqrt(t^2+1);s=c*t;elsec=1;s=0;endR=[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 算法（稀疏矩陣）

目標：求解稀疏矩陣的部分極端特征值。
步驟：

1. 生成 Krylov 子空間的正交基底。
2. 投影到三對角矩陣 $T_k$：
   $$T_k=V_k^{?}A{V}_k$$
3. 對 $T_k$ 應用 QR 算法求特征值。
   特點：

• 僅需矩陣-向量乘法，適合大規模稀疏矩陣。

---

6. 分治法

目標：高效求解對稱三對角矩陣的所有特征值。
步驟：

1. 將矩陣分解為子矩陣。
2. 遞歸求解子矩陣特征值。
3. 合并子問題解并修正。

特點：

• 復雜度 $O(n^2)$，適合大規模三對角矩陣。

---

7. 方法選擇指南

場景	推薦方法
中小規模稠密矩陣	QR 算法
對稱矩陣	Jacobi 或 QR 算法
稀疏矩陣的極端特征值	Lanczos/Arnoldi 迭代
最小/靠近 $\mu$ 的特征值	逆冪迭代法 + 位移
工程問題中的部分特征值	子空間迭代法

---

8. 關鍵公式與說明

特征方程
矩陣 $A$ 的特征值滿足：
$$\det (A?\lambda I)=0$$

• 2x2 矩陣：
  $${\lambda }^2?\text{tr}(A)\lambda +\det (A)=0$$
• n 階矩陣：
  $$P(\lambda )=(?1{)}^n{\lambda }^n+?+\det (A)$$

---

提示：

• 高階矩陣避免解析法，優先使用數值庫（如 LAPACK、ARPACK）。
• 對稱矩陣的特征向量可正交化，提升計算穩定性。

參考文獻

[1] *數值計算day5-特征值與特征向量
[2] 數值計算方法 Chapter7. 計算矩陣的特征值和特征向量
[3] 數值線性代數：Arnoldi求解特征值/特征向量
[4] 使用Matlab實現：冪法、反冪法(原點位移)
[5] MATLAB求解矩陣特征值的六種方法