在科学计算、图像重建、机器学习等领域，我们经常需要求解大型线性系统 Ax ≈ b，尤其是最小二乘问题 min ||Ax - b||₂。这种问题往往出现在系数矩阵 A 规模巨大、稀疏的情况下，直接使用 SVD 或 QR 分解会因内存和时间成本过高而不可行。这时，迭代方法就成了首选，而 LSQR（Least Squares QR）算法正是其中最经典、最稳定的方法之一。
LSQR 由 Paige 和 Saunders 于 1982 年提出，本质上是将共轭梯度法（CG）应用于正规方程 AᵀAx = Aᵀb，但通过巧妙的双对角化过程避免了直接形成 AᵀA，从而保持数值稳定性和对稀疏性的充分利用。它特别适合求解超定系统（m > n）、欠定系统或带阻尼的正则化问题 min ||[A; damp·I]x - [b; 0]||₂。
今天，我们来详细剖析一个 MATLAB 实现的 LSQR 函数。它支持直接传入矩阵 A 和右端向量 B（可多列），无需用户提供矩阵-向量乘法函数，简化了使用门槛。同时支持阻尼参数、迭代上限、收敛容差和迭代日志显示，非常适合实际工程应用。
LSQR 算法核心思想
LSQR 通过 Lanczos 双对角化过程逐步构建下双对角矩阵，同时使用平面旋转（Givens 旋转）将其转化为上双对角形式。在此过程中，它维护两个正交基：左向量 u 和右向量 v，满足：
β u = A v - α u
α v = Aᵀ u - β v
通过一系列旋转，算法逐步更新解 X、残差范数、A 的条件数估计等关键量，并在每步检查多种收敛条件，确保在达到精度或检测
