ESL 4: Linear Methods for Classification

4.2 Linear Regression of an Indicator Matrix

對于分類問題屉栓，我們可以把每個分類用 indicator variable 來編碼。

例如，如果有 $K$ 個類型挺尿，我們可以定義：

$\textbf{y} = [y_1, ..., y_K]^T$

當樣本是第 $i$ 類時监署， $y_i = 1$ 辈末，其它值為 0亡驰。

當我們有 N 個樣本時锯七，可以寫成一個 $N \times K$ 矩陣（indicator response matrix）：

$\textbf{Y} = [\textbf{y}_1, ..., \textbf{y}_N]^T$

這個矩陣每一行僅有一個 1源织，其余值為 0翩伪。這個 1 代表該行樣本的類型。

這實際上就是 one-hot 編碼谈息，這種編碼相對于直接使用整數 $1, ..., K$ 來表示的有優(yōu)勢在于 類型間的距離都是一樣的缘屹。

利用第三章的 Linear Regression，我們可以得出系數矩陣 $\hat{\textbf{B}}$ ：

$\hat{\textbf{B}} = \textbf{X}(\textbf{X}^T \textbf{X})^{-1} \textbf{X}^T \textbf{Y}$

$\hat{\textbf{B}}$ 是一個 $(p+1) \times K$ 的矩陣侠仇。 $\textbf{X}$ 是 $N \times (p+1)$ 矩陣轻姿。

得到估計值：

$\hat{\textbf{Y}} = \textbf{X}(\textbf{X}^T \textbf{X})^{-1} \textbf{X}^T \textbf{Y}$

該方法的步驟是：

對類型進行 one-hot 編碼得到 $\textbf{Y}$
將問題視作一個多變量的線性回歸問題，對 one-hot 編碼后的每個 bit 進行擬合得到 $\hat{\textbf{B}}$
在判斷類別時逻炊，選擇 $\textbf{X}\hat{\textbf{B}}$ 每行的最大值所表示的類型

使用 Linear Regression 解決分類問題的最大問題在于互亮，當類別數 $K \leq 3$ 時，可能會出現某個類被掩蓋（mask）的情況余素。其本質原因是 Linear Regression 的“剛性”特質豹休，即它的分界面是不夠靈活的。

舉個簡單的例子桨吊，對于以下三個 1 維正態(tài)分布：

class1: $x \sim N(1, 1)$
class2: $x \sim N(5, 1)$
class3: $x \sim N(9, 1)$

分布如下：

masking1.png

如果我們用 Linear Regression 擬合威根，那我們可以得到 3 組 $\beta_0, \beta_1$ ，分別對應第 $i$ 組视乐，可以用來計算 $y_i$ 的值：

$y_i = \beta_0 + \beta_1 x$

masking2.png

可以看到洛搀， $y_1$ （對應class2）從來不是最大值。也就是說我們的分類結果中只有 class1 和 class3 了炊林，class2 被 mask 了姥卢。可以通過結果驗證：

pd.DataFrame(X*B).T.idxmax().value_counts()
# 2    1505
# 0    1495
# dtype: int64

代碼：

import numpy as np
import pandas as pd

def generate_nd_sample(name, mu_array, sigma, N):
    xx = {}
    for i in range(1, len(mu_array) + 1):
        xi = np.random.normal(mu_array[i-1], sigma, N)
        xx[f"x{i}"] = xi
    xx["name"] = name
    return pd.DataFrame(xx).astype({"name":"category"})

s1 = generate_nd_sample("class1", [1], 1, 1000)
s2 = generate_nd_sample("class2", [5], 1, 1000)
s3 = generate_nd_sample("class3", [9], 1, 1000)

s = pd.concat([s1, s2, s3]).reset_index(drop=True)

# Linear Regression

tmp = s.copy()
tmp.insert(0, "ones", np.ones(s.shape[0]))
X = np.matrix(tmp.drop("name", axis=1).to_numpy().T).T

Y = np.matrix(pd.get_dummies(s.name))

def LR_beta(X, Y):
    # class count
    K = Y.shape[1]
    return (X.T * X).I * X.T * Y

B = LR_beta(X, Y)

# Plot

from bokeh.io import output_notebook
output_notebook()
from bokeh.palettes import viridis
from bokeh.plotting import figure, show
import itertools

def create_histogram_figure(sample_data):
    # sample data must be like:
    # | x | category |
    assert sample_data.shape[1] == 2, "input must be of shape (N, 2)"

    x = sample_data.iloc[:, 0]
    categories = sample_data.iloc[:, 1].unique()

    fig = figure(x_axis_label=x.name, y_axis_label="counts")

    color_gen = itertools.cycle(viridis(len(categories)))
    for (category, color) in zip(categories, color_gen):
        data = sample_data[sample_data.iloc[:, 1] == category]
        counts, bins = np.histogram(data.iloc[:, 0], bins='auto')
        fig.quad(top=counts, bottom=0, left=bins[:-1], right=bins[1:],
                   alpha=0.5, color=color, legend_label=str(category))

    return fig

