[Scikit-learn] 1.1 Generalized Linear Models - Logistic regression & Softmax

時間 2019-11-18

標籤 scikit learn 1.1 generalized linear models logistic regression softmax 欄目應用數學简体版

原文原文鏈接

二分類：Logistic regressionphp

多分類：Softmax分類函數html

對於損失函數，咱們求其最小值，git

對於似然函數，咱們求其最大值。算法

Logistic是loss function，即：shell

　　在邏輯迴歸中，選擇了「對數似然損失函數」，L(Y,P(Y|X)) = -logP(Y|X)。app

　　對似然函數求最大值，其實就是對對數似然損失函數求最小值。dom

Logistic regression, despite its name, is a linear model for classification rather than regression. 機器學習

（但非迴歸，實際上是分類，經過後驗機率比較，那麼如何求MLE就是收斂的核心問題）ide

邏輯迴歸與正則化

As an optimization problem, binary class L2 penalized logistic regression minimizes the following cost function:函數

Similarly, L1 regularized logistic regression solves the following optimization problem

Optimal Solvers

邏輯迴歸提供的四種「漸近算法」

The solvers implemented in the class LogisticRegression are 「liblinear」, 「newton-cg」, 「lbfgs」 and 「sag」: （收斂算法，例如梯度降低等）

「lbfgs」和「newton-cg」只支持L2罰項，而且對於一些高維數據收斂很是快。L1罰項產生稀疏預測的權重。

「liblinear」使用了基於Liblinear的座標降低法(CD)。

　　對於L1罰項， sklearn.svm.l1_min_c 容許計算C的下界以得到一個非」null」的模型（全部特徵權重爲0）。

　　這依賴於很是棒的一個庫 LIBLINEAR library，用在scikit-learn中。然而，CD算法在liblinear中的實現沒法學習一個真正的多維（多類）的模型；反而，最優問題被分解爲「one-vs-rest」多個二分類問題來解決多分類。

　　因爲底層是這樣實現的，因此使用了該庫的 LogisticRegression 類就能夠做爲多類分類器了。

特色分析

LogisticRegression 使用「lbfgs」或者「newton-cg」程序來設置 multi_class 爲「」，則該類學習了一個真正的多類邏輯迴歸模型，也就是說這種機率估計應該比默認「one-vs-rest」設置要更加準確。

可是「lbfgs」, 「newton-cg」和「sag」程序沒法優化含L1罰項的模型，因此」multinomial」的設置沒法學習稀疏模型。

「sag」程序使用了隨機平均梯度降低（ Stochastic Average Gradient descent）。它沒法解決多分類問題，並且對於含L2罰項的模型有侷限性。然而在超大數據集下計算要比其餘程序快不少，當樣本數量和特徵數量都很是大的時候。

In a nutshell, one may choose the solver with the following rules:

Case	Solver
Small dataset or L1 penalty	「liblinear」
Multinomial loss or large dataset	「lbfgs」, 「sag」 or 「newton-cg」
Very Large dataset	「sag」

注：

For large dataset, you may also consider using SGDClassifier(Stochastic Gradient Descent) with ‘log’ loss.

Goto: [Scikit-learn] 1.5 Generalized Linear Models - Stochastic Gradient Descent

數學原理

原理大放送

Ref: 機器學習算法與Python實踐之（七）邏輯迴歸（Logistic Regression）

Ref: 對於logistic函數的交叉熵損失函數

(a) 機率模型 - Logistic Regression

將「參數」與「變量」的關係轉化爲了機率關係；σ是sigmoid函數。

本質就是：轉化爲機率的思惟模式；轉化的方式有不少種，這種相對最好。

LogisticRegression 就是一個被logistic方程歸一化後的線性迴歸。

(b) 代價函數 - MLE

LogisticRegression最基本的學習算法是最大似然。

似然做爲 "cost function"

【交叉熵】

這個似然的角度得出的上述公式結論，其實也就是交叉熵：goto 簡單的交叉熵損失函數，你真的懂了嗎？

(c) 優化似然

Log 似然做爲 "cost function" 會容易計算一些。

(d) 漸近求解

導數 = 0：但貌似很差求，便有了 solver (「liblinear」, 「newton-cg」, 「lbfgs」 and 「sag」)。

