Chapter 4
1. 最小二乘和正規方程
1.1 最小二乘的兩種視角
從數值計算視角看最小二乘法
我們在學習數值線性代數時，學習了當方程的解存在時，如何找到\(\textbf{A}\bm{x}=\bm{b}\)的解。但是當解不存在的時候該怎麼辦呢？當方程不一致(無解)時，有可能方程的個數超過未知變數的個數，我們需要找到第二可能好的解，即最小二乘近似。這就是最小二乘法的數值計算視角。
從統計視角看最小二乘法
我們在數值計算中學習過如何找出多項式精確擬合資料點(即插值)，但是如果有大量的資料點，或者採集的資料點具有一定的誤差，在這種情況下我們一般不採用高階多項式進行擬合，而是採用簡單模型近似資料。這就是最小二乘法的統計學視角。
1.2 數值計算的角度推導正規方程
對於線性方程組\(\textbf{A}\bm{x}=\bm{b}\)，\(\textbf{A}\)為\(n\times m\)的矩陣，我們將\(\textbf{A}\)看做\(\textbf{F}^m \rightarrow \textbf{F}^n\)的線性對映。對於對映\(\textbf{A}\)，方程組無解的情況對應於\(\bm{b}∉\text{range}(\textbf{A})\)，即\(\bm{b}\)不在對映\(\textbf{A}\)的值域內，即\(b∉ \{\textbf{A}\bm{x} | \bm{x}∈\textbf{R}^{m}\}\)。(而這等價於\(\bm{b}\) 不在矩陣\(\textbf{A}\)的列向量 的張成空間中)。我們對於
將\(\left( \begin{matrix} 1\\ 1\\ 1 \end{matrix}\right)\)視為列向量\(\bm{a}_1\)，\(\left(\begin{matrix}1\\-1\\1\end{matrix}\right)\)視為列向量\(\bm{a}_2\)，\(\left(\begin{matrix}2\\1\\3\end{matrix}\right)\)視為列向量\(\bm{b}\)，則該線性系統如下圖所示：
雖然解\(\bm{x}\)不存在，但我們仍然想要找到所有候選點構成的平面\(\textbf{A}\bm{x}\)中存在的與\(\bm{b}\)最接近的點\(\overline{\bm{x}}\)(讀作\(\bm{x}\)_bar)。
觀察圖１，我們發現特殊向量\(\textbf{A}\)具有以下特點：餘項\(\bm{b}-\textbf{A}\overline{\bm{x}}\)和平面\(\{\textbf{A}\bm{x} | \bm{x} \in \textbf{R}^n \}\)垂直。即對於\(\textbf{R}^n\)上的所有\(\bm{x}\)，都有向量\(\textbf{A}\bm{x}\)和向量\(\bm{b}-\textbf{A}\overline{\bm{x}}\)正交，即\(\langle \textbf{A}\bm{x},\bm{b}-\textbf{A}\overline{\bm{x}}\rangle \geqslant0\)，寫成矩陣乘法的形式，我們有：
基於矩陣轉置的性質，我們有：
因為對\(\bm{x}\)取任意值該式子都滿足，那我們必須有：
然後我們就得到了正規方程(因為我們這裡用法線推匯出來的，也可稱為法線方程)：
該方程的解\(\overline{\bm{x}}\)是方程組\(\textbf{A}\bm{x} = \bm{b}\)的最小二乘解。如果我們採用數值優化的形式而不是解析的形式，這個最小二乘解也可以看做優化式子\(\underset{x}{\text{argmin}} \frac{1}{2} ||\textbf{A}\bm{x} − \bm{b}||^2\)的解。
1.3 統計的角度推導正規方程
統計裡面有一個常用的模型叫做線性迴歸，將資料矩陣\(\textbf{X}∈\textbf{R}^{m\times n}\)（\(m\) 個樣本， \(n\) 個特徵維度）做為輸入，預測的向量\(\bm{y}∈\textbf{R}^m\)(\(\bm{m}\)個樣本，每個樣本一個標籤)做為輸出。線性迴歸的輸出是輸入的線性函式，令\(\hat{\bm{y}}\)表示模型預測\(\bm{y}\)應該取的值。 我們定義
其中\(\bm{w}∈\textbf{R}^n\)是引數(parameter)向量。
接下來我們採用機器學習的視角。《機器學習：一種人工智慧方法》中Tom.Mitchell 對機器學習的定義是：“假設用\(P\)來評估計算機程式在某任務類\(T\)上的效能，若一個程式通過利用經驗\(E\)在任務\(T\)上獲得了效能改善，則我們說關於\(T\)和\(P\)，該任務對\(E\)進行了學習。”我們將樣本劃分為訓練資料\(\{\bm{x}_{train}^{(i)}, y_{train}^{(i)}\}\)和測試資料\(\{\bm{x}_{test}^{(i)}, y_{test}^{(i)}\}\)採用訓練資料\(\{\bm{x}_{train}^{(i)}, y_{train}^{(i)}\}\)做為經驗，按照一定策略去訓練模型，得到\(w\)向量的值，然後 在測試集上採用同樣策略對模型的效能進行評估。
我們的模型採用的評估策略是評估測試集上的均方誤差(MSE,mean squared error)：\(\text{MSE}_{test} = \frac{1}{m}\sum_{i}(\hat{y}_{test}^{(i)} − y_{test}^{(i)})^2\)。我們通過觀察訓練集\(\{\bm{x}_{train}^{(i)}, y_{train}^{(i)}\}\)獲得經驗， 降低訓練集上的均方誤\(\text{MSE}_{test} = \frac{1}{m}\sum_{i}(\hat{y}_{test}^{(i)} − y_{test}^{(i)})^2\)以改進權重$bm{w}$(注意，機器學習和傳統的最優化不同，最優化單純優化一個目標函式，而 機器學習是降低訓練集上的誤差以期求降低測試集上的誤差，即追求泛化性。而 這往往面臨欠擬合和過擬合的問題，需要後面提到的正則化技巧來解決)，最小化這個均方誤差式子實際上和最 大化對數似然的期望 \(\mathbb{E}_{\bm{x}\sim \hat{p}_{data}}\text{log}_{p_{model}}(y_{train}|\bm{x}_{train}, \bm{θ})\)的得到的引數是相同的。這裡\(\hat{p}_{data}\)是樣本的經驗分佈。然而。我們要最小化的這個均方誤差(或稱為損失函式)是凸函式，也就是說可以通過我們前面學習的數值優化演算法(如牛頓法，梯度下降法，參見《一階和二階優化演算法》)可以找到唯一的最優解。然而，最小二乘問題有個美妙之處就是它可以找到解析解而無需數值迭代。計算其解析解的步驟如下：
我們這裡簡單地採用一階導數條件判別法即可，即令\(∇_\bm{w}\text{MSE}_{train}\)，即：
這樣就得到了統計視角的正規方程，等同於我們前面從數值計算角度推導的正規方程。
值得注意的是，我們除了權重\(\textbf{W}\)之外，經常會加一個偏置\(b\)，這樣\(\bm{y} = \textbf{X}\bm{w} + b\)。
我們想繼續上面推導的正規方程求解，就需要將這個\(b\)併入權重向量\(\bm{w}\)中（即下面的\(w_1\)）。 對於\(\textbf{X}\in \mathbb{R}^{m\times n}\)，我們給\(\textbf{X}\)擴充一列，有：
然後我們繼續按照我們前面的正規方程來求即可。
我們還有最小二乘的多項式擬合形式，允許我們用多項式曲線去擬合資料 點。(關於給定向量\(\bm{w}\)的函式是非線性的，但是關於\(\bm{w}\)的函式是線性的，我們的正規方程解法仍然有效)不過這要求我們先對資料矩陣進行預處理，得到：
然後求解方法與我們之前所述相同。
1.4 矩陣方程的近似解和矩陣的 Moore-Penrose 偽擬(廣義逆)
我們前面得到的基於數值計算視角的正規方程是：\(\textbf{A}^T\textbf{A}\hat{\bm{x}}=\textbf{A}^T\bm{b}\)
基於概率統計視角的正規方程式：\(\textbf{X}^T\textbf{X}\bm{w}=\textbf{X}^T\bm{y}\)
很顯然，這兩個式是等效的，只不過是字母命名不太一樣而已。一般而言， 我們會根據上下文采取字母命名風格。如果上下文是數值計算，那麼我們會採用 第一種寫法，如果上下文是概率統計，我們會採取第二種寫法。
我們現在繼續數值計算的符號命名風格，但是採用統計和機器學習的思想對 我們的理解也是很有幫助的。對於\(\textbf{A}^T\textbf{A}\bm{x}=\textbf{A}^T\bm{b}\)，我們將\(\textbf{A}^T\textbf{A}\)移項，得到：
當線性系統\(\textbf{A}\bm{x}=\bm{b}\)無解時，即因\(\textbf{A}\)不可逆而導致無法按照\(\bm{x}=\textbf{A}^{-1}\bm{b}\)求解時，我們可以將等式右邊的\((\textbf{A}^T\textbf{A})^{-1}\textbf{A}^T\)做為線性系統\(\textbf{A}\bm{x}=\bm{b}\)中\(\textbf{A}\)的偽逆。
由前面的推導知道這個偽逆是由最小二乘問題的正規方程得到的。最小二乘還可以描述成最優化問題的形式進行數值迭代求解，所求解的最優解\(\bm{x}\)為
這裡的\(\frac{1}{2}\)為為了方人工求導配湊的。而在機器學習中，我們為了避免模型 過擬合，我們會對目標函式進行正則化(regularization)。以最優化的視角，即添 加懲罰項以避免權重過大。我們將正則化後的最優化式子寫作：
人工求導可知，這個最小二乘問題對應的正規方程是：
從而得到新的矩陣偽逆形式為：
我們將這個更“全面”的偽逆形式稱為 Moore-Penrose 偽逆（Moore-Penrose pseudoinverse），常常記做\(\textbf{A}^+\)。注意和另外一個在擬牛頓法中做為對Hessian矩陣近似的Shaann偽逆做區別。
偽逆的作用常常用於當\(\textbf{A}\)不可逆（奇異）對矩陣方程\(\textbf{A}\bm{x}=\bm{b}\)的近似求解，在數值計算領域有重要的應用。
2. Ridge迴歸和Lasso
前面我們給出了未對目標函式施以正則化時的最優解\(\bm{x}\)為
因為我們下面要討論機器學習的正則化問題，於是我們採用機器學習視角， 將\(\textbf{A}\)做為資料矩陣\(\textbf{X}\)，\(\bm{x}\)做為優化向量 \(\bm{w}\)，捨去為了方便人工求導配湊的\(\frac{1}{2}\)， 並配湊一個\(\frac{1}{m}\)，\(m\)為樣本個數。則我們就得到均方誤差(MSE, mean squared error)
\(\text{MSE} = \frac{1}{m}||\textbf{X}\bm{w} − \bm{y}||_2^2\)，在機器學習領域一般稱為均方誤差損失函式(MSE Loss)，這樣最優解寫作：
新增正則項，寫作：
其中正則化引數\(λ>0\)。
這裡，因為引入了正則項\(||\bm{w}||_2^2\)，我們稱使用了\(L_2\)正則化。相應的迴歸問題稱為ridge regression(嶺迴歸/極迴歸)
我們定義\(p\)範數，正則項為\(||\bm{w}||_p^p\)的正則化為\(L_p\)正則化。
有類特殊的正則化，叫 Lasso(讀作“拉索”)，它對應的是\(L_1\)正則化：
其中正則化引數\(λ>0\)。
\(L_1\)範數和\(L_2\)範數正則化都有助於降低過擬合的風險，但前者還會帶來一個額外的好處：它比後者更易於獲得”稀疏“(sparse)解，即它求得的\(w\)會有更少的非零分量。（事實上，對\(w\)施加“稀疏約束”（即希望\(w\)的非零分量儘可能少） 最自然的是使用\(L_0\)範數（\(L_0\)範數可以理解為向量所有維度取絕對值，然後取最大的維度）。但\(L_0\)範數不連續，不是凸的，其優化是 NP-hard 的，因此常使用\(L_1\)範數來近似。\(L_1\)範數的優化是凸優化問題，在多項式時間內有數值最優解）
為了理解這一點，我們來看一個直觀的例子：假定任意一個樣本\(\bm{x}\)只有兩個屬性，於是無論是\(L_1\)正則化還是\(L_2\)正則化解出的\(\bm{w}\)都只有兩個分量，即\(w_1\)和\(w_2\)。我們將其做為兩個座標軸，然後在圖中繪製出式\((1)\)和式\((2)\)第一項的”等值線“，即在\((w_1,w_2)\)空間中平方誤差項取值相同的點的連線，在分別繪製出\(L_1\)範數和\(L_2\)範數的等值線，即在\((w_1,w_2)\)空間中\(L_1\)範數取值相同的點的連線，以及\(L_2\)範數取值相同的點的連線，如下圖所示。（圖片來自《西瓜書》253 頁）
式\((1)\)和式\((2)\)的解要在平方誤差項與正則化項之間折中，即出現在圖中平方誤差項等值線與正則化項等值線相交處。如圖可以看出，採用\(L_1\)範數時平方誤差項等值線與正則化項等值線的交點常出現在座標軸上，即\(w_1\)和\(w_2\)為0，而在採用\(L_2\)範數時，兩者的交點出現在某個象限中國，即\(w_1\)和\(w_2\)均非0；換言之，採用\(L_1\)範數比\(L_2\)範數更易於得到稀疏解。
L2正則化還可以有從二次規劃的角度來理解，參見《統計學習3：線性支援向量機(Pytorch實現)》
以下是最小二乘解析解的演算法描述：
例 1 最小二乘的解析解（待優化變數\(\bm{c}=(w,b)\)，\(x\)為資料點的橫座標，不過需要先進行擴充常數的預處理，\(y\)為資料點的縱座標）
import numpy as np
if __name__ == '__main__':
x = np.array(
[
[-1],
[0],
[1],
[2]
]
)
y = [1, 0, 0, -2]
# 對資料矩陣x預處理，即擴充常數1的列
# 此處A一共兩列，最高次數只有1次
A = np.concatenate([np.ones([x.shape[0], 1]), x], axis=1)
AT_A = A.T.dot(A)
AT_y = A.T.dot(y)
c_bar = np.linalg.solve(AT_A, AT_y)
print("最小二乘估計得到的引數:", c_bar)
# 條件數
print("條件數:", np.linalg.cond(AT_A))
這裡因為涉及到拼接操作，所以向量\(\bm{x}\)寫成矩陣(列向量)形式。
解析解和條件數分別為：
最小二乘估計得到的引數: [ 0.2 -0.9]
條件數: 2.6180339887498953
例2 最小二乘的多項式擬合（待優化變數\(\bm{c}=(w_{1},w_{2},b)\)，資料點仍然為二維平面上的點，\(x\)仍然為資料點的橫座標，除了擴充外還需要進行計算次方項的預處理。\(y\)仍然為資料點的縱座標）
import numpy as np
if __name__ == '__main__':
x = np.array(
[
[-1],
[0],
[1],
[2]
]
)
y = [1, 0, 0, -2]
# 對資料矩陣x預處理，從第三列開始依次計算出第二列的次方值(還是擬合平面上的點，不過擴充了)
# 此處A一共三列，最高次數有2次，即拋物線
A = np.concatenate([np.ones([x.shape[0], 1]), x, x**2], axis=1)
AT_A = A.T.dot(A)
AT_y = A.T.dot(y)
c_bar = np.linalg.solve(AT_A, AT_y) # 該API AT_y是一維/二維的都行
print("最小二乘估計得到的引數:", c_bar)
# 條件數
print("條件數:", np.linalg.cond(AT_A))
解析解和A的條件數分別為：
最小二乘估計得到的引數: [ 0.45 -0.65 -0.25]
條件數: 20.608278259652856
最小二乘也可以寫成迭代解的形式，以下是用梯度下降法求解最小二乘。
例2 最小二乘的迭代解
import numpy as np
import torch
def mse_loss(y_pred, y, w, lam, reg):
m = y.shape[0]
return 1/m*torch.square(y-y_pred).sum()
def linear_f(X, w):
return torch.matmul(X, w)
# 之前實現的梯度下降法，做了一些小修改
def gradient_descent(X, w, y, n_iter, eta, lam, reg, loss_func, f):
# 初始化計算圖引數，注意：這裡是建立新物件，非引數引用
w = torch.tensor(w, requires_grad=True)
X = torch.tensor(X)
y = torch.tensor(y)
for i in range(1, n_iter+1):
y_pred = f(X, w)
loss_v = loss_func(y_pred, y, w, lam, reg)
loss_v.backward()
with torch.no_grad():
w.sub_(eta*w.grad)
w.grad.zero_()
w_star = w.detach().numpy()
return w_star
# 本模型按照多分類架構設計
def linear_model(
X, y, n_iter=200, eta=0.001, loss_func=mse_loss, optimizer=gradient_descent):
# 初始化模型引數
# 我們使w和b融合，X後面新增一維
X = np.concatenate([np.ones([X.shape[0], 1]), X], axis=1)
w = np.zeros((X.shape[1],), dtype=np.float64)
# 呼叫梯度下降法對函式進行優化
# 這裡採用單次迭代對所有樣本進行估計，後面我們會介紹小批量法減少時間複雜度
w_star = optimizer(X, w, y, n_iter, eta, mse_loss, linear_f)
return w_star
if __name__ == '__main__':
X = np.array(
[
[-1],
[0],
[1],
[2]
], dtype=np.float32
)
y = np.array([1, 0, 0, -2], dtype=np.float32)
n_iter, eta = 200, 0.1, 0.1
w_star = linear_model(X, y, n_iter, eta, mse_loss, gradient_descent)
print("最小二乘估計得到的引數:", w_star)
可以看出經過足夠的迭代次數(200)次，最終收斂到最優解(等於我們之前求的解析解(0.2, -0.9))：
最小二乘估計得到的引數: [ 0.2 -0.9]
最小二乘也可以用正則化項對權重進行約束，以下是用梯度下降法求解待正則項的最小二乘。
例3 最小二乘的迭代解(帶正則項)
import numpy as np
import torch
def mse_loss(y_pred, y, w, lam, reg):
# 這裡y是建立新的物件，這裡將y變為(batch_size. 1)形式
m = y.shape[0]
if reg == 'L2':
reg_item = lam*torch.square(w).sum()
elif reg == 'L1':
reg_item = lam*torch.norm(w, p=1).sum()
else:
reg_item = 0
return 1/m*torch.square(y-y_pred).sum() + reg_item
def linear_f(X, w):
return torch.matmul(X, w)
# 之前實現的梯度下降法，做了一些小修改
def gradient_descent(X, w, y, n_iter, eta, lam, reg, loss_func, f):
# 初始化計算圖引數，注意：這裡是建立新物件，非引數引用
w = torch.tensor(w, requires_grad=True)
X = torch.tensor(X)
y = torch.tensor(y)
for i in range(1, n_iter+1):
y_pred = f(X, w)
loss_v = loss_func(y_pred, y, w, lam, reg)
loss_v.backward()
with torch.no_grad():
w.sub_(eta*w.grad)
w.grad.zero_()
w_star = w.detach().numpy()
return w_star
# 本模型按照多分類架構設計
def linear_model(
X, y, n_iter=200, eta=0.001, lam=0.01, reg="L2", loss_func=mse_loss, optimizer=gradient_descent):
# 初始化模型引數
# 我們使w和b融合，X後面新增一維
X = np.concatenate([np.ones([X.shape[0], 1]), X], axis=1)
w = np.zeros((X.shape[1],), dtype=np.float64)
# 呼叫梯度下降法對函式進行優化
# 這裡採用單次迭代對所有樣本進行估計，後面我們會介紹小批量法減少時間複雜度
w_star = optimizer(X, w, y, n_iter, eta, lam, reg, mse_loss, linear_f)
return w_star
if __name__ == '__main__':
X = np.array(
[
[-1],
[0],
[1],
[2]
], dtype=np.float32
)
y = np.array([1, 0, 0, -2], dtype=np.float32)
n_iter, eta, lam = 200, 0.1, 0.1
reg = "L1"
w_star = linear_model(X, y, n_iter, eta, lam, reg, mse_loss, gradient_descent)
print("最小二乘估計得到的引數:", w_star)
不採用正則化，即令正則化係數\(\lambda\)為0，迭代次數和學習率不變，可以看到最終估計的結果和上一個例子一樣：
最小二乘估計得到的引數: [ 0.2 -0.9]
採用\(L_2\)範數正則化，正則化係數為0.1，迭代次數和學習率不變，可以看到引數向量整體的範數(可以理解成大小)變小了，可見正則化對引數的範數進行了限制：
最小二乘估計得到的引數: [ 0.14900662 -0.82781457]
注意這裡採用的梯度下降法不是隨機演算法，具有可復現性，演算法執行多次也還是同樣的結果。
引用
- [1] Timothy sauer. 數值分析(第2版)[M].機械工業出版社, 2018.
- [2] Golub, Van Loan. 矩陣計算[M]. 人民郵電出版社, 2020.
- [3] 周志華. 機器學習[M]. 清華大學出版社, 2016.
- [4] Ian Goodfellow，Yoshua Bengio等.深度學習[M].人民郵電出版社, 2017.