梯度下降法解逻辑斯蒂回归

本文是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的参数做正则化，但是即使做了影响也不大：

有关正则化的详细讨论可以见我爱公开课的笔记。