PRML讀書筆記——線性回歸模型（上）

本章開始學習第一個有監督學習模型——線性回歸模型。"線性"在這里的含義僅限定了模型必須是參數的線性函數。而正如我們接下來要看到的，線性回歸模型可以是輸入變量\(x\)的非線性函數。

書中首先對回歸問題給出了一個簡短的不那么正式的定義：

Given a training data set comprising \(N\) observations \(\{x_n\}\), where \(n = 1, ... , N\), together with corresponding target values \(\{t_n\}\), the goal is to predict the value of \(t\) for a new value of \(x\).

而解決這個問題有兩種思路：一種是直接對\(x \rightarrow t\)的映射\(y(x)\)進行建模，類似判別式模型的建模思路；另一種則從概率的角度出發，首先建模預測變量\(t\)的條件概率分布\(p(t|x)\)，然后通過最小化一個損失函數——對於實數變量通常是均方誤差（squared loss）——給出\(t\)的一個最優估計：\(y(x)=E_t[t|x]\)。（注意區分這里的損失函數和用於參數回歸的損失函數。）

線性基函數模型（Linear Basis Function Models）

一個最簡單的線性回歸模型是將預測變量表示為輸入變量的一組線性組合：

\[y(\mathbf{x}, \mathbf{w}) = w_o + w_1 x_1 + ... + w_D x_D \]

其中，\(\mathbf{x} = (x_1, ..., x_D)^T\)。

這個模型的特點是，它不僅是參數\(w\)，也是輸入變量\(x\)的線性函數，因此也有着極大的局限性。

因此，我們將這個簡單模型擴展為對輸入變量進行非線性變換后的線性組合：

\[y(\mathbf{x}, \mathbf{w}) = w_0 + \sum_{j=1}^{M-1}{w_j \phi_j(\mathbf{x})} \]

這個模型被稱為線性基函數模型。其中，\(\phi_j(\mathbf{x})\)被稱為模型的基函數，是一個\(\mathbb{R}^D \rightarrow \mathbb{R}\)的映射，也可以被理解為是對原始輸入特征\(\mathbf{x}\)的預處理函數。參數\(w_0\)被稱為模型的bias（區別於統計學意義上的bias），代表數據中的"偏移量"。

如果我們定義\(\phi_0(\mathbf{x})=1\)，則可以用一個緊湊的矢量運算來表示這個模型：

\[y(\mathbf{x}, \mathbf{w}) = \sum_{j=0}^{M-1}{w_j \phi_j(\mathbf{x})} = \mathbf{w}^T \mathbf{\phi(x)} \]

顯然，線性基函數模型可以是輸入變量\(\mathbf{x}\)的非線性函數（如果我們選定一組非線性的基函數）。而我們之所以仍稱其為"線性"的模型，是因為它始終是參數\(\mathbf{w}\)的線性函數。

對於單變量的輸入\(x\)，常見的基函數有四種：

1. Power basis：\(\phi_j(x) = x^j\)
2. Gaussian basis：\(\phi_j(x) = \exp\{-\frac{(x-\mu_j)^2}{2s^2}\}\)
3. Sigmoidal basis：\(\phi_j(x) = \sigma(\frac{x-\mu_j}{s})\)，其中：\(\sigma(x) = \frac{1}{1+\exp(-x)}\)
4. Fourier basis：類似於信號處理里的小波變換

最大似然和最小二乘法

我們可以從概率角度出發進一步理解這個模型。

假定目標變量\(t\)服從以\(y(\mathbf{x}, \mathbf{w})\)為均值的高斯分布：

\[t \sim p(t|\mathbf{x}, \mathbf{w}, \beta)=\mathcal{N}(t|y(\mathbf{x}, \mathbf{w}), \beta^{-1}) \]

根據決策理論（decision theory），當我們選擇最小化均方誤差函數：

\[E[L]=\int\int{\{y(x)-t\}^2p(x, t)dxdt} \]

做為優化目標時，得到\(t\)的最優估計量為：

\[y(\mathbf{x})=E[t|\mathbf{x}]=\int{tp(t|\mathbf{x})}dt=y(\mathbf{x}, \mathbf{w}) \]

那么，如何得到函數\(y(\mathbf{x}, \mathbf{w})\)的具體形式呢？換言之，我們該如何計算出\(t\)的條件概率分布參數\(\mathbf{w}\)和\(\beta\)呢？頻率學派的最大似然法給出了一個可行解。

