統計學習：邏輯回歸與交叉熵損失（Pytorch實現）

1 Logistic 分布和對率回歸

監督學習的模型可以是概率模型或非概率模型，由條件概率分布\(P(Y|\bm{X})\)或決 策函數(decision function)\(Y=f(\bm{X})\)表示，隨具體學習方法而定。對具體的輸入\(\bm{x}\)進行相應的輸出預測並得到某個結果時，寫作\(P(y|\bm{x})\)或\(y=f(\bm{x})\)。

我們這里的 Logistic 分類模型是概率模型，模型\(P(Y|\bm{X})\)表示給定隨機向量\(\bm{X}\)下，輸出類別\(Y\)的條件概率分布。這里我們只討論二分類問題，后面我們介紹多層感知機的時候會介紹多分類問題，道理是一樣的。

首先，我們先來介紹 Logistic 分布。

Logistic 的分布函數為：

\[F(x) = P(X<=x) = \frac{1}{1+e^{-(x-u)/\gamma}} \]

該分布函數的圖像為：

該函數以點\((u,1/2)\)中心對稱，在中心附近增長速度較快，在兩端增長速度較慢。

我們后面對於未歸一化的數\(z\)，采用\(\text{sigmoid}\)函數

\[σ(z) = \frac{1}{1 + \text{exp}(−z)} \]

將其“壓”在(0, 1)之間。這樣就可以使\(z\)歸一化，一般我們使歸一化后的值表示置信度（belief），可以理解成概率，但是與概率有着細微的差別，更體現 “屬於××類的把握”的意思。

對於二分類，輸出\(Y=0\)或\(Y=1\)，我們將\(P(Y=0| \bm{x})\)和\(P(Y=1|\bm{x})\)表示如下，也就定義了二項邏輯回歸模型：

\[\begin{aligned} P(Y=1|\bm{x}) &= \sigma(\bm{w}^T\bm{x}+b) \\ &= \frac{1}{1+\text{exp}(-(\bm{w}^T\bm{x}+b))} \end{aligned}\space (1) \\ \begin{aligned} P(Y=0|\bm{x}) &= 1-\sigma(\bm{w}^T\bm{x}+b) \\ &= \frac{1}{1+\text{exp}(\bm{w}^T\bm{x}+b)} \end{aligned} \qquad (2) \]

現在考察 Logistic 回歸模型的特點，一個事件的幾率(odds)是指該事件發生的概率與不發生的概率的比值。如果事件發生的概率是\(p\)，那么該事件的幾率是\(p/(1-p)\)，該事件的對數幾率(log odds，簡稱對率)或 logit 函數是

\[\text{logit}(p) = \text{log}\space p/(1-p) \]

對 Logistic 回歸而言，由式子\((1)\)和式子\((2)\)得到:

\[\text{log}\frac{P(Y = 1|\bm{x})}{1 − P(Y = 1|\bm{x})} = \bm{w}^T\bm{x} + b \]

這玩意在統計學里面稱之為“對率回歸”，其實就是“Logistic regression 名稱”的由來。這里的 Logistic 和“邏輯”沒有任何關系，和對率才是有關系的。 可以看出，輸出\(Y=1\)的對數幾率是由輸入\(\bm{x}\)的線性函數表示的模型，即 Logistic 回歸模型。

2 Logistic 回歸模型的參數估計

我們在《統計推斷：極大似然估計、貝葉斯估計與方差偏差分解》)中說過，從概率模型的視角來看，機器模型學習的過程就是概率分布參數估計的過程，這里就是估計條件概率分布\(P(Y|\bm{X})\)的參數\(\bm{θ}\)。對於給定的訓練數據集\(T=\{\{\bm{x}_1, y_1\},\{\bm{x}_2, y_2\}, \dots, \{\bm{x}_N, y_N\}\}\)，其中，\(\bm{x}_i∈\mathbb{R}^n\)，\(y_i∈\{0, 1\}\)，可以應用極大似然估計法估計模型參數，從而得到Logistic回歸模型。

