379      阿裏雲  技術社區[雲棲]

邏輯回歸那些事—使用牛頓法解決實際問題

更多深度文章，請關注雲計算頻道：https://yq.aliyun.com/cloud

前言

在本篇博客中，我們要介紹的是牛頓法的原理，並且將之應用到實際的邏輯回歸問題中。邏輯回歸的主要知識點包括伯努利分布的對數似然和用來平滑的sigmoid函數。
我們還要介紹Hession，這是一個二階偏導的方陣。看完了本片博客，您就知道如何使用Hession結合梯度來實現牛頓法。
和之前的博客一樣，我們這篇也將從牛頓法的整體概述、數學推導以及編程實現幾個方麵展開。最終將理論和實踐的結合，靈活運用牛頓法解決邏輯回歸問題。

推薦的先驗知識

• 導數和鏈式法則（微積分）
• 偏導數和梯度（多元微積分）
• 基本的向量操作（線性代數）
• 對Numpy的基本了解
• 獨立概率的基本認識

數據

我們的本次實驗采用的數據集是南波士頓的房產數據集。數據集中的每一行包含了一座房屋的價值，並且用一個布爾類型的數值表示是否有超過兩個的浴室。

$$\widehat{x} =HomeValue = <550000.00,600000.00,...578000.00>^{T}$$
$$\widehat{y} = MoreThan2Bathrooms = <1,0,0,...1>^{T}$$

south_boston_bedroom_logreg_raw

數學模型

首先，我們需要掌握基本的線性回歸模型。當我們給出房屋的價格的時候，該模型可以用二元分類器的方式為我們輸出是否該房屋有兩個以上的浴室。
接下來我們需要解決的是特征和權重的線性組合問題。我們需要將這個線性組合封裝在一個函數中，而這個函數必須是平滑的，且值域是[0,1]。
按照慣例，我們選擇滿足以上要求的Logistic函數，它的公式如下所示：

$$h(x)=1/(1+e^{-z}),$$

$$z=\theta_{1} x+ \theta_{2} $$

Note:我們增加了\(\theta_{2}\)參數來增加模型的靈活度，雖然我們的數據隻有一維，但是我們卻用了一個二維的模型。

sigmoid_function

在線性回歸中，我們需要定義一個平方和的目標函數，邏輯回歸也是它相似。我們用我們的假設函數h(x)來生成似然函數，最終使用似然函數來求得合適的權重。以下是具體的數學推導。

數學原理：定義似然函數

首先我們要定義一個概率質量函數：
$$P(y=1|x;\theta)=h_{\theta}(x)$$

$$P(y=0|x;\theta)=1-h_{\theta}(x)$$

Note:
上邊的第一個公式的等式左邊讀作“在已知參數\(\theta\)和特征向量x的條件下，y=1的概率”，我們使用右邊的假設函數來計算這個概率。

這兩個公式也可以合並成一個：