在頻率學派眼里，模型的參數\(\mathbf{w}\)和\(\beta\)是一個確定值，而數據集\(\mathcal{D}\)是一個隨機變量。具體地，我們假定數據集\(\mathcal{D}\)里的每一對樣本\((\mathbf{x}, t)\)都是獨立同分布的隨機變量。而對於一般的有監督學習問題，我們通常不考慮輸入變量\(\mathbf{x}\)的分布。因此，我們可以直接用樣本里目標變量\(t\)的條件概率分布的乘積來表示這個樣本集合\(\mathcal{D}\)的概率分布：

\[p(\mathbf{t}|\mathbf{w}, \beta) = \prod_{n=1}^N{\mathcal{N}(t_n|\mathbf{w}^T\mathbf{\phi(x_n)}, \beta^{-1})} \]

對於一個給定的數據集\(\mathcal{D}\)和對應的目標變量\(\mathbf{t} = (t_1, ..., t_N)^T\)，上面這個概率分布又被稱為這個數據集的似然函數（注意，與以\(\mathbf{t}\)為變量的概率分布函數不同，有着相同表達式的似然函數是以模型參數\(\mathbf{w}\)和\(\beta\)為變量的函數）。通常，我們取其對數形式：

\[\begin{align*} \ln{p(\mathbf{t}|\mathbf{w}, \beta)} & = \sum_{n=1}^N{\ln\mathcal{N}(t_n|\mathbf{w}^T\mathbf{\phi(x_n)}, \beta^{-1})} \\ & = \frac{N}{2}\ln\beta - \frac{N}{2}\ln(2\pi) - \beta E_D(\mathbf{w}) \end{align*}\]

其中，\(E_D(\mathbf{w})=\frac{1}{2}\sum_{n=1}^N{\{t_n-\mathbf{w}^T\mathbf{\phi(x_n)}\}^2}\)又被稱為平方誤差和（sum-of-squares error）。

最大似然法假定對\(\mathbf{w}\)和\(\beta\)的點估計應使得在當前數據集上的似然函數最大。注意到對於高斯分布的似然函數，求解\(\mathbf{w}\)的過程等價於最小化平方誤差和\(E_D(\mathbf{w})\)。

我們可以通過令似然函數對\(\mathbf{w}\)和\(\beta\)的偏導為\(0\)直接得到參數估計的解析解：

\[\begin{align*} \mathbf{w}_{ML} & = (\mathbf{\Phi}^T\mathbf{\Phi})^{-1}\mathbf{\Phi}^T\mathbf{t} \\ \beta_{ML}^{-1} & = \frac{1}{N}\sum_{n=1}^N{\{t_n-\mathbf{w}_{ML}^T\mathbf{\phi(x_n)}\}^2} \end{align*}\]

其中，\(\mathbf{\Phi}\)（\(\mathbf{\Phi}_{nj}=\phi_j(\mathbf{x}_n)\)）被稱為設計矩陣（design matrix），\(\mathbf{\Phi}^\dagger = (\mathbf{\Phi}^T\mathbf{\Phi})^{-1}\mathbf{\Phi}^T\)是矩陣\(\mathbf{\Phi}\)的偽逆（Moore-Penrose pseudo-inverse）。這兩個等式又被稱為是線性回歸模型的正則方程解。

最小二乘法的幾何解釋

對於一個大小為\(N\)的數據集，我們可以將其目標變量\(\mathbf{t}=(t_1,...,t_N)^T\)看作是\(N\)維空間里的一個向量。同時，模型的每一個基函數作用在這個數據集上得到的\(N\)個輸入變換也可以被視為是一個\(N\)維向量。假設我們有\(M\)個基函數，則對應着\(M\)個\(N\)維向量。這些向量構成了\(N\)維空間里的一個\(M\)維子空間\(\mathcal{S}\)。我們對目標變量\(\mathbf{t}\)作出的預估\(\mathbf{y}\)則對應着這個\(M\)維子空間上的一個向量。

最小二乘法告訴我們，我們應選擇使得目標向量\(\mathbf{t}\)和預估向量\(\mathbf{y}\)的（歐式）距離盡可能短的\(\mathbf{y}\)。這等價於將\(\mathbf{y}\)設置為\(\mathbf{t}\)在\(M\)維子空間\(\mathcal{S}\)上的投影向量。這就是最小二乘法的幾何意義。

