标题：基于Givens旋转的QR分解在求解最小二乘问题中的数值稳定实现
字数 4743 2025-12-24 03:04:41

好的，我将从线性代数算法领域中，选取一个你尚未学习过的经典算法进行详细讲解。

标题：基于Givens旋转的QR分解在求解最小二乘问题中的数值稳定实现

题目描述

我们考虑一个超定线性最小二乘问题：给定一个矩阵 \(A \in \mathbb{R}^{m \times n}\)（其中 \(m > n\)）和一个向量 \(b \in \mathbb{R}^{m}\)，目标是找到向量 \(x \in \mathbb{R}^{n}\)，使得残差 \(\|Ax - b\|_2\) 最小化。

直接求解正规方程 \(A^T A x = A^T b\) 在数值上可能不稳定，尤其当 \(A\) 是病态（接近秩亏）时，因为计算 \(A^T A\) 会使其条件数平方，放大舍入误差。一种更稳定、更精确的求解方法是直接对系数矩阵 \(A\) 进行QR分解。

具体任务是：详细讲解如何利用一系列Givens旋转，以高度稳定且适用于结构化矩阵（尤其是稀疏矩阵）的方式，对矩阵 \(A\) 进行QR分解，并利用分解结果高效求解最小二乘问题 \(x = \arg\min \|Ax - b\|_2\)。这包括分解过程、求解过程，以及为什么Givens旋转比Householder反射在某些场景下更具优势。

解题过程详解

第一步：理解核心思想——QR分解与最小二乘

  1. QR分解定理：任何实矩阵 \(A_{m \times n}\)（\(m \ge n\)）都可以分解为 \(A = QR\)。
    • \(Q\) 是一个 \(m \times m\) 的正交矩阵（\(Q^T Q = I\)）。
    • \(R\) 是一个 \(m \times n\) 的上三角矩阵（当 \(m > n\) 时，其底部是零行）。
  2. 应用到最小二乘问题：
    • 将 \(A = QR\) 代入目标函数：\(\|Ax - b\|_2^2 = \|QRx - b\|_2^2\)。
    • 利用正交变换不改变向量的2-范数：\(\|QRx - b\|_2 = \|Q(Rx - Q^T b)\|_2 = \|Rx - Q^T b\|_2\)。
    • 令 \(c = Q^T b\)，我们将原问题转化为求解 \(\min \|Rx - c\|_2\)。
    • 由于 \(R\) 是上三角矩阵（我们通常只使用其前 \(n\) 行非零部分，记作 \(\tilde{R}_{n \times n}\)），对应的 \(c\) 也取其前 \(n\) 个分量记作 \(\tilde{c}\)，则该问题等价于求解一个简单的上三角方程组 \(\tilde{R}x = \tilde{c}\)，可以通过回代法（Back Substitution） 快速求解。

所以，关键在于如何稳定地得到 \(Q\) 和 \(R\)，特别是 \(R\) 和 \(Q^T b\)。

第二步：引入构建工具——Givens旋转矩阵

Givens旋转是用于在二维平面内进行旋转的正交矩阵。在更高维空间中，它用于消去矩阵特定位置（例如，\(a_{ji}\)）的元素。

  1. 定义：一个Givens旋转矩阵 \(G(i, j, \theta)\) 是一个单位矩阵，仅在第 \(i\) 行 \(i\) 列、\(i\) 行 \(j\) 列、\(j\) 行 \(i\) 列、\(j\) 行 \(j\) 列处被修改：

\[ G_{ii} = c, \quad G_{ij} = s, \quad G_{ji} = -s, \quad G_{jj} = c, \quad \text{其余位置与单位矩阵相同}。 \]

其中 $ c = \cos\theta $，$ s = \sin\theta $，且满足 $ c^2 + s^2 = 1 $。
  1. 几何作用：左乘 \(G^T\)（等价于右乘 \(G\) 的逆，因为正交）到向量 \(v\) 上，其效果是在 \((i, j)\) 坐标平面内将 \(v\) 旋转角度 \(\theta\)，使得旋转后向量的第 \(j\) 个分量变为0。
  2. 数值计算：我们不需要显式计算角度 \(\theta\)。给定两个数 \(a\) 和 \(b\)（不全为0），我们希望构造 \(c\) 和 \(s\)，使得：

