Lanczos 迭代算法概述

Lanczos 迭代算法概述

关键词：Krylov 子空间、三项递推、三对角矩阵、Ritz 值、对称稀疏矩阵、特征值/特征向量、CG、Lanczos‑Bidiagonalization

Lanczos 迭代（Lanczos iteration）是 针对实对称（或复埃尔米特）矩阵 的一种 Krylov 子空间方法
它在每一步只需要一次矩阵‑向量乘（A·v）, 并通过 三项递推 产生一组正交基, 使得原矩阵在这组基下被压缩为 三对角矩阵
三对角矩阵的特征值（Ritz 值）快速逼近原矩阵的极端特征值, 从而实现 大规模稀疏特征值、线性方程组、奇异值 等问题的求解

下面从原理、算法细节、实现要点以及典型应用四个层面展开说明

1. 原理

其核心思想是通过正交投影将大规模对称矩阵的特征值问题转化为更小的三对角矩阵问题, 从而降低计算复杂度

1.1 Krylov 子空间

给定对称矩阵 $A \in \mathbb R^{n \times n}$（或 Hermitian）和非零向量 $v_{1}$, 定义 Krylov 子空间

$$ \mathcal Km(A,v{1}) = \operatorname{span}{v*{1},Av*{1},A^{2}v*{1},\dots,A^{m-1}v*{1}}. $$

Lanczos 迭代的目标是 在 $\mathcal K_m$ 中构造一组正交基 ${q*{1},q*{2},\dots,q*{m}}$
（记为矩阵 $Q_m=[q*{1},\dots,q_{m}]$）, 使得

$$ Q_m^{!T} A Q_m = T_m, $$

其中 $T_m$ 为 实对称三对角矩阵

$$ Tm= \begin{pmatrix} \alpha{1} & \beta*{2} & & & 0 \ \beta*{2} & \alpha*{2} & \beta*{3} & & \ & \beta*{3} & \ddots & \ddots & \ & & \ddots & \alpha*{m-1} & \beta*{m} \ 0 & & & \beta*{m} & \alpha_{m} \end{pmatrix}. $$

Ritz 值（$T_m$ 的特征值）和 Ritz 向量（$Q_m y$, 其中 $y$ 为对应特征向量）分别近似原矩阵 $A$ 的特征值和特征向量

1.2 三项递推（核心公式）

因为 $A$ 对称, Krylov 基可以通过 三项递推 产生——这正是 Lanczos 方法的精髓
设已知 $q*{k-1},q*{k}$（其中 $q_{0}=0$）, 则

$$ \begin{aligned} w &= A q*{k} - \beta*{k} q*{k-1}, \ \alpha*{k} &= q*{k}^{!T} w,\ w &\leftarrow w - \alpha*{k} q*{k},\ \beta*{k+1} &= |w|{2},\ q{k+1} &= \begin{cases} \displaystyle \frac{w}{\beta*{k+1}}, & \beta*{k+1}\neq0,\[12pt] \text{停止（Krylov 子空间已达满秩）}, & \beta_{k+1}=0. \end{cases} \end{aligned} $$

这三个系数 $\alpha*{k},\beta*{k+1}$ 正好是三对角矩阵 $T_m$ 的对角元和次对角元
整个过程只涉及一次矩阵‑向量乘（A·q_k）和向量的点乘、向量加减、标量乘除, 计算量极低

1.3 正交性与失去正交性

在 无限精度 下, 上述递推产生的向量 ${q_k}$ 完全正交
但在 有限精度（浮点数）中, 正交性会逐渐丧失, 导致 伪特征值（ghost eigenvalues） 出现。常用的对策：

• 正交性维护：通过 Gram-Schmidt 正交化确保向量 $ {v_i} $ 的正交性, 但数值计算中可能因舍入误差导致正交性丢失, 需采用重新正交化策略。
• 对称性利用：作为 Arnoldi 算法的对称特例, Lanczos 算法仅需存储三个向量（当前、前一步、后一步）, 显著降低内存需求。

2. 算法细节（伪代码）

下面给出最简版的 Lanczos 迭代（不含重新正交化）：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19

Input : symmetric matrix A (or linear operator),
        initial vector v (‖v‖=1), max iterations m
Output: tridiagonal T_m (α,β), orthonormal basis Q_m

β_0 ← 0
q_0 ← 0
q_1 ← v

for k = 1,…,m
    w ← A * q_k - β_{k-1} * q_{k-1}
    α_k ← q_kᵀ * w
    w ← w - α_k * q_k
    β_k ← ‖w‖₂
    if β_k == 0 then break
    q_{k+1} ← w / β_k
end for