序列學習算法

在實際應用中，由於對整個數據集進行矩陣運算的代價太大，正則方程解並不可行。事實上，我們通常采用一種序列學習的算法，例如SGD，批量地處理數據並學習到模型的參數。

SGD應用在線性回歸模型時的更新法則如下：

\[\mathbf{w}^{(\tau+1)}=\mathbf{w}^{(\tau)}+\eta(t_n-\mathbf{w}^{(\tau)T}\mathbf{\phi_n})\mathbf{\phi_n} \]

正則化的最小二乘法

通過在誤差函數中引入正則項，正則化提供了一種控制模型過擬合程度的能力。最常見的正則項為L2正則：

\[E_W(\mathbf{w}) = \frac{1}{2}\mathbf{w^Tw} \]

對於線性回歸模型，采用L2正則修正后的誤差函數為：

\[E(\mathbf{w})=\frac{1}{2}\sum_{n=1}^N\{t_n-\mathbf{w}^T\mathbf{\phi(x_n)}\}^2+\frac{\lambda}{2}\mathbf{w^Tw} \]

其中\(\lambda\)為正則系數，用於控制正則項在計算誤差時的重要性，進而調節模型的復雜度和過擬合的程度。

L2正則又被稱為weight decay技術，這是因為如果將上式對\(\mathbf{w}\)求導，我們將得到一個與\(\mathbf{w}\)成正比的衰減項。更一般的正則項具有如下的形式：

\[E_W(\mathbf{w}) = \frac{1}{2}\sum_{j=1}^M|w_j|^q \]

例如，當\(q=1\)時，誤差函數對\(\mathbf{w}\)的梯度存在一個常數衰減項，使得\(\mathbf{w}\)更加稀疏（存在大量\(0\)元素）。該正則項被稱為L1正則，在一些統計學文獻里也被稱為lasso技術。

我們可以用一個帶約束的優化問題來理解正則項對模型參數的影響。此時，正則化項等價於約束條件：\(\sum_{j=1}^M|w_j|^q \leq \eta\)。下圖給出了二維空間里的一個圖例。

正則化項的本質是通過限制模型的有效復雜度使得那些表面上看起來很復雜的模型可以在小數據集上也取得不錯的泛化效果。然而，模型過擬合的問題並沒有在這個框架下得到根本性的解決——我們只不過是將控制模型復雜度的參數從基函數的個數\(M\)轉移到了正則化系數\(\lambda\)上。我們依然需要Validation Set幫助我們找到合適的超參數\(\lambda\)。

Bias-Variance分解

在這一節里，我們將從頻率學派的角度分析模型過擬合的成因。

我們討論的出發點依然是目標變量的條件概率分布\(p(t|\mathbf{x})\)。假設我們知道了這個概率分布的具體形式，我們該如何對一個新的輸入變量\(\mathbf{x}\)給出預估的\(t\)呢？

決策理論告訴我們，通過最小化一個均方損失函數（區分於最大似然法里的均方誤差函數）\(E[L]\)，我們可以對任意輸入\(\mathbf{x}\)給出一個最優的估計值\(y(\mathbf{x})\)：

\[E[L]=\int\int{\{y(\mathbf{x})-t\}^2p(\mathbf{x}, t)d\mathbf{x}dt} \]

進一步地對\(E[L]\)進行分解，我們得到：

\[E[L] = \int{\{y(\mathbf{x})-h(\mathbf{x})\}^2p(\mathbf{x})}d\mathbf{x} + \int{\{h(\mathbf{x})-t\}^2p(\mathbf{x}, t)}d\mathbf{x}dt \]

其中，\(h(\mathbf{x})\)是目標變量\(t\)的條件概率期望，也被稱為回歸函數：

\[h(\mathbf{x})=E[t|\mathbf{x}]=\int{tp(t|\mathbf{x})}dt \]

上式右邊的第二項來源於數據分布內在的噪聲，我們無法改變。因此，對\(E[L]\)的最小化等價於對上式的第一項：\(\{y(\mathbf{x})-h(\mathbf{x})\}^2\)的最小化。顯然，當\(y(\mathbf{x})=h(\mathbf{x})\)時\(E[L]\)達到最小值。

然而，我們並不知道回歸函數\(h(\mathbf{x})\)的具體形式。如果我們有足夠多的數據以及計算資源，理論上，我們可以以任意精度逼近回歸函數\(h(\mathbf{x})\)（非參估計），從而給出\(y(\mathbf{x})\)的一個最佳估計。遺憾的是，在實際中我們往往只有一個有限大小的數據集。