我們設\(P(Y=1|x)=π(\bm{x})\)，\(P(Y=0|\bm{x})=1-π(\bm{x})\)，很顯然\(P(Y|\bm{x})\)服從一個伯努利分 布，不過這個伯努利分布很復雜，它的參數\(p\)是一個函數\(π(\bm{x})\)。我們這里采用一個抽象的觀念，將函數\(π(\bm{x})\)看做一個整體，這樣\(P(Y|\bm{x})\)服從二項分布，可以方便我們理解。我們的參數\(\bm{θ}=(w, b)\)，不過我們還是采用之前的技巧，將權值向量和輸入向量加以擴充，直接記做\(\bm{θ}=\bm{w}\)。

這樣，對於\(N\)個樣本，定義似然函數：

\[L(\bm{w}|\bm{x}_1,\bm{x}_2,...,\bm{x}_N) = \Pi_{i=1}^N[\pi(\bm{x}_i)]^{y_i}[1-\pi(\bm{x}_i)]^{1-y_i} \]

因為函數比較復雜，我們采用數值優化進行求解。然而這個函數是非凸的， 如果直接用數值迭代算法作用於次函數，可能陷入局部最優解不能保證到收斂到全局最優解。我們為了將其變為負對數似然函數（想一想為什么要取負號），也就是一個凸函數，方便用梯度下降法、牛頓法進行求解：

\[-L(\bm{w}|\bm{x}_1,\bm{x}_2,...,\bm{x}_N) = -\sum_{i=1}^N[y_i\text{log}\space \pi(\bm{x}_i)+(1-y_i)\text{log}(1-\pi(\bm{x}_i))] \]

配湊 \(1/N\)轉換為關於數據經驗分布期望的無偏估計 ，即\(\mathbb{E}_{\bm{x}\sim \hat{p}_{data}}\text{log}\space p_{model}(\bm{x}_i| \bm{\theta})\)而不改變目標函數的優化方向。

\[-\frac{1}{N}\sum_{i=1}^N[y_i\text{log}\space \pi(\bm{x}_i)+(1-y_i)\text{log}(1-\pi(\bm{x}_i))] \]

上面這個式子，就是我們在深度學習里面常用的交叉熵損失函數(cross entropy loss function)的特殊情況。（交叉熵損失函數一般是多分類，這里是二分類）。后面我們會學習，在深度學習領域中，我們用\(f(\cdot)\)表示神經網絡對輸入\(\bm{x}\)加以的 映射（這里沒有激活函數，是線性的映射），\(f(\bm{x})\)就是輸出的條件概率分布\(P(Y|\bm{X}=\bm{x})\)。 設類別總數為\(K\)，\(\bm{x}_i\)為第\(i\)個樣本的特征向量(特征維度為\(D\))，\(f(\bm{x}_i)\)輸出維度也為\(K\)。\(f(\bm{x}_i)_k\)表示\(P(y_k | \bm{x}_i)\)。因為是多分類，\(P(y_k | \bm{x}_i)\)的計算方法就不能通過\(\text{sigmoid}\)函數來得到了，取而代之是更通用的\(\text{softmax}\)函數。我們給定輸入樣本\(\bm{x}\)，設其被預測為第\(k\)類的logit \(z_k=\sum_{j=1}^Dw_{kj}x_{j}\)，其中\(\textbf{W} ∈\mathbb{R}^{K\times D}\)神經網絡的權重矩陣\(D\)為特征向量維度，\(K\)為類別總數。則我們有：

\[f(\bm{x})_k = \text{softmax}(\bm{z})_k=\frac{\text{exp}(z_k)}{\sum_{k^{\prime}}\text{exp}(z_{k^{\prime}})} \]

我們規定第\(i\)個樣本的訓練標簽\(\bm{y}_i\)為\(K\)維one-hot向量。則交叉熵函數表述如下：

\[-\frac{1}{N}\sum_{i=1}^N\sum_{k=1}^K y_{ik}\text{log}f(\bm{x}_i)_k \]

因為\(\bm{y}_i\)為 one-hot 向量，只有一個維度為 1，設這個維度為\(c\)（即表示類別\(c\)），則交叉熵函數可以化簡為：

