數值分析：最小二乘與嶺回歸（Pytorch實現）

算法的完整實現代碼我已經上傳到了GitHub倉庫：NumericalAnalysis-Python（包括其它數值分析算法），感興趣的童鞋可以前往查看。

1. 最小二乘和正規方程

1.1 最小二乘的兩種視角

從數值計算視角看最小二乘法
我們在學習數值線性代數時，學習了當方程的解存在時，如何找到\(\textbf{A}\boldsymbol{x}=\boldsymbol{b}\)的解。但是當解不存在的時候該怎么辦呢？當方程不一致(無解)時，有可能方程的個數超過未知變量的個數，我們需要找到第二可能好的解，即最小二乘近似。這就是最小二乘法的數值計算視角。

從統計視角看最小二乘法
我們在數值計算中學習過如何找出多項式精確擬合數據點(即插值)，但是如果有大量的數據點，或者采集的數據點具有一定的誤差，在這種情況下我們一般不采用高階多項式進行擬合，而是采用簡單模型近似數據。這就是最小二乘法的統計學視角。

1.2 數值計算的角度推導正規方程

對於線性方程組\(\textbf{A}\boldsymbol{x}=\boldsymbol{b}\)，\(\textbf{A}\)為\(n\times m\)的矩陣，我們將\(\textbf{A}\)看做\(\textbf{F}^m \rightarrow \textbf{F}^n\)的線性映射。對於映射\(\textbf{A}\)，方程組無解的情況對應於\(\boldsymbol{b}∉\text{range}(\textbf{A})\)，即\(\boldsymbol{b}\)不在映射\(\textbf{A}\)的值域內，即\(b∉ \{\textbf{A}\boldsymbol{x} | \boldsymbol{x}∈\textbf{R}^{m}\}\)。(而這等價於\(\boldsymbol{b}\) 不在矩陣\(\textbf{A}\)的列向量 的張成空間中)。我們對於

\[\left ( \begin{matrix} 1 & 1 \\ 1 & -1 \\ 1 & 1 \end{matrix} \right) \left( \begin{matrix} x_1 \\ x_2 \end{matrix} \right)= \left( \begin{matrix} 2 \\ 1 \\ 3 \end{matrix} \right) \]

將\(\left( \begin{matrix} 1\\ 1\\ 1 \end{matrix}\right)\)視為列向量\(\boldsymbol{a}_1\)，\(\left(\begin{matrix}1\\-1\\1\end{matrix}\right)\)視為列向量\(\boldsymbol{a}_2\)，\(\left(\begin{matrix}2\\1\\3\end{matrix}\right)\)視為列向量\(\boldsymbol{b}\)，則該線性系統如下圖所示：

雖然解\(\boldsymbol{x}\)不存在，但我們仍然想要找到所有候選點構成的平面\(\textbf{A}\boldsymbol{x}\)中存在的與\(\boldsymbol{b}\)最接近的點\(\overline{\boldsymbol{x}}\)(讀作\(\boldsymbol{x}\)_bar)。

觀察圖１，我們發現特殊向量\(\textbf{A}\)具有以下特點：余項\(\boldsymbol{b}-\textbf{A}\overline{\boldsymbol{x}}\)和平面\(\{\textbf{A}\boldsymbol{x} | \boldsymbol{x} \in \textbf{R}^n \}\)垂直。即對於\(\textbf{R}^n\)上的所有\(\boldsymbol{x}\)，都有向量\(\textbf{A}\boldsymbol{x}\)和向量\(\boldsymbol{b}-\textbf{A}\overline{\boldsymbol{x}}\)正交，即\(\langle \textbf{A}\boldsymbol{x},\boldsymbol{b}-\textbf{A}\overline{\boldsymbol{x}}\rangle \geqslant0\)，寫成矩陣乘法的形式，我們有：

\[(\textbf{A}\boldsymbol{x})^T(\boldsymbol{b}-\textbf{A}\overline{\boldsymbol{x}})=0 \]

基於矩陣轉置的性質，我們有：

\[\boldsymbol{x}^T\textbf{A}^T(\boldsymbol{b}-\textbf{A}\overline{\boldsymbol{x}})=0 \]

因為對\(\boldsymbol{x}\)取任意值該式子都滿足，那我們必須有：

\[\textbf{A}^T(\boldsymbol{b}-\textbf{A}\overline{\boldsymbol{x}})=0 \]

然后我們就得到了正規方程(因為我們這里用法線推導出來的，也可稱為法線方程)：

\[\textbf{A}^T\textbf{A}\overline{\boldsymbol{x}}=\textbf{A}^T\boldsymbol{b} \]