def create_line_figure(lines, x_start, x_end):
    fig = figure(x_axis_label="x", y_axis_label="y")
    color_gen = itertools.cycle(viridis(lines.shape[0]))
    for i in range(lines.shape[0]):
        b0, b1 = lines[i,0], lines[i, 1]
        fig.line(x=[x_start, x_end], y = [b0 + b1*x_start, b0 + b1*x_end],
                 line_width=2, alpha=0.5, color=next(color_gen),
                 legend_label=f"y{i} = {b0:.6f} + {b1:.6f}x")
    return fig

hist_fig = create_histogram_figure(s)
line_fig = create_line_figure(B.T, s.x1.min(), s.x1.max())
from bokeh import layouts
show(layouts.column(hist_fig, line_fig))

4.3 Linear discriminant analysis

為了獲得最優(yōu)的分類結果，我們需要知道后驗概率（ $X = x$ 時屬于第 $k$ 類的概率）：

$\text{Pr}(G=k | X=x) = \frac{f_k(x) \pi_k}{\sum_{i=1}^K f_i(x) \pi_i}$

因為本質上独榴，我們是在找到一個 k 使得后驗概率最大僧叉，即：

$\begin{align} \hat{k} &= \mathop{\arg \max}_{k} \text{Pr}(G=k | X=x) \\ &= \mathop{\arg \max}_{k} f_k(x) \pi_k \\ &= \mathop{\arg \max}_{k} [\ln(f_k(x)) + \ln(\pi_k)] \end{align}$

這被稱為 判別函數（discriminant function），其中：

$f_i(x)$ 是第 i 類樣本取 x 的概率
$\pi_i$ 是屬于第 i 類的先驗概率

這里的難點在于確定 $f_i(x)$ 棺榔，顯然 $\pi_i$ 的估計是可以通過樣本數據直接得到的瓶堕。

線性判別分析（Linear Discriminant Analysis, LDA）假設變量 X 服從多維高斯分布（X 包含多維）：

$f_k(x) = \frac{1}{(2 \pi)^{p/2} |\mathbf{\Sigma}_k|^{1/2}} e^{-\frac{1}{2}(x - \mu_k)^T \mathbf{\Sigma}_k^{-1} (x - \mu_k)}$

帶入最優(yōu)分類的式子, 逐步去掉與 $k$ 無關的部分：

$\begin{align} \hat{k} &= \mathop{\arg \max}_{k} [\ln(f_k(x)) + \ln(\pi_k)] \\ &= \mathop{\arg \max}_{k} [- \ln((2 \pi)^{p/2} |\mathbf{ \Sigma }_k|^{1/2}) - \frac{1}{2}(x - \mu_k)^T \mathbf{ \Sigma }_k^{-1} (x - \mu_k) + \ln(\pi_k)] \\ &= \mathop{\arg \max}_{k} [- \frac{1}{2} \ln |\mathbf{ \Sigma }_k| - \frac{1}{2}(x - \mu_k)^T \mathbf{\Sigma}_k^{-1} (x - \mu_k) + \ln(\pi_k)] \\ \end{align}$

此時，判別函數為：

$\delta_k(x) = - \frac{1}{2} \ln |\mathbf{ \Sigma }_k| - \frac{1}{2}(x - \mu_k)^T \mathbf{\Sigma}_k^{-1} (x - \mu_k) + \ln(\pi_k)$

是 $x$ 的二次函數症歇。因此稱為二次判別分析(Quadratic Discriminant Analysis, QDA)郎笆。

我們再 假設每個類中變量 X 分布的方差是相等的，則 $\mathbf{\Sigma}$ 也與 $k$ 無關了忘晤，可以進步一化簡判別函數為：

$\delta_k(x) = x^T\mathbf{\Sigma}^{-1}\mu_k - \frac{1}{2}\mu_k^T\mathbf{\Sigma}^{-1} \mu_k + \ln(\pi_k)$

我們可以看出宛蚓，化簡后判別函數對于 $x$ 是線性的。這說明兩個類的分界面（即判別函數相等時）也是線性的设塔。因此叫做線性判別分析(Linear Discriminant Analysis, LDA)凄吏。

實際中，我們可以通過樣本估計高斯分布的參數：

