基于Householder反射的QR分解算法在求解最小二乘问题中的应用
1. 题目描述
本题目探讨基于Householder反射的QR分解算法如何应用于求解最小二乘问题。
给定一个矩阵 \(A \in \mathbb{R}^{m \times n}\)（通常 \(m \ge n\)）和一个向量 \(b \in \mathbb{R}^{m}\)，最小二乘问题可表述为：
\[\min_{x \in \mathbb{R}^{n}} \|Ax - b\|_2^2 \]
其中 \(\| \cdot \|_2\) 表示向量的2-范数。该问题的核心是找到 \(x\) 使得残差向量的平方和最小。
核心思想：利用Householder反射对矩阵 \(A\) 进行QR分解，将 \(A\) 化为上三角矩阵 \(R\)，同时相应地变换右端项 \(b\)，从而将原问题转化为一个易于求解的上三角线性方程组。
2. 解题过程
2.1 背景知识
- QR分解：将矩阵 \(A\) 分解为正交矩阵 \(Q\) 和上三角矩阵 \(R\) 的乘积，即 \(A = QR\)。
- Householder反射：是一种构造正交矩阵（反射变换）的方法，通过对向量进行镜像反射，可将向量的部分分量化为零。每个Householder反射矩阵 \(H\) 满足 \(H = I - 2uu^T / (u^Tu)\)，其中 \(u\) 是反射面的法向量。
2.2 算法步骤
假设 \(A\) 是列满秩的（即 \(\text{rank}(A) = n\)），算法分为三个阶段：
阶段一：对矩阵 \(A\) 进行QR分解
目标：构造一系列Householder反射矩阵 \(H_1, H_2, \dots, H_n\)，使得
\[H_n H_{n-1} \cdots H_1 A = \begin{bmatrix} R \\ 0 \end{bmatrix} \]
其中 \(R \in \mathbb{R}^{n \times n}\) 是上三角矩阵，\(0 \in \mathbb{R}^{(m-n) \times n}\) 是零矩阵。
每一步的详细操作（以第 \(k\) 步为例，\(1 \le k \le n\)）：
- 从当前矩阵的第 \(k\) 列中，取第 \(k\) 行到第 \(m\) 行的子向量，记为 \(v_k = A[k:m, k]\)。
- 计算反射向量 \(u_k\)：
- 计算 \(v_k\) 的2-范数：\(\sigma_k = \|v_k\|_2\)。
- 令 \(u_k = v_k + \text{sign}(v_{k1}) \sigma_k e_1\)，其中 \(e_1 = [1, 0, \dots, 0]^T\) 是单位向量，\(v_{k1}\) 是 \(v_k\) 的第一个元素，\(\text{sign}\) 是符号函数（若 \(v_{k1}=0\)，则取 +1）。
- 归一化：\(u_k = u_k / \|u_k\|_2\)。
- 构造Householder反射矩阵 \(H_k = I - 2u_k u_k^T\)。
- 对 \(A\) 的第 \(k\) 列及右侧所有列（即子矩阵 \(A[k:m, k:n]\)）应用反射：
\[ A[k:m, k:n] \leftarrow (I - 2u_k u_k^T) \cdot A[k:m, k:n] \]
由于反射的性质，这一步后 \(A[k+1:m, k]\) 的元素全部变为零。
5. 同时，为了后续计算方便，通常将 \(u_k\) 存储在 \(A[k:m, k]\) 中（覆盖原列），并将 \(A[k,k]\) 替换为 \(-\text{sign}(v_{k1}) \sigma_k\)。
最终，矩阵 \(A\) 的上三角部分（包括对角线）存储了 \(R\)，而下三角部分（严格下三角）存储了各步的反射向量 \(u_k\)（不包含归一化因子）。
阶段二：变换右端项 \(b\)
同样的Householder反射也应用于 \(b\)：
\[\hat{b} = H_n H_{n-1} \cdots H_1 b \]
由于 \(H_k\) 是正交矩阵，变换不改变2-范数。计算时不必显式构造 \(H_k\)，而是用存储的反射向量 \(u_k\) 逐次更新 \(b\)：
对 \(k = 1\) 到 \(n\)：
\[b[k:m] \leftarrow b[k:m] - 2u_k (u_k^T b[k:m]) \]
最后得到变换后的向量：
\[\hat{b} = \begin{bmatrix} c \\ d \end{bmatrix}, \quad c \in \mathbb{R}^{n}, \quad d \in \mathbb{R}^{m-n} \]
阶段三：回代求解
由于正交变换保持范数，原最小二乘问题等价于：
\[\min_x \left\| \begin{bmatrix} R \\ 0 \end{bmatrix} x - \begin{bmatrix} c \\ d \end{bmatrix} \right\|_2^2 = \min_x \left( \|Rx - c\|_2^2 + \|d\|_2^2 \right) \]
第二项 \(\|d\|_2^2\) 是常数，因此最小化只需令第一项为零，即解上三角方程组：
\[R x = c \]
通过回代法（从第 \(n\) 行到第 \(1\) 行）求解 \(x\)：
\[x_n = c_n / R_{nn}, \quad x_i = \left( c_i - \sum_{j=i+1}^{n} R_{ij} x_j \right) / R_{ii}, \quad i=n-1, \dots, 1 \]
最终得到的 \(x\) 即为最小二乘解。
2.3 算法复杂度与稳定性
- 时间复杂度：QR分解需 \(O(mn^2 - n^3/3)\) 次浮点运算，变换 \(b\) 需 \(O(mn)\) 次，回代需 \(O(n^2)\) 次。整体为 \(O(mn^2)\)。
- 数值稳定性：Householder反射是数值稳定的正交变换，不会放大舍入误差，因此该算法是求解最小二乘问题的标准稳定方法。
3. 总结
基于Householder反射的QR分解算法通过正交变换将最小二乘问题转化为上三角方程组，避免了正规方程 \(A^T A x = A^T b\) 可能出现的条件数平方问题。其核心在于逐步构造反射向量，同时更新矩阵和右端项，最后回代求解。该算法是数值线性代数中求解最小二乘问题的基础且高效的方法。