該方程的解\(\overline{\boldsymbol{x}}\)是方程組\(\textbf{A}\boldsymbol{x} = \boldsymbol{b}\)的最小二乘解。如果我們采用數值優化的形式而不是解析的形式，這個最小二乘解也可以看做優化式子\(\underset{x}{\text{argmin}} \frac{1}{2} ||\textbf{A}\boldsymbol{x} − \boldsymbol{b}||^2\)的解。

1.3 統計的角度推導正規方程

統計里面有一個常用的模型叫做線性回歸，將數據矩陣\(\textbf{X}∈\textbf{R}^{m\times n}\)（\(m\) 個樣本， \(n\) 個特征維度）做為輸入，預測的向量\(\boldsymbol{y}∈\textbf{R}^m\)(\(\boldsymbol{m}\)個樣本，每個樣本一個標簽)做為輸出。線性回歸的輸出是輸入的線性函數，令\(\hat{\boldsymbol{y}}\)表示模型預測\(\boldsymbol{y}\)應該取的值。 我們定義

\[ \hat{\boldsymbol{y}}=\textbf{X}\boldsymbol{w} \]

其中\(\boldsymbol{w}∈\textbf{R}^n\)是參數(parameter)向量。

接下來我們采用機器學習的視角。《機器學習：一種人工智能方法》中Tom.Mitchell 對機器學習的定義是：“假設用\(P\)來評估計算機程序在某任務類\(T\)上的性能，若一個程序通過利用經驗\(E\)在任務\(T\)上獲得了性能改善，則我們說關於\(T\)和\(P\)，該任務對\(E\)進行了學習。”我們將樣本划分為訓練數據\(\{\boldsymbol{x}_{train}^{(i)}, y_{train}^{(i)}\}\)和測試數據\(\{\boldsymbol{x}_{test}^{(i)}, y_{test}^{(i)}\}\)采用訓練數據\(\{\boldsymbol{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\)。我們通過觀察訓練集\(\{\boldsymbol{x}_{train}^{(i)}, y_{train}^{(i)}\}\)獲得經驗， 降低訓練集上的均方誤\(\text{MSE}_{test} = \frac{1}{m}\sum_{i}(\hat{y}_{test}^{(i)} − y_{test}^{(i)})^2\)以改進權重\(\boldsymbol{w}\)(注意，機器學習和傳統的最優化不同，最優化單純優化一個目標函數，而 機器學習是降低訓練集上的誤差以期求降低測試集上的誤差，即追求泛化性。而 這往往面臨欠擬合和過擬合的問題，需要后面提到的正則化技巧來解決)，最小化這個均方誤差式子實際上和最 大化對數似然的期望 \(\mathbb{E}_{\boldsymbol{x}\sim \hat{p}_{data}}\text{log}_{p_{model}}(y_{train}|\boldsymbol{x}_{train}, \boldsymbol{θ})\)的得到的參數是相同的。這里\(\hat{p}_{data}\)是樣本的經驗分布。然而。我們要最小化的這個均方誤差(或稱為損失函數)是凸函數，也就是說可以通過我們前面學習的數值優化算法(如牛頓法，梯度下降法，參見我的博客《一階和二階優化算法》)可以找到唯一的最優解。然而，最小二乘問題有個美妙之處就是它可以找到解析解而無需數值迭代。計算其解析解的步驟如下：

我們這里簡單地采用一階導數條件判別法即可，即令\(∇_{\boldsymbol{w}}\text{MSE}_{train}\)，即：

\[\nabla_{\boldsymbol{w}}\frac{1}{m}||\hat{\boldsymbol{y}}_{train}-\boldsymbol{y}_{train}||_2^2 = 0\\ \frac{1}{m}\nabla_{\boldsymbol{w}}||\textbf{X}_{train}\boldsymbol{w}-\boldsymbol{y}_{train}||_2^2=0 \\ \nabla_{\boldsymbol{w}}(\textbf{X}_{train}\boldsymbol{w}-\boldsymbol{y}_{train})^T(\textbf{X}_{train}\boldsymbol{w}-\boldsymbol{y}_{train})= 0 \\ \nabla_{\boldsymbol{w}}(\boldsymbol{w}^T\textbf{X}_{train}^T\textbf{X}_{train}\boldsymbol{w}-2\boldsymbol{w}^T\textbf{X}_{train}^T\boldsymbol{y}_{train}+\boldsymbol{y}_{train}^T\boldsymbol{y}_{train})=0\\ 2\textbf{X}_{train}^T\textbf{X}_{train}\boldsymbol{w} - 2\textbf{X}_{train}^T\boldsymbol{y}_{train} = 0\\ \textbf{X}_{train}^T\textbf{X}_{train}\boldsymbol{w} = \textbf{X}_{train}^T\boldsymbol{y}_{train} \]

這樣就得到了統計視角的正規方程，等同於我們前面從數值計算角度推導的正規方程。

值得注意的是，我們除了權重\(\textbf{W}\)之外，經常會加一個偏置\(b\)，這樣\(\boldsymbol{y} = \textbf{X}\boldsymbol{w} + b\)。

我們想繼續上面推導的正規方程求解，就需要將這個\(b\)並入權重向量\(\boldsymbol{w}\)中（即下面的\(w_1\)）。 對於\(\textbf{X}\in \mathbb{R}^{m\times n}\)，我們給\(\textbf{X}\)擴充一列，有：

\[\left( \begin{matrix} 1 &X_{1,1}&\dots & X_{1,n} \\ 1 &\vdots & \ddots & \vdots \\ 1 &X_{m,1} &\dots & X_{m,n} \end{matrix} \right) \left( \begin{matrix} w_1\\ \vdots\\ w_{n+1} \end{matrix} \right) = \left( \begin{matrix} y_1\\ \vdots\\ y_m \end{matrix} \right) \]

然后我們繼續按照我們前面的正規方程來求即可。
我們還有最小二乘的多項式擬合形式，允許我們用多項式曲線去擬合數據 點。(關於給定向量\(\boldsymbol{w}\)的函數是非線性的，但是關於\(\boldsymbol{w}\)的函數是線性的，我們的正規方程解法仍然有效)不過這要求我們先對數據矩陣進行預處理，得到：

\[\left( \begin{matrix} 1 &X_{1,1}&\dots & X_{1,n}^n \\ 1 &\vdots & \ddots & \vdots \\ 1 &X_{m,1} &\dots & X_{m,n}^n \end{matrix} \right) \]

然后求解方法與我們之前所述相同。

1.4 矩陣方程的近似解和矩陣的 Moore-Penrose 偽擬(廣義逆)

我們前面得到的基於數值計算視角的正規方程是：\(\textbf{A}^T\textbf{A}\hat{\boldsymbol{x}}=\textbf{A}^T\boldsymbol{b}\)

基於概率統計視角的正規方程式：\(\textbf{X}^T\textbf{X}\boldsymbol{w}=\textbf{X}^T\boldsymbol{y}\)

很顯然，這兩個式是等效的，只不過是字母命名不太一樣而已。一般而言， 我們會根據上下文采取字母命名風格。如果上下文是數值計算，那么我們會采用 第一種寫法，如果上下文是概率統計，我們會采取第二種寫法。

我們現在繼續數值計算的符號命名風格，但是采用統計和機器學習的思想對 我們的理解也是很有幫助的。對於\(\textbf{A}^T\textbf{A}\boldsymbol{x}=\textbf{A}^T\boldsymbol{b}\)，我們將\(\textbf{A}^T\textbf{A}\)移項，得到：

\[\overline{\boldsymbol{x}} = (\textbf{A}^T\textbf{A})^{-1}\textbf{A}^T\boldsymbol{b} \]

當線性系統\(\textbf{A}\boldsymbol{x}=\boldsymbol{b}\)無解時，即因\(\textbf{A}\)不可逆而導致無法按照\(\boldsymbol{x}=\textbf{A}^{-1}\boldsymbol{b}\)求解時，我們可以將等式右邊的\((\textbf{A}^T\textbf{A})^{-1}\textbf{A}^T\)做為線性系統\(\textbf{A}\boldsymbol{x}=\boldsymbol{b}\)中\(\textbf{A}\)的偽逆。

由前面的推導知道這個偽逆是由最小二乘問題的正規方程得到的。最小二乘還可以描述成最優化問題的形式進行數值迭代求解，所求解的最優解\(\boldsymbol{x}\)為

\[\underset{\boldsymbol{x}}{\text{argmin}}\frac{1}{2} ||\textbf{A}\boldsymbol{x} - \boldsymbol{b} ||_2^2 \]

這里的\(\frac{1}{2}\)為為了方人工求導配湊的。而在機器學習中，我們為了避免模型 過擬合，我們會對目標函數進行正則化(regularization)。以最優化的視角，即添 加懲罰項以避免權重過大。我們將正則化后的最優化式子寫作：

\[\underset{\boldsymbol{x}}{\text{argmin} }\frac{1}{2}||\textbf{A}\boldsymbol{x} − \boldsymbol{b}||_2^2+c||\boldsymbol{x}||_2^2 \]

人工求導可知，這個最小二乘問題對應的正規方程是：

\[(\textbf{A}^T\textbf{A}+c\textbf{I})\overline{\boldsymbol{x}} = \textbf{A}^T\boldsymbol{b} \]

從而得到新的矩陣偽逆形式為：

\[(\textbf{A}^T\textbf{A}+c\textbf{I})^{-1}\textbf{A}^T \]

我們將這個更“全面”的偽逆形式稱為 Moore-Penrose 偽逆（Moore-Penrose pseudoinverse），常常記做\(\textbf{A}^+\)。注意和另外一個在擬牛頓法中做為對Hessian矩陣近似的Shaann偽逆做區別。

偽逆的作用常常用於當\(\textbf{A}\)不可逆（奇異）對矩陣方程\(\textbf{A}\boldsymbol{x}=\boldsymbol{b}\)的近似求解，在數值計算領域有重要的應用。

2 Ridge回歸和Lasso

前面我們給出了未對目標函數施以正則化時的最優解\(\boldsymbol{x}\)為

\[ \underset{\boldsymbol{x}}{\text{argmin}} \frac{1}{2} ||\textbf{A}\boldsymbol{x} − \boldsymbol{b}||_2^2 \]

因為我們下面要討論機器學習的正則化問題，於是我們采用機器學習視角， 將\(\textbf{A}\)做為數據矩陣\(\textbf{X}\)，\(\boldsymbol{x}\)做為優化向量 \(\boldsymbol{w}\)，舍去為了方便人工求導配湊的\(\frac{1}{2}\)， 並配湊一個\(\frac{1}{m}\)，\(m\)為樣本個數。則我們就得到均方誤差(MSE, mean squared error)

\(\text{MSE} = \frac{1}{m}||\textbf{X}\boldsymbol{w} − \boldsymbol{y}||_2^2\)，在機器學習領域一般稱為均方誤差損失函數(MSE Loss)，這樣最優解寫作：

\[\underset{w}{\text{argmin}} \frac{1}{m}||\textbf{X}\boldsymbol{w} − \boldsymbol{y}||_2^2 \]

添加正則項，寫作：

\[\underset{\boldsymbol{w}}{\text{argmin}}\frac{1}{m}||\textbf{X}\boldsymbol{w} − \boldsymbol{y}||_2^2+λ||\boldsymbol{w}||_2^2 \tag{1} \]

其中正則化參數\(λ>0\)。

這里，因為引入了正則項\(||\boldsymbol{w}||_2^2\)，我們稱使用了\(L_2\)正則化。相應的回歸問題稱為ridge regression(嶺回歸/極回歸)

我們定義\(p\)范數，正則項為\(||\boldsymbol{w}||_p^p\)的正則化為\(L_p\)正則化。

有類特殊的正則化，叫 Lasso(讀作“拉索”)，它對應的是\(L_1\)正則化：

\[\underset{\boldsymbol{x}}{\text{argmin}} \frac{1}{m}||\textbf{X}\boldsymbol{w} − \boldsymbol{y}||_2^2+λ||\boldsymbol{w}||_1^1 \tag{2} \]

其中正則化參數\(λ>0\)。

\(L_1\)范數和\(L_2\)范數正則化都有助於降低過擬合的風險，但前者還會帶來一個額外的好處：它比后者更易於獲得”稀疏“(sparse)解，即它求得的\(w\)會有更少的非零分量。（事實上，對\(w\)施加“稀疏約束”（即希望\(w\)的非零分量盡可能少） 最自然的是使用\(L_0\)范數（\(L_0\)范數可以理解為向量所有維度取絕對值，然后取最大的維度）。但\(L_0\)范數不連續，不是凸的，其優化是 NP-hard 的，因此常使用\(L_1\)范數來近似。\(L_1\)范數的優化是凸優化問題，在多項式時間內有數值最優解）

為了理解這一點，我們來看一個直觀的例子：假定任意一個樣本\(\boldsymbol{x}\)只有兩個屬性，於是無論是\(L_1\)正則化還是\(L_2\)正則化解出的\(\boldsymbol{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　最小二乘的解析解（待優化變量\(\boldsymbol{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))

這里因為涉及到拼接操作，所以向量\(\boldsymbol{x}\)寫成矩陣(列向量)形式。

解析解和條件數分別為：

最小二乘估計得到的參數: [ 0.2 -0.9]
條件數: 2.6180339887498953

例2 最小二乘的多項式擬合（待優化變量\(\boldsymbol{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.