miniyuan

特征值和特征向量的计算（QR 分解，SVD 分解）

QR 分解

我们想要对 A∈Rn×nA \in \mathbb{R}^{n \times n} 进行 QR 分解，也即：

A=QRA = QR

其中 QQ 为正交阵，RR 为上三角阵。

Householder 变换

对单位向量 w∈Rnw \in \mathbb{R}^n，Householder 矩阵定义为：

H=I−2wwTH = I - 2ww^T

性质：

  • 对称：HT=HH^T = H
  • 正交：HTH=I⇒H2=IH^T H = I \Rightarrow H^2 = I

几何意义：将向量 xx 关于与 ww 正交的超平面 {y∣wTy=0}\{y \mid w^T y = 0\} 镜像反射。

定理：

对非零向量 x∈Rnx \in \mathbb{R}^n，存在 Householder 矩阵 HH，使得：

Hx=−σe1,σ=sign(x1)∥x∥2Hx = -\sigma e_1, \quad \sigma = \text{sign}(x_1)\|x\|_2

其中 e1=[1,0,…,0]Te_1 = [1,0,\dots,0]^T。

证明：

只需证明如下引理：对非零向量 x,y∈Rnx, y \in \mathbb{R}^n 且 ∥x∥2=∥y∥2\|x\|_2 = \|y\|_2，存在 Householder 矩阵 HH，使得：

Hx=yHx = y

由几何直观，构造 w=x−y∥x−y∥w = \dfrac{x - y}{\|x - y\|} 即可。

Householder QR 分解

输入：A∈Rn×nA \in \mathbb{R}^{n \times n}

输出：正交矩阵 QQ 和上三角矩阵 RR，使得 A=QRA = QR

  1. 初始化 A(1)=AA^{(1)} = A，Q(1)=InQ^{(1)} = I_n
  2. 对 k=1,2,…,n−1k = 1, 2, \dots, n-1：
    1. 取 x(k)=A(k)[k:n,k]∈R(n−k+1)x^{(k)} = A^{(k)}[k:n, k] \in \mathbb{R}^{(n-k+1)}
    2. 构造 Householder 矩阵 H(k)∈R(n−k+1)×(n−k+1)H^{(k)} \in \mathbb{R}^{(n-k+1) \times (n-k+1)} 使 H(k)x(k)=sign(x1(k))∥x(k)∥2e1H^{(k)} x^{(k)} = \text{sign}(x_1^{(k)})\|x^{(k)}\|_2 e_1
    3. 令 H~(k)=[Ik−100H(k)]\tilde{H}^{(k)} = \begin{bmatrix} I_{k-1} & 0 \\ 0 & H^{(k)} \end{bmatrix}，注意其仍为 Householder 矩阵
    4. A(k+1)=H~(k)A(k)A^{(k+1)} = \tilde{H}^{(k)} A^{(k)}
    5. Q(k+1)=Q(k)(H~(k))TQ^{(k+1)} = Q^{(k)} (\tilde{H}^{(k)})^T
  3. 令 Q=Q(n)Q = Q^{(n)}，R=A(n)R = A^{(n)}，输出 Q,RQ, R 即为所求

证明：

下归纳证明第 k−1k-1 步后，A(k)A^{(k)} 具有以下结构：

A(k)=[R(k)C(k)0D(k)]A^{(k)} = \begin{bmatrix} R^{(k)} & C^{(k)} \\ 0 & D^{(k)} \end{bmatrix}

其中 R(k)∈R(k−1)×(k−1)R^{(k)} \in \mathbb{R}^{(k-1) \times (k-1)} 是上三角矩阵，D(k)∈R(n−k+1)×(n−k+1)D^{(k)} \in \mathbb{R}^{(n-k+1) \times (n-k+1)}。

k=1k=1 时，A(1)=AA^{(1)} = A，R(1)R^{(1)} 为空，成立。假设 k−1k-1 时成立，下证 kk 时成立。

设 x(k)=D(k)[:,1]x^{(k)} = D^{(k)}[:, 1]。 构造 H(k)H^{(k)} 使 H(k)x(k)=∥x(k)∥2e1H^{(k)} x^{(k)} = \|x^{(k)}\|_2 e_1。 令 H~(k)=diag(Ik−1,H(k))\tilde{H}^{(k)} = \text{diag}(I_{k-1}, H^{(k)})，则：

