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

字数 4469 2025-12-24 11:29:15

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

题目描述

在科学计算和数据分析中，最小二乘问题 是求解方程 \(A\mathbf{x} = \mathbf{b}\) 的最小二乘解，其中 \(A \in \mathbb{R}^{m \times n} (m \ge n)\) 是列满秩的矩阵，\(\mathbf{b} \in \mathbb{R}^m\) 是已知向量。该问题等价于求解法方程 \(A^T A \mathbf{x} = A^T \mathbf{b}\)。然而，直接通过法方程求解可能因 \(A^T A\) 的病态导致数值不稳定。QR分解 提供了稳定的求解方法。其中，利用 Givens旋转 实现QR分解，尤其适用于处理稀疏矩阵、带状结构矩阵或需要增量更新分解的情形。本题目讲解基于Givens旋转的QR分解，并应用于求解最小二乘问题，重点分析其数值稳定性和逐步计算过程。

解题步骤

第一步：理解最小二乘问题与QR分解的基本关系

给定最小二乘问题：

\[ \min_{\mathbf{x}} \|A\mathbf{x} - \mathbf{b}\|_2^2 \]

对 \(A\) 进行 QR分解：\(A = QR\)，其中：
- \(Q \in \mathbb{R}^{m \times m}\) 是正交矩阵（\(Q^T Q = I\)），
- \(R \in \mathbb{R}^{m \times n}\) 是上三角矩阵（其前 \(n\) 行构成上三角矩阵 \(R_1 \in \mathbb{R}^{n \times n}\)，后 \(m-n\) 行为零）。
将分解代入目标函数：

\[ \|A\mathbf{x} - \mathbf{b}\|_2^2 = \|QR\mathbf{x} - \mathbf{b}\|_2^2 = \|R\mathbf{x} - Q^T \mathbf{b}\|_2^2 \]

记 \(Q^T \mathbf{b} = \begin{pmatrix} \mathbf{c} \\ \mathbf{d} \end{pmatrix}\)，其中 \(\mathbf{c} \in \mathbb{R}^n\)，\(\mathbf{d} \in \mathbb{R}^{m-n}\)，则最小化问题简化为：

\[ \min_{\mathbf{x}} \|R_1 \mathbf{x} - \mathbf{c}\|_2^2 + \|\mathbf{d}\|_2^2 \]

最小二乘解通过回代求解上三角线性系统得到：

\[ R_1 \mathbf{x} = \mathbf{c} \]

残差平方和为 \(\|\mathbf{d}\|_2^2\)。

第二步：Givens旋转的原理

Givens旋转是一种正交变换，用于在向量中零化特定元素。定义 Givens矩阵 \(G(i,j,\theta) \in \mathbb{R}^{m \times m}\)：

在二维子空间 \((i, j)\) 上，形式为：

\[ G(i,j,\theta) = \begin{pmatrix} 1 & \cdots & 0 & \cdots & 0 & \cdots & 0 \\ \vdots & \ddots & \vdots & & \vdots & & \vdots \\ 0 & \cdots & c & \cdots & s & \cdots & 0 \\ \vdots & & \vdots & \ddots & \vdots & & \vdots \\ 0 & \cdots & -s & \cdots & c & \cdots & 0 \\ \vdots & & \vdots & & \vdots & \ddots & \vdots \\ 0 & \cdots & 0 & \cdots & 0 & \cdots & 1 \end{pmatrix} \]

其中 \(c = \cos\theta\)，\(s = \sin\theta\)，非零元素位于 \((i,i)\)、\((i,j)\)、\((j,i)\)、\((j,j)\) 位置。

作用在向量 \(\mathbf{y} \in \mathbb{R}^m\) 上时，通过选择 \(c, s\) 可使 \(y_j\) 变为零。具体地，若 \(y_i = a\)，\(y_j = b\)，则取：

\[ c = \frac{a}{\sqrt{a^2 + b^2}}, \quad s = \frac{b}{\sqrt{a^2 + b^2}} \]

变换后，新向量在位置 \(i\) 的分量为 \(\sqrt{a^2 + b^2}\)，位置 \(j\) 的分量为 0。

第三步：用Givens旋转对矩阵 \(A\) 进行QR分解

Givens旋转通过逐元素消元将 \(A\) 化为上三角矩阵 \(R\)，同时累积所有旋转得到正交矩阵 \(Q^T\)。

算法步骤如下：

初始化：令 \(R^{(0)} = A\)，\(Q^T = I\)。
对 \(j = 1, 2, \dots, n\)（列索引）执行：
- 对 \(i = m, m-1, \dots, j+1\)（从底部向上扫描行）执行：
  - 若 \(R^{(k)}_{i,j} = 0\) 则跳过（其中 \(k\) 为当前步骤）。
  - 否则，构造Givens旋转 \(G(i-1, i, \theta)\) 使得旋转后元素 \(R^{(k)}_{i,j} = 0\)。即：