\[ \begin{bmatrix} c & s \\ -s & c \end{bmatrix}^T \begin{bmatrix} a \\ b \end{bmatrix} = \begin{bmatrix} r \\ 0 \end{bmatrix} \]

解得：
*   $ r = \sqrt{a^2 + b^2} $
*   为避免溢出和下溢，常用稳定计算：若 $ |b| > |a| $，则令 $ \tau = -a/b, s = 1/\sqrt{1+\tau^2}, c = s\tau $；否则令 $ \tau = -b/a, c = 1/\sqrt{1+\tau^2}, s = c\tau $。
*   这样，我们直接用原始数据 $ a, b $ 算出了 $ c, s $，实现了将 $ (a, b) $ 变为 $ (r, 0) $ 的“旋转变换”。

第三步：使用Givens旋转进行QR分解

我们的目标是通过一系列Givens旋转，将矩阵 \(A\) 逐步化为上三角矩阵 \(R\)。同时，这些旋转的累积效应就构成了正交矩阵 \(Q^T\)。

  1. 过程（逐列消元）：
    • 从第一列开始。对于第一列，我们需要消去 \(a_{21}, a_{31}, ..., a_{m1}\) 这些对角线以下的元素。
    • 以消去 \(a_{21}\) 为例：构造Givens旋转 \(G_1 = G(1, 2, \theta_1)\)，其中 \(c, s\) 根据 \(a_{11}\) 和 \(a_{21}\) 计算，使得 \(G_1^T\) 左乘 \(A\) 后，新的 \(a‘_{21} = 0\)。同时，\(b\) 向量也需要左乘 \(G_1^T\)：\(b := G_1^T b\)。
    • 接下来用 \(G(1, 3, \theta_2)\) 消去新的第一列第三行元素，依此类推，直到第一列对角线以下全为0。记这一步的所有旋转乘积为 \(Q_1^T\)。
    • 关键点：每个Givens旋转只改变它作用的两行。这在处理稀疏矩阵时优势巨大，因为可以只更新相关的少数元素，避免了Householder反射对一整列的稠密更新。
    • 然后处理第二列。注意，此时第一行已被“冻结”。从第二列的第二行元素 \(a_{22}\) 开始，消去其下方的 \(a_{32}, ..., a_{m2}\)。构造的Givens旋转矩阵维度仍是 \(m \times m\)，但只作用于第2行及以下的行。
    • 重复此过程，直到前 \(n\) 列都被处理完毕。最终我们得到：

\[ Q^T A = R, \quad Q^T b = c \]

    其中 $ Q^T = G_k^T ... G_2^T G_1^T $，$ R $ 是上三角矩阵，$ c $ 是变换后的右端项。

第四步：求解最小二乘解

  1. 经过上一步，我们得到了等价的上三角最小二乘问题 \(\min \|Rx - c\|_2\)。
  2. 将 \(R\) 和 \(c\) 进行分块：

\[ R = \begin{bmatrix} \tilde{R}_{n \times n} \\ 0_{(m-n) \times n} \end{bmatrix}, \quad c = \begin{bmatrix} \tilde{c}_{n \times 1} \\ d_{(m-n) \times 1} \end{bmatrix} \]

其中 $ \tilde{R} $ 是可逆的上三角矩阵（假设 $ A $ 列满秩）。
  1. 目标函数变为：

\[ \|Rx - c\|_2^2 = \|\tilde{R}x - \tilde{c}\|_2^2 + \|d\|_2^2 \]

$ \|d\|_2^2 $ 是一个与 $ x $ 无关的常数，即**残差平方和的最小值**。
  1. 因此，使目标函数最小化的 \(x\) 就是上三角方程组 \(\tilde{R} x = \tilde{c}\) 的解。
  2. 使用回代法求解 \(x\)：
    • 从最后一个方程开始：\(\tilde{r}_{nn} x_n = \tilde{c}_n\)，解得 \(x_n = \tilde{c}_n / \tilde{r}_{nn}\)。
    • 代入倒数第二个方程：\(\tilde{r}_{n-1, n-1} x_{n-1} + \tilde{r}_{n-1, n} x_n = \tilde{c}_{n-1}\)，解得 \(x_{n-1} = (\tilde{c}_{n-1} - \tilde{r}_{n-1, n} x_n) / \tilde{r}_{n-1, n-1}\)。
    • 以此类推，直到求出 \(x_1\)。

