梯度下降法解邏輯斯蒂回歸

本文是Andrew Ng在Coursera的機器學習課程的筆記。

Logistic回歸屬於分類模型。回顧線性回歸，輸出的是連續的實數，而Logistic回歸輸出的是[0,1]區間的概率值，通過概率值來判斷因變量應該是1還是0。因此，雖然名字中帶着“回歸”（輸出范圍常為連續實數），但Logistic回歸屬於分類模型（輸出范圍為一組離散值構成的集合）。

整體步驟

假如我們的自變量是“數學課和英語課的成績”，x={x1,x2}，因變量是“能否被哥大錄取”，y∈{0,1}。我們要通過這兩個自變量的分布，來預測因變量的值。Logistic回歸的步驟為：

1. 設定擬合函數（hypothesis function）：h_θ(x)，其意義是給定參數θ，根據輸入x，給出輸出h_θ(x)，當輸出值大於0.5時預測錄取，否則預測被拒。
2. 設定代價函數（cost function）：J(θ)，其意義是累加所有樣本的 預測結果h_θ(x) 與 真實結果y 之間的差距。
3. 利用梯度下降法，來調整參數θ，使得代價函數J(θ)的值最小。

比較線性回歸與Logistic回歸，可以看出二者非常相似，但是Logistic回歸的擬合函數（步驟一）和代價函數（步驟二）的定義方法與線性回歸有所不同。

Step 1：擬合函數

線性回歸的擬合函數為：h_θ(x) = θ^Tx，輸出范圍為所有實數，而其因變量的取值范圍也確實屬於所有實數。但是Logistic回歸的最終輸出要么是0，要么是1，我們不能直接套用線性回歸的擬合函數。對於Logistic回歸，我們需要控制輸出在[0,1]之間，因此借助函數g:

\( g(z)=\frac{1}{1+e^{-z}} \)

函數g為S型函數（Sigmoid function），也稱為Logistic function，“Logistic回歸”就是得名於此。最終的擬合函數為：

\( h_\theta(x)=g(\theta^{T}x)=\frac{1}{1+e^{-\theta^{T}x}} \)

這個擬合函數的輸出范圍在[0,1]之間，表示分類結果為1的可能性。例如，我輸入我的成績，得到的擬合函數輸出值為0.7，就表示我有70%的概率被哥大錄取（30%的概率被拒）。當輸出值超過0.5，我們將其分類為1（這表示模型最終預測我會被哥大錄取）。值為0.5的線稱為“Decision Boundary”（可以是曲線）。

想象一個三維坐標系(x1,x2,y)，對於任意的地面坐標(x1,x2)，都有唯一的y值與之對應。首先計算 θ^Tx，值可正可負且沒有邊界。然后將其作為S型函數g的輸入，得到的輸出固定在[0,1]之間。當 θ^Tx≥0時，h≥0.5，預測為1分類，否則為0分類。擬合函數的意義就在於將值固定在0到1之間。

Step 2：代價函數

如果直接套用線性回歸的代價函數，那么得到的代價函數將非凸（non-convex），利用梯度下降我們可能停留在局部最小值，而不是我們想要的全局最小值。因此我們需要重新定義代價函數。

對於一個樣本，我們新的代價函數定義如下：

公式可以進一步化簡為：

\( Cost(h_\theta(x),y) = -y log(h_\theta(x)) - (1-y) log(1-h_\theta(x)) \)

下圖是代價函數曲線，橫坐標為h的值，縱坐標為代價。可以看出，在y=1的前提下，如果預測值越接近1，那么相應代價就越小，反之則越大。

將所有m個樣本的代價累加並平均，我們有最終的代價函數：

這個代價函數滿足convex的條件，所以有全局最小值。其來源與最大似然估計有關，此處略去細節。

Step 3：梯度下降

我們采用與線性回歸中一樣的梯度下降法來確定θ的值，即設置一個合適的學習率α之后，同步更新所有j=1 to n：

 重復更新步驟，直到代價函數的值收斂為止。

高級操作

我們在第三步使用的梯度下降法雖然可行，但是收斂速度比較慢。有不少高級梯度下降算法已經被提出，包括 Conjugate gradient、BFGS、L-BFGS 等等。這些算法的優點是不需要手動挑選學習率，速度也較快，但是缺點就是比較復雜，難以手動實現。不過，借助matlab我們就可以利用這些算法來計算了。我們可以利用matlab中的 fminunc函數(Find minimum of unconstrained multivariable function) 來實現高級操作。