\[ \begin{pmatrix} c & s \\ -s & c \end{pmatrix} \begin{pmatrix} R^{(k)}_{i-1,j} \\ R^{(k)}_{i,j} \end{pmatrix} = \begin{pmatrix} \sqrt{R^{(k)}_{i-1,j}^2 + R^{(k)}_{i,j}^2} \\ 0 \end{pmatrix} \]

   - 计算：

\[ r = \sqrt{R^{(k)}_{i-1,j}^2 + R^{(k)}_{i,j}^2}, \quad c = \frac{R^{(k)}_{i-1,j}}{r}, \quad s = \frac{R^{(k)}_{i,j}}{r} \]

   - 将旋转应用于当前矩阵 $R^{(k)}$ 的第 $i-1$ 和 $i$ 行，仅需更新第 $j$ 列至第 $n$ 列的元素，因为前 $j-1$ 列已为零。对 $l = j, \dots, n$：

\[ \begin{aligned} t_1 &= c \cdot R^{(k)}_{i-1,l} + s \cdot R^{(k)}_{i,l} \\ t_2 &= -s \cdot R^{(k)}_{i-1,l} + c \cdot R^{(k)}_{i,l} \\ R^{(k+1)}_{i-1,l} &= t_1, \quad R^{(k+1)}_{i,l} = t_2 \end{aligned} \]

   - 同时更新 $Q^T$：将同一旋转应用于 $Q^T$ 的第 $i-1$ 和 $i$ 行（或记录旋转参数，最后计算 $Q^T \mathbf{b}$ 时直接应用）。

最终得到上三角矩阵 \(R = R^{(final)}\)，且累积变换满足 \(Q^T A = R\)，即 \(A = QR\)。

注意：实际计算中，通常不显式形成完整 \(Q\) 矩阵，而是在处理向量 \(\mathbf{b}\) 时应用相同的旋转序列。

第四步：用Givens QR求解最小二乘问题

对 \(A\) 和向量 \(\mathbf{b}\) 同时进行变换：
- 初始化：令 \(R = A\)，\(\mathbf{z} = \mathbf{b}\)。
- 对 \(A\) 执行第三步的Givens旋转消元，同时对 \(\mathbf{z}\) 应用相同的旋转。即，每当用旋转 \(G(i-1,i,\theta)\) 更新 \(R\) 的行时，更新：

\[ \begin{pmatrix} z_{i-1} \\ z_i \end{pmatrix} \leftarrow \begin{pmatrix} c & s \\ -s & c \end{pmatrix} \begin{pmatrix} z_{i-1} \\ z_i \end{pmatrix} \]

消元完成后，得到上三角矩阵 \(R_1\)（前 \(n\) 行）和变换后的向量 \(\mathbf{z} = \begin{pmatrix} \mathbf{c} \\ \mathbf{d} \end{pmatrix}\)，其中 \(\mathbf{c} \in \mathbb{R}^n\) 对应前 \(n\) 个分量。

回代求解：
- 解上三角系统 \(R_1 \mathbf{x} = \mathbf{c}\)：

\[ x_n = c_n / R_{n,n}, \quad x_i = \left(c_i - \sum_{j=i+1}^{n} R_{i,j} x_j \right) / R_{i,i}, \quad i = n-1, \dots, 1 \]

残差计算：
- 最小二乘残差范数平方为 \(\|\mathbf{d}\|_2^2 = \sum_{i=n+1}^{m} z_i^2\)。

第五步：数值稳定性分析

Givens旋转的稳定性：由于旋转矩阵是正交的，变换不放大误差。计算 \(c, s\) 时采用稳定的数值公式：

\[ \text{if } |b| > |a|: \quad \tau = -a/b, \quad s = 1/\sqrt{1+\tau^2}, \quad c = s\tau \\ \text{else}: \quad \tau = -b/a, \quad c = 1/\sqrt{1+\tau^2}, \quad s = c\tau \]

避免溢出，并通过符号调整确保 \(r \ge 0\)。

相比Householder变换，Givens旋转更适合处理稀疏结构或增量更新问题，但计算量通常更大（约两倍浮点运算量）。
与法方程相比，QR分解不显式计算 \(A^T A\)，避免条件数平方导致的数值不稳定。

总结

基于Givens旋转的QR分解，通过一系列正交旋转逐步将矩阵 \(A\) 化为上三角形式，并在变换中同时处理右端项 \(\mathbf{b}\)，最终通过回代得到最小二乘解。该方法数值稳定性高，适用于结构化或稀疏矩阵，是求解最小二乘问题的可靠算法。