$\hat{\pi}_k = N_k / N$ 闰蛔，即第 k 類的樣本數占總樣本數的比例
$\hat{\mu}_k = \sum_{g_i = k} x_i / N_k$ 痕钢，即第 k 類樣本 X 的平均值
$\hat{\mathbf{\Sigma}} = \sum_{k=1}^K \sum_{g_i = k} (x_i - \hat{\mu}_k)(x_i - \hat{\mu}_k)^T / (N - K)$ ，對協方差矩陣的無偏估計序六，證明在 ESL3 中

有了判別函數的表達式 $\delta_k(x)$ 任连，我們只需要依次帶入 $k = 1, ..., K$ , 當得到的 $\delta_k(x)$ 最大時的 $k$ 即為最佳分類。

4.3.2 Computation of LDA

協方差矩陣 $\mathbf{\Sigma}$ 是一個對稱矩陣例诀，可以進行特征值分解：

$\mathbf{\Sigma} = \mathbf{Q}\mathbf{\Lambda}\mathbf{Q}^T$

其中： $\mathbf{Q}$ 是單位正交矩陣随抠， $\mathbf{\Lambda}$ 是對角陣。帶入判別函數有：

$\begin{align} \delta_k(x) &= x^T\mathbf{\Sigma}^{-1}\mu_k - \frac{1}{2}\mu_k^T\mathbf{\Sigma}^{-1} \mu_k + \ln(\pi_k) \\ &= x^T\mathbf{Q}^{T}\mathbf{\Lambda}^{-1}\mathbf{Q}\mu_k - \frac{1}{2}\mu_k^T\mathbf{Q}^{T}\mathbf{\Lambda}^{-1}\mathbf{Q}\mu_k + \ln(\pi_k) \end{align}$

令：

$x^{*} = \mathbf{\Lambda}^{-\frac{1}{2}}\mathbf{Q}x$

$\mu^{*} = \mathbf{\Lambda}^{-\frac{1}{2}}\mathbf{Q} \mu$

有：

$\delta_k(x^{*}) = x^{* T}\mu_k^{*} - \frac{1}{2}\mu_k^{* T} \mu_k^{*} + \ln(\pi_k)$

當我們判斷某個樣本 $x_1$ 屬于 m 和 n 中的哪一個類時余佃，可以比較其判別函數暮刃，我們判斷它是 m 類如果滿足：

$\delta_m(x_1^*) > \delta_n(x_1^*)$

帶入表達式有：

$x^{*T} (\mu^*_m - \mu^*_n) > \frac{1}{2} (\mu^*_m + \mu^*_n)^T(\mu^*_m - \mu^*_n) - \ln(\pi_m/\pi_n)$

這樣看起來就非常直觀了。 LDA 是將樣本投影在兩個類中心的連線上爆土，并且比較它更靠近哪一邊椭懊，以此決定它屬于哪個類。當然步势，這個過程還考慮了兩個類的先驗概率（ $\ln(\pi_m/\pi_n)$ 項）氧猬。

4.3.3 Reduced-Rank Linear Discriminant Analysis

LDA 也是一種降維的手段。假設我們有 $p$ 維特征坏瘩， $K$ 個類別盅抚。根據 4.3.2 中介紹的計算方式，我們一共有 $K$ 個類中心點倔矾。他們一定在一個最高 $K-1$ 維的空間里妄均。

例：紅酒分類

例如柱锹，對于 2 個類的分類問題，無論特征是多少維丰包，我們只有 2 個類中心點禁熏。他們必定在一條直線（1維）上。同理邑彪，對于 3 個類的分類問題瞧毙，我們只有 3 個類中心點。他們必定在一個平面（2維）內寄症，如果特征維度大于等于 2宙彪。

因此，經過 LDA有巧，原始數據總能被投影到一個超平面上释漆，其維度為（對應 sklearn LDA 方法中的 n_components 參數）：

$L = \min(p, K-1)$

這說明，在 $p \gg K$ 時剪决，使用 LDA 可以將一個 $p$ 維的輸入降維到 $K-1$ 維灵汪。

我們以 sklearn 中的 wine 數據集為例。它具有 13 維特征柑潦，3 個類別。我們使用 LDA 可以將這些數據投影到一個 2 維的平面上峻凫。

lda.png

代碼：

import pandas as pd
import numpy as np
from sklearn import datasets

wine = datasets.load_wine()
X = pd.DataFrame(wine.data, columns = wine.feature_names)
y = wine.target

# LDA projection
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis
model = LinearDiscriminantAnalysis(n_components=2).fit(X, y)
model.transform(X)

# plot
data_to_plot = pd.DataFrame(
    np.insert(model.transform(X), 2, y, axis=1),
    columns=["x1", "x2", "class"])
show(create_scatter_figure("LDA projection for wine data", data_to_plot))

