矩阵的QR分解

QR分解

定义:对任意实矩阵ARm×nA \in \mathbb{R}^{m \times n},存在正交矩阵QQ和上三角矩阵RR,使得

A=QR.A = QR.

  • 唯一性:若AA非奇异且为方阵(m=nm = n),且规定RR的对角元全为正,则分解唯一。
  • 经济分解:当m>nm > n时,QQm×nm \times n列正交矩阵,RRn×nn \times n上三角阵;当m<nm < n,分解形式需调整。

核心思想:通过正交变换(Householder反射或Givens旋转)逐步将矩阵AA转化为上三角阵RR,同时记录变换过程得到正交阵QQ

Householder变换

定义5.8:给定单位向量wRnw \in \mathbb{R}^n,Householder矩阵定义为:

H(w)=I2wwT,H(w) = I - 2ww^T,

其性质包括对称性(HT=HH^T = H)和正交性(HTH=IH^TH = I)。

几何意义:
-HxHx表示向量xx关于以ww为法向量的超平面的镜像反射。

  • 保持向量2-范数不变,即Hx2=x2\|Hx\|_2 = \|x\|_2

关键定理:

  • 定理5.18:若x,yRnx, y \in \mathbb{R}^n满足x2=y2\|x\|_2 = \|y\|_2,则存在Householder矩阵HH,使得Hx=yHx = y
  • 定理5.19:可构造HH将任意向量xx变换为±x2e1\pm \|x\|_2 e_1,即

Hx=[σ00],σ=sign(x1)x2. Hx = \begin{bmatrix} -\sigma \\ 0 \\ \vdots \\ 0 \end{bmatrix}, \quad \sigma = \text{sign}(x_1)\|x\|_2.

构造方法(取xxσe1-\sigma e_1的连线作为法向量,即沿着xxσe1-\sigma e_1的垂直平分面对xx做镜像):

  1. v=x+σe1v = x + \sigma e_1,其中σ=sign(x1)x2\sigma = \text{sign}(x_1)\|x\|_2
  2. 归一化:w=v/v2w = v / \|v\|_2
  3. 定义H=I2wwTH = I - 2ww^T

示例:
对向量a=[212]a = \begin{bmatrix} 2 \\ 1 \\ 2 \end{bmatrix},构造HH使其变换为[300]\begin{bmatrix} -3 \\ 0 \\ 0 \end{bmatrix}

  1. 计算σ=3\sigma = 3,则v=a+σe1=[512]v = a + \sigma e_1 = \begin{bmatrix} 5 \\ 1 \\ 2 \end{bmatrix}
    2.w=v/v2w = v / \|v\|_2,得H=I2wwTH = I - 2ww^T
  2. 验证:Ha=3e1Ha = -3e_1

Givens旋转变换

定义:Givens矩阵是平面旋转矩阵的推广,用于在特定平面内旋转向量以消元。

二维形式:

G=[cssc],c=cosθ,s=sinθ,G = \begin{bmatrix} c & s \\ -s & c \end{bmatrix}, \quad c = \cos\theta, \, s = \sin\theta,

满足GTG=IG^TG = I

消元原理:对向量[xixj]\begin{bmatrix} x_i \\ x_j \end{bmatrix},选择c=xi/αc = x_i / \alphas=xj/αs = x_j / \alpha,其中α=xi2+xj2\alpha = \sqrt{x_i^2 + x_j^2},使得

G[xixj]=[α0].G \begin{bmatrix} x_i \\ x_j \end{bmatrix} = \begin{bmatrix} \alpha \\ 0 \end{bmatrix}.

高维扩展:将二维Givens矩阵嵌入n维单位阵,仅修改两行两列。例如,对向量x=[2012]x = \begin{bmatrix} 2 \\ 0 \\ 1 \\ 2 \end{bmatrix}

  1. 构造G1G_1消去第三分量:

G1=[2/501/5001001/502/500001]    G1x=[5002]. G_1 = \begin{bmatrix} 2/\sqrt{5} & 0 & 1/\sqrt{5} & 0 \\ 0 & 1 & 0 & 0 \\ -1/\sqrt{5} & 0 & 2/\sqrt{5} & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix} \implies G_1x = \begin{bmatrix} \sqrt{5} \\ 0 \\ 0 \\ 2 \end{bmatrix}.

  1. 构造G2G_2消去第四分量:

G2=[5/3002/3010000102/3005/3]    G2G1x=[3000]. G_2 = \begin{bmatrix} \sqrt{5}/3 & 0 & 0 & 2/3 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ -2/3 & 0 & 0 & \sqrt{5}/3 \end{bmatrix} \implies G_2G_1x = \begin{bmatrix} 3 \\ 0 \\ 0 \\ 0 \end{bmatrix}.

优势:每次旋转仅影响两个元素,适合稀疏矩阵的局部消元。

QR分解的算法实现

Householder方法:

  1. 列消元:从第一列开始,逐列构造Householder矩阵HkH_k,消去当前列下方的非零元素。
  2. 矩阵更新:

HnH2H1A=R    A=H1H2HnR=QR, H_n \cdots H_2H_1A = R \quad \implies \quad A = H_1H_2 \cdots H_n R = QR,

其中Q=H1H2HnQ = H_1H_2 \cdots H_n
3. 经济存储:对于ARm×nA \in \mathbb{R}^{m \times n}mnm \geq n),仅需存储前nn个Householder向量。