T_m ← tridiagonal(α_1,…,α_m ; β_2,…,β_m)
Q_m ← [q_1,…,q_m]

• 矩阵‑向量乘 A * q_k 是唯一的高成本操作, 若 A 为稀疏矩阵或只提供乘法子程序（如在量子 Monte‑Carlo 中）, 该算法仍然有效。
• 若需要 特征向量, 只需对 T_m 求特征分解 $T_m = Y\Lambda Y^{!T}$, 再回代得到近似特征向量 $x_i = Q_m Y e_i$。

2.1 计算 Ritz 值的步骤

1. 完成 m 步迭代得到 T_m。
2. 对 T_m（尺寸仅为 m×m）使用标准的 QR 或 分裂对称三对角求特征值 方法, 得到特征值 ${\theta_i}$（Ritz 值）和特征向量 ${y_i}$。
3. 近似原矩阵的特征向量为 $x_i = Q_m y_i$。

因为 T_m 规模很小（通常 m 只取几百）, 这一步几乎不耗时。

3. 与其它 Krylov 方法的关系

隐式重启（IRL）是 ARPACK、MATLAB eigs、SciPy eigsh 等库的核心技术, 它在每次 QR 迭代后把不感兴趣的 Ritz 向量“压回”子空间, 从而保持子空间维数固定, 适合 大规模稀疏特征值 计算。

4. 典型应用

4.1 大规模稀疏特征值（极值）

• 求最大/最小特征值：Lanczos 只需几百次迭代即可得到前几个极端特征值, 广泛用于结构动力学（模态分析）、量子化学（哈密顿量对角化）等。
• 谱聚类 / 图划分：对图拉普拉斯矩阵 $L$（对称、半正定）使用 Lanczos 求前几个最小非零特征值及对应特征向量, 得到嵌入空间用于聚类。

4.2 线性方程组（CG）

• 对称正定系统 $Ax=b$ 的 共轭梯度法 本质上是 Lanczos 在 Krylov 子空间上做最小二乘（或最小残差）求解。
• 在求解离散偏微分方程（有限元、有限差分）的大型稀疏系统时, CG 是首选迭代求解器。

4.3 奇异值分解（Lanczos Bidiagonalization）

• Golub–Kahan 双对角化（又称 Lanczos‑Bidiagonalization）把一般矩阵 $A$ 投影到 Krylov 子空间 $\mathcal K_m(A^{!T}A, v)$ 中, 得到上/下双对角矩阵 $B_m$。
• 用于 大型稀疏 SVD（如信息检索中的潜在语义分析、推荐系统的矩阵分解）以及 LSQR / LSMR 线性最小二乘求解。

4.4 计算矩阵函数（如指数、对数）

• 通过 Lanczos 生成的三对角矩阵 $T_m$, 可以使用 Gauss–Quadrature（Lanczos‑quadrature）近似 $\phi(A)v$（$\phi$ 为任意标量函数）, 常用于 量子动力学（$e^{-iHt}v$）和 热力学（$\exp(-\beta A)v$）的时间演化。

4.5 密度状态（DOS）与谱密度估计

• 采用 随机向量 与 Lanczos 迭代结合（Stochastic Lanczos Quadrature）, 可高效估计大矩阵的 谱密度 $\rho(\lambda)$, 在材料科学、电子结构计算中用于计算态密度。

4.6 模型降阶（Model Order Reduction）

• 在控制系统与电路仿真中, Lanczos 用于构造 Krylov 子空间投影, 得到低阶近似模型, 保持系统的瞬态和频域特性。

4.7 机器学习与大数据

• 核主成分分析（KPCA）、大规模图卷积网络 中, 需要求解稀疏相似度矩阵的前几大特征值/特征向量, Lanczos 提供了可伸缩的求解方案。

5. 实现示例（Python）

下面演示 使用 SciPy 的 eigsh（内部基于 Lanczos） 计算稀疏对称矩阵的前 5 个最小特征值。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19

import numpy as np
import scipy.sparse as sp
import scipy.sparse.linalg as spla

# 生成一个 10000×10000 的 5‑对角稀疏矩阵（离散 Laplacian）
n = 10000
diagonals = [
    -1*np.ones(n-1),   # 下对角
    2*np.ones(n),      # 主对角
    -1*np.ones(n-1)    # 上对角
]
A = sp.diags(diagonals, offsets=[-1, 0, 1], format='csr')

# 计算前 5 个最小特征值（不包括零特征值）
k = 5
eigvals, eigvecs = spla.eigsh(A, k=k, which='SM')  # SM: smallest magnitude

print("前 {} 个最小特征值:".format(k))
print(eigvals)

