最小二乘问题的解法

最小二乘问题

定义 3.1.1 给定矩阵 $A\in\mathbb{R}^{m\times n}$ 及向量 $b\in\mathbb{R}^m$ ，确定 $x\in\mathbb{R}^n$ ，使得

\[\|b-Ax\|_2 = \|r(x)\|_2 = \min_{y\in\mathbb{R}^n}\|Ay-b\|_2 \]

称为最小二乘问题，简称为 LS 问题，其中 $r(x)$ 称为残向量.

若 $r$ 线性依赖于 $x$ ，则称为线性最小二乘问题；否则称为非线性最小二乘问题.

最小二乘问题的解 $x$ 又可称作线性方程组

\[Ax = b,\quad A\in\mathbb{R}^{m\times n} \]

的最小二乘解，当 $m>n$ 时，称为超定方程组或矛盾方程组；当 \(m 时，称为欠定方程组.

设 $A\in\mathbb{R}^{m\times n}$ ， $A$ 的值域定义为

\[\mathcal{R}(A) = \{y\in\mathbb{R}^m:y=Ax,\ x\in\mathbb{R}^n\} \]

容易发现 $\mathcal{R}(A)$ 是由 $A$ 的列向量张成的线性空间； $A$ 的零空间定义为

\[\mathcal{N}(A) = \{x\in\mathbb{R}^n:Ax=0\} \]

其维数记为 $\mathrm{null}(A)$ .

一个子空间 $\mathcal{S}\subset\mathbb{R}^n$ 的正交补定义为

\[\mathcal{S}^{\perp} = \{y\in\mathbb{R}^n:y^Tx=0,\ \forall x\in\mathcal{S}\} \]

定理 3.1.1 线性方程组 $Ax = b,\ A\in\mathbb{R}^{m\times n}$ 解存在当且仅当

\[\mathrm{rank}(A) = \mathrm{rank}([A,b]) \]

即 $A$ 与其增广矩阵的秩相同.

证明这是非常显然的，如果存在解，那么 $b$ 一定在 $A$ 列向量张成的线性空间中.

定理 3.1.2 若上述解存在，且 $x$ 是一个给定的解，则全部解的集合是 $x+\mathcal{N}(A)$ .

推论 3.1.1 上述解存在唯一当且仅当 $\mathcal{N}(A) = \{0\}$ .

定理 3.1.3 线性最小二乘问题的解总是存在的，并且解唯一当且仅当 $\mathcal{N}(A)=\{0\}$ .

证明将 $b$ 分为 $b_1+b_2$ ，其中 $b_1$ 是在 $\mathcal{R}(A)$ 中的分量，简单分析得到

\[\|r(x)\|_2^2 = \|b-Ax\|_2^2 = \|b_1-Ax\|_2^2 + \|b_2\|_2^2 \]

显然 $Ax=b_1$ 有解，由推论即证.

记最小二乘问题的解集为 $\chi_{LS}$ ，则显然它非空，用 $x_{LS}$ 表示其中 $2$ 范数最小的解，称为最小 $2$ 范数解。下面我们介绍正则化方法：

定理 3.1.4 $x\in\chi_{LS}$ 当且仅当

\[A^TAx = A^Tb \]

证明对上式做简单变换 $A^T(Ax-b) = 0$ ，容易看出，如果 $ x\in\chi_{LS} $ ，则残向量 $r$ 必然与 $A^T$ 正交，等式成立；反之，满足等式的 $x$ ，做任意扰动 $y$ 有

\[\|b-A(x+y)\|_2^2 = \|b-Ax\|_2^2 - 2y^TA^T(b-Ax) + \|Ay\|_2^2 = \|b-Ax\|_2^2 + \|Ay\|_2^2 \ge \|b-Ax\|_2^2 \]

则 $x$ 处就是最小值；事实上， $(b-Ax)^T(b-Ax) = b^Tb - 2x^TA^Tb + x^TA^TAx$ ，该函数有梯度 $A^TAx-A^Tb=0$ ，并且二阶导 $A^TA>0$ ，因此在此处取得极小值.

定理 3.1.5 考虑 $b$ 的扰动对 $x$ 的影响，设 $b_1$ 和 $\widetilde{b}_1$ 分别是 $b$ 和 $\widetilde{b}$ 在 $\mathcal{R}(A)$ 上的正交投影，若 $b_1\neq0$ ，则

\[\dfrac{\|\delta x\|_2}{\|x\|_2} \le \kappa_2(A)\dfrac{\|b_1-\widetilde{b}_1\|_2}{\|b_1\|_2} \]

其中 $\kappa_2(A) = \|A\|_2\|A^{\dagger}\|_2,\ A^{\dagger} = \left(A^TA\right)^{-1}A^T$ ；这里 $A^{\dagger}$ 其实是 $A$ 的广义逆，

\[A^TAx = A^Tb\ \Rightarrow\ x = \left(A^TA\right)^{-1}A^Tb \]

并且恰为 Moore-Penrose 广义逆.

证明利用 $\delta x = A^{\dagger}(b-\widetilde{b}) = A^{\dagger}(b_1-\widetilde{b}_1),\ b_1 = Ax $ ，取范数后由不等式关系即证.

定理 3.1.6 设 $A$ 的列向量线性无关，则 $\kappa_2(A)^2 = \kappa_2(A^TA)$ .

初等正交变换

Householder 变换

定义 3.2.1 设 $\omega\in\mathbb{R}^{n},\ \|\omega\|_2=1$ ，定义

\[H = I - 2\omega\omega^T \]

则称 $H$ 为 Householder 变换，也叫做初等反射矩阵或镜像变换.

定理 3.2.1 设 $H$ 是 Householder 变换，则有性质

对称性： $H^T = H$
正交性： $H^TH = I$
对合性： $H^2 = I$
反射性：对任意的 $x\in\mathbb{R}^n$ ， $Hx$ 是 $x$ 关于 $\omega$ 的垂直超平面 $\mathrm{span}\{\omega\}^{\perp}$ 的镜像反射

定理 3.2.2 设 $n\neq x\in\mathbb{R}^n$ ，则可构造单位向量 $\omega\in\mathbb{R}^n$ ，使得 Householder 变换 $H$ 满足

\[Hx = \alpha e_1,\quad \alpha = \pm\|x\|_2 \]

证明注意到反射性，即 Householder 变换不改变模长，可以确定 $\alpha = \pm\|x\|_2$ ，因此我们只需要找到方向

\[Hx = \left(I-2\omega\omega^T\right)x = x - 2\left(\omega^Tx\right)\omega,\quad 2\left(\omega^Tx\right)\omega = x - Hx = x - \alpha e_1 \]

这样 $\omega$ 与 $x-\alpha e_1$ 方向相同，只需

\[\omega = \dfrac{x-\alpha e_1}{\|x-\alpha e_1\|_2} \]

可以验证这是成立的.

为了让 $\alpha>0$ ，我们取 $v = x - \|x\|_2e_1$ ，但是如果 $x$ 非常接近 $\alpha e_1$ ，则可能导致巨大的舍入误差，因此在 $x_1>0$ 时我们计算

\[v_1 = x_1 - \|x\|_2 = \dfrac{x_1^2-\|x\|_2^2}{x_1+\|x\|_2} = \dfrac{-(x_2^2+\cdots+x_n^2)}{x_1+\|x\|_2} \]

避免误差；其次，注意到

\[H = I - 2\omega\omega^T = I - \dfrac{2}{v^Tv}vv^T = I - \beta vv^T,\quad \beta = \dfrac{2}{v^Tv} \]

因此我们不需要求出 $\omega$ ，只需求出 $\beta$ 和 $v$ 即可；在实际计算时，我们将 $v$ 规格化为第一个分量是 $1$ 的向量，这样恰好将 $v$ 存放在 $x$ 的后 $n-1$ 个化为 $0$ 的分量中.

代码实现

double householder(vector &x, vector &v)
{
	double model = inf_norm(x); // 存放 x 的无穷范数

	// protect x
	if (model == 0)
	{
		return 0;
	}

	x = x / model;	// 防止溢出
	double sumOfx = x * x;
	double sumOfx1 = sumOfx - x[0] * x[0];
	v = x, v[0] = 0;

	double beta;
	if (sumOfx1 == 0)
	{
		beta = 0;
	}
	else
	{
		if (x[0] <= 0)
		{
			v[0] = x[0] - sqrt(sumOfx);
		}
		// 如果是正的分量，则用特殊算法减小舍入误差
		else
		{
			v[0] = -1 * sumOfx1 / (x[0] + sqrt(sumOfx));
		}
		// 对 beta 乘 v[0] * v[0] ，抵消规格化为 1 的影响
		beta = 2 * v[0] * v[0] / (sumOfx1 + v[0] * v[0]);
		v = v / v[0];
	}
    return beta;
}

如果要对 $x\in\mathbb{R}^n$ 从 $k+1$ 到 $j$ 分量引入 $0$ ，考虑 $Hx = (x_1,\cdots,x_k,0,\cdots,0,x_{j+1},\cdots,x_n)$ ，这显然是不能保证成立的，因为 $Hx$ 与 $x$ 的模相同，因此我们把一个分量 $x_k$ 替换为未知量 $\alpha$ ，并比较变换前后的模有

\[Hx = (x_1,\cdots,x_{k-1},\alpha,0,\cdots,0,x_{j+1},\cdots,x_n),\quad \alpha^2 = \sum_{i=k}^jx_i^2 \]

则 $\alpha$ 已知，代入得到 $Hx$ ，由

\[\omega = \dfrac{x-Hx}{\|x-Hx\|_2} \]

得到 $\omega$ 即可.

Givens 变换

有时我们希望将向量中一个分量化为 $0$ ，可以考虑变换

\[G = \left( \begin{matrix} 1 & & \vdots & & \vdots\\ & \ddots & \vdots & & \vdots\\ \cdots & \cdots & c & \cdots & s & \cdots & \cdots\\ & & \vdots & & \vdots\\ \cdots & \cdots & -s & \cdots & c & \cdots & \cdots\\ & & \vdots & & \vdots & \ddots\\ & & \vdots & & \vdots & & 1\\ \end{matrix} \right) \]

这种变换只改变两个分量，只考虑其作用的二维分量 $i,k$ 则有

\[\left( \begin{matrix} y_i\\ y_k \end{matrix} \right) = \left( \begin{matrix} c & s\\ -s & c \end{matrix} \right) \left( \begin{matrix} x_i\\ x_k \end{matrix} \right) \]

我们希望有 $y_k = 0$ ，因此可取

\[c = \dfrac{x_i}{\sqrt{x_i^2+x_k^2}},\quad s = \dfrac{x_k}{\sqrt{x_i^2+x_k^2}} \]

则有

\[y_i = \sqrt{x_i^2+x_k^2},\quad y_k = 0 \]

在几何上，变换 $G$ 是对 $x$ 在 $(i,k)$ 坐标平面顺时针旋转 $\theta$ 角，事实上有

\[\cos\theta = \dfrac{x_i}{\sqrt{x_i^2+x_k^2}},\quad \sin\theta = \dfrac{x_k}{\sqrt{x_i^2+x_k^2}} \]

所以称 Givens 变换为平面旋转变换.

代码实现

void givens(double x[], double cs[])
{
	double sqrt_x = sqrt(x[0] * x[0] + x[1] * x[1]);
	cs[0] = x[0] / sqrt_x;
	cs[1] = x[1] / sqrt_x;
}

正交变换法

由于 $2$ 范数具有正交不变性，因此我们可以对最小二乘法进行正交变换来简化求解

定理 3.3.1（QR 分解定理） 设 $A\in\mathbb{R}^{m\times n},\ m\ge n$ ，则 $A$ 有 QR 分解

\[A = Q\left( \begin{matrix} R\\ 0 \end{matrix} \right) \]

其中 $Q\in\mathbb{R}^{m\times m}$ 为正交阵， $R\in\mathbb{R}^{n\times n}$ 是具有非负对角元的上三角阵；若 $m=n$ 且 $A$ 可逆，该分解唯一.

证明应用 Householder 变换和对列 $n$ 的数学归纳法，先对第一列做 Householder 变换有

\[H_1A = \left( \begin{matrix} \|a_1\|_2 & v^T\\ 0 & A_1 \end{matrix} \right) \]

其中 $A_1$ 是 $n-1$ 的矩列阵，由归纳假设，存在分解

\[A_1 = Q_2\left( \begin{matrix} R_2\\ 0 \end{matrix} \right) \]

只需令

\[Q = H_1\left( \begin{matrix} 1 & 0\\ 0 & Q_2 \end{matrix} \right),\quad R = \left( \begin{matrix} \|a_1\|_2 & v^T\\ 0 & R_2\\ 0 & 0 \end{matrix} \right) \]

故分解完成；如果 $m=n$ 且 $A$ 可逆，设有 $A = Q_1R_1 = Q_2R_2$ ，则 $Q_2^TQ_1 = R_2R_1^{-1}$ ，对比得到右边是正交且对角元为正的上三角阵，则只能是单位阵，唯一性即证.

利用 QR 分解，我们得到正交变换法

\[\|Ax-b\|_2^2 = \left\|Q^TAx-Q^Tb\right\|_2^2 = \|Rx-c_1\|_2^2 + \|c_2\|_2^2 \]

这样只需求解上三角阵方程 $Rx = c_1$ ；通过 Householder 变换进行 QR 分解，注意到

\[H_1A = \left( \begin{matrix} \|a_1\|_2 & v^T\\ 0 & A_1\\ \end{matrix} \right),\quad \widetilde{H}_2 = \left( \begin{matrix} 1 & \\ & H_2\\ \end{matrix} \right) \]

只需每次对右下角的矩阵的第一列做变换 $H_k$ ，然后补全阶数得到 $\widetilde{H}_k$ 即可.

?完成分解后，我们不需要存放整个 $Q$ 矩阵，而是保存每次变换的 $H_k$ ，而对于 $H_k$ 又只需要保存 $v_k,\beta_k$ ，例如

\[A = \left( \begin{matrix} r_{11} & r_{12} & r_{13}\\ v_2^{(1)} & r_{22} & r_{23} \\ v_3^{(1)} & v_3^{(2)} & r_{33}\\ v_4^{(1)} & v_4^{(2)} & v_4^{(3)} \end{matrix} \right),\quad d = \left( \begin{matrix} \beta_1\\ \beta_2\\ \beta_3 \end{matrix} \right) \]

其中每次变换的 $v_k$ 第一个分量为 $1$ ，因此可以省略.

代码实现

vector QR(vector &A)
{
    int row = A.size();
    int column = A[0].size();
    double beta;
    
	vector d(column);
	for (int i = 0; i < column; i++)
	{
		vector x(row - i);
		vector v(row - i);
		for (int j = i; j < row; j++)
		{
			x[j - i] = A[j][i];
		}
		beta = householder(x, v);
		vector w(column - i);
		// get w = beta * AT * v
		for (int k = i; k < column; k++)
		{
			double sum = 0;
			for (int p = i; p < row; p++)
			{
				sum += A[p][k] * v[p - i];
			}
			// note: it is w[k-i], not w[k]
			w[k - i] = beta * sum;
		}
		// get HA = A - v * wT
		for (int k = i; k < row; k++)
		{
			for (int p = i; p < column; p++)
			{
				if (p == i && k > i)
				{
					A[k][p] = v[k - i];
				}
				else
				{
					A[k][p] -= v[k - i] * w[p - i];
				}
			}
		}
		d[i] = beta;
	}
	return d;
}

数值代数

最小二乘问题的解法