第五步：算法总结与优势分析

算法流程：

  1. 初始化：令 \(R = A\)，\(\tilde{b} = b\)。
  2. Givens QR分解：
    • for j = 1 to n (遍历列)
      • for i = m downto j+1 (从该列底部向上消元，顺序很重要！)
        • 若 R[i, j] 已为0（在稀疏矩阵中常见），则跳过。
        • 否则，针对元素 R[j, j] 和 R[i, j]，计算Givens旋转参数 \(c, s\)。
        • 将旋转应用到 R 的第 j 行和第 i 行（仅需更新从第 j 列到第 n 列的元素）。
        • 同样将旋转应用到 \(\tilde{b}\) 的第 j 个和第 i 个分量上。
  3. 回代求解：
    • 提取 \(\tilde{R}\)（即 R 的前 n 行前 n 列部分）和 \(\tilde{c}\)（即 \tilde{b} 的前 n 个分量）。
    • 调用回代法求解 \(\tilde{R} x = \tilde{c}\)，得到解向量 \(x\)。

为什么选择Givens旋转？

  1. 数值稳定性：与Householder反射一样，Givens旋转基于正交变换，非常稳定，不会放大误差。
  2. 适用于稀疏/结构化矩阵：这是其主要优势。对于带状矩阵、Hessenberg矩阵或一般稀疏矩阵，Householder反射会破坏稀疏性（产生“填入”），而Givens旋转可以精确地只消去需要为零的元素，并能利用现有的零元结构，保持矩阵的稀疏性，从而大幅节省存储和计算量。
  3. 并行性：多个不相交的Givens旋转（例如，消去奇数行和消去偶数行的操作）可以同时进行，适合并行计算。

通过以上五步，我们完成了基于Givens旋转的QR分解及其在求解最小二乘问题中的应用全流程。这种方法在需要高数值稳定性且矩阵具有特殊稀疏结构时，是首选的高效算法。