A(k+1)=H~(k)A(k)=[R(k)C(k)0H(k)D(k)]A^{(k+1)} = \tilde{H}^{(k)} A^{(k)} = \begin{bmatrix} R^{(k)} & C^{(k)} \\ 0 & H^{(k)} D^{(k)} \end{bmatrix}

由于 H(k)D(k)H^{(k)} D^{(k)} 的第一列变为 ∥x(k)∥2e1\|x^{(k)}\|_2 e_1，故 A(k+1)A^{(k+1)} 的左上角 k×kk \times k 块

[R(k)C:,1(k)0∥x(k)∥2]\begin{bmatrix} R^{(k)} & C^{(k)}_{:,1} \\ 0 & \|x^{(k)}\|_2 \end{bmatrix}

是上三角矩阵。记此块为 R(k+1)R^{(k+1)}，则：

A(k+1)=[R(k+1)C(k+1)0D(k+1)]A^{(k+1)} = \begin{bmatrix} R^{(k+1)} & C^{(k+1)} \\ 0 & D^{(k+1)} \end{bmatrix}

满足归纳假设。

复杂度：O(2n3/3)\mathcal{O}(2n^3/3) 次浮点运算，优于 Givens。但稀疏时 Givens 更优。

Givens 旋转变换

Givens 旋转用于有选择地消去矩阵中的特定元素，每次只影响两行（列）。

2×2 Givens 旋转：

对向量 [ab]\begin{bmatrix} a \\ b \end{bmatrix}，构造旋转矩阵：

G=[cs−sc],c=aa2+b2,s=ba2+b2G = \begin{bmatrix} c & s \\ -s & c \end{bmatrix}, \quad c = \frac{a}{\sqrt{a^2+b^2}}, \quad s = \frac{b}{\sqrt{a^2+b^2}}

则：

G[ab]=[a2+b20]G \begin{bmatrix} a \\ b \end{bmatrix} = \begin{bmatrix} \sqrt{a^2+b^2} \\ 0 \end{bmatrix}

也即将向量 (a,b)(a,b) 旋转到 xx 轴正方向。

n×n Givens 旋转：

构造 Gij(θ)∈Rn×nG_{ij}(\theta) \in \mathbb{R}^{n \times n}，仅在以下位置非单位元：

  • Gii=c,Gjj=cG_{ii} = c,\quad G_{jj} = c
  • Gij=s,Gji=−sG_{ij} = s,\quad G_{ji} = -s
  • 其余对角元为 1，其余非对角元为 0

从而：

  • 左乘 GijTG_{ij}^T：只修改第 ii 行和第 jj 行
  • 右乘 GijG_{ij}：只修改第 ii 列和第 jj 列

注： 显然之前的 Jacobi 法求特征值也是利用了 Givens 旋转。

Givens QR 分解

输入：A∈Rn×nA \in \mathbb{R}^{n \times n}

输出：正交矩阵 QQ 和上三角矩阵 RR，使得 A=QRA = QR

  1. 初始化 A(1)=AA^{(1)} = A，Q(1)=InQ^{(1)} = I_n
  2. 对 k=1,2,…,n−1k = 1, 2, \dots, n-1， 对 i=k+1,…,ni = k+1, \dots, n：
    1. 取 a=A(k)[k,k]a = A^{(k)}[k, k]，b=A(k)[i,k]b = A^{(k)}[i, k]
    2. 计算 c=aa2+b2c = \frac{a}{\sqrt{a^2+b^2}}，s=ba2+b2s = \frac{b}{\sqrt{a^2+b^2}}
    3. 构造 Givens 矩阵 Gik∈Rn×nG_{ik} \in \mathbb{R}^{n \times n}： Gik=[1⋱c⋯s⋮⋱⋮−s⋯c1]G_{ik} = \begin{bmatrix} 1 & & & & & \\ & \ddots & & & & \\ & & c & \cdots & s & \\ & & \vdots & \ddots & \vdots & \\ & & -s & \cdots & c & \\ & & & & & 1 \end{bmatrix} 其中非零元位于 (k,k),(k,i),(i,k),(i,i)(k,k), (k,i), (i,k), (i,i) 位置
    4. A(k+1)=GikTA(k)A^{(k+1)} = G_{ik}^T A^{(k)}
    5. Q(k+1)=Q(k)GikQ^{(k+1)} = Q^{(k)} G_{ik}
  3. 令 Q=Q(n)Q = Q^{(n)}，R=A(n)R = A^{(n)}，输出 Q,RQ, R 即为所求