因此，我們只能用一個含參的函數猜想\(y(\mathbf{x}, \mathbf{w})\)對\(h(\mathbf{x})\)進行建模，然后在給定的數據集\(\mathcal{D}\)上用最大似然法找到參數\(\mathbf{w}\)的點估計。最后，我們用這個擬合模型\(y(\mathbf{x}, \mathcal{D})\)近似真實的條件期望\(h(\mathbf{x})\)，作為\(y(\mathbf{x})\)的一個最佳估計。這一切，正是造成模型存在泛化誤差的幕后根源。下面給出更深入的分析。

按照頻率學派的觀點，數據集\(\mathcal{D}\)是一個隨機變量，因此，我們在這個數據集上對\(h(\mathbf{x})\)的近似\(y(\mathbf{x}, \mathcal{D})\)也是一個隨機變量——采樣自同一分布的不同數據集將得到不同參數的回歸函數和回歸誤差。也因此，對於上式中\(E[L]\)分解的第一項：\(\int{\{y(\mathbf{x})-h(\mathbf{x})\}^2p(\mathbf{x})}d\mathbf{x}\)，我們可以用一個期望：\(E_{\mathcal{D}}[\{y(\mathbf{x}, \mathcal{D})-h(\mathbf{x})\}^2]\)替代。

進一步對\(E_{\mathcal{D}}[\{y(\mathbf{x}, \mathcal{D})-h(\mathbf{x})\}^2]\)進行分解，我們得到一個十分重要的公式，也就是本小節的主角——Bias-Variance分解：

\[E_{\mathcal{D}}[\{y(\mathbf{x}, \mathcal{D})-h(\mathbf{x})\}^2] = \{E_{\mathcal{D}}[y(\mathbf{x}, \mathcal{D})]-h(\mathbf{x})\}^2 + E_{\mathcal{D}}[\{y(\mathbf{x}, \mathcal{D})-E_{\mathcal{D}}[y(\mathbf{x}, \mathcal{D})]\}^2] \]

這個等式的第一項是模型的Bias——它衡量了我們的回歸模型\(y(\mathbf{x}, \mathbf{w})\)與真實的回歸函數\(h(\mathbf{x})\)之間的差異；第二項是模型的Variance——它衡量了在不同數據集上模型擬合的方差。

直覺上來看，模型的Bias和Variance均與模型的復雜度有着密切的聯系。顯然，簡單的模型可能會過於偏離真實的\(h(\mathbf{x})\)而使得模型的bias過高；而在同一份數據集上，復雜的模型的擬合往往具有更高的波動性。進一步的實驗證實了這一猜想。我們使用正則系數\(\lambda\)控制回歸模型的有效復雜度，並在抽樣自同一分布的多個數據集上擬合模型，得到的結果如下圖：

圖中自上而下模型的復雜度由低到高；左邊一列繪出了每個數據集上的擬合曲線；右邊一列則匯出了所有擬合曲線的均值（\(E_{\mathcal{D}}[y(\mathbf{x}, \mathcal{D})]\)的無偏估計）；綠色曲線是真實的分布。可見，隨着模型復雜度的提高，模型的Variance也越來越高；而對於較低的模型復雜度，雖然模型有着較低的Variance，但也有着極高的Bias。

我們也可以直接匯出模型Bias和Variance隨着正則系數\(\lambda\)的變化曲線：

可見，模型的Bias和Variance存在一個trade-off——我們無法使得一個模型有着低Bias的同時也有着較低的Variance。在實際的應用中，我們往往需要通過一個驗證集來調節模型的復雜度（例如調節正則化系數\(\lambda\)）從而尋求Bias和Variance的某種平衡。事實上，高Bias正對應着一個欠擬合的模型；而高Variance則對應着一個過擬合的模型。

最后，在上面的實驗中，我們也看到對在不同數據集上擬合的復雜模型取平均，可以得到一個最佳的擬合曲線。這種平均模型的思想在后面提到的Bayesian方法和在神經網絡上風靡的Dropout正則化方法均有所體現。然而在實際應用中，我們很少會對同一個模型在不同數據集上擬合的結果取平均——如果我們有那么多的數據集，一個更自然的想法是將這些數據集合並為一個更大的數據集來訓練一個模型。