好的，我将从线性代数算法领域中，选取一个你尚未学习过的经典算法进行详细讲解。 标题：基于Givens旋转的QR分解在求解最小二乘问题中的数值稳定实现 题目描述 我们考虑一个超定线性最小二乘问题：给定一个矩阵 \( A \in \mathbb{R}^{m \times n} \)（其中 \( m > n \)）和一个向量 \( b \in \mathbb{R}^{m} \)，目标是找到向量 \( x \in \mathbb{R}^{n} \)，使得残差 \( \|Ax - b\|_ 2 \) 最小化。 直接求解正规方程 \( A^T A x = A^T b \) 在数值上可能不稳定，尤其当 \( A \) 是病态（接近秩亏）时，因为计算 \( A^T A \) 会使其条件数平方，放大舍入误差。一种更稳定、更精确的求解方法是直接对系数矩阵 \( A \) 进行QR分解。 具体任务是：详细讲解如何利用一系列 Givens旋转 ，以高度稳定且适用于结构化矩阵（尤其是稀疏矩阵）的方式，对矩阵 \( A \) 进行QR分解，并利用分解结果高效求解最小二乘问题 \( x = \arg\min \|Ax - b\|_ 2 \)。这包括分解过程、求解过程，以及为什么Givens旋转比Householder反射在某些场景下更具优势。 解题过程详解 第一步：理解核心思想——QR分解与最小二乘 QR分解定理 ：任何实矩阵 \( A_ {m \times n} \)（\( m \ge n \)）都可以分解为 \( A = QR \)。 \( Q \) 是一个 \( m \times m \) 的正交矩阵（\( Q^T Q = I \)）。 \( R \) 是一个 \( m \times n \) 的上三角矩阵（当 \( m > n \) 时，其底部是零行）。 应用到最小二乘问题 ： 将 \( A = QR \) 代入目标函数：\(\|Ax - b\|_ 2^2 = \|QRx - b\|_ 2^2\)。 利用正交变换不改变向量的2-范数：\(\|QRx - b\|_ 2 = \|Q(Rx - Q^T b)\|_ 2 = \|Rx - Q^T b\|_ 2\)。 令 \( c = Q^T b \)，我们将原问题转化为求解 \(\min \|Rx - c\|_ 2\)。 由于 \( R \) 是上三角矩阵（我们通常只使用其前 \( n \) 行非零部分，记作 \( \tilde{R}_ {n \times n} \)），对应的 \( c \) 也取其前 \( n \) 个分量记作 \( \tilde{c} \)，则该问题等价于求解一个简单的上三角方程组 \( \tilde{R}x = \tilde{c} \)，可以通过 回代法（Back Substitution） 快速求解。 所以，关键在于如何稳定地得到 \( Q \) 和 \( R \)，特别是 \( R \) 和 \( Q^T b \)。 第二步：引入构建工具——Givens旋转矩阵 Givens旋转是用于在二维平面内进行旋转的正交矩阵。在更高维空间中，它用于消去矩阵特定位置（例如，\( a_ {ji} \)）的元素。 定义 ：一个Givens旋转矩阵 \( G(i, j, \theta) \) 是一个单位矩阵，仅在第 \( i \) 行 \( i \) 列、\( i \) 行 \( j \) 列、\( j \) 行 \( i \) 列、\( j \) 行 \( j \) 列处被修改： \[ G_ {ii} = c, \quad G_ {ij} = s, \quad G_ {ji} = -s, \quad G_ {jj} = c, \quad \text{其余位置与单位矩阵相同}。 \] 其中 \( c = \cos\theta \)，\( s = \sin\theta \)，且满足 \( c^2 + s^2 = 1 \)。 几何作用 ：左乘 \( G^T \)（等价于右乘 \( G \) 的逆，因为正交）到向量 \( v \) 上，其效果是在 \( (i, j) \) 坐标平面内将 \( v \) 旋转角度 \( \theta \)，使得旋转后向量的第 \( j \) 个分量变为0。 数值计算 ：我们不需要显式计算角度 \( \theta \)。给定两个数 \( a \) 和 \( b \)（不全为0），我们希望构造 \( c \) 和 \( s \)，使得： \[ \begin{bmatrix} c & s \\ -s & c \end{bmatrix}^T \begin{bmatrix} a \\ b \end{bmatrix} = \begin{bmatrix} r \\ 0 \end{bmatrix} \] 解得： \( r = \sqrt{a^2 + b^2} \) 为避免溢出和下溢，常用稳定计算：若 \( |b| > |a| \)，则令 \( \tau = -a/b, s = 1/\sqrt{1+\tau^2}, c = s\tau \)；否则令 \( \tau = -b/a, c = 1/\sqrt{1+\tau^2}, s = c\tau \)。 这样，我们直接用原始数据 \( a, b \) 算出了 \( c, s \)，实现了将 \( (a, b) \) 变为 \( (r, 0) \) 的“旋转变换”。 第三步：使用Givens旋转进行QR分解 我们的目标是通过一系列Givens旋转，将矩阵 \( A \) 逐步化为上三角矩阵 \( R \)。同时，这些旋转的累积效应就构成了正交矩阵 \( Q^T \)。 过程（逐列消元） ： 从第一列开始。对于第一列，我们需要消去 \( a_ {21}, a_ {31}, ..., a_ {m1} \) 这些对角线以下的元素。 以消去 \( a_ {21} \) 为例：构造Givens旋转 \( G_ 1 = G(1, 2, \theta_ 1) \)，其中 \( c, s \) 根据 \( a_ {11} \) 和 \( a_ {21} \) 计算，使得 \( G_ 1^T \) 左乘 \( A \) 后，新的 \( a‘_ {21} = 0 \)。同时，\( b \) 向量也需要左乘 \( G_ 1^T \)：\( b := G_ 1^T b \)。 接下来用 \( G(1, 3, \theta_ 2) \) 消去新的第一列第三行元素，依此类推，直到第一列对角线以下全为0。记这一步的所有旋转乘积为 \( Q_ 1^T \)。 关键点 ：每个Givens旋转只改变它作用的两行。这在处理 稀疏矩阵 时优势巨大，因为可以只更新相关的少数元素，避免了Householder反射对一整列的稠密更新。 然后处理第二列。注意，此时第一行已被“冻结”。从第二列的第二行元素 \( a_ {22} \) 开始，消去其下方的 \( a_ {32}, ..., a_ {m2} \)。构造的Givens旋转矩阵维度仍是 \( m \times m \)，但只作用于第2行及以下的行。 重复此过程，直到前 \( n \) 列都被处理完毕。最终我们得到： \[ Q^T A = R, \quad Q^T b = c \] 其中 \( Q^T = G_ k^T ... G_ 2^T G_ 1^T \)，\( R \) 是上三角矩阵，\( c \) 是变换后的右端项。 第四步：求解最小二乘解 经过上一步，我们得到了等价的上三角最小二乘问题 \( \min \|Rx - c\|_ 2 \)。 将 \( R \) 和 \( c \) 进行分块： \[ R = \begin{bmatrix} \tilde{R} {n \times n} \\ 0 {(m-n) \times n} \end{bmatrix}, \quad c = \begin{bmatrix} \tilde{c} {n \times 1} \\ d {(m-n) \times 1} \end{bmatrix} \] 其中 \( \tilde{R} \) 是可逆的上三角矩阵（假设 \( A \) 列满秩）。 目标函数变为： \[ \|Rx - c\|_ 2^2 = \|\tilde{R}x - \tilde{c}\|_ 2^2 + \|d\|_ 2^2 \] \( \|d\|_ 2^2 \) 是一个与 \( x \) 无关的常数，即 残差平方和的最小值 。 因此，使目标函数最小化的 \( x \) 就是 上三角方程组 \( \tilde{R} x = \tilde{c} \) 的解。 使用 回代法 求解 \( x \)： 从最后一个方程开始：\( \tilde{r}_ {nn} x_ n = \tilde{c}_ n \)，解得 \( x_ n = \tilde{c} n / \tilde{r} {nn} \)。 代入倒数第二个方程：\( \tilde{r} {n-1, n-1} x {n-1} + \tilde{r} {n-1, n} x_ n = \tilde{c} {n-1} \)，解得 \( x_ {n-1} = (\tilde{c} {n-1} - \tilde{r} {n-1, n} x_ n) / \tilde{r}_ {n-1, n-1} \)。 以此类推，直到求出 \( x_ 1 \)。 第五步：算法总结与优势分析 算法流程 ： 初始化 ：令 \( R = A \)，\( \tilde{b} = b \)。 Givens QR分解 ： for j = 1 to n (遍历列) for i = m downto j+1 (从该列底部向上消元，顺序很重要！) 若 R[i, j] 已为0（在稀疏矩阵中常见），则跳过。 否则，针对元素 R[j, j] 和 R[i, j] ，计算Givens旋转参数 \( c, s \)。 将旋转应用到 R 的第 j 行和第 i 行（仅需更新从第 j 列到第 n 列的元素）。 同样将旋转应用到 \( \tilde{b} \) 的第 j 个和第 i 个分量上。 回代求解 ： 提取 \( \tilde{R} \)（即 R 的前 n 行前 n 列部分）和 \( \tilde{c} \)（即 \tilde{b} 的前 n 个分量）。 调用回代法求解 \( \tilde{R} x = \tilde{c} \)，得到解向量 \( x \)。 为什么选择Givens旋转？ 数值稳定性 ：与Householder反射一样，Givens旋转基于正交变换，非常稳定，不会放大误差。 适用于稀疏/结构化矩阵 ：这是其主要优势。对于带状矩阵、Hessenberg矩阵或一般稀疏矩阵，Householder反射会破坏稀疏性（产生“填入”），而Givens旋转可以精确地只消去需要为零的元素，并能利用现有的零元结构，保持矩阵的稀疏性，从而大幅节省存储和计算量。 并行性 ：多个不相交的Givens旋转（例如，消去奇数行和消去偶数行的操作）可以同时进行，适合并行计算。 通过以上五步，我们完成了基于Givens旋转的QR分解及其在求解最小二乘问题中的应用全流程。这种方法在需要高数值稳定性且矩阵具有特殊稀疏结构时，是首选的高效算法。