Givens方法:

  1. 元素消元:从矩阵左下角开始,逐行消元。对每个非零元素aija_{ij}i>ji > j),构造Givens矩阵GijG_{ij}将其变为零。
  2. 累积变换:

GkG2G1A=R    A=G1TG2TGkTR=QR. G_{k} \cdots G_2G_1A = R \quad \implies \quad A = G_1^TG_2^T \cdots G_k^T R = QR.

计算矩阵所有特征值的方法

1. 实Schur分解

实Schur分解定理(Th5.21):
对任意实矩阵ARn×nA \in \mathbb{R}^{n \times n},存在正交矩阵QQ和拟上三角矩阵SS(实Schur型),使得

QTAQ=S,Q^T A Q = S,

其中SS的形式为分块上三角矩阵,其对角块为1阶或2阶实矩阵:

  • 1阶对角块:对应实特征值。
  • 2阶对角块:对应一对复共轭特征值(例如α±βi\alpha \pm \beta i)。

实Schur型的结构:
SS的形式为:

S=[S11S12S1m0S22S2m00Smm],S = \begin{bmatrix} S_{11} & S_{12} & \cdots & S_{1m} \\ 0 & S_{22} & \cdots & S_{2m} \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & S_{mm} \end{bmatrix},

其中每个SiiS_{ii}是1阶或2阶实矩阵:

  • SiiS_{ii}为1阶块,则其唯一元素即为AA的实特征值。
  • SiiS_{ii}为2阶块[abcd]\begin{bmatrix} a & b \\ c & d \end{bmatrix},其特征值为复数α±βi\alpha \pm \beta i,满足α=a+d2\alpha = \frac{a+d}{2}β=bc(ad)24\beta = \sqrt{bc - \frac{(a-d)^2}{4}}

存在性与唯一性:

  • 存在性:定理保证实Schur分解对任意实矩阵均存在。
  • 唯一性:若矩阵AA的特征值均为实数且互异,则实Schur型退化为对角矩阵(特征值按任意顺序排列)。对于重复特征值或复特征值,分解不唯一,但分块结构由特征值的代数重数决定。

优势与应用:

  • 避免复数运算:通过2阶对角块表示复共轭特征值对,避免直接处理复数矩阵。
  • 直接读取特征值:从分块结构中可直接提取所有实数和复共轭特征值对。
  • 数值稳定性:正交相似变换(Householder或Givens)保范性强,适合数值计算。

2. QR迭代算法

理论基础:
QR迭代通过正交相似变换序列{Ak}\{A_k\},逐步逼近实Schur型矩阵SS。迭代公式为:

Ak+1=RkQk,其中 Ak=QkRk(QR分解).A_{k+1} = R_k Q_k, \quad \text{其中 } A_k = Q_k R_k \text{(QR分解)}.

收敛性(Th5.22):

  • AA的等模特征值为实重根或复共轭对,则{Ak}\{A_k\}收敛到拟上三角阵(实Schur型)。
  • 对称矩阵的QR迭代收敛到对角阵,特征值为实数。

算法实现:

  1. 预处理为上Hessenberg型:
    通过Householder变换将AA化为上Hessenberg矩阵(次对角线以下全零):

[000]\begin{bmatrix}* & * & * & * \\ * & * & * & * \\ 0 & * & * & * \\ 0 & 0 & * & * \end{bmatrix}

优势:QR分解的计算量从O(n3)O(n^3)降至O(n2)O(n^2)

  1. Givens旋转的应用:
    • 对上Hessenberg矩阵执行QR分解时,逐行消元仅需n1n-1次Givens旋转。
    • 例:对AkA_k,构造Givens矩阵P1,,Pn1P_1, \dots, P_{n-1},使得Pn1P1Ak=RkP_{n-1} \cdots P_1 A_k = R_k,则Qk=(Pn1P1)TQ_k = (P_{n-1} \cdots P_1)^T,更新Ak+1=RkQkA_{k+1} = R_k Q_k

结构保持性:

  • 每步迭代后Ak+1A_{k+1}仍保持上Hessenberg型,确保算法高效性。

3. 收缩技术(Deflation)

核心步骤:

  1. 已知特征向量处理:若AA有特征值λ1\lambda_1及对应特征向量x1x_1,构造Householder变换HH,使Hx1=σe1Hx_1 = \sigma e_1
  2. 正交相似变换:计算HAHTHAH^T,其分块形式为:

HAHT=[λ1r1T0A1], HAH^T = \begin{bmatrix} \lambda_1 & r_1^T \\ 0 & A_1 \end{bmatrix},

其中A1A_1为降阶后的子矩阵,其特征值为AA的剩余特征值。

示例:
对矩阵A=[311201112]A = \begin{bmatrix} 3 & -1 & 1 \\ 2 & 0 & 1 \\ 1 & -1 & 2 \end{bmatrix},已知特征值λ1=2\lambda_1 = 2,特征向量x1=[1,1,0]Tx_1 = [1, 1, 0]^T

  1. 构造Householder变换HH,将x1x_1反射到σe1\sigma e_1
  2. 计算HAHTHAH^T,子矩阵A1A_1的特征值为1122

局限性:误差可能累积,且效率较低。


矩阵的QR分解
https://xiao-ao-jiang-hu.github.io/2025/05/26/numerical-analysis/numerical-analysis2/
作者
wst
发布于
2025年5月27日
许可协议