1. 引言

復習上一篇文章《最小二乘問題詳解1：線性最小二乘》中的知識，對于一個線性問題模型：

那么線性最小二乘問題可以表達為求一組待定值\(\theta\)，使得殘差的平方和最小：

本質上是求解超定線性方程組：

具體的線性最小二乘解是：

2. 求解

2.1 問題

雖然線性最小二乘解已經給出，但是并不意味著在實際的數值計算中就能按照式(1)來進行求解。一個典型的問題就是求逆矩陣：在工程實踐和數值計算中，直接求解逆矩陣通常是一個性能消耗大且可能不精確的操作，應該盡量避免。舉例來說，我們按照大學本科《線性代數》課程中的方法寫程序來求解一個逆矩陣，假設使用伴隨矩陣法：

其中：

• \(\det(A)\) 是矩陣 \(A\) 的行列式。
• \(\text{adj}(A)\) 是 \(A\) 的伴隨矩陣。

為了求解伴隨矩陣 \(\text{adj}(A)\)：

1. 求代數余子式 (Cofactor)：對于矩陣 \(A\) 中的每一個元素 \(a_{ij}\)，計算其代數余子式 \(C_{ij}\)。

   • 代數余子式 \(C_{ij} = (-1)^{i+j} \cdot M_{ij}\)
   • \(M_{ij}\) 是刪去 \(A\) 的第 \(i\) 行和第 \(j\) 列后得到的子矩陣的行列式（稱為余子式）。

2. 構造余子式矩陣：將所有代數余子式 \(C_{ij}\) 按照原來的位置排列，形成一個新矩陣 \(C\)（稱為余子式矩陣）。

3. 轉置：將余子式矩陣 \(C\) 進行轉置，得到的矩陣就是伴隨矩陣 \(\text{adj}(A)\)。

   \[\text{adj}(A) = C^T \]

4. 代入公式：將 \(\det(A)\) 和 \(\text{adj}(A)\) 代入公式 \(A^{-1} = \frac{1}{\det(A)} \cdot \text{adj}(A)\) 即可。

這里我們大概能估算，使用伴隨矩陣法求逆矩陣的理論復雜度是\(O(n!)\)，這是一個階乘級的增長，算法效率非常低。《線性代數》中介紹的另外一種算法高斯消元法也只能達到\(O(n^3)\)，呈指數級增加。其實效率只是一方面的問題，使用計算機求解的另外一個問題是舍入誤差累積：在計算機中，浮點數運算存在固有的舍入誤差；求逆過程涉及大量的除法和減法運算，這些誤差會在計算過程中不斷累積和傳播。總而言之，使用通解求解逆矩陣，可能存在不精確且性能消耗大的問題。

2.2 QR分解

那么不使用逆矩陣怎么辦呢？我們需要注意的是，最小二乘問題的本質是求解，而不是求逆矩陣，因此關鍵是要求解正規方程：

對矩陣\(A\)作QR分解：

其中:

• \(Q_1\in\mathbb R^{m\times n}\) 列正交，滿足\(Q_1^T Q_1 = I_n\)；
• \(R\in\mathbb R^{n\times n}\)是上三角矩陣，如果\(A\)列滿秩，則\(R\)的對角元均非零，可逆。

那么把\(A=Q_1R\)代入正規方程，得到：

左邊整理，因為\(Q_1^T Q_1 = I_n\)：

右邊為

因此正規方程等價于

若\(R\)可逆（即\(A\)滿秩，\(\operatorname{rank}(A)=n\)），則\(R^T\)也可逆。左右兩邊左乘\((R^T)^{-1}\)，得到：

令\(c = Q_1^\top b\)（這是一個長度為\(n\)的向量），于是我們得到一個簡單的上三角線性系統：

這就是QR方法把正規方程化簡得到的核心結果：只需解上三角方程\(R x = Q_1^T b\)。

以上只是對\(A\)列滿秩的情況做了推導，如果\(A\)列滿秩，那么QR分解可以表示為\(x = R^{-1}Q_1^\top b\)；如果\(A\)列不滿秩（\(R\)奇異），需要使用列主元QR方法對\(R^T R x = R^T (Q_1^T b)\)進行求解，或者干脆使用下面要介紹的SVD分解（奇異值分解）法。

2.3 SVD分解

另外一種求解的方法是SVD分解。對任意矩陣\(A\)，存在奇異值分解：

其中:

• \(U\in\mathbb R^{m\times m}\)為正交（列為左奇異向量），

• \(V\in\mathbb R^{n\times n}\)為正交（列為右奇異向量），

• \(\Sigma\in\mathbb R^{m\times n}\)為“對角塊”矩陣，通常寫成

  \[\Sigma=\begin{bmatrix}\Sigma_r & 0 \\ 0 & 0\end{bmatrix} \]

  其中\(\Sigma_r=\operatorname{diag}(\sigma_1,\dots,\sigma_r)\)，\((\sigma_1\ge\sigma_2\ge\cdots\ge\sigma_r>0)\)，\(r=\operatorname{rank}(A)\)。

將SVD代入正規方程，先計算\(A^\top A\)：

注意\(U^T U=I\)。而\(\Sigma^T\Sigma\)是\(n\times n\)的對角塊矩陣，其非零對角元就是\(\sigma_i^2(i=1..r)\)，其余為零。

同樣的，計算\(A^T b\)：

于是正規方程變為：

兩邊左乘\(V^T\)，因為\(V\)正交，\(V^TV=I\)，得到：

把\(y=V^T x\)與\(c=U^T b\)代入，得到更簡單的對角方程：