Find the optimal subspace

在實際中渗鬼，如果 $K$ 也很大，那還需要進一步降低維度荧琼。假設我們目標是降低到 $L$ 維（ $L \ll K-1$ )譬胎，即尋找超平面 $H_{K-1}$ 的最優(yōu)子空間 $H_L$ 。

Fisher 將這個問題提煉為：

找到線性組合 $Z = a^T X$ 使得類間的方差相對于類內方差最大

$X$ 的類內(Within)方差為：

$\textbf{W} = \sum_{k=1}^K \sum_{g_i=k} (x_i - \mu_k)(x_i - \mu_k)^T$

$X$ 的類間(Between)方差為：

$\textbf{B} = \sum_{k=1}^K (\mu - \mu_k)(\mu - \mu_k)^T$

根據向量的統(tǒng)計特性命锄，有 $Z$ 的類內方差為 $a^T \textbf{W} a$ 堰乔，類間方差為 $a^T \textbf{B} a$ 。

于是脐恩，Fisher 實際上是在解決這個優(yōu)化問題：

$\mathop{\arg \max}_a \frac{a^T \textbf{B} a}{a^T \textbf{W} a}$

由于我們總可以通過調節(jié)求得的 $a$ 的系數使得 $a^T \textbf{W} a = 1$ 镐侯，我們可以將其改寫為：

$\begin{align} \mathop{\arg \max}_a & ~ a^T \textbf{B} a \\ \text{s.t.} & ~ a^T \textbf{W} a = 1 \end{align}$

假設 $\mathbf{W}$ 可逆，令 $u = \mathbf{W}^{\frac{1}{2}} a$ 驶冒，由于 $\textbf{W}$ 是對稱矩陣苟翻，有：

$\begin{align} \mathop{\arg \max}_u & ~ u^T \mathbf{W}^{-\frac{1}{2}}\textbf{B}\mathbf{W}^{-\frac{1}{2}} u \\ \text{s.t.} & ~ u^Tu = 1 \end{align}$

$\mathbf{W}^{-\frac{1}{2}}\textbf{B}\mathbf{W}^{-\frac{1}{2}}$ 也是對稱矩陣，必定存在特征值分解：

$\mathbf{W}^{-\frac{1}{2}}\textbf{B}\mathbf{W}^{-\frac{1}{2}} = \mathbf{Q} \mathbf{\Lambda} \mathbf{Q}^T$

因此化簡為：

$\begin{align} \mathop{\arg \max}_u & ~ u^T \mathbf{Q} \mathbf{\Lambda} \mathbf{Q}^T u \\ \text{s.t.} & ~ u^Tu = 1 \end{align}$

再令 $v = \mathbf{Q}^T u$ 骗污，由于 $\mathbf{Q}$ 是單位正交矩陣崇猫， $v^T v = 1$ 依然成立：

$\begin{align} \mathop{\arg \max}_v & ~ v^T \mathbf{\Lambda} v \\ \text{s.t.} & ~ v^Tv = 1 \end{align}$

$\mathbf{\Lambda}$ 是對角矩陣，所以 該優(yōu)化問題本質是求最大特征值需忿。假設 $\lambda_i$ 最大诅炉，顯然在 $v_i = 1$ 時取得最大值 $\lambda_i^2$ 蜡歹。

由于 $\mathbf{\Lambda}$ 是 $\mathbf{W}^{-\frac{1}{2}}\textbf{B}\mathbf{W}^{-\frac{1}{2}}$ 的特征值，為了簡化求解涕烧，我們利用定理：

$\mathbf{AB}$ 與 $\mathbf{BA}$ 具有同樣的特征值月而，如果 $x$ 是 $\mathbf{AB}$ 的某個特征向量，則對應的 $\mathbf{BA}$ 的特征向量是 $y = \mathbf{B}x$

可以得到 $\mathbf{\Lambda}$ 也是 $\mathbf{W}^{-1}\textbf{B}$ 的特征值澈魄。

假設 $\mathbf{W}^{-1}\textbf{B}$ 的最大特征值為 $\lambda_i$ 景鼠，對應的特征向量為 $\xi$ ，則所求的線性變換為：

$a = \mathbf{W}^{-\frac{1}{2}} \xi$

這樣就找到了 $H_L$ 的 1 個維度痹扇，同理铛漓，我們選取 top L 個維度，即得到了 $H_L$ 鲫构。

定理證明

假設 $\lambda$ 是 $\mathbf{AB}$ 的任意特征值浓恶。

$\mathbf{AB}x = \lambda x$

令 $y = \mathbf{B}x$ ，則有：

$\mathbf{A}y = \lambda x$

同時左乘 $\mathbf{B}$ ：