\[-\frac{1}{N}\sum_{i=1}^N \text{log}f(\bm{x}_i)_c \]

這里\(f(\bm{x}_i)_c=\frac{\text{exp}(z_c)}{\sum_{k}\text{exp}(z_k)}\)。

很明顯，如果\(f(\bm{x}_i)_c\)越大，表示神經網絡預測\(\bm{x}_i\)樣本類別為第\(c\)類的概率越大。這也是目標函數優化的方向——即對於類別為\(c\)的樣本\(\bm{x}_i\)，優化器盡量使樣本\(\bm{x}_i\)被預測為第\(c\)類的logits（也即\(z_c\)）更大，而使其被預測為其它類別的logits（也即\(z_{k\neq c}\)）更小。如果你對這點有疑惑，可以進一步將\(f(\bm{x}_i)_c\)寫成\(\frac{1}{1 + \sum_{k \neq c}\text{exp}(z_{k} - z_c)}\)，這樣看的更明顯一些，這就是糖水不等式\(\frac{b + c}{a + c} > \frac{b}{a} (a>b>0, c>0)\)的原理233。

注 關於sigmoid函數和softmax函數的聯系和區別這里額外強調一下。首先softmax函數是sigmoid函數在多分類情況下的擴展。不過，即使在二分類的情況下，二者也有一些區別。比如給定一個正類樣本\(\bm{x}\)（其訓練標簽\(y=1\)），設\(z = \bm{w}^T\bm{x}+b\in\mathbb{R}\)，sigmoid函數會使得下列概率最大化：

\[P(Y=1|\bm{x}) = σ(z) = \frac{1}{1 + \text{exp}(−{z})} \]

這里本質上是在使得\(z\)最大化。

而對softmax函數而言，設\(\bm{z} = \bm{W}\bm{x}+b\in\mathbb{R}^2\)。則對於正類樣本（其訓練標簽為一個one-hot向量\(\bm{y}\)，滿足\(y_1=1, y_2 = 0\)），softmax函數會使得下列概率最大化：

\[P(Y=1|\bm{x}) = \text{softmax}(\bm{z})_1=\frac{\text{exp}(z_1)}{\text{exp}(z_1) + \text{exp}(z_2)} =\frac{1}{1 + \text{exp}(z_2 - z_1)} \]

這里不但使得\(z_1\)最大化，而且還使得\(z_2\)最小化。
因此我們可以認為對softmax而言，參與其中的單元之間形成了競爭，這和神經科學中的現象很吻合。事實上，softmax的輸出總滿足和為1，所以一個單元的值增加必然對應着其他單元值的減少。這與皮質中相鄰神經元間的側抑制類似。在極端情況下，它變成了贏者通吃（winner-take-all） 的形式（其中一個輸出接近1，其它的接近0）。

如下為使用梯度下降法求解二分類邏輯回歸問題（並采用breast_cancer數據集進行測試，該數據集是二分類數據集）：

from locale import normalize
import numpy as np
import random
import torch
from sklearn.datasets import load_breast_cancer
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score
# batch_size表示單批次用於參數估計的樣本個數
# y_pred大小為(batch_size, 1)
# y大小為(batch_size, )，為類別型變量
def cross_entropy(y_pred, y):
    y = y.reshape(-1, 1)
    return -(torch.mul(y, torch.log(y_pred)) + torch.mul(1-y, torch.log(1-y_pred))).sum()/y.shape[0]

# 前向函數
def logistic_f(X, w): 
    z = torch.matmul(X, w).reshape(-1, 1) 
    return 1/(torch.exp(-z) + 1)

# 之前實現的梯度下降法，做了一些小修改
def gradient_descent(X, w, y, n_iter, eta, loss_func, f):
    for i in range(1, n_iter+1):
        y_pred = f(X, w)
        #print(y_pred)
        loss_v = loss_func(y_pred, y)
        #print(loss_v)
        loss_v.backward() 
        with torch.no_grad(): 
            w -= eta*w.grad
        w.grad.zero_()  
    w_star = w.detach()
    return w_star 