接下來按奇異值分塊展開對角方程，先寫出\(\Sigma\)相關的形狀：

對\(y\)和\(c\)也做相應分塊：

其中\(y_1,c_1\in\mathbb R^r\)對應非零奇異值，\(y_2,c_2\)對應奇異值為0的部分（維度 \(n-r\)）,代入得到分塊方程:

即等價于兩組方程：

由于\(\Sigma_r\)為對角且可逆，第一式等價于

而\(y_2\)（對應零奇異值的分量）在正規方程中不受約束——這反映了在列秩不足時普通最小二乘解不是唯一的（可以在零空間方向任意加解）。為得到最小范數解（慣常的選擇），取 \(y_2=0\)。

最后回到\(x\)的求解，對于\(y\)有：

將\(c_1\)與\(c=U^\top b\)關系代回：

由于\(y=V^T x\)，于是：

定義\(\Sigma^+\)為將非零奇異值取倒數后轉置得到的偽逆矩陣（對角塊為\(\Sigma_r^{-1}\)，其余為0)，則

這就是 最小二乘的 Moore–Penrose 偽逆解：

• 若\(A\)列滿秩，則為唯一最小二乘解，由于那么\(\Sigma^+=\Sigma^{-1}\)，SVD求解公式退化為常見的\(x = V\Sigma^{-1}U^T b\)
• 若秩虧，它給出 在所有最小二乘解中范數最小的那個（minimum-norm solution）。

2.4 比較

從以上論述可以看到，SVD分解穩定且能處理秩虧的情況，但比QR分解慢，復雜度高，通常\(O(mn^2)\)；QR分解在列滿秩、條件數不是太差時更快；若需要判定秩或求最小范數解，SVD是首選。

3. 補充

在最后補充一些基礎知識，也是筆者很感興趣的一點，那就是為什么一個矩陣A可以進行QR分解或者SVD分解呢？

3.1 QR分解

QR分解其實非常好理解，它的本質其實就是大學本科《線性代數》課程中提到的施密特（Gram–Schmidt）正交化。我們先復習一下施密特正交化相關的知識。

設有一組線性無關向量

我們想把它們變成一組正交（再歸一化后變成標準正交）的向量\(q_1, q_2, \dots, q_n\)。具體步驟如下：

1. 取第一個向量，歸一化：

   \[q_1 = \frac{a_1}{|a_1|} \]

2. 對第 2 個向量，先減去在\(q_1\)上的投影：

   \[\tilde{q}_2 = a_2 - (q_1^T a_2) q_1 \]

   然后歸一化：

   \[q_2 = \frac{\tilde{q}_2}{|\tilde{q}_2|} \]

3. 對第 3 個向量，減去它在前兩個正交向量上的投影：

   \[\tilde{q}_3 = a_3 - (q_1^T a_3) q_1 - (q_2^T a_3) q_2 \]

   然后歸一化：

   \[q_3 = \frac{\tilde{q}_3}{|\tilde{q}_3|} \]

4. 一般地，對第\(j\)個向量：

   \[\tilde{q}_j = a_j - {\sum_{i=1}^{j-1}} (q_i^T a_j) q_i, \quad q_j = \frac{\tilde{q}_j}{|\tilde{q}_j|} \]

這樣得到的\({q_i}\)就是標準正交基，且每個\(q_j\)只用到了前\(j-1\)個。

現在把矩陣\(A\)看成由列向量組成：

把施密特正交化寫成矩陣形式，我們得到一組正交向量：

同時，原向量可以寫成：

其中：

把這些關系拼成矩陣形式：

其中：

• \(R = (r_{ij})\)是\(n \times n\)上三角矩陣，因為第\(j\)列只用到前\(j\)個\(q_i\)。
• \(Q_1\)的列正交，所以\(Q_1^T Q_1 = I\)。

3.2 SVD分解

SVD分解其實也非常有意思，同樣也可以順著《線性代數》中基礎知識來進行推導。首先復習一下特征值和特征向量。對于一個方陣 $ A \in \mathbb{R}^{n \times n} $，如果存在一個非零向量 $ \mathbf{v} \in \mathbb{R}^n $ 和一個實數 $ \lambda $，使得：

那么：

• $ \lambda $ 稱為 特征值（eigenvalue）
• $ \mathbf{v} $ 稱為對應于 $ \lambda $ 的 特征向量（eigenvector）

接下來復習一下什么叫做對角化。如果一個\(n \times n\)矩陣\(A\)可以寫成：

其中：

• $ D $ 是一個對角矩陣（只有對角線上有元素）
• $ P $ 是一個可逆矩陣

我們就說 \(A\) 是 可對角化的。

而且通常：

• $ D $ 的對角元素是 $ A $ 的特征值：$ D = \operatorname{diag}(\lambda_1, \dots, \lambda_n) $
• $ P $ 的列是對應的特征向量

即：

對角化非常重要，因為對角矩陣計算非常簡單，比如計算\(D^k\)只需把對角元各自取\(k\)次方即可。對角化的本質就是把復雜的線性變換，變成旋轉 → 拉伸 → 逆旋轉的過程。注意，不是所有矩陣都能對角化，只有當矩陣有\(n\)個線性無關的特征向量時，才能對角化。但是，所有對稱矩陣（如 $ A^T A $）都可以對角化，而且可以使用正交矩陣對角化。

也就是說，存在正交矩陣\(V\)，使得：

然后根據這個對角化公式，構造\(U\)和\(\Sigma\)，最終得到SVD：

這里具體構造\(U\)和\(\Sigma\)的過程還是有點繁瑣的，這里就不進一步推導了，免得離題太遠。