options = optimset(‘GradObj’, ‘on’, ‘MaxIter’, ‘100’);
initialTheta = zeros(2,1);
[optTheta, functionVal, exitFlag] = fminunc(@costFunction, initialTheta, options);

fminunc函數 的：

• 第一個參數是一個指向代價函數的指針，此代價函數需要返回代價值與各個參數的偏導數；
• 第二個參數是θ的初始值；
• 第三個參數是一些開關選項。

matlab代碼

輸入是ex2data1.txt文件，前兩列是自變量，最后一列是因變量，用逗號分隔，格式類似於：

34.62365962451697,78.0246928153624,0

30.28671076822607,43.89499752400101,0

……

這個代碼用兩種方法來進行梯度下降，一種是常規方法，一種是高級方法。實驗中發現常規方法經過了很長時間（大約好幾分鍾）還是沒有收斂，且需要考慮學習率的大小；而高級方法只需要幾秒鍾就能收斂，並且不需要考慮學習率。

function logisticRegressionDemo(  )
%A demo of logistic regression
%   Logistic regression using traditional gradient descent and advanced methods

%% Load Data
data = load('ex2data1.txt');
X = data(:, [1, 2]); y = data(:, 3);

%% ============ Part 1: Compute Cost and Gradient ============
%  Setup the data matrix appropriately, and add ones for the intercept term
[m, n] = size(X);

% Add intercept term to x and X_test
X = [ones(m, 1) X];

% Initialize fitting parameters
initial_theta = zeros(n + 1, 1);

% Compute and display initial cost and gradient
[cost, grad] = costFunction(initial_theta, X, y);

theta = initial_theta;

iterNum = 50000;
costList = zeros(iterNum,1);

% this method fails because it doesn't converge even I set iteration number to 50000, which takes quite a long time
for i=1:iterNum
    theta = theta - 0.001 * grad;
    [cost, grad] = costFunction(theta, X, y);
    costList(i,1)= cost;
end;
plot(costList);

% ======================= Predict  ============================
prob = sigmoid([1 45 85] * theta);
fprintf(['For a student with scores 45 and 85, we predict an admission ' ...
         'probability of %f\n\n'], prob);

%% ============= Part 2: Optimizing using fminunc  =============
%  Set options for fminunc
options = optimset('GradObj', 'on', 'MaxIter', 400);

%  Run fminunc to obtain the optimal theta
%  This function will return theta and the cost 
[theta, cost] = fminunc(@(t)(costFunction(t, X, y)), initial_theta, options);

% ======================= Predict  ============================
prob = sigmoid([1 45 85] * theta);
fprintf(['For a student with scores 45 and 85, we predict an admission ' ...
         'probability of %f\n\n'], prob);

end

function g = sigmoid(z)
%SIGMOID Compute sigmoid functoon
%   J = SIGMOID(z) computes the sigmoid of z.

g = 1 ./ (1 + exp(-z));
end

function [J, grad] = costFunction(theta, X, y)
%COSTFUNCTION Compute cost and gradient for logistic regression
%   J = COSTFUNCTION(theta, X, y) computes the cost of using theta as the
%   parameter for logistic regression and the gradient of the cost
%   w.r.t. to the parameters.

% Initialize some useful values
m = length(y); % number of training examples

J = 0;
grad = zeros(size(theta));

for i = 1:m
    J = J + (-y(i)) * log(sigmoid(theta' * X(i,:)')) - (1-y(i)) * log(1-sigmoid(theta' * X(i,:)'));
    grad = grad + (sigmoid(theta' * X(i,:)') - y(i)) * X(i,:)';
end
J = J / m;
grad = grad ./ m;

end

正則化：防止過擬合

我們采用正則化（Regularization）的方法，通過修改代價函數來防止過擬合。

過擬合問題一般歸咎於過多的特征類型。有兩種方法來減少過擬合：

1. 丟掉一些特征，不過這樣也丟失了一些信息；
2. 正則化，修改代價函數，來限制參數值的大小

正則化除了能防止過擬合之外，還有一個好處就是可以避免利用normal equation一步到位地求得參數值中需要考慮的矩陣可逆問題，因為加入正則參數之后的矩陣總是可逆的。

修改之后的代價函數和梯度方程見下，出於習慣我們不對j=0的參數做正則化，但是即使做了影響也不大：

有關正則化的詳細討論可以見我愛公開課的筆記。