基于Givens旋转的QR分解在求解最小二乘问题中的数值稳定实现题目描述在科学计算和数据分析中，最小二乘问题是求解方程 \(A\mathbf{x} = \mathbf{b}\) 的最小二乘解，其中 \(A \in \mathbb{R}^{m \times n} (m \ge n)\) 是列满秩的矩阵，\(\mathbf{b} \in \mathbb{R}^m\) 是已知向量。该问题等价于求解法方程 \(A^T A \mathbf{x} = A^T \mathbf{b}\)。然而，直接通过法方程求解可能因 \(A^T A\) 的病态导致数值不稳定。 QR分解提供了稳定的求解方法。其中，利用 Givens旋转实现QR分解，尤其适用于处理稀疏矩阵、带状结构矩阵或需要增量更新分解的情形。本题目讲解基于Givens旋转的QR分解，并应用于求解最小二乘问题，重点分析其数值稳定性和逐步计算过程。解题步骤第一步：理解最小二乘问题与QR分解的基本关系给定最小二乘问题： \[ \min_ {\mathbf{x}} \|A\mathbf{x} - \mathbf{b}\|_ 2^2 \] 对 \(A\) 进行 QR分解：\(A = QR\)，其中： \(Q \in \mathbb{R}^{m \times m}\) 是正交矩阵（\(Q^T Q = I\)）， \(R \in \mathbb{R}^{m \times n}\) 是上三角矩阵（其前 \(n\) 行构成上三角矩阵 \(R_ 1 \in \mathbb{R}^{n \times n}\)，后 \(m-n\) 行为零）。将分解代入目标函数： \[ \|A\mathbf{x} - \mathbf{b}\|_ 2^2 = \|QR\mathbf{x} - \mathbf{b}\|_ 2^2 = \|R\mathbf{x} - Q^T \mathbf{b}\|_ 2^2 \] 记 \(Q^T \mathbf{b} = \begin{pmatrix} \mathbf{c} \\ \mathbf{d} \end{pmatrix}\)，其中 \(\mathbf{c} \in \mathbb{R}^n\)，\(\mathbf{d} \in \mathbb{R}^{m-n}\)，则最小化问题简化为： \[ \min_ {\mathbf{x}} \|R_ 1 \mathbf{x} - \mathbf{c}\|_ 2^2 + \|\mathbf{d}\|_ 2^2 \] 最小二乘解通过回代求解上三角线性系统得到： \[ R_ 1 \mathbf{x} = \mathbf{c} \] 残差平方和为 \(\|\mathbf{d}\|_ 2^2\)。第二步：Givens旋转的原理 Givens旋转是一种正交变换，用于在向量中零化特定元素。定义 Givens矩阵 \(G(i,j,\theta) \in \mathbb{R}^{m \times m}\)：在二维子空间 \((i, j)\) 上，形式为： \[ G(i,j,\theta) = \begin{pmatrix} 1 & \cdots & 0 & \cdots & 0 & \cdots & 0 \\ \vdots & \ddots & \vdots & & \vdots & & \vdots \\ 0 & \cdots & c & \cdots & s & \cdots & 0 \\ \vdots & & \vdots & \ddots & \vdots & & \vdots \\ 0 & \cdots & -s & \cdots & c & \cdots & 0 \\ \vdots & & \vdots & & \vdots & \ddots & \vdots \\ 0 & \cdots & 0 & \cdots & 0 & \cdots & 1 \end{pmatrix} \] 其中 \(c = \cos\theta\)，\(s = \sin\theta\)，非零元素位于 \((i,i)\)、\((i,j)\)、\((j,i)\)、\((j,j)\) 位置。作用在向量 \(\mathbf{y} \in \mathbb{R}^m\) 上时，通过选择 \(c, s\) 可使 \(y_ j\) 变为零。具体地，若 \(y_ i = a\)，\(y_ j = b\)，则取： \[ c = \frac{a}{\sqrt{a^2 + b^2}}, \quad s = \frac{b}{\sqrt{a^2 + b^2}} \] 变换后，新向量在位置 \(i\) 的分量为 \(\sqrt{a^2 + b^2}\)，位置 \(j\) 的分量为 0。第三步：用Givens旋转对矩阵 \(A\) 进行QR分解 Givens旋转通过逐元素消元将 \(A\) 化为上三角矩阵 \(R\)，同时累积所有旋转得到正交矩阵 \(Q^T\)。算法步骤如下：初始化：令 \(R^{(0)} = A\)，\(Q^T = I\)。对 \(j = 1, 2, \dots, n\)（列索引）执行：对 \(i = m, m-1, \dots, j+1\)（从底部向上扫描行）执行：若 \(R^{(k)}_ {i,j} = 0\) 则跳过（其中 \(k\) 为当前步骤）。否则，构造Givens旋转 \(G(i-1, i, \theta)\) 使得旋转后元素 \(R^{(k)} {i,j} = 0\)。即： \[ \begin{pmatrix} c & s \\ -s & c \end{pmatrix} \begin{pmatrix} R^{(k)} {i-1,j} \\ R^{(k)} {i,j} \end{pmatrix} = \begin{pmatrix} \sqrt{R^{(k)} {i-1,j}^2 + R^{(k)}_ {i,j}^2} \\ 0 \end{pmatrix} \] 计算： \[ r = \sqrt{R^{(k)} {i-1,j}^2 + R^{(k)} {i,j}^2}, \quad c = \frac{R^{(k)} {i-1,j}}{r}, \quad s = \frac{R^{(k)} {i,j}}{r} \] 将旋转应用于当前矩阵 \(R^{(k)}\) 的第 \(i-1\) 和 \(i\) 行，仅需更新第 \(j\) 列至第 \(n\) 列的元素，因为前 \(j-1\) 列已为零。对 \(l = j, \dots, n\)： \[ \begin{aligned} t_ 1 &= c \cdot R^{(k)} {i-1,l} + s \cdot R^{(k)} {i,l} \\ t_ 2 &= -s \cdot R^{(k)} {i-1,l} + c \cdot R^{(k)} {i,l} \\ R^{(k+1)} {i-1,l} &= t_ 1, \quad R^{(k+1)} {i,l} = t_ 2 \end{aligned} \] 同时更新 \(Q^T\)：将同一旋转应用于 \(Q^T\) 的第 \(i-1\) 和 \(i\) 行（或记录旋转参数，最后计算 \(Q^T \mathbf{b}\) 时直接应用）。最终得到上三角矩阵 \(R = R^{(final)}\)，且累积变换满足 \(Q^T A = R\)，即 \(A = QR\)。注意：实际计算中，通常不显式形成完整 \(Q\) 矩阵，而是在处理向量 \(\mathbf{b}\) 时应用相同的旋转序列。第四步：用Givens QR求解最小二乘问题对 \(A\) 和向量 \(\mathbf{b}\) 同时进行变换：初始化：令 \(R = A\)，\(\mathbf{z} = \mathbf{b}\)。对 \(A\) 执行第三步的Givens旋转消元，同时对 \(\mathbf{z}\) 应用相同的旋转。即，每当用旋转 \(G(i-1,i,\theta)\) 更新 \(R\) 的行时，更新： \[ \begin{pmatrix} z_ {i-1} \\ z_ i \end{pmatrix} \leftarrow \begin{pmatrix} c & s \\ -s & c \end{pmatrix} \begin{pmatrix} z_ {i-1} \\ z_ i \end{pmatrix} \] 消元完成后，得到上三角矩阵 \(R_ 1\)（前 \(n\) 行）和变换后的向量 \(\mathbf{z} = \begin{pmatrix} \mathbf{c} \\ \mathbf{d} \end{pmatrix}\)，其中 \(\mathbf{c} \in \mathbb{R}^n\) 对应前 \(n\) 个分量。回代求解：解上三角系统 \(R_ 1 \mathbf{x} = \mathbf{c}\)： \[ x_ n = c_ n / R_ {n,n}, \quad x_ i = \left(c_ i - \sum_ {j=i+1}^{n} R_ {i,j} x_ j \right) / R_ {i,i}, \quad i = n-1, \dots, 1 \] 残差计算：最小二乘残差范数平方为 \(\|\mathbf{d}\| 2^2 = \sum {i=n+1}^{m} z_ i^2\)。第五步：数值稳定性分析 Givens旋转的稳定性：由于旋转矩阵是正交的，变换不放大误差。计算 \(c, s\) 时采用稳定的数值公式： \[ \text{if } |b| > |a|: \quad \tau = -a/b, \quad s = 1/\sqrt{1+\tau^2}, \quad c = s\tau \\ \text{else}: \quad \tau = -b/a, \quad c = 1/\sqrt{1+\tau^2}, \quad s = c\tau \] 避免溢出，并通过符号调整确保 \(r \ge 0\)。相比Householder变换，Givens旋转更适合处理稀疏结构或增量更新问题，但计算量通常更大（约两倍浮点运算量）。与法方程相比，QR分解不显式计算 \(A^T A\)，避免条件数平方导致的数值不稳定。总结基于Givens旋转的QR分解，通过一系列正交旋转逐步将矩阵 \(A\) 化为上三角形式，并在变换中同时处理右端项 \(\mathbf{b}\)，最终通过回代得到最小二乘解。该方法数值稳定性高，适用于结构化或稀疏矩阵，是求解最小二乘问题的可靠算法。