证明：

第 kk 步固定，对 i=k+1,…,ni = k+1, \dots, n 依次消去 A[i,k]A[i, k]。

设消去 A[i,k]A[i, k] 前，当前矩阵为 AA。取 a=A[k,k]a = A[k, k]，b=A[i,k]b = A[i, k]，构造 GikG_{ik} 使得：

GikT[ab]=[a2+b20]G_{ik}^T \begin{bmatrix} a \\ b \end{bmatrix} = \begin{bmatrix} \sqrt{a^2+b^2} \\ 0 \end{bmatrix}

左乘 GikTG_{ik}^T 只修改第 kk 行和第 ii 行，故：

  • 新 (i,k)(i, k) 位置变为 0
  • 第 kk 行第 kk 列变为 a2+b2\sqrt{a^2+b^2}
  • 其他已消零的位置不受影响

对 k=1,…,n−1k = 1, \dots, n-1 依次执行，最终 A(n)A^{(n)} 下三角全为 0，即为上三角矩阵 RR。

且正交阵乘积 Q=G12G13⋯Gn−1,nQ = G_{12} G_{13} \cdots G_{n-1,n} 仍为正交阵。

复杂度：O(4n3/3)\mathcal{O}(4n^3/3) 次浮点运算，比 Householder 多约一倍。但 Givens 每次只改两行，适合稀疏矩阵和并行计算。

两种 QR 分解对比

方法每次操作适用场景
Householder消去一整列稠密矩阵
Givens消去单个元素稀疏矩阵、并行计算

QR 算法

我们知道：

  • 容易求特征值的矩阵：对角阵、上三角阵、分块上三角阵（对角块小）
  • 保持特征值的变换：相似变换，正交相似变换（数值稳定）

我们利用好的变换把原矩阵化为好的矩阵，这就是 QR 算法。

实 Schur 分解

定义拟上三角矩阵（quasi-upper triangular）为对角块为 1 阶或 2 阶的分块上三角矩阵。

实 Schur 分解：

对任意 A∈Rn×nA \in \mathbb{R}^{n \times n}，存在正交矩阵 Q∈Rn×nQ \in \mathbb{R}^{n \times n}，使：

QTAQ=SQ^T A Q = S

其中 SS 为拟上三角矩阵，且：

  • 1 阶对角块对应一个实特征值
  • 2 阶对角块对应一对共轭复特征值

证明：略

QR 算法步骤

迭代步骤：

  1. QR 分解：A(k)=Q(k)R(k)A^{(k)} = Q^{(k)} R^{(k)}
  2. 更新：A(k+1)=R(k)Q(k)=(Q(k))TA(k)Q(k)A^{(k+1)} = R^{(k)} Q^{(k)} = (Q^{(k)})^T A^{(k)} Q^{(k)}

显然 A(k+1)A^{(k+1)} 与 A(k)A^{(k)} 正交相似，故特征值不变。

收敛性：

设 AA 是 n×nn \times n 实矩阵，特征值满足：

∣λ1∣≥∣λ2∣≥⋯≥∣λn∣|\lambda_1| \ge |\lambda_2| \ge \cdots \ge |\lambda_n|

且等号仅出现在共轭复特征值对（即 λ=a±bi,b≠0\lambda = a \pm bi, b \ne 0）的情形。

则 QR 迭代产生的 A(k)A^{(k)} 收敛到拟上三角矩阵（实 Schur 标准型）。

证明：略

注：实际中常加入位移（shift）加速收敛。

QR 算法实现

k = 0

while k < N:
    k += 1

    # QR 分解
    Q, R = qr_decomposition(A)
    # 进行迭代
    A = R @ Q
    
    # 检查收敛
    if np.all(np.abs(np.diag(A, k=-1)) < EPS): # 次对角线
        break

3×3 实对称矩阵的特征对