其中：solver = 「sag」時，參考：[Scikit-learn] 1.5 Generalized Linear Models - SGD for Classification中的

clf = SGDClassifier(loss="log", alpha=0.01, n_iter=200, fit_intercept=True)

一個 loss function 對應一堆 solvers；

一個 solver 對應一堆 loss functions；

逼近的過程形式化爲迭代公式：

Logistic 代碼實現

邏輯實現

from __future__ import print_function import numpy as np import matplotlib.pyplot as plt # Define the logistic function
def logistic(z): return 1 / (1 + np.exp(-z))　　` # Plot the logistic function
z = np.linspace(-6,6,100) plt.plot(z, logistic(z), 'b-')　　// <-- 畫圖 plt.xlabel('$z$', fontsize=15) plt.ylabel('$\sigma(z)$', fontsize=15) plt.title('logistic function') plt.grid() plt.show()

Result:

求導實現

# Define the logistic function
def logistic_derivative(z): return logistic(z) * (1 - logistic(z))




 # Plot the derivative of the logistic function
z = np.linspace(-6,6,100) plt.plot(z, logistic_derivative(z), 'r-')　　// <-- 畫圖 plt.xlabel('$z$', fontsize=15) plt.ylabel('$\\frac{\\partial \\sigma(z)}{\\partial z}$', fontsize=15) plt.title('derivative of the logistic function') plt.grid() plt.show()

Result:

線性多分類

Ref: LogisticRegression（參數解密）

Logistic Regression (aka logit, MaxEnt) classifier.

多分類策略

‘multi_class’ option :

- ‘ovr’: one-vs-rest (OvR)
- 'multinomial': cross- entropy loss【it is supported only by the ‘lbfgs’, ‘sag’ and ‘newton-cg’ solvers.】

正則化的參數

參數 penalty

【str, ‘l1’, ‘l2’, ‘elasticnet’ or ‘none’, optional (default=’l2’)】

This class implements regularized logistic regression using the ‘liblinear’ library, ‘newton-cg’, ‘sag’ and ‘lbfgs’ solvers.

參數 solver

【str, {‘newton-cg’, ‘lbfgs’, ‘liblinear’, ‘sag’, ‘saga’}, optional (default=’liblinear’)】

It can handle both dense and sparse input. Use C-ordered arrays or CSR matrices containing 64-bit floats for optimal performance; any other input format will be converted (and copied).

- The ‘newton-cg’, ‘sag’, and ‘lbfgs’ solvers support only L2 regularization with primal formulation.
- The ‘liblinear’ solver supports both L1 and L2 regularization, with a dual formulation only for the L2 penalty.

參數 C

【float, default: 1.0】

Inverse of regularization strength; must be a positive float. Like in support vector machines, smaller values specify stronger regularization.

Answer: 參數越小則代表越強的正則化：正則項係數的倒數

參數 multi_class

【str, {‘ovr’, ‘multinomial’}, default: ‘ovr’】

Multiclass option can be either ‘ovr’ or ‘multinomial’. If the option chosen is ‘ovr’, then a binary problem is fit for each label. Else the loss minimised is the multinomial loss fit across the entire probability distribution. Works only for the ‘newton-cg’, ‘sag’ and ‘lbfgs’ solver.

New in version 0.18: Stochastic Average Gradient descent solver for ‘multinomial’ case.

正則化邏輯迴歸

例一：L1 Penalty and Sparsity in Logistic Regression

# 有點僅一層的logit NN的意思

""" ============================================== L1 Penalty and Sparsity in Logistic Regression ============================================== Comparison of the sparsity (percentage of zero coefficients) of solutions when L1 and L2 penalty are used for different values of C. We can see that large values of C give more freedom to the model. Conversely, smaller values of C constrain the model more. In the L1 penalty case, this leads to sparser solutions. We classify 8x8 images of digits into two classes: 0-4 against 5-9. The visualization shows coefficients of the models for varying C. """

print(__doc__) # Authors: Alexandre Gramfort <alexandre.gramfort@inria.fr> # Mathieu Blondel <mathieu@mblondel.org> # Andreas Mueller <amueller@ais.uni-bonn.de> # License: BSD 3 clause

import numpy as np import matplotlib.pyplot as plt from sklearn.linear_model import LogisticRegression from sklearn import datasets from sklearn.preprocessing import StandardScaler digits = datasets.load_digits() X, y = digits.data, digits.target X = StandardScaler().fit_transform(X) # 四舍五如
y = (y > 4).astype(np.int) # Set regularization parameter
for i, C in enumerate((100, 1, 0.01)): # turn down tolerance for short training time
    clf_l1_LR = LogisticRegression(C=C, penalty='l1', tol=0.01) clf_l2_LR = LogisticRegression(C=C, penalty='l2', tol=0.01) clf_l1_LR.fit(X, y) clf_l2_LR.fit(X, y) coef_l1_LR = clf_l1_LR.coef_.ravel() coef_l2_LR = clf_l2_LR.coef_.ravel() # coef_l1_LR contains zeros due to the
    # L1 sparsity inducing norm
 sparsity_l1_LR = np.mean(coef_l1_LR == 0) * 100 sparsity_l2_LR = np.mean(coef_l2_LR == 0) * 100

    print("C=%.2f" % C) print("Sparsity with L1 penalty: %.2f%%" % sparsity_l1_LR) print("score with L1 penalty: %.4f" % clf_l1_LR.score(X, y)) print("Sparsity with L2 penalty: %.2f%%" % sparsity_l2_LR) print("score with L2 penalty: %.4f" % clf_l2_LR.score(X, y)) l1_plot = plt.subplot(3, 2, 2 * i + 1) l2_plot = plt.subplot(3, 2, 2 * (i + 1)) if i == 0: l1_plot.set_title("L1 penalty") l2_plot.set_title("L2 penalty") l1_plot.imshow(np.abs(coef_l1_LR.reshape(8, 8)), interpolation='nearest', cmap='binary', vmax=1, vmin=0) l2_plot.imshow(np.abs(coef_l2_LR.reshape(8, 8)), interpolation='nearest', cmap='binary', vmax=1, vmin=0) plt.text(-8, 3, "C = %.2f" % C) l1_plot.set_xticks(()) l1_plot.set_yticks(()) l2_plot.set_xticks(()) l2_plot.set_yticks(()) plt.show()

Result:

C=100.00 Sparsity with L1 penalty: 6.25% score with L1 penalty: 0.9104 Sparsity with L2 penalty: 4.69% score with L2 penalty: 0.9098 C=1.00 Sparsity with L1 penalty: 9.38% score with L1 penalty: 0.9098 Sparsity with L2 penalty: 4.69% score with L2 penalty: 0.9093 C=0.01 Sparsity with L1 penalty: 85.94% score with L1 penalty: 0.8598 Sparsity with L2 penalty: 4.69% score with L2 penalty: 0.8915

Jeff：如上，權重的歸零化顯而易見。

例二：Path with L1- Logistic Regression

print(__doc__) # Author: Alexandre Gramfort <alexandre.gramfort@inria.fr> # License: BSD 3 clause

from datetime import datetime import numpy as np import matplotlib.pyplot as plt from sklearn import linear_model from sklearn import datasets from sklearn.svm import l1_min_c iris = datasets.load_iris() X = iris.data y = iris.target X = X[y != 2] y = y[y != 2] X -= np.mean(X, 0) ############################################################################### # Demo path functions cs = l1_min_c(X, y, loss='log') * np.logspace(0, 3) # 對正則化的強度取不一樣值，實驗各自效果 print("Computing regularization path ...") start = datetime.now() clf = linear_model.LogisticRegression(C=1.0, penalty='l1', tol=1e-6) coefs_ = [] for c in cs: clf.set_params(C=c) clf.fit(X, y) coefs_.append(clf.coef_.ravel().copy()) print("This took ", datetime.now() - start) coefs_ = np.array(coefs_) plt.plot(np.log10(cs), coefs_) ymin, ymax = plt.ylim() plt.xlabel('log(C)') plt.ylabel('Coefficients') plt.title('Logistic Regression Path') plt.axis('tight') plt.show()

係數的變化

coefs_ Out[30]: array([[ 0. , 0. , 0. , 0. ], [ 0. , 0. , 0.17813685, 0. ], [ 0. , 0. , 0.3384641 , 0. ], ..., [ 0. , -1.29412716,  5.71841258, 0. ], [ 0. , -1.41538897,  5.82221819, 0. ], [ 0. , -1.53535948,  5.92878532,  0.        ]])

Result:

Jeff：四個係數對應四根線，相似於trace plot，其實仍是在說L1, L2的東西。

Softmax 多分類

Ref: Softmax函數解密

Ref: http://www.cnblogs.com/daniel-D/archive/2013/05/30/3109276.html

Ref: http://ufldl.stanford.edu/wiki/index.php/Softmax%E5%9B%9E%E5%BD%92

立體 sigmoid

先來個炫酷的立體sigmoid做爲開場：

import numpy as np import matplotlib.pyplot as plt from matplotlib import cm # Define the softmax function
def softmax(z): return np.exp(z) / np.sum(np.exp(z)) #-------------------------------------------------------------------------- # Plot the softmax output for 2 dimensions for both classes # Plot the output in function of the weights # Define a vector of weights for which we want to plot the ooutput
nb_of_zs = 200 zs = np.linspace(-10, 10, num=nb_of_zs) # input 
zs_1, zs_2 = np.meshgrid(zs, zs) # generate grid
y = np.zeros((nb_of_zs, nb_of_zs, 2)) # initialize output
 # Fill the output matrix for each combination of input z's
for i in range(nb_of_zs): for j in range(nb_of_zs): y[i,j,:] = softmax(np.asarray([zs_1[i,j], zs_2[i,j]]))
 # Plot the cost function surfaces for both classes
fig = plt.figure() # Plot the cost function surface for t=1
ax = fig.gca(projection='3d') surf = ax.plot_surface(zs_1, zs_2, y[:,:,0], linewidth=0, cmap=cm.coolwarm) ax.view_init(elev=30, azim=70) cbar = fig.colorbar(surf) ax.set_xlabel('$z_1$', fontsize=15) ax.set_ylabel('$z_2$', fontsize=15) ax.set_zlabel('$y_1$', fontsize=15) ax.set_title ('$P(t=1|\mathbf{z})$') cbar.ax.set_ylabel('$P(t=1|\mathbf{z})$', fontsize=15) plt.grid() plt.show()

Result:

基本原理

在數字手寫識別中，咱們要識別的是十個類別，每次從輸入層輸入 28×28 個像素，輸出層就能夠獲得本次輸入可能爲 0, 1, 2… 的機率。

簡易版本的，看起來更直觀:

OK, 左邊是輸入層，輸入的 x 經過中間的黑線 w (包含了 bias 項)做用下，獲得 w.x, 到達右邊節點，右邊節點經過紅色的函數將這個值映射成一個機率，預測值就是輸入機率最大的節點，這裏可能的值是 {0, 1, 2}。

在實現Softmax迴歸時，將 $\textstyle \theta$ 用一個 $\textstyle k \times(n+1)$ 的矩陣來表示會很方便，該矩陣是將 $\theta_1, \theta_2, \ldots, \theta_k$ 按行羅列起來獲得的，以下所示： (k是類別數，n是輸入channel數)

$\theta = \begin{bmatrix}\mbox{---} \theta_1^T \mbox{---} \\\mbox{---} \theta_2^T \mbox{---} \\\vdots \\\mbox{---} \theta_k^T \mbox{---} \\\end{bmatrix}$

在 softmax regression 中，輸入的樣本屬於第 j 類的機率能夠寫成：

注意到，這個迴歸的參數向量減去一個常數向量，會有什麼結果：

沒有變化！(指數函數的無記憶性)

若是某一個向量是代價函數的極小值點，那麼這個向量在減去一個任意的常數向量也是極小值點，這是由於 softmax 模型被過分參數化了。

進一步而言，若是參數 $\textstyle (\theta_1, \theta_2,\ldots, \theta_k)$ 是代價函數 $\textstyle J(\theta)$ 的極小值點，那麼 $\textstyle (\theta_1 - \psi, \theta_2 - \psi,\ldots,\theta_k - \psi)$ 一樣也是它的極小值點，其中 $\textstyle \psi$ 能夠爲任意向量。所以使 $\textstyle J(\theta)$ 最小化的解不是惟一的。

（有趣的是，因爲 $\textstyle J(\theta)$ 仍然是一個凸函數，所以梯度降低時不會遇到局部最優解的問題。可是 Hessian 矩陣是奇異的/不可逆的，這會直接致使採用牛頓法優化就遇到數值計算的問題）

注意，當 $\textstyle \psi = \theta_1$ 時，咱們老是能夠將 $\textstyle \theta_1$ 替換爲 $\textstyle \theta_1 - \psi = \vec{0}$ （即替換爲全零向量），而且這種變換不會影響假設函數。所以咱們能夠去掉參數向量 $\textstyle \theta_1$ （或者其餘 $\textstyle \theta_j$ 中的任意一個）而不影響假設函數的表達能力。實際上，與其優化所有的 $\textstyle k\times(n+1)$ 個參數 $\textstyle (\theta_1, \theta_2,\ldots, \theta_k)$ （其中 $\textstyle \theta_j \in \Re^{n+1}$ ），咱們能夠令 $\textstyle \theta_1 =\vec{0}$ ，只優化剩餘的 $\textstyle (k-1)\times(n+1)$ 個參數，這樣算法依然可以正常工做。

在實際應用中，爲了使算法實現更簡單清楚，每每保留全部參數 $\textstyle (\theta_1, \theta_2,\ldots, \theta_n)$ ，而不任意地將某一參數設置爲 0。但此時咱們須要對代價函數作一個改動：加入權重衰減。權重衰減能夠解決 softmax 迴歸的參數冗餘所帶來的數值問題。

既然模型被過分參數化了，咱們就事先肯定一個參數，好比將 w1 替換成全零向量，將 w1.x = 0 帶入 binomial softmax regression ，獲得了咱們最開始的二項 logistic regression (能夠動手算一算), 用圖就能夠表示爲：

(注：虛線表示爲 0 的權重，在第一張圖中沒有畫出來，能夠看到 logistic regression 就是 softmax regression 的一種特殊狀況)

權重衰減 - 正則化

咱們經過添加一個權重衰減項 $\textstyle \frac{\lambda}{2} \sum_{i=1}^k \sum_{j=0}^{n} \theta_{ij}^2$ 來修改代價函數（L2原理？），這個衰減項會懲罰過大的參數值，如今咱們的代價函數變爲：

$\begin{align}J(\theta) = - \frac{1}{m} \left[ \sum_{i=1}^{m} \sum_{j=1}^{k} 1\left\{y^{(i)} = j\right\} \log \frac{e^{\theta_j^T x^{(i)}}}{\sum_{l=1}^k e^{ \theta_l^T x^{(i)} }} \right] + \frac{\lambda}{2} \sum_{i=1}^k \sum_{j=0}^n \theta_{ij}^2\end{align}$

有了這個權重衰減項之後 ( $\textstyle \lambda > 0$ )，代價函數就變成了嚴格的凸函數，這樣就能夠保證獲得惟一的解了。此時的 Hessian矩陣變爲可逆矩陣，而且由於 $\textstyle J(\theta)$ 是凸函數，梯度降低法和 L-BFGS 等算法能夠保證收斂到全局最優解。

爲了使用優化算法，咱們須要求得這個新函數 $\textstyle J(\theta)$ 的導數，以下：

$\begin{align}\nabla_{\theta_j} J(\theta) = - \frac{1}{m} \sum_{i=1}^{m}{ \left[ x^{(i)} ( 1\{ y^{(i)} = j\} - p(y^{(i)} = j | x^{(i)}; \theta) ) \right] } + \lambda \theta_j\end{align}$

經過最小化 $\textstyle J(\theta)$ ，咱們就能實現一個可用的 softmax 迴歸模型。

Softmax 迴歸 vs. k 個二元分類器

若是你在開發一個音樂分類的應用，須要對k種類型的音樂進行識別，那麼是選擇使用 softmax 分類器呢，仍是使用 logistic 迴歸算法創建 k 個獨立的二元分類器呢？

這一選擇取決於你的類別之間是否互斥，

- 例如，若是你有四個類別的音樂，分別爲：古典音樂、鄉村音樂、搖滾樂和爵士樂，那麼你能夠假設每一個訓練樣本只會被打上一個標籤（即：一首歌只能屬於這四種音樂類型的其中一種），此時你應該使用類別數 $k = 4 的softmax迴歸。（若是在你的數據集中，有的歌曲不屬於以上四類的其中任何一類，那麼你能夠添加一個「其餘類」，並將類別數 k 設爲5。）$
- 若是你的四個類別以下：人聲音樂、舞曲、影視原聲、流行歌曲，那麼這些類別之間並非互斥的。例如：一首歌曲能夠來源於影視原聲，同時也包含人聲。這種狀況下，使用4個二分類的 logistic 迴歸分類器更爲合適。這樣，對於每一個新的音樂做品，咱們的算法能夠分別判斷它是否屬於各個類別。

如今咱們來看一個計算視覺領域的例子，你的任務是將圖像分到三個不一樣類別中。

- (i) 假設這三個類別分別是：室內場景、戶外城區場景、戶外荒野場景。你會使用sofmax迴歸仍是 3個logistic 迴歸分類器呢？ —— 三個類別是互斥的，所以更適於選擇softmax迴歸分類器。
- (ii) 如今假設這三個類別分別是室內場景、黑白圖片、包含人物的圖片，你又會選擇 softmax 迴歸仍是多個 logistic 迴歸分類器呢？—— 創建三個獨立的 logistic迴歸分類器更加合適。

例三： Plot multinomial and One-vs-Rest Logistic Regression

""" ==================================================== Plot multinomial and One-vs-Rest Logistic Regression ==================================================== Plot decision surface of multinomial and One-vs-Rest Logistic Regression. The hyperplanes corresponding to the three One-vs-Rest (OVR) classifiers are represented by the dashed lines. """
print(__doc__) # Authors: Tom Dupre la Tour <tom.dupre-la-tour@m4x.org> # License: BSD 3 clause

import numpy as np import matplotlib.pyplot as plt from sklearn.datasets import make_blobs from sklearn.linear_model import LogisticRegression # make 3-class dataset for classification
centers = [[-5, 0], [0, 1.5], [5, -1]] X, y = make_blobs(n_samples=1000, centers=centers, random_state=40) transformation = [[0.4, 0.2], [-0.4, 1.2]] X = np.dot(X, transformation) for multi_class in ('', 'ovr'): clf = LogisticRegression(solver='sag', max_iter=100, random_state=42, multi_class=multi_class).fit(X, y) # print the training scores
    print("training score : %.3f (%s)" % (clf.score(X, y), multi_class)) # create a mesh to plot in
    h = .02  # step size in the mesh
    x_min, x_max = X[:, 0].min() - 1, X[:, 0].max() + 1 y_min, y_max = X[:, 1].min() - 1, X[:, 1].max() + 1 xx, yy = np.meshgrid(np.arange(x_min, x_max, h), np.arange(y_min, y_max, h)) # Plot the decision boundary. For that, we will assign a color to each
    # point in the mesh [x_min, x_max]x[y_min, y_max].
    Z = clf.predict(np.c_[xx.ravel(), yy.ravel()]) # Put the result into a color plot
    Z = Z.reshape(xx.shape) plt.figure() plt.contourf(xx, yy, Z, cmap=plt.cm.Paired) plt.title("Decision surface of LogisticRegression (%s)" % multi_class) plt.axis('tight') # Plot also the training points
    colors = "bry"
    for i, color in zip(clf.classes_, colors): idx = np.where(y == i) plt.scatter(X[idx, 0], X[idx, 1], c=color, cmap=plt.cm.Paired) # Plot the three one-against-all classifiers
    xmin, xmax = plt.xlim() ymin, ymax = plt.ylim() coef = clf.coef_ intercept = clf.intercept_ def plot_hyperplane(c, color): def line(x0): return (-(x0 * coef[c, 0]) - intercept[c]) / coef[c, 1] plt.plot([xmin, xmax], [line(xmin), line(xmax)], ls="--", color=color) for i, color in zip(clf.classes_, colors): plot_hyperplane(i, color) plt.show()