# 迭代次數
n_iter = 600

# 學習率
eta = 10 # 因為我們前面loss除了樣本總數，這里調大些

if __name__ == '__main__':

    X, y = load_breast_cancer(return_X_y=True)

    D = X.shape[1]
    X_train, X_test, y_train, y_test = train_test_split(
        X, y, test_size=0.3, random_state=0)
    n_train, n_test = X_train.shape[0], X_test.shape[0]

    # 初始化權重（含偏置）
    w = torch.tensor(2 * np.random.ranf(size=D + 1) - 1, requires_grad=True)
    print(w)
    
    # 擴充1維
    X_train = torch.tensor(np.concatenate([X_train, np.ones([X_train.shape[0], 1])], axis=1))
    X_train = torch.nn.functional.normalize(X_train) # 先歸一化，否則后面sigmoid部分輸出無窮趨近於0，導致log()得-inf
    X_test = torch.tensor(np.concatenate([X_test, np.ones([X_test.shape[0], 1])], axis=1))
    X_test = torch.nn.functional.normalize(X_test) # 先歸一化，否則后面sigmoid部分輸出無窮趨近於0，導致log()得-inf

 
    y_train, y_test = torch.tensor(y_train), torch.tensor(y_test)

    # 調用梯度下降法對函數進行優化
    # 這里采用單次迭代對所有樣本進行估計，后面我們會介紹小批量法減少時間復雜度
    w_star = gradient_descent(X_train, w, y_train, n_iter, eta, cross_entropy, logistic_f)
    print(w_star)
    
    y_pred = logistic_f(X_test, w_star)
    pred_label = torch.where(y_pred < 0.5, 0, 1)
    acc = accuracy_score(y_test, pred_label)
    print("accuracy: %f" % (acc))

注意輸入數據要先歸一化到(0, 1)，否則后面\(|Xw + b|\)過大，以致\(Xw + b<0\)時使sigmoid函數的輸出無窮趨近於0，最終導致log函數得-inf。我們來看一下算法的輸出結果。初始權重向量、最終得到的權重向量和模型精度如下：

Initial w: tensor([-0.3109, -0.4637, -0.0326,  0.0081, -0.0108, -0.9996, -0.2989,  0.2311,
        -0.9197,  0.7227,  0.2692,  0.4146,  0.6281, -0.8207,  0.3153,  0.9207,
        -0.4418,  0.6558, -0.1684,  0.9044, -0.7847, -0.1003, -0.2369, -0.0608,
        -0.3722,  0.3982, -0.7417, -0.8493,  0.0326,  0.1718, -0.3145],
       dtype=torch.float64, requires_grad=True)
Final w: tensor([  3.0186,   5.7167,  19.8512,  18.5805,   0.0232,  -1.0004,  -0.3366,
          0.2130,  -0.8558,   0.7494,   0.2977,   0.9121,   0.7568,  -6.9429,
          0.3183,   0.9233,  -0.4391,   0.6570,  -0.1606,   0.9056,   2.5354,
          7.5881,  19.1409, -19.0107,  -0.3286,   0.3691,  -0.8191,  -0.8705,
          0.1200,   0.1985,   0.1274], dtype=torch.float64)
Final acc: 0.918129

可見模型收斂正常。

如下是多分類邏輯回歸問題，我們需要將\(\text{sigmoid}\)函數修改為\(\text{softmax}\)函數，損失函數修改為更為通用的交叉熵函數，根據輸出維度變化將權重向量修改為權重矩陣等。我們采用iris數據集進行測試，該數據集是一個3分類數據集。

from locale import normalize
import numpy as np
import random
import torch
from sklearn.datasets import load_iris
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score
# batch_size表示單批次用於參數估計的樣本個數
# n_feature為特征向量維度
# n_class為類別個數
# y_pred大小為(batch_size, n_class)
# y大小為(batch_size, )，為類別型變量
def cross_entropy(y_pred, y):
    # 這里y是創建新的對象，這里將y變為(batch_size. 1)形式
    y = y.reshape(-1, 1)
    return -torch.log(torch.gather(y_pred, 1, y)).sum()/ y_pred.shape[0]