3×3 实对称矩阵的特征值

对实对称矩阵 A=[aij]∈R3×3A = [a_{ij}] \in \mathbb{R}^{3\times 3}，特征方程为：

P(λ)=λ3−I1λ2+I2λ−I3=0P(\lambda) = \lambda^3 - I_1 \lambda^2 + I_2 \lambda - I_3 = 0

其中：

  • I1=tr(A)=a11+a22+a33I_1 = \text{tr}(A) = a_{11}+a_{22}+a_{33}
  • I2=12[(trA)2−tr(A2)]=a11a22+a22a33+a33a11−a122−a232−a132I_2 = \frac{1}{2}[(\text{tr} A)^2 - \text{tr}(A^2)] = a_{11}a_{22} + a_{22}a_{33} + a_{33}a_{11} - a_{12}^2 - a_{23}^2 - a_{13}^2
  • I3=det⁡AI_3 = \det A

令：

p=I12−3I2q=272I3+I13−92I1I2ϕ=13arctan⁡(p3−q2q)\begin{aligned} p &= I_1^2 - 3I_2 \\ q &= \frac{27}{2}I_3 + I_1^3 - \frac{9}{2}I_1 I_2 \\ \phi &= \frac{1}{3} \arctan\left(\frac{\sqrt{p^3 - q^2}}{q}\right) \end{aligned}

则三个实特征值为：

λ1=2p3cos⁡ϕ+I13λ2=2p3cos⁡(ϕ+2π3)+I13λ3=2p3cos⁡(ϕ+4π3)+I13\begin{aligned} \lambda_1 &= \frac{2\sqrt{p}}{3} \cos\phi + \frac{I_1}{3} \\ \lambda_2 &= \frac{2\sqrt{p}}{3} \cos\left(\phi + \frac{2\pi}{3}\right) + \frac{I_1}{3} \\ \lambda_3 &= \frac{2\sqrt{p}}{3} \cos\left(\phi + \frac{4\pi}{3}\right) + \frac{I_1}{3} \end{aligned}

注：对于实对称阵而言，p3−q2≥0p^3 - q^2 \ge 0 恒成立。

3×3 实对称矩阵的特征向量

对特征值 λi\lambda_i，解 (A−λiI)vi=0(A - \lambda_i I)v_i = 0。

取转置：

viT(A−λiI)=0v_i^T (A - \lambda_i I) = 0

右乘 e1,e2e_1, e_2 得：

viT(a1−λie1)=0viT(a2−λie2)=0\begin{aligned} v_i^T (a_1 - \lambda_i e_1) &= 0 \\ v_i^T (a_2 - \lambda_i e_2) &= 0 \end{aligned}

当 a1−λie1,a1−λie1a_1 - \lambda_i e_1, a_1 - \lambda_i e_1 线性无关时，取：

vi=(a1−λie1)×(a2−λie2)v_i = (a_1 - \lambda_i e_1) \times (a_2 - \lambda_i e_2)

即可。

注：当两向量线性相关时（如重特征值），需改用其他列组合或 SVD 求解

SVD 分解

SVD 分解的数学原理

对 A∈Rm×nA \in \mathbb{R}^{m \times n}，若存在 σ≥0, u∈Rm, v∈Rn\sigma \ge 0,\ u \in \mathbb{R}^m,\ v \in \mathbb{R}^n 满足：

Av=σu,ATu=σvAv = \sigma u, \quad A^T u = \sigma v

且 ∥u∥2=∥v∥2=1\|u\|_2 = \|v\|_2 = 1，则称 σ\sigma 为奇异值，u/vu/v 为对应的左/右奇异向量。

SVD 分解：

任意 A∈Rm×nA \in \mathbb{R}^{m \times n} 存在正交矩阵 U∈Rm×m, V∈Rn×nU \in \mathbb{R}^{m\times m},\ V \in \mathbb{R}^{n\times n} 和

Σ=[diag(σ1,…,σr)000]∈Rm×n\Sigma = \begin{bmatrix} \mathrm{diag}(\sigma_1,\dots,\sigma_r) & 0 \\ 0 & 0 \end{bmatrix} \in \mathbb{R}^{m\times n}

使得

A=UΣVTA = U \Sigma V^T