• eigsh 默认使用 隐式重启 Lanczos (IRL), 只需 O(k·nnz(A)) 的矩阵‑向量乘次数

• 若想自己手写最原始的 Lanczos（不带重启）, 可以参考下面的极简实现（仅用于教学）：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

def lanczos_basic(A, v, m):
    """A: scipy.sparse matrix or LinearOperator
       v: initial vector (must be normalized)
       m: number of Lanczos steps
       Returns α, β, Q (orthonormal basis)"""
    n = A.shape[0]
    Q = np.zeros((n, m+1))
    α = np.zeros(m)
    β = np.zeros(m+1)

    Q[:,0] = v
    q_prev = np.zeros(n)          # q_0 = 0
    for k in range(m):
        w = A.dot(Q[:,k]) - β[k] * q_prev
        α[k] = np.dot(Q[:,k], w)
        w = w - α[k] * Q[:,k]
        β[k+1] = np.linalg.norm(w)
        if β[k+1] < 1e-12:        # Krylov 子空间已满秩
            Q = Q[:,:k+1]
            α = α[:k+1]
            β = β[:k+2]
            break
        q_prev = Q[:,k]
        Q[:,k+1] = w / β[k+1]
    return α, β[1:], Q[:,:len(α)]

注意：上述 lanczos_basic 没有重新正交化, 实际大规模计算时请使用成熟库（ARPACK、SLEPc、PRIMME）或自行加入 选择性重新正交化。

6. 收敛性与误差分析

• 极端特征值（最大/最小）收敛最快；中间特征值往往需要更多迭代或使用 阻塞（deflation） 技术。
• Ritz 值 的误差可用 残差向量 $r_i = A x_i - \theta_i x_i$ 衡量, 残差范数 $|r_i|$ 随迭代指数衰减（在无重数、间隔较大的情况下）。
• 间隙（spectral gap）对收敛速率至关重要：若 $\lambda*{k}$ 与 $\lambda*{k+1}$ 之间的间隙大, 则前 $k$ 个 Ritz 值会迅速收敛。
• 误差界（如 Saad’s bound） $$ |\lambdai - \theta_i| \le |r_i| \frac{|p{i-1}|}{|p’{i-1}(\lambda_i)|}, $$ 其中 $p{i-1}$ 为对应的 Lanczos 多项式。

7. 常见变体与扩展

8. 算法局限性与改进方向

1. 正交性丢失问题

   • 数值误差可能导致迭代向量正交性退化, 影响结果精度。改进策略包括：
     • 完全重新正交化：每次迭代后对所有向量重新正交化（计算量较大）。
     • 局部重新正交化：仅对可能失去正交性的向量进行修正。

2. 起始向量选择敏感

   • 起始向量 $ v_1 $ 的选择影响收敛速度。可采用随机向量或基于问题先验知识的向量（如模态分析中加载模式向量）。

3. 并行化与分布式优化

   • 现代高性能计算环境下, 需结合并行计算技术（如 MPI、CUDA）加速大规模矩阵运算, 进一步拓展 Lanczos 算法的应用边界。

9. 小结

• Lanczos 迭代是针对对称（Hermitian）矩阵的高效 Krylov 子空间方法, 核心是 三项递推, 只需一次矩阵‑向量乘即可生成正交基, 使原矩阵在该基下压缩为 三对角矩阵。
• Ritz 值/向量提供了原矩阵特征值/特征向量的快速近似, 特别适合 极端特征值 的求解。
• 通过 重新正交化、隐式重启、块迭代 等技巧, Lanczos 能在 数十万至数百万维 的稀疏矩阵上保持数值稳定性。
• 与 CG、Arnoldi、Golub–Kahan 等 Krylov 方法相互关联, 构成了现代数值线性代数的核心工具箱。
• 典型应用包括：大规模稀疏特征值/奇异值计算、共轭梯度求解线性系统、谱聚类与图划分、密度状态估计、矩阵函数逼近、模型降阶、量子物理中的哈密顿量对角化等。

参考文献（供进一步阅读）

1. C. Lanczos, An Iteration Method for the Solution of the Eigenvalue Problem, 1950.
2. Y. Saad, Iterative Methods for Sparse Linear Systems, 2nd ed., SIAM, 2003.
3. G. H. Golub & C. F. Van Loan, Matrix Computations, 4th ed., Johns Hopkins, 2013.
4. D. C. Sorensen, Implicitly Restarted Arnoldi/Lanczos Methods for Large Scale Eigenvalue Problems, 1992.
5. J. H. Wilkinson, The Algebraic Eigenvalue Problem, Oxford, 1965.
6. E. Anderson et al., LAPACK Users’ Guide, 3rd ed., SIAM, 1999.