# 前向函數，注意，這里的sigmoid改為多分類的softmax函數
def logistic_f(X, W): 
    z = torch.matmul(X, W)
    return torch.exp(z)/torch.exp(z).sum()

# 之前實現的梯度下降法，做了一些小修改
def gradient_descent(X, W, y, n_iter, eta, loss_func, f):
    for i in range(1, n_iter+1):
        y_pred = f(X, W)
        loss_v = loss_func(y_pred, y)
        loss_v.backward() 
        with torch.no_grad(): 
            W -= eta*W.grad
        W.grad.zero_()  
    W_star = W.detach()
    return W_star 

# 迭代次數
n_iter = 1200

# 學習率
eta = 1 # 因為我們前面loss除了樣本總數，這里也要調大些

if __name__ == '__main__':

    X, y = load_iris(return_X_y=True)

    n_classes = y.max() - y.min() # 分類類別總數
    D = X.shape[1]
    X_train, X_test, y_train, y_test = train_test_split(
        X, y, test_size=0.3, random_state=0)
    n_train, n_test = X_train.shape[0], X_test.shape[0]

    # 初始化權重（含偏置）
    W = torch.tensor(2 * np.random.ranf(size=(D + 1, n_classes)), requires_grad=True)
    print("Initial W: %s" % W)
    
    # 擴充1維
    X_train = torch.tensor(np.concatenate([X_train, np.ones([X_train.shape[0], 1])], axis=1))
    X_train = torch.nn.functional.normalize(X_train) # 先歸一化，否則后面sigmoid輸出無窮趨近於0，導致log()得-inf
    X_test = torch.tensor(np.concatenate([X_test, np.ones([X_test.shape[0], 1])], axis=1))
    X_test = torch.nn.functional.normalize(X_test) # 先歸一化，否則后面sigmoid輸出無窮趨近於0，導致log()得-inf

    print(y_test)
    y_train, y_test = torch.tensor(y_train), torch.tensor(y_test)

    # 調用梯度下降法對函數進行優化
    # 這里采用單次迭代對所有樣本進行估計，后面我們會介紹小批量法減少時間復雜度
    W_star = gradient_descent(X_train, W, y_train, n_iter, eta, cross_entropy, logistic_f)
    print("Final w: %s" % W_star)
    
    y_pred = logistic_f(X_test, W_star)
    pred_label = torch.argmax(y_pred, axis=1)

    acc = accuracy_score(y_test, pred_label)
    print("Final acc: %f" % (acc))

同樣需要注意輸入數據要先歸一化到(0, 1)。我們來看一下算法的輸出結果。初始權重向量、最終得到的權重向量和模型精度如下（同樣，權重最后一行為偏置向量\(\bm{b}\)。因為是多分類，每一個類別維度\(k\)都會對應一個偏置\(b_k\)）：

Initial W: tensor([[0.2088, 0.7352, 1.3713],
        [0.6877, 0.1925, 0.7291],
        [0.4245, 1.7149, 1.8838],
        [0.5134, 1.7090, 0.8319],
        [1.8033, 1.8256, 0.1155]], dtype=torch.float64, requires_grad=True)
Final w: tensor([[ 3.8652,  4.2337, -2.7884],
        [ 5.9122, -3.3024, -3.2518],
        [-7.2260,  3.5287, 10.9562],
        [-3.4692, -1.2400,  7.3598],
        [ 3.5502,  3.5274, -2.3503]], dtype=torch.float64)
Final acc: 0.933333

可見模型收斂正常。

參考

• [1] 李航. 統計學習方法(第2版)[M]. 清華大學出版社, 2019.
• [2] 周志華. 機器學習[M]. 清華大學出版社, 2016.
• [3] Goodfellow I, Bengio Y, Courville A. Deep learning[M]. MIT press, 2016.
• [4] Pytorch官網： torch.nn.functional.normalize 函數
• [5] Scikit Learn 文檔： Toy datasets