最小二乘問題
定義 3.1.1 給定矩陣 \(A\in\mathbb{R}^{m\times n}\) 及向量 \(b\in\mathbb{R}^m\) ，確定 \(x\in\mathbb{R}^n\) ，使得
稱為最小二乘問題，簡稱為 LS 問題，其中 \(r(x)\) 稱為殘向量.
若 \(r\) 線性依賴於 \(x\) ，則稱為線性最小二乘問題；否則稱為非線性最小二乘問題.
最小二乘問題的解 \(x\) 又可稱作線性方程組
的最小二乘解，當 \(m>n\) 時，稱為超定方程組或矛盾方程組；當 \(m<n\) 時，稱為欠定方程組.
設 \(A\in\mathbb{R}^{m\times n}\) ， \(A\) 的值域定義為
容易發現 \(\mathcal{R}(A)\) 是由 \(A\) 的列向量張成的線性空間； \(A\) 的零空間定義為
其維數記為 \(\mathrm{null}(A)\) .
一個子空間 \(\mathcal{S}\subset\mathbb{R}^n\) 的正交補定義為
定理 3.1.1 線性方程組 \(Ax = b,\ A\in\mathbb{R}^{m\times n}\) 解存在當且僅當
即 \(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)\) 中的分量，簡單分析得到
顯然 \(Ax=b_1\) 有解，由推論即證.
記最小二乘問題的解集為 \(\chi_{LS}\) ，則顯然它非空，用 \(x_{LS}\) 表示其中 \(2\) 范數最小的解，稱為最小 \(2\) 范數解。下面我們介紹正則化方法：
定理 3.1.4 \(x\in\chi_{LS}\) 當且僅當
證明 對上式做簡單變換 \(A^T(Ax-b) = 0\) ，容易看出，如果 $ x\in\chi_{LS} $ ，則殘向量 \(r\) 必然與 \(A^T\) 正交，等式成立；反之，滿足等式的 \(x\) ，做任意擾動 \(y\) 有
則 \(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\) ，則
其中 \(\kappa_2(A) = \|A\|_2\|A^{\dagger}\|_2,\ A^{\dagger} = \left(A^TA\right)^{-1}A^T\) ；這里 \(A^{\dagger}\) 其實是 \(A\) 的廣義逆，
並且恰為 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\) 為 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\) 滿足
證明 注意到反射性，即 Householder 變換不改變模長，可以確定 \(\alpha = \pm\|x\|_2\) ，因此我們只需要找到方向
這樣 \(\omega\) 與 \(x-\alpha e_1\) 方向相同，只需
可以驗證這是成立的.
為了讓 \(\alpha>0\) ，我們取 \(v = x - \|x\|_2e_1\) ，但是如果 \(x\) 非常接近 \(\alpha e_1\) ，則可能導致巨大的舍入誤差，因此在 \(x_1>0\) 時我們計算
避免誤差；其次，注意到
因此我們不需要求出 \(\omega\) ，只需求出 \(\beta\) 和 \(v\) 即可；在實際計算時，我們將 \(v\) 規格化為第一個分量是 \(1\) 的向量，這樣恰好將 \(v\) 存放在 \(x\) 的后 \(n-1\) 個化為 \(0\) 的分量中.
代碼實現
double householder(vector<double> &x, vector<double> &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\) ，並比較變換前后的模有
則 \(\alpha\) 已知，代入得到 \(Hx\) ，由
得到 \(\omega\) 即可.
Givens 變換
有時我們希望將向量中一個分量化為 \(0\) ，可以考慮變換
這種變換只改變兩個分量，只考慮其作用的二維分量 \(i,k\) 則有
我們希望有 \(y_k = 0\) ，因此可取
則有
在幾何上，變換 \(G\) 是對 \(x\) 在 \((i,k)\) 坐標平面順時針旋轉 \(\theta\) 角，事實上有
所以稱 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 分解
其中 \(Q\in\mathbb{R}^{m\times m}\) 為正交陣， \(R\in\mathbb{R}^{n\times n}\) 是具有非負對角元的上三角陣；若 \(m=n\) 且 \(A\) 可逆，該分解唯一.
證明 應用 Householder 變換和對列 \(n\) 的數學歸納法，先對第一列做 Householder 變換有
其中 \(A_1\) 是 \(n-1\) 的矩列陣，由歸納假設，存在分解
只需令
故分解完成；如果 \(m=n\) 且 \(A\) 可逆，設有 \(A = Q_1R_1 = Q_2R_2\) ，則 \(Q_2^TQ_1 = R_2R_1^{-1}\) ，對比得到右邊是正交且對角元為正的上三角陣，則只能是單位陣，唯一性即證.
利用 QR 分解，我們得到正交變換法
這樣只需求解上三角陣方程 \(Rx = c_1\) ；通過 Householder 變換進行 QR 分解，注意到
只需每次對右下角的矩陣的第一列做變換 \(H_k\) ，然后補全階數得到 \(\widetilde{H}_k\) 即可.
完成分解后，我們不需要存放整個 \(Q\) 矩陣，而是保存每次變換的 \(H_k\) ，而對於 \(H_k\) 又只需要保存 \(v_k,\beta_k\) ，例如
其中每次變換的 \(v_k\) 第一個分量為 \(1\) ，因此可以省略.
代碼實現
vector<double> QR(vector<vector<double> &A)
{
int row = A.size();
int column = A[0].size();
double beta;
vector<double> d(column);
for (int i = 0; i < column; i++)
{
vector<double> x(row - i);
vector<double> v(row - i);
for (int j = i; j < row; j++)
{
x[j - i] = A[j][i];
}
beta = householder(x, v);
vector<double> 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;
}