$$P(y|x;\theta)=h_{\theta}(x)^{y}(1-h_{\theta}(x)^{1-y}$$

下表顯示了我們假設函數的輸出與預測概率的關係。

很明顯，我們需要最大化等式的右邊，因此我們需要引申出似然函數的概念。
以我的理解，似然函數就是表示給定特征向量\(\widehat{x}\)能正確預測出y的可能性。當然，從統計的角度講，似然和概率是不一樣的，這裏給出了詳細的解釋。
現在我們可以詳細的談一談似然函數的事情啦！在這個項目中，我們應用在訓練集中的似然函數很簡單，我們將每一種情況的可能性乘起來，這樣就得到了我們模型正確預測的累積似然，公式如下。
$$L(\theta) = \prod_{n}^{i=1}p(y_{i}|x_{i};\theta)$$
更通用的，可以用以下的公式表示：
$$L(\theta) = h_{\theta}(x_{i})^{y_{i}}(1-h_{\theta}(x_{i}))^{1-y_{i}}$$
顯然，我們將n個概率相乘了，而相乘的結果必然是處於[0,1]區間的。n表示樣本的個數，所以最終的概率可能是個\(10^{-n}\)數量級的一個概率。這簡直就是一個糟糕的消息，因為計算的過程中可能會損失精度，而且如果用python來計算的話，我們最終看到的結果很有可能被四舍五入成為0。
我們的解決辦法是對似然函數求對數，這樣一來，連乘的形式就變成了連加，公式如下：
$$\ell (\theta)=log(L(\theta))$$
$$\ell (\theta)=\sum_{i=1}^{n}y_{i}log(h_{\theta}(x_{i}))+(1-y_{i})log(1-h_{\theta})$$

Note:
這裏給出對數函數的基本運算法則
$$log(xy)=log(x)+log(y)$$
$$log(x)^{n}=nlog(x)$$

這個函數被稱為我們假設函數的對數似然函數，因為對數函數具有單調性，因此不影響我們下麵的計算。
記住，當數據非常小的時候，我們之前假設的函數會變得非常不精確（因為數值太小，連乘運算之後誤差很大），而對數似然就顯現出良好的精度，對數似然的圖如下所示，而接下來我們要做的就是最大化對數似然函數。

log_likelihood_normalized

數學原理：單變量牛頓法

在我們求對數似然函數的最大值之前，我們先來介紹一下牛頓法

牛頓法使用迭代的思想，它是求多項式函數根的一種算法。在簡單的單變量情況下，牛頓法的實現步驟如下：

1. 在\(p(x_{n},y_{n})\)處求出f(x)的切線\($y={f}'(x_{n})(x-x_{n})+f(x_{n})$\)
2. 計算出切線在x軸上的截距 $$0={f}'(x_{n})(x-x_{n})+f(x_{n})$$ $$-f(x_{n})={f}'(x_{n})(x-x_{n})$$ $$x_{n+1}=x_{n}-\frac{f(n)}{{f}'(n)}$$
3. 求出該位置的y值 $$y_{n+1}=f(x_{n+1})$$
4. 如果\(y_{n+1}\approx y{n}\)恭喜你，我們的函數收斂了
5. 否則繼續更新\(p(x_{n},y_{n})\) \(x=x_{n+1},y=y_{n+1}\),返回第一步迭代 下麵我們使用維基百科上的gif圖來更直觀的理解這一過程：

NewtonIteration_Ani

如果你理解了以上的步驟，我們可以把這一過程歸納為：
$$x_{n+1}=x_{n}-\frac{f(n)}{{f}'(n)}$$ （直到\(x_{n}-x_{n+1}\approx 0\)）

任何高中學習過微積分的人都能理解上麵的內容。但是我們如何推廣到多元的，也就是“n維”的情況呢？

數學原理：n維牛頓法

在多維的情況下，我們使用偏導數的向量，也就是梯度來代替一維情況下的導數。
如果你對梯度的概念有疑問，那麼請看這篇帖子複習一下。

相應的，我們的更新規則也發生了改變，具體的公式如下所示：
$$x_{n+1}=x_{n}-\frac{f(n)}{\triangledown {f}'(n)}$$

Note從一維推廣到n維，其實就是將導數換成了梯度，也就是將\({f}'(n)$變成\triangledown {f}'(n)\)

數學原理：用牛頓法最大化對數似然

還記得我們之前提到的，我們的目的是最大化對數似然函數，也就是最大化\(\ell(\theta)\)函數，因為log函數是單調上漲的，因此，隻要保證對數似然函數是最大的，那麼我們假設的函數\(H_{\theta}(x)\)也就達到了最大值。

通常，求最大值的方式是對函數求偏導數，並讓其等於0，求解\(\theta_{1}和\theta_{2}\)來算出極值,極值也就是我們要找的最大值了。
Note:我們的log函數是遞增的，因此我們隻有一個極值點，而本方法就是求這個極值點的唯一方法。
這個聽起來是不是非常的熟悉，接下來，我們就要求讓偏導數為0的\(\theta_{1},\theta_{2}\)兩個參數的值了，這正是牛頓法解決的問題，讓我們一起再來回顧下：
$$x_{n+1}=x_{n}-\frac{f(n)}{\triangledown {f}'(n)}$$

我們用梯度來替換$f(x_{n})$,$\triangledown \l(\theta)$，這樣我們的公式就會變成
$$x_{n+1}=x_{n}-\frac{\triangledown \ell(\theta)}{?}$$

?表示什麼？直覺告訴我們，我們需要取梯度向量的導數，就像我們以前取f(x)的導數一樣。
接下來我們就需要Hessian矩陣

數學原理：Hessian矩陣

根據以前學的多元微積分的知識，我們應該知道，為了找到函數的二階偏導數，我們取每一階偏導數的偏導數。如果我們有n個參數，然後我們就會得到n×n階導數。
因此，Hessian矩陣是n*n階偏導數的平方矩陣。
在我們的案例中，我們隻有兩個參數\(\theta_{1},\theta_{2}\)，我們的Hessian矩陣就是這樣:
$$
H_{\ell(\theta)}=\begin{bmatrix} \frac{\partial^{2} \ell }{\partial {\theta_{1}}^{2}}& \frac{\partial^{2} \ell }{\partial \theta_{1}\partial \theta_{2}} \ \frac{\partial^{2} \ell }{\partial \theta_{2}\partial \theta_{1}} & \frac{\partial^{2} \ell }{\partial {\theta_{2}}^{2}} \end{bmatrix}\quad
$$

數學原理：整合

在我們的通過將Hessian矩陣替換進牛頓法的更新步驟，那麼我們的公式就變成：
$$\theta_{n+1}=\theta_{n}+H_{\ell(\theta)}^{-1}\bigtriangledown \ell(\theta)$$
Note在這裏我們取Hessian矩陣的逆，而不是取它的倒數，因為它是一個矩陣。

為了簡潔起見，這篇文章省去了梯度和Hessian矩陣的數學推導。如果你想了解下麵推導的資源可以在：
1. 對數似然函數的梯度的推導，參考Andrew Ng的講義17-18頁。
2. Hessian矩陣的推導是非常直白的，當你計算出梯度的時候，求Sigmoid函數導數也可以參考Andrew Ng的講義
\(\ell(\theta)\)的梯度是
$$
\begin{bmatrix}
\sum_{i=1}^{n}(y_{i} - h_{\theta}(x_{i}))x_{i} \
\sum_{i=1}^{n}(y_{i} - h_{\theta}(x_{i}))
\end{bmatrix}
$$
$\ell(\theta)$的Hessian矩陣是
$$H_{\ell(\theta)}=\begin{bmatrix}
\sum_{i=1}^{n}h_{\theta}(x_{i})(1-h_{\theta}(x_{i}))\theta_{1}\theta_{1} & \sum_{i=1}^{n}h_{\theta}(x_{i})(1-h_{\theta}(x_{i}))\theta_{1} \ \sum_{i=1}^{n}h_{\theta}(x_{i})(1-h_{\theta}(x_{i}))\theta_{1} & \sum_{i=1}^{n}h_{\theta}(x_{i})(1-h_{\theta}(x_{i})
\end{bmatrix}\quad
$$
其中\(h_{\theta}(x) = \frac{1}{1 + e^{-z}} z = \theta_{1}x + \theta_{2}\)

牛頓法實現

我們首先定義我們的假設函數為Sigmoid函數。

def sigmoid(x, Θ_1, Θ_2):                                                        
    z = (Θ_1*x + Θ_2).astype("float_")                                              
    return 1.0 / (1.0 + np.exp(-z))      

然後我們定義\(\ell(\theta)\),也就是我們的對數似然函數

def log_likelihood(x, y, Θ_1, Θ_2):                                                                
    sigmoid_probs = sigmoid(x, Θ_1, Θ_2)                                        
    return np.sum(y * np.log(sigmoid_probs)
                  + (1 - y) * np.log(1 - sigmoid_probs))

最後，我們實現了對數似然函數的梯度和Hessian矩陣。

def gradient(x, y, Θ_1, Θ_2):                                                         
    sigmoid_probs = sigmoid(x, Θ_1, Θ_2)                                        
    return np.array([[np.sum((y - sigmoid_probs) * x),                          
                     np.sum((y - sigmoid_probs) * 1)]])                         

def hessian(x, y, Θ_1, Θ_2):                                                          
    sigmoid_probs = sigmoid(x, Θ_1, Θ_2)                                        
    d1 = np.sum((sigmoid_probs * (1 - sigmoid_probs)) * x * x)                  
    d2 = np.sum((sigmoid_probs * (1 - sigmoid_probs)) * x * 1)                  
    d3 = np.sum((sigmoid_probs * (1 - sigmoid_probs)) * 1 * 1)                  
    H = np.array([[d1, d2],[d2, d3]])                                           
    return H

我們在這4個函數中實現了所有的數學公式，我們創建了一個While循環，並使用牛頓的方法迭代，直到收斂到最大值為止。

def newtons_method(x, y):                                                             
    """
    :param x (np.array(float)): Vector of Boston House Values in dollars
    :param y (np.array(boolean)): Vector of Bools indicting if house has > 2 bedrooms:
    :returns: np.array of logreg's parameters after convergence, [Θ_1, Θ_2]
    """

    # Initialize log_likelihood & parameters                                                                   
    Θ_1 = 15.1                                                                     
    Θ_2 = -.4 # The intercept term                                                                 
    Δl = np.Infinity                                                                
    l = log_likelihood(x, y, Θ_1, Θ_2)                                                                 
    # Convergence Conditions                                                        
    δ = .0000000001                                                                 
    max_iterations = 15                                                            
    i = 0                                                                           
    while abs(Δl) > δ and i < max_iterations:                                       
        i += 1                                                                      
        g = gradient(x, y, Θ_1, Θ_2)                                                      
        hess = hessian(x, y, Θ_1, Θ_2)                                                 
        H_inv = np.linalg.inv(hess)                                                 
        # @ is syntactic sugar for np.dot(H_inv, g.T)¹
        Δ = H_inv @ g.T                                                             
        ΔΘ_1 = Δ[0][0]                                                              
        ΔΘ_2 = Δ[1][0]                                                              

        # Perform our update step                                                    
        Θ_1 += ΔΘ_1                                                                 
        Θ_2 += ΔΘ_2                                                                 

        # Update the log-likelihood at each iteration                                     
        l_new = log_likelihood(x, y, Θ_1, Θ_2)                                                      
        Δl = l - l_new                                                           
        l = l_new                                                                
    return np.array([Θ_1, Θ_2]) 

可視化：牛頓法

讓我們用可視化的方式看一看，對數似然函數的表麵上繪製牛頓法的每一次迭代會發生什麼：

newtons_method_ll

Note:第一次迭代我們用紅色來標識，第二次迭代我們用橘黃色來標識...最後一次迭代的顏色是紫色

newtons_method_max

這個圖形可以很直觀的看到我們的“紫色”點是最大值的點，這就說明，我們方法已經收斂了！

可視化：我們的方案

我們擁有的數據集其實是一個一維數據集，通常來說，一維數據的可視化是將這些點放置在線上，然後在線上做一些劃分的邊界。可問題是我們數據集裏的點太多，如果畫在一條線上，數據點會串在一起。
所以，我就將這些點沿著y軸給拉了出來，並用不同的顏色表示不同的類別。另外，我還在圖中做了三條輔助線來按百分比劃分這些點。

south_boston_bedroom_logreg

結論

本篇已經涵蓋了許多新的主題，比如Hessian矩陣、對數似然函數、Sigmoid函數。結合這些方法，我們已經能夠實現牛頓的方法來解決logistic回歸。
盡管我們的方法非常的具體，但是需要擔心的是我們的牛頓法也有發散的可能，當然，這個超出了本文的討論範圍，如果感興趣，可以查閱一些相關的資料。

作者介紹: Sean Harrington 全棧工程師、數據科學家，熟悉Python, Javascript, C, C++, Bash等語言

個人主頁：https://thelaziestprogrammer.com/author/sharrington/

以上為譯文

本文由北郵@愛可可-愛生活 老師推薦，阿裏雲雲棲社區組織翻譯。
文章原標題《Solving Logistic Regression with Newton's Method》，作者：Sean Harrington，譯者：愛小乖，審校：主題曲。

文章為簡譯，更為詳細的內容，請查看原文。

本文由用戶為個人學習及研究之目的自行翻譯發表，如發現侵犯原作者的版權，請與社區聯係處理yqgroup@service.aliyun.com

最後更新：2017-09-12 23:32:49

  上一篇：  Logtail技術分享(一) : Polling + Inotify 組合下的日誌保序采集方案

  下一篇：  log4j2配置文件log4j2.xml解析