其中 σ1≥σ2≥⋯≥σr>0\sigma_1 \ge \sigma_2 \ge \cdots \ge \sigma_r > 0。

证明：略

性质：

  1. σi2\sigma_i^2 是 ATAA^T A（以及 AATAA^T）的非零特征值。

    证明：

    ATAvi=AT(σiui)=σiATui=σi(σivi)=σi2viA^T A v_i = A^T (\sigma_i u_i) = \sigma_i A^T u_i = \sigma_i (\sigma_i v_i) = \sigma_i^2 v_i

    类似：

    AATui=A(σivi)=σiAvi=σi2uiAA^T u_i = A (\sigma_i v_i) = \sigma_i A v_i = \sigma_i^2 u_i

    所以 σi2\sigma_i^2 是特征值，vi, uiv_i,\ u_i 是对应特征向量。

  2. viv_i 是 ATAA^T A 的特征向量，uiu_i 是 AATAA^T 的特征向量。

    证明：由上一步直接得出。

  3. ui=1σiAviu_i = \frac{1}{\sigma_i} A v_i

    证明： 从 Avi=σiuiA v_i = \sigma_i u_i，两边除以 σi>0\sigma_i > 0 即得。

  4. 对 i≠ji \ne j：

    uiTuj=0,viTvj=0u_i^T u_j = 0,\quad v_i^T v_j = 0

    证明： 取 i≠ji \ne j：

    σjuiTuj=uiTAvj=(ATui)Tvj=σiviTvj\sigma_j u_i^T u_j = u_i^T A v_j = (A^T u_i)^T v_j = \sigma_i v_i^T v_j

    若 σi≠σj\sigma_i \ne \sigma_j 或 i≠ji \ne j 对应不同特征值，由对称矩阵 ATAA^T A 的特征向量正交性得 viTvj=0v_i^T v_j = 0，进而 uiTuj=0u_i^T u_j = 0。 若特征值相同，可通过 Gram–Schmidt 正交化选正交基。

  5. 设 r=rank(A)r = \mathrm{rank}(A)，则：

    • {v1,…,vr}\{v_1,\dots,v_r\} 张成行空间 C(AT)C(A^T)
    • {vr+1,…,vn}\{v_{r+1},\dots,v_n\} 张成零空间 N(A)N(A)
    • {u1,…,ur}\{u_1,\dots,u_r\} 张成列空间 C(A)C(A)
    • {ur+1,…,um}\{u_{r+1},\dots,u_m\} 张成左零空间 N(AT)N(A^T)

    证明：

    • Avi=σiui≠0  ⟹  vi⊥N(A)Av_i = \sigma_i u_i \ne 0 \implies v_i \perp N(A)，且 dim⁡C(AT)=r\dim C(A^T) = r，因此 {v1,…,vr}\{v_1,\dots,v_r\} 是 C(AT)C(A^T) 的一组正交基。
    • vr+1,…,vnv_{r+1},\dots,v_n 对应 σ=0\sigma=0，即 Avi=0Av_i=0，故属于 N(A)N(A)，维数 n−rn-r。
    • 类似得 uiu_i 的结论。

  6. rank(A)=#{σi>0}\mathrm{rank}(A) = \#\{\sigma_i > 0\}

    证明： Avi=σiuiA v_i = \sigma_i u_i，前 rr 个 uiu_i 线性无关且属于 C(A)C(A)，后 n−rn-r 个 viv_i 映射到 0。因此 dim⁡C(A)=r\dim C(A) = r。

  7. ∥A∥2=σ1\|A\|_2 = \sigma_1

    证明： 对任意单位 x∈Rnx \in \mathbb{R}^n，x=∑αivix = \sum \alpha_i v_i，则

    ∥Ax∥2=∑σi2αi2≤σ12\|Ax\|^2 = \sum \sigma_i^2 \alpha_i^2 \le \sigma_1^2

    且 x=v1x=v_1 时等号成立。

Eckart-Young-Mirsky 定理

设矩阵的谱分解为 A=∑i=1rσiuiviTA = \sum_{i=1}^r \sigma_i u_i v_i^T，对 k<rk < r，定义截断 SVD：

Ak=∑i=1kσiuiviTA_k = \sum_{i=1}^k \sigma_i u_i v_i^T