$\mathbf{BA}y = \lambda \mathbf{B}x = \lambda y$

例：手寫數字分類

手寫數字分類是通過 8x8 的手寫數字圖片判斷是 0-9 中的哪個數字结笨。顯然包晰，這是一個具有 $p = 8 \times 8 = 64$ 維特征， $K = 10$ 個類別的分類任務炕吸。我們可以用 sklearn 庫的 LDA 得到下面的結果（降至 2 維 plot）：

sklearn_lda.png

import pandas as pd
import numpy as np
import scipy
from sklearn import datasets

digits = datasets.load_digits()
X = pd.DataFrame(digits.data, columns = digits.feature_names)
y = digits.target

trainX = X[:len(X)//2]
trainY = y[:len(y)//2]

from sklearn.discriminant_analysis import LinearDiscriminantAnalysis

model = LinearDiscriminantAnalysis(n_components=2, solver="eigen", shrinkage=0.01).fit(trainX, trainY)

# reduce rank to 2 dimension data
reduceRankX = model.transform(trainX)

顯然伐憾，直接調包隱藏了太多細節(jié)，為了加深理解赫模，我們根據公式推導自己手動實現一個：

def shrink(cov, shrinkage):
    dimensions = cov.shape[0]
    return (1-shrinkage) * cov  + shrinkage * np.trace(cov) / dimensions * np.identity(dimensions)

def my_lda(X, y, n_components, shrinkage):
    T = X.cov()
    W = X.groupby(y).apply(lambda g: g.cov() * (len(g)-1)
                          ).groupby(level=1).sum() / X.shape[0]

    shrunkW = shrink(W, shrinkage)
    invShrunkW = np.linalg.inv(shrunkW)
    shrunkB = shrink(T, shrinkage) - shrunkW

    # eigen values is ascending
    eigenvalues, eigenvectors = np.linalg.eigh(invShrunkW.dot(shrunkB))

    # eigen vectors with greatest eigen values
    xi = eigenvectors[:,::-1][:,:n_components]


    # L = np.linalg.cholesky(invShrunkW)
    # return L.dot(xi)

    # W^{-1/2}, not Cholesky
    return scipy.linalg.sqrtm(invShrunkW).dot(xi)

其中树肃，shrink 函數是為了處理當 $\mathbf{W}$ 是 奇異矩陣 (不可逆)的情形。我們可以使用 shrunk covariance 來使其變?yōu)榭赡婢仃囋偬幚砥俾蕖＿@在輸入 $X$ 是稀疏矩陣時很常見胸嘴，例如手寫數字分類。

$\mathbf{\Sigma}_{\text{shrunk}} = (1 - \alpha) \hat{\mathbf{\Sigma}} + \alpha \frac{\mathop{tr} \hat{\mathbf{\Sigma}}}{p} \mathbf{I}$

其中 $\alpha$ 是 shrinkage斩祭， $p$ 是特征維度劣像。

另外還需要注意區(qū)分矩陣平方根 sqrtm 與 Cholesky Decomposition。

矩陣平方根指 $\mathbf{B} \mathbf{B} = \mathbf{B}^T \mathbf{B} = \mathbf{A}$ 摧玫，要求 $\mathbf{B}$ 是對稱矩陣耳奕。

Cholesky Decomposition 指 $\mathbf{L}^T \mathbf{L} = \mathbf{A}$ ，要求 $\mathbf{L}$ 是三角矩陣席赂。

使用自己實現的 LDA 得到與 sklearn 實現類似的結果：

my_lda.png

可以看出吮铭，在降到 2 維后，LDA 還是能夠清楚區(qū)分出數字 0, 1, 2, 3, 4, 6颅停，但是存在一些數字的類別重疊在一起的情況谓晌。此時可以 增加維度 $L$ 解決。

下面是第3癞揉、4 維：

extra_dimension_lda.png

可以看出纸肉，它較好的補足了 1溺欧、2 維無法正確分類的數字。對于數字 5柏肪、7 有了清楚的分類姐刁。

4.4 Logistic regression

邏輯回歸希望用輸入 $x$ 的線性組合建模屬于 k 類的后驗概率，并且要求所有概率之和為 1烦味。

以最后一個類（K 類）的后驗概率 $\text{Pr} (G = K | X = x)$ 為比較基準聂使，有：

$\ln \frac{\text{Pr} (G = k | X = x)}{\text{Pr} (G = K | X = x)} = \beta_{k0} + \beta_{k1}^T x$

可以寫為：

$\text{Pr} (G = k | X = x) = \begin{cases} \dfrac{e^{\beta_{k0} + \beta_{k1}^Tx}}{1 + \sum_{l=1}^{K-1} e^{\beta_{l0} + \beta_{l1}^Tx}}, & k = 1, 2, ..., K-1 \\ \dfrac{1}{1 + \sum_{l=1}^{K-1} e^{\beta_{l0} + \beta_{l1}^Tx}}, & k = K \end{cases}$

4.4.1 Fitting Logistic Regression Models

對于第 k 類的某個樣本 i，我們希望找到一組模型參數 $\theta$ 谬俄，使模型認為該樣本屬于 k 的概率（即 正確分類概率）盡量大柏靶。我們將這個概率計作：

$\text{Pr} (G = k | X = x) = p_{k} (x; \theta)$

則該優(yōu)化問題是：

$\begin{align} \hat{\theta} = \mathop{\arg \max}_{\theta} ~ p_{k} (x; \theta) \end{align}$

對于所有樣本，我們目標是令他們 被正確分類的概率和最大：

$\ell (\theta) = \sum_{i=1}^N \ln p_{g_i} (x_i; \theta)$

為了簡化這個問題溃论，我們只考慮而 $K=2$ 即 二分類 問題情況屎蜓。此時樣本要么屬于 1 類要么屬于 2 類。

為了簡單表示：

$p_{k} (x; \theta) = \begin{cases} p_1 (x; \theta), & k = 1 \\ 1 - p_1(x; \theta), & k = 2 \end{cases}$

我們可以假設樣本觀測值 $y_i = 1$ 代表屬于第 1 類钥勋， $y_i = 0$ 代表 不屬于 第 1 類（那么肯定屬于第 2 類）炬转。則

$p_{k} (x_i; \theta) = y_i p_1 (x_i; \theta) + (1-y_i)(1-p_1(x_i; \theta))$

則有：

$\begin{align} \ell (\theta) &= \sum_{i=1}^N \ln p_{g_i} (x_i; \theta) \\ &= \sum_{i=1}^N {y_i \ln p_1 (x_i; \theta) + (1-y_i) \ln (1-p_1(x_i; \theta))} \\ &= \sum_{i=1}^N {\ln(1 - p_1(x_i; \theta)) + y_i \ln \frac{ p_1 (x_i; \theta)}{1 - p_1 (x_i; \theta)} } \\ &= \sum_{i=1}^N {\ln(1 - p_1(x_i; \theta)) + y_i \mathbf{\beta}^Tx_i } \\ &= \sum_{i=1}^N {y_i \mathbf{\beta}^Tx_i - \ln(1 + e^{\mathbf{\beta}^Tx_i}) } \end{align}$

其中 $\mathbf{\beta} = [ \beta_{10}, \beta_{11}, ... \beta_{1p} ]$ ，對應 $x = [1, x_1, x_2, ... x_p]$ 算灸。

求解這個優(yōu)化問題可以令其對 $\beta$ 的導數為 0:

$\frac{\partial \ell(\beta)}{\partial \beta} = \sum_{i=1}^N x_i(y_i - p_1(x_i; \beta)) = 0,$

為了求解 $\dfrac{\partial \ell(\beta)}{\partial \beta} = 0$ 扼劈，我們可以用凸優(yōu)化中的 Newton-Raphson 法。即選取一個任意初始 $x_n$ （不是初始解，因為要求的就是解）烦绳，并取其切線與 x 軸交點 $x_{n+1}$ ：

$\beta^{\text{new}} = \beta - (\frac{\partial^2 \ell(\beta)}{\partial \beta \partial \beta^T})^{-1} \frac{\partial \ell(\beta)}{\partial \beta}$

其中 $(p+1) \times (p+1)$ 維 Hessian 矩陣：

$\begin{align} \frac{\partial^2 \ell(\beta)}{\partial \beta \partial \beta^T} &= - \sum_{i=1}^N x_i \frac{\partial p_1(x_i; \beta)}{\partial \beta^T} \\ &= - \sum_{i=1}^N x_i \dfrac{\partial \dfrac{e^{\beta^T x_i}}{1 + e^{\beta^T x_i}} }{\partial \beta^T} \\ &= - \sum_{i=1}^N x_i \frac{x_i^T e^{\beta^T x_i}}{(1 + e^{\beta^T x_i})^2} \\ &= -\sum_{i=1}^N x_i x_i^T p_1(x_i; \beta) (1 - p_1(x_i; \beta)) \end{align}$

我們將他們寫成矩陣形式：

$\frac{\partial \ell(\beta)}{\partial \beta} = \mathbf{X}^T(\mathbf{y} - \mathbf{p})$

$\frac{\partial^2 \ell(\beta)}{\partial \beta \partial \beta^T} = \mathbf{X}^T \mathbf{W} \mathbf{X}$

其中 $\mathbf{X}$ 是 $N \times (p+1)$ 矩陣， $\mathbf{y} - \mathbf{p}$ 是 $N$ 維列向量， $\mathbf{W}$ 是 $N \times N$ 對角矩陣檐嚣，第 i 個元素是 $p_1(x_i; \beta) (1 - p_1(x_i; \beta))$ 。

帶入 Newton-Raphson 公式：

$\begin{align} \beta^{\text{new}} &= \beta - (\frac{\partial^2 \ell(\beta)}{\partial \beta \partial \beta^T})^{-1} \frac{\partial \ell(\beta)}{\partial \beta} \\ &= \beta + (\mathbf{X}^T \mathbf{W} \mathbf{X})^{-1} \mathbf{X}^T(\mathbf{y} - \mathbf{p}) \\ &= (\mathbf{X}^T \mathbf{W} \mathbf{X})^{-1}\mathbf{X}^T \mathbf{W}(\mathbf{X}\beta + \mathbf{W}^{-1}(\mathbf{y} - \mathbf{p})) \\ &= (\mathbf{X}^T \mathbf{W} \mathbf{X})^{-1}\mathbf{X}^T \mathbf{W}\mathbf{z} \end{align}$

這就是 迭代更新 的計算公式蘸鲸。由于 $\mathbf{p}$ 和 $\mathbf{W}$ 都隨 $\beta$ 變化噪伊，因此我們需要在每輪迭代重新計算他們的值。

上式中磁携，我們定義了：

$\mathbf{z} = \mathbf{X}\beta + \mathbf{W}^{-1}(\mathbf{y} - \mathbf{p})$

回想 最小二乘法 的公式：

$\hat{\beta} = (\mathbf{X}^T \mathbf{X})^{-1} \mathbf{X}^T \mathbf{y}$

可以發(fā)現兩者非常相似褒侧，區(qū)別在于邏輯回歸中多了一個權重對角矩陣 $\mathbf{W}$ 。因此這個算法也稱為 iteratively reweighted least squares（重新加權迭代最小二乘法）谊迄。

接下來我們還需要一個 $\beta$ 的 初始值闷供。一般情況下可以選 $\beta = 0$ 。該算法 不保證收斂统诺。

例：乳腺癌診斷

該案例通過 30 維特征來判斷乳腺癌是惡性還是良性歪脏，是一個典型的 2 分類問題。

import pandas as pd
import numpy as np
from sklearn import datasets, model_selection

def load_sklearn_data(sk_data):
    X = pd.DataFrame(sk_data.data, columns = sk_data.feature_names)
    y = pd.Series(sk_data.target, name="target").apply(
        lambda index: sk_data.target_names[index]).astype("category")
    return X, y

X, y = load_sklearn_data(datasets.load_breast_cancer())
X_train, X_test, y_train, y_test = model_selection.train_test_split(
    X, y, test_size=0.2)

from sklearn.metrics import accuracy_score
from sklearn.linear_model import LogisticRegression

model = LogisticRegression(solver="newton-cg").fit(X_train, y_train)
accuracy_score(model.predict(X_test), y_test)

結果是 0.9385964912280702粮呢。我們再用自己實現的 LR 進行判斷：

def sigmoid(x, beta):
    tmp = np.exp(x.dot(beta))
    return tmp / (1 + tmp)

class MyLogisticRegression:
    def __init__(self, max_iter, tolerance = 0.01):
        self.max_iter_ = max_iter
        self.tolerance_ = tolerance
        self.beta = None

    def fit(self, train_X, train_y):
        N, p = train_X.shape
        X = train_X.to_numpy()
        X = np.concatenate((np.atleast_2d(np.ones(N)).T,
                            train_X.to_numpy()), axis=1)
        self.beta = np.zeros(p+1)
        beta = self.beta
        for i in range(self.max_iter_):
            prob = np.apply_along_axis(sigmoid, 1, X, beta)
            W = np.diag(prob * (1 - prob))
            z = X.dot(beta) + np.linalg.inv(W).dot(train_y.cat.codes - prob)
            new_beta = np.linalg.inv(X.T @ W @ X) @ X.T @ W.dot(z)
            beta = new_beta

        self.beta = beta


    def predict(self, test_X):
        X = test_X.copy()
        X.insert(0, "x_0", np.ones(test_X.shape[0]))
        prob = X.apply(lambda row: sigmoid(row, self.beta), axis=1)
        return prob.apply(lambda x: 1 if x >= 0.5 else 0)

根據公式寫代碼婿失，注意 $\beta$ 的迭代更新即可钞艇。最大的難度是終止條件，因為迭代數次之后會遇到 Hessian 矩陣 $\mathbf{X}^T \mathbf{W} \mathbf{X}$ 成了奇異矩陣的問題豪硅。

運行的結果表明分類正確率可以達到 0.9473684210526315哩照。

TODO: how to find a stopping criteria?

4.4.5 Logistic Regression or LDA?

在 4.3 中我們介紹了 LDA，回憶其判斷某個樣本屬于 $k$ 還是 $l$ 的函數（當該式大于 0 時屬于 $k$ 類）：

$\begin{align} \ln \frac{\text{Pr}(G=k | X=x)}{\text{Pr}(G=l | X=x)} &= \ln \frac{f_k(x) \pi_k}{ f_l(x) \pi_l} \\ &= \ln \frac{\pi_k}{\pi_l} - \frac{1}{2}(\mu_k + \mu_l) \mathbf{\Sigma}^{-1} (\mu_k - \mu_l) + x^T \mathbf{\Sigma}^{-1} (\mu_k - \mu_l) \\ &= \alpha_{k0} + \alpha_{k1}^T x \end{align}$

可以看出其形式與 logistic regression 一致：

$\ln \frac{\text{Pr}(G=k | X=x)}{\text{Pr}(G=l | X=x)} = \beta_{k0} + \beta_{k1}^T x$

那么他們倆是不是一樣的方法呢懒浮？并不是飘弧。他們的主要區(qū)別是：

假設不同。LDA 假設了 $x$ 服從正態(tài)分布砚著，且各個類的協方差矩陣相同次伶。LR 沒有限定 $x$ 的分布。由于 LDA 假設了正態(tài)分布赖草，它對于極端的 outlier 魯棒性差（會影響分布函數)学少。
優(yōu)化目標不同。LDA 最大化 full log-likelihood秧骑，需要考慮 $x$ 的邊緣分布函數版确。LR 最大化 conditional likelihood，不需要考慮 $x$ 的邊緣分布函數乎折。