则 AkA_k 是所有秩 ≤k\le k 矩阵中，在谱范数（2-范数）下对 AA 的最佳逼近：

∥A−Ak∥2=min⁡rank(B)≤k∥A−B∥2=σk+1\|A - A_k\|_2 = \min_{\text{rank}(B) \le k} \|A - B\|_2 = \sigma_{k+1}

证明：略

注：也即能量集中在前几个奇异值。

应用1：图像压缩

  • 灰度图等价于矩阵 A∈Rm×nA \in \mathbb{R}^{m \times n}
  • 作低秩近似，计算 Ak=UkΣkVkTA_k = U_k \Sigma_k V_k^T
  • 压缩率：原存储 mnmn，压缩后 k(m+n+1)k(m+n+1)（Uk,VkU_k, V_k 各 kk 列，Σk\Sigma_k 有 kk 个值）
  • 质量：由 σk+1\sigma_{k+1} 决定，kk 越大越清晰

应用2：Netflix 推荐系统

  • 用户-电影评分矩阵 A∈RN×MA \in \mathbb{R}^{N \times M}（稀疏）
  • 用SVD分解 A≈UkΣkVkTA \approx U_k \Sigma_k V_k^T
  • UkU_k 的行：用户隐因子向量 VkV_k 的列：电影隐因子向量
  • 预测用户 ii 对电影 jj 的评分：(UkΣkVkT)ij(U_k \Sigma_k V_k^T)_{ij}

应用3：LLM 微调中的 LoRA

  • 全参数微调：Wfine=W0+ΔWW_{\text{fine}} = W_0 + \Delta W，ΔW\Delta W 存储成本高

  • 关键观察：ΔW\Delta W 的奇异值衰减极快，故可作低秩近似

  • LoRA 公式：

    Wfine=W0+BAW_{\text{fine}} = W_0 + B A

    其中：

    • W0∈Rd×kW_0 \in \mathbb{R}^{d \times k} 为预训练权重矩阵，不参与微调，冻结。
    • B∈Rd×lB \in \mathbb{R}^{d \times l}，A∈Rl×kA \in \mathbb{R}^{l \times k}，l≪min⁡(d,k)l \ll \min(d,k) 为人为选定维数，训练时仅更新 A,BA, B。

  • ΔW\Delta W 的最优低秩近似为 ΔW≈UlΣlVlT\Delta W \approx U_l \Sigma_l V_l^T。 在数学上，我们可取 B=UlΣl1/2B = U_l \Sigma_l^{1/2} 和 A=Σl1/2VlTA = \Sigma_l^{1/2} V_l^T，使得 BA=UlΣlVlTBA = U_l \Sigma_l V_l^T。 这表明，最优低秩近似确实能被分解为 LoRA 所采用的 BABA 形式。

应用4：潜在语义分析（LSA）

  • 文档-词项矩阵 A∈Rm×nA \in \mathbb{R}^{m \times n}（稀疏，元素为 TF-IDF）
  • 用 SVD 分解 A≈UkΣkVkTA \approx U_k \Sigma_k V_k^T
  • UkΣkU_k \Sigma_k 的行：文档的语义向量（m×km \times k）
  • VkΣkV_k \Sigma_k 的行（或 ΣkVkT\Sigma_k V_k^T 的列）：词的语义向量（n×kn \times k）
  • 查询 qq（词的集合）的语义向量： qsem=qTVkΣk−1q_{\text{sem}} = q^T V_k \Sigma_k^{-1} 这样 qsemq_{\text{sem}} 与文档语义向量在同一空间，可计算余弦相似度。

复杂度与方法对比总结

方法适用问题时间复杂度空间复杂度
QR 算法一般矩阵特征值O(n3)O(n^3) / 迭代O(n2)O(n^2)
Householder QRQR 分解（稠密）O(2n3/3)O(2n^3/3)O(n2)O(n^2)
Givens QRQR 分解（稀疏/并行）O(n3)O(n^3)（系数更大）O(n2)O(n^2)
3×3 解析法3×3 实对称矩阵O(1)O(1)O(1)O(1)
SVD低秩近似、非方阵O(mn2)O(mn^2) (m≥nm \ge n)O(mn)O(mn)