比較難理解的是第二點绒疗。

首先，根據 Bayes 公式骂澄，在已知 B 發(fā)生時吓蘑，A 的條件概率為：

$P(A|B) = \frac{P( A \cap B ) }{P(B)}$

因此對于分類問題，條件概率(conditional likelihood) 是 在已知 $X = x$ 時坟冲，樣本類別 $G = k$ 類的概率磨镶。它與 全概率(full likelihood) 之間存在關系：

$\begin{align} \text{Pr}(G=k | X=x) &= \frac{\text{Pr}(X, G=k)}{\text{Pr}(X)} \\ &= \frac{\text{Pr}(X, G=k)}{\sum_{l=1}^K \text{Pr}(X, G=l)} \\ \end{align}$

LDA 的目標是找到一個類別 $k$ 使得后驗概率最大，即：

$\begin{align} \hat{k} &= \mathop{\arg \max}_{k} \text{Pr}(G=k | X=x) \\ &= \mathop{\arg \max}_{k} \frac{\text{Pr}(X, G=k)}{\sum_{l=1}^K \text{Pr}(X, G=l)} \end{align}$

LDA 選擇 對全概率建模健提，需要知道 $x$ 的在每個類的密度函數 $f_i(x)$ 琳猫，LDA 假設其為正態(tài)分布，且協方差相等私痹。全概率為：

$\text{Pr}(X, G=k) = f_k(x) \pi_k = \phi(x; \mu_k, \mathbf{\Sigma}) \pi_k$

所以說它考慮的是 full log-likelihood脐嫂。

對于 LR，它直接 對條件概率建模紊遵。即假設存在一組 $\beta$ 能夠使條件概率的 log ratio 和 $x$ 有如下線性關系：

$\ln \frac{\text{Pr}(G=k | X=x)}{\text{Pr}(G=l | X=x)} = \beta_{k0} + \beta_{k1}^T x$

并且账千，這一組 $\beta$ 使 所有樣本正確分類的概率和 最大：

$\begin{align} \hat{\theta} &= \mathop{\arg \max}_{\theta} ~ \ell (\theta) \\ &= \mathop{\arg \max}_{\theta} ~ \sum_{i=1}^N \ln p_{g_i} (x_i; \theta) \\ &= \sum_{i=1}^N y_i \mathbf{\beta}^Tx_i - \ln(1 + e^{\mathbf{\beta}^Tx_i}) \end{align}$

其中不包含任何與邊緣密度函數 $\text{Pr}(X)$ 相關的部分，也沒有假設 $x$ 的分布暗膜。因此說它只考慮 conditional likelihood匀奏。

Reference

masking in linear regression for multiple classes
Linear discriminant analysis, explained
shrunk covariance
between-class scatter matrix calculation
邏輯回歸
Implementing logistic regression from scratch in Python