逻辑回归（Logistics Regression）

一、逻辑回归的概念

逻辑回归又称logistic回归分析，是一种广义的线性回归分析模型，常用于数据挖掘，疾病自动诊断，经济预测等领域。逻辑回归从本质来说属于二分类问题。
二分类问题是指预测的结果只有两种独立的情况，即$y$值只能取0或者1代表不同的结果，例如，一个人要么是左撇子，要么是右撇子，那么用0来表示左撇子，1来表示右撇子，我们就可以构造一个映射$X\to Y$，其中 X = { 所有人 } X = \{ \text{所有人} \} X={所有人}， Y = { 0 , 1 } Y = \{ 0, 1 \} Y={0,1}，这就是一个二分类问题。

1. 回顾线性回归

我们之前学过的线性回归模型是对于给定的特征$\mathbf{x}$，给定目标值$\mathbf{y}$，构造待求解的目标函数$h_{\theta }(\mathbf{x})={\mathbf{\theta }}^T\mathbf{x}$，其中$\mathbf{\theta }={\left(\begin{array}{cccc}{\theta }_0& {\theta }_1& \cdots & {\theta }_n\end{array}\right)}^T$，$\mathbf{x}={\left(\begin{array}{cccc}1& x_1& \cdots & x_n\end{array}\right)}^T$，通过各种方式（例如梯度下降法、共轭梯度法、牛顿法、$\cdots$）使误差函数$J(\mathbf{\theta })=\frac{1}{2}\sum _{i=1}^n{\left(h_{\theta }({\mathbf{x}}^{(i)})-y^{(i)}\right)}^2$取到最小值，从而使目标函数$h_{\theta }(\mathbf{x})$能够很好地刻画目标值$\mathbf{y}$，从而得到线性回归模型$h_{\theta }(\mathbf{x})$.

2. 为什么要用逻辑回归

那么既然已经有线性回归了，为什么我们还要研究逻辑回归？

我们来看这样一个例子：下图是肿瘤的大小与肿瘤是否是恶性的之间的关系。其中我们用$0$去表示良性肿瘤，用$1$去表示恶性肿瘤。

那么这时候我们如果用线性回归去建模，就会得到如下图所示的一条直线，其中当$h_{\theta }(\mathrm{t}\mathrm{u}\mathrm{m}\mathrm{o}\mathrm{r}\mathrm{S}\mathrm{i}\mathrm{z}\mathrm{e})>0.5$时，我们预测$\mathrm{m}\mathrm{a}\mathrm{l}\mathrm{i}\mathrm{g}\mathrm{n}\mathrm{a}\mathrm{n}\mathrm{t}=1$；当$h_{\theta }(\mathrm{t}\mathrm{u}\mathrm{m}\mathrm{o}\mathrm{r}\mathrm{S}\mathrm{i}\mathrm{z}\mathrm{e})<0.5$时，我们预测$\mathrm{m}\mathrm{a}\mathrm{l}\mathrm{i}\mathrm{g}\mathrm{n}\mathrm{a}\mathrm{n}\mathrm{t}=0$。

但是如果我们得到的训练集是如下图所示的，那么就会训练出来一条斜率略小的直线。

这时我们画出$\mathrm{m}\mathrm{a}\mathrm{l}\mathrm{i}\mathrm{g}\mathrm{n}\mathrm{a}\mathrm{n}\mathrm{t}=0.5$这条直线，从下图中我们可以有相当一部分恶性肿瘤如果按照肿瘤大小去划分，则被分为了良性肿瘤，导致了划分的不准确。

所以为了划分集的准确性，我们引入了逻辑回归。

二、逻辑回归

那么对于上面那个问题，线性回归会导致分类效果很差，所以我们需要一个可以很好拟合上面训练集的模型。

考虑Sigmoid函数$g(z)=\frac{1}{1+e^{-z}}$，其函数图形如下图所示。

我们发现在做二分类问题时，Sigmoid函数开始趋于平缓，然后急速上升，最后又趋于平缓，即Sigmoid处在趋近于0和趋近于1的时候比较多，能将$y=0$与$y=1$较为准确地划分。

例如，还是用按照肿瘤大小去判断是否为恶性肿瘤这个例子，如果我们用Sigmoid函数去拟合，就会得到如下图所示的结果。

这时如果我们再规定当$h_{\theta }(\mathrm{t}\mathrm{u}\mathrm{m}\mathrm{o}\mathrm{r}\mathrm{S}\mathrm{i}\mathrm{z}\mathrm{e})>0.5$时，我们预测$\mathrm{m}\mathrm{a}\mathrm{l}\mathrm{i}\mathrm{g}\mathrm{n}\mathrm{a}\mathrm{n}\mathrm{t}=1$；当$h_{\theta }(\mathrm{t}\mathrm{u}\mathrm{m}\mathrm{o}\mathrm{r}\mathrm{S}\mathrm{i}\mathrm{z}\mathrm{e})<0.5$时，我们预测$\mathrm{m}\mathrm{a}\mathrm{l}\mathrm{i}\mathrm{g}\mathrm{n}\mathrm{a}\mathrm{n}\mathrm{t}=0$，就会得到分类效果更好的模型。

我们把这种用Sigmoid函数去预测二分类问题的模型叫做逻辑回归模型。

1. “回归”的基本原理

1. 找一个合适的预测函数（hypothesis），一般表示为$h$函数，该函数就是我们需要找的分类函数，它用来预测输入数据的判断结果。这个过程时非常关键的，需要对数据有一定的了解或分析，知道或者猜测预测函数的“大概”形式，比如是线性函数还是非线性函数。

2. 构造一个$Cost$函数（损失函数），该函数表示预测的输出$h(\mathbf{x})$与训练数据类别$\mathbf{y}$之间的偏差，可以是二者之间的方差$(h(\mathbf{x})-\mathbf{y}{)}^2$或者是其他的形式。综合考虑所有训练数据的“损失”，将$Cost$求和或者求平均，记为$J(\mathbf{\theta })$函数，表示所有训练数据预测值与实际类别的偏差。

3. 显然，$J(\mathbf{\theta })$函数的值越小表示预测函数越准确（即$h(\mathbf{x})$对$\mathbf{y}$的预测越准确），所以这一步需要做的是找到$J(\mathbf{\theta })$函数的最小值。找函数的最小值有很多不同的方法，其中逻辑回归实现时梯度下降法、共轭梯度法和牛顿法居多。

2. 具体过程

（1）引言

下面我们来考察如下图这样的一个分类问题。现在我想将红色的叉与蓝色的圆圈划分开。

这时我们很容易就能想到用下图中的一条直线就可以将红色的叉与蓝色的圆圈划分开。

但是现在的问题是，我们如何找到这条直线呢？可能有多条直线都能将这两种点分开，哪条直线在面对接下来有可能还要加进来的数据的效果会更好呢？

于是就引入了逻辑回归算法。

（2）逻辑回归模型的构建

现在我们结合之前讲过的预测肿瘤为良性还是恶性的例子，我们发现之所以肿瘤那个例子那么好预测的原因是我们只需要根据0和1来判断肿瘤是良性的还是恶性的，所以现在面对这个划分点的问题，我们是不是也可以像肿瘤那样处理呢？

现在我们令红色的叉就是$y=1$，蓝色的圆圈就是$y=0$，如下图所示，进行划分。

那么对应过去的横坐标又应该是什么呢？我们想一下肿瘤那个例子的横坐标，是肿瘤的大小，换言之，是肿瘤的特征，那么在这个划分叉和圆圈的问题中，叉和圆圈的特征是什么呢？我想答案显而易见，就是这些叉和圆圈所在整个坐标轴中的位置。

所以我们现在想，如何刻画这些圆圈和叉的位置？答案应该是通过叉和圆圈的横纵坐标$x_1$和$x_2$，在这个问题里我们可以清楚地看到当$-3+x_1+x_2>0$时坐标平面里分布的全是红色的叉；当$-3+x_1+x_2<0$时坐标平面里分布的全是蓝色的圆圈。所以我们以$-3+x_1+x_2=0$为分界线，$-3+x_1+x_2$越小，分布就越是蓝色的圆圈，$y$值就越趋近于$0$，$-3+x_1+x_2$越大，分布就越是蓝色的圆圈，$y$值就越趋近于$1$，这样我们就可以得到一个如下图一样的映射。

于是我们就可以使用Sigmoid函数进行二分类问题的划分。

（3）构造预测函数

对于标准Sigmoid函数$g(z)=\frac{1}{1+e^{-z}}$，我们可以看到是以$z=0$为分界线的，所以我们现在令

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

其中

$$\mathbf{\theta }={\left(\begin{array}{cccc}{\theta }_0& {\theta }_1& \cdots & {\theta }_n\end{array}\right)}^T,\mathbf{x}={\left(\begin{array}{cccc}1& x_1& \cdots & x_n\end{array}\right)}^T$$

我们不妨设$x_0=1$，则

$${\mathbf{\theta }}^T\mathbf{x}=\sum _{i=0}^n{\theta }_i{x}_i={\theta }_0+{\theta }_1{x}_1+\cdots +{\theta }_n{x}_n$$

这样的话，我们就成功实现了以${\mathbf{\theta }}^T\mathbf{x}=0$为分界的模型，其概率表达为

$$\begin{array}{rl}P(y=1\mid \mathbf{x};\mathbf{\theta })& =h_{\theta }(\mathbf{x})\\ P(y=0\mid \mathbf{x};\mathbf{\theta })& =1-h_{\theta }(\mathbf{x})\end{array}$$

即对于真实标记是$1$的样本，我们希望预测值越接近于$1$，损失越小；对于真实标记是$0$的样本，我们希望预测值越接近于$0$，损失越小。

（4）构造损失函数$J(\mathbf{\theta })$

我们都知道，评价一个模型好坏用的是损失函数，损失函数越小，就能说明模型越好。

最直观的想法是希望预测值与实际值尽可能接近，即看预测值与实际值之间的均方误差是否最小，所以定义线性回归损失函数为

$$J(\mathbf{\theta })=\frac{1}{2m}\sum _{i=1}^m{\left(h_{\theta }({\mathbf{x}}^{(i)})-y^{(i)}\right)}^2=\frac{1}{2m}{\Vert h_{\theta }(\mathbf{x})-\mathbf{y}\Vert }_2^2$$

那么这个问题就变成了一个优化问题，即要找到令损失函数$J(\mathbf{\theta })$最小的$\mathbf{\theta }$，很自然的，我们需要使用梯度下降法。

但是很不巧的是，这个函数不是凸函数，它的函数图形如下图所示。

也就是说，这个函数我们如果使用梯度下降法，则很有可能在某一个局部凸函数的最小点停止，而不是整个函数的最小值点。所以我们不能使用均方差函数作为损失函数。

那么我们换个思路，跟据前面的推导，$y$只能取$0$或$1$，服从二项分布，则后验概率为

$$P(y\mid \mathbf{x};\mathbf{\theta })=h_{\theta }(\mathbf{x}{)}^y{\left(1-h_{\theta }(\mathbf{x})\right)}^{1-y}$$

在这个式子中，无论$y=0$还是$y=1$，$P(y\mid \mathbf{x};\mathbf{\theta })$的结果都是越趋近于1，模型预测能力越强。

那么我们就可以使用极大似然估计去求得$P(y\mid \mathbf{x};\mathbf{\theta })$的最大值，从而得到预测能力最好的模型。

对于$m$个独立同分布的训练样本$\mathbf{x}$，其似然函数$L(\mathbf{\theta })$为

$$\begin{array}{rl}L(\mathbf{\theta })& =\prod _{i=1}^{m}P(y^{(i)}\mid {\mathbf{x}}^{(i)};\mathbf{\theta })\\ & =\prod _{i=1}^m{h}_{\theta }({\mathbf{x}}^{(i)}{)}^{y^{(i)}}{\left(1-h_{\theta }({\mathbf{x}}^{(i)})\right)}^{1-y^{(i)}}\end{array}$$

取对数，则对数似然函数$l(\mathbf{\theta })$为

$$\begin{array}{rl}l(\mathbf{\theta })& =\log L(\mathbf{\theta })\\ & =\log \prod _{i=1}^{m}P(y^{(i)}\mid {\mathbf{x}}^{(i)};\mathbf{\theta })\\ & =\sum _{i=1}^m\log P(y^{(i)}\mid {\mathbf{x}}^{(i)};\mathbf{\theta })\\ & =\sum _{i=1}^m\log \left[h_{\theta }({\mathbf{x}}^{(i)}{)}^{y^{(i)}}{\left(1-h_{\theta }({\mathbf{x}}^{(i)})\right)}^{1-y^{(i)}}\right]\\ & =\sum _{i=1}^m\left[y^{(i)}\log h_{\theta }({\mathbf{x}}^{(i)})+\left(1-y^{(i)}\right)\log \left(1-h_{\theta }({\mathbf{x}}^{(i)})\right)\right]\end{array}$$

但是我们可以看出来，随着$h_{\theta }(\mathbf{x})$逐渐变得准确，$P(y\mid \mathbf{x};\mathbf{\theta })$的值越来越趋近于$1$，那么$l(\mathbf{\theta })$的值从$-\infty$逐渐变为0，为单调递增函数。

那么定义损失函数为

$$\begin{array}{rl}J(\mathbf{\theta })& =-\frac{1}{m}l(\mathbf{\theta })\\ & =-\frac{1}{m}\sum _{i=1}^m\left[y^{(i)}\log h_{\theta }({\mathbf{x}}^{(i)})+\left(1-y^{(i)}\right)\log \left(1-h_{\theta }({\mathbf{x}}^{(i)})\right)\right]\end{array}$$

就可以实现损失函数是凸函数（见下图），从而可以使用梯度下降法去求$\arg \underset{\mathbf{\theta }}{\min }J(\mathbf{\theta })$.

事实上，如果我们定义$Cost$函数为

$$Cost(h_{\theta }(x),y)=\{\begin{array}{ll}-\log h_{\theta }(x)& ,\text{if y=1}\\ -\log \left(1-h_{\theta }(x)\right)& ,\text{if y=0}\end{array}$$

这个函数的图形如下图所示

可以看到，随着$h(x)$的增大，$-\log h(x)$逐渐由$\infty$降到$0$；$-\log \left(1-h(x)\right)$逐渐由$0$增到$\infty$，我们把这个$Cost$函数叫做交叉熵。

那么就有

$$\begin{array}{rl}J(\mathbf{\theta })& =\frac{1}{m}\sum _{i=1}^{m}Cost(h_{\theta }({\mathbf{x}}^{(i)}),y^{(i)})\\ & =-\frac{1}{m}\sum _{i=1}^m\left[y^{(i)}\log h_{\theta }({\mathbf{x}}^{(i)})+\left(1-y^{(i)}\right)\log \left(1-h_{\theta }({\mathbf{x}}^{(i)})\right)\right]\end{array}$$

但是如果我们的模型选取的太复杂的话，会造成过拟合的后果，所以我们还需要加上惩罚项来避免过拟合，则损失函数变为

$$J(\mathbf{\theta })=-\frac{1}{m}\sum _{i=1}^m\left[y^{(i)}\log h_{\theta }({\mathbf{x}}^{(i)})+\left(1-y^{(i)}\right)\log \left(1-h_{\theta }({\mathbf{x}}^{(i)})\right)\right]+\frac{\lambda }{2m}{\Vert \mathbf{\theta }\Vert }_2^2$$

（5）梯度下降法找最小值

因为损失函数$J(\mathbf{\theta })$满足凸函数的性质，所以可以与线性回归一样，使用梯度下降法来求解$\arg \underset{\mathbf{\theta }}{\min }J(\mathbf{\theta })$.

即

$$\begin{array}{rl}{\theta }_j^{(n+1)}& ={\theta }_j^{(n)}-\alpha \frac{\partial J(\mathbf{\theta })}{\partial {\theta }_j^{(n)}}\\ & ={\theta }_j^{(n)}-\alpha \left[\frac{1}{m}\sum _{i=1}^m\left(h_{\theta }({\mathbf{x}}^{(i)})-y^{(i)}\right)x_j^{(i)}+\frac{\lambda }{m}{\theta }_j\right]\end{array}$$

其中

$$\begin{array}{rl}\frac{\partial J(\mathbf{\theta })}{\partial {\theta }_j}& =-\frac{1}{m}\sum _{i=1}^m\left(y^{(i)}\frac{1}{h_{\theta }({\mathbf{x}}^{(i)})}\frac{\partial h_{\theta }({\mathbf{x}}^{(i)})}{\partial {\theta }_j}-\left(1-y^{(i)}\right)\frac{1}{1-h_{\theta }({\mathbf{x}}^{(i)})}\frac{\partial h_{\theta }({\mathbf{x}}^{(i)})}{\partial {\theta }_j}\right)+\frac{\lambda }{m}{\theta }_j\\ & =-\frac{1}{m}\sum _{i=1}^m\left(y^{(i)}\frac{1}{g({\mathbf{\theta }}^T{\mathbf{x}}^{(i)})}-\left(1-y^{(i)}\right)\frac{1}{1-g({\mathbf{\theta }}^T{\mathbf{x}}^{(i)})}\right)\frac{\partial g({\mathbf{\theta }}^T{\mathbf{x}}^{(i)})}{\partial {\theta }_j}+\frac{\lambda }{m}{\theta }_j\\ & =-\frac{1}{m}\sum _{i=1}^m\left(y^{(i)}\frac{1}{g({\mathbf{\theta }}^T{\mathbf{x}}^{(i)})}-\left(1-y^{(i)}\right)\frac{1}{1-g({\mathbf{\theta }}^T{\mathbf{x}}^{(i)})}\right)g({\mathbf{\theta }}^T{\mathbf{x}}^{(i)})\left(1-g({\mathbf{\theta }}^T{\mathbf{x}}^{(i)})\right)\frac{\partial \left({\mathbf{\theta }}^T{\mathbf{x}}^{(i)}\right)}{\partial {\theta }_j}+\frac{\lambda }{m}{\theta }_j\\ & =-\frac{1}{m}\sum _{i=1}^m\left(y^{(i)}\left(1-g({\mathbf{\theta }}^T{\mathbf{x}}^{(i)})\right)-\left(1-y^{(i)}\right)g({\mathbf{\theta }}^T{\mathbf{x}}^{(i)})\right)x_j^{(i)}+\frac{\lambda }{m}{\theta }_j\\ & =-\frac{1}{m}\sum _{i=1}^m\left(y^{(i)}-g({\mathbf{\theta }}^T{\mathbf{x}}^{(i)})\right)x_j^{(i)}+\frac{\lambda }{m}{\theta }_j\\ & =-\frac{1}{m}\sum _{i=1}^m\left(y^{(i)}-h_{\theta }({\mathbf{x}}^{(i)})\right)x_j^{(i)}+\frac{\lambda }{m}{\theta }_j\\ & =\frac{1}{m}\sum _{i=1}^m\left(h_{\theta }({\mathbf{x}}^{(i)})-y^{(i)}\right)x_j^{(i)}+\frac{\lambda }{m}{\theta }_j\end{array}$$

这里用到复合函数积分，若$f(x)=\frac{1}{1+e^{-g(x)}}$，则

$$\begin{array}{rl}\frac{\mathrm{d}f(x)}{\mathrm{d}x}& =\frac{1}{{\left(1+e^{-g(x)}\right)}^2}{e}^{-g(x)}\frac{\mathrm{d}g(x)}{\mathrm{d}x}\\ & =\frac{1}{1+e^{-g(x)}}\frac{e^{-g(x)}}{1+e^{-g(x)}}\frac{\mathrm{d}g(x)}{\mathrm{d}x}\\ & =f(x)\left(1-f(x)\right)\frac{\mathrm{d}g(x)}{\mathrm{d}x}\end{array}$$

（6）向量化结果表示

上面已经推导出了梯度下降法的公式，现在我们对公式进行向量化。

记

$$\begin{array}{rl}\mathbf{x}& ={\left(\begin{array}{cccc}{\mathbf{x}}^{(1)}& {\mathbf{x}}^{(2)}& \cdots & {\mathbf{x}}^{(m)}\end{array}\right)}^T\\ \mathbf{y}& ={\left(\begin{array}{cccc}{y}^{(1)}& y^{(2)}& \cdots & y^{(m)}\end{array}\right)}^T\end{array}$$

其中

$${\mathbf{x}}^{(i)}={\left(\begin{array}{cccc}{x}_0^{(i)}& x_1^{(i)}& \cdots & x_n^{(i)}\end{array}\right)}^T,i=1,2,\cdots ,m$$

则

$$\begin{array}{rl}\frac{\partial J(\mathbf{\theta })}{\partial \mathbf{\theta }}& =\left(\begin{array}{c}\frac{\partial J(\mathbf{\theta })}{\partial {\theta }_0}\\ \frac{\partial J(\mathbf{\theta })}{\partial {\theta }_1}\\ ⋮\\ \frac{\partial J(\mathbf{\theta })}{\partial {\theta }_n}\end{array}\right)\\ & =\left(\begin{array}{c}\frac{1}{m}\sum _{i=1}^m\left(h_{\theta }({\mathbf{x}}^{(i)})-y^{(i)}\right)x_0^{(i)}+\frac{\lambda }{m}{\theta }_0\\ \frac{1}{m}\sum _{i=1}^m\left(h_{\theta }({\mathbf{x}}^{(i)})-y^{(i)}\right)x_1^{(i)}+\frac{\lambda }{m}{\theta }_1\\ ⋮\\ \frac{1}{m}\sum _{i=1}^m\left(h_{\theta }({\mathbf{x}}^{(i)})-y^{(i)}\right)x_n^{(i)}+\frac{\lambda }{m}{\theta }_n\end{array}\right)\\ & =\frac{1}{m}\left(\begin{array}{c}\sum _{i=1}^m\left(h_{\theta }({\mathbf{x}}^{(i)})-y^{(i)}\right)x_0^{(i)}\\ \sum _{i=1}^m\left(h_{\theta }({\mathbf{x}}^{(i)})-y^{(i)}\right)x_1^{(i)}\\ ⋮\\ \sum _{i=1}^m\left(h_{\theta }({\mathbf{x}}^{(i)})-y^{(i)}\right)x_n^{(i)}\end{array}\right)+\frac{\lambda }{m}\left(\begin{array}{c}{\theta }_0\\ {\theta }_1\\ ⋮\\ {\theta }_n\end{array}\right)\\ & =\frac{1}{m}\left(\begin{array}{c}\left(\begin{array}{cccc}{x}_0^{(1)}& x_0^{(2)}& \cdots & x_0^{(m)}\end{array}\right)\left(\begin{array}{c}{h}_{\theta }({\mathbf{x}}^{(1)})-y^{(1)}\\ h_{\theta }({\mathbf{x}}^{(2)})-y^{(2)}\\ ⋮\\ h_{\theta }({\mathbf{x}}^{(m)})-y^{(m)}\end{array}\right)\\ \left(\begin{array}{cccc}{x}_1^{(1)}& x_1^{(2)}& \cdots & x_1^{(m)}\end{array}\right)\left(\begin{array}{c}{h}_{\theta }({\mathbf{x}}^{(1)})-y^{(1)}\\ h_{\theta }({\mathbf{x}}^{(2)})-y^{(2)}\\ ⋮\\ h_{\theta }({\mathbf{x}}^{(m)})-y^{(m)}\end{array}\right)\\ ⋮\\ \left(\begin{array}{cccc}{x}_n^{(1)}& x_n^{(2)}& \cdots & x_n^{(m)}\end{array}\right)\left(\begin{array}{c}{h}_{\theta }({\mathbf{x}}^{(1)})-y^{(1)}\\ h_{\theta }({\mathbf{x}}^{(2)})-y^{(2)}\\ ⋮\\ h_{\theta }({\mathbf{x}}^{(m)})-y^{(m)}\end{array}\right)\end{array}\right)+\frac{\lambda }{m}\mathbf{\theta }\\ & =\frac{1}{m}\left(\begin{array}{cccc}{x}_0^{(1)}& x_0^{(2)}& \cdots & x_0^{(m)}\\ x_1^{(1)}& x_1^{(2)}& \cdots & x_1^{(m)}\\ ⋮& ⋮& \ddots & ⋮\\ x_n^{(1)}& x_n^{(2)}& \cdots & x_n^{(m)}\end{array}\right)\left(\begin{array}{c}{h}_{\theta }({\mathbf{x}}^{(1)})-y^{(1)}\\ h_{\theta }({\mathbf{x}}^{(2)})-y^{(2)}\\ ⋮\\ h_{\theta }({\mathbf{x}}^{(m)})-y^{(m)}\end{array}\right)+\frac{\lambda }{m}\mathbf{\theta }\\ & =\frac{1}{m}{\mathbf{x}}^T\left(\begin{array}{c}{\left(\begin{array}{c}{h}_{\theta }({\mathbf{x}}^T)\end{array}\right)}^T-\mathbf{y}\end{array}\right)+\frac{\lambda }{m}\mathbf{\theta }\end{array}$$

在这里，如果我们记

$$\begin{array}{rl}\mathbf{A}& =\mathbf{x}\cdot \mathbf{\theta }\\ & =\left(\begin{array}{cccc}{x}_0^{(1)}& x_1^{(1)}& \cdots & x_n^{(1)}\\ x_0^{(2)}& x_1^{(2)}& \cdots & x_n^{(2)}\\ ⋮& ⋮& \ddots & ⋮\\ x_0^{(m)}& x_1^{(m)}& \cdots & x_n^{(m)}\end{array}\right)\cdot \left(\begin{array}{c}{\theta }_0\\ {\theta }_1\\ ⋮\\ {\theta }_n\end{array}\right)\end{array}$$

则很容易地，我们可以做代换

$${\left(\begin{array}{c}{h}_{\theta }({\mathbf{x}}^T)\end{array}\right)}^T=h_{\theta }(\mathbf{A})$$

那么

$$\frac{\partial J(\mathbf{\theta })}{\partial \mathbf{\theta }}=\frac{1}{m}{\mathbf{x}}^T\left(\begin{array}{c}{h}_{\theta }(\mathbf{A})-\mathbf{y}\end{array}\right)+\frac{\lambda }{m}\mathbf{\theta }$$

所以在迭代的时候，我们有

$$\begin{array}{rl}{\mathbf{\theta }}^{(p+1)}& ={\mathbf{\theta }}^{(p)}-\alpha \frac{\partial J({\mathbf{\theta }}^{(p)})}{\partial {\mathbf{\theta }}^{(p)}}\\ & ={\mathbf{\theta }}^{(p)}-\frac{\alpha }{m}\left[{\mathbf{x}}^T\left(\begin{array}{c}{h}_{\theta }(\mathbf{A})-\mathbf{y}\end{array}\right)+\lambda {\mathbf{\theta }}^{(p)}\right]\end{array}$$

如此，变完成了梯度下降法的向量转化。

三、代码实现

1. 构造数据

本次构造数据采用二维高斯分布构造数据，构造满足朴素贝叶斯的数据，即协方差矩阵是对角矩阵的数据。

构造出来的数据如下图所示（满足朴素贝叶斯假设）。

2. 实现逻辑回归分类

按照上面所讲的算法进行分类，最终得到如下图所示的分类效果，$J(\mathbf{\theta })$经过训练，最终降到$0.000489596368502944$，$\mathbf{\theta }$的值为$\left(\begin{array}{ccc}-18.09239552& 2.61418602& 0.66697504\end{array}\right)$

3. 代码（Python实现）

import numpy as np
import matplotlib.pyplot as plt


# 构造符合独立性的二维高斯分布的数据点
def create_data(mu1, mu2, sigma1, sigma2, size):
    mean = np.array([mu1, mu2])
    cov = np.diag([sigma1, sigma2])
    f1 = np.random.multivariate_normal(mean, cov, size)
    X1 = f1[:, 0]
    X2 = f1[:, 1]
    return X1, X2


# Sigmoid函数
def sigmoid(z):
    return 1.0 / (1 + np.exp(-z))


# 计算损失函数
def loss(X, Y, theta, lamda, size):
    X = np.array(X)
    Y = np.array(Y)
    theta = np.array(theta)

    res = 0
    for i in range(size):
        z = np.dot(theta.T, X[i])
        res = res + Y[i] * np.log(sigmoid(z)) + (1 - Y[i]) * np.log(1 - sigmoid(z))
    res = -res / size
    res = res + lamda / (2 * size) * np.dot(theta.T, theta)
    return res


# 梯度下降法
def gradient_method(X, Y, theta, lamda, alpha, size):
    X = np.array(X)
    Y = np.array(Y)
    theta = np.array(theta)

    error1 = loss(X, Y, theta, lamda, size)
    times = 0
    while True:
        tmp = []
        for i in range(size):
            tmp.append(sigmoid(np.dot(theta.T, X[i]))-Y[i])
        tmp = np.array(tmp)
        gra = np.dot(X.T, tmp) + lamda * theta
        gra = gra / size
        theta = theta - alpha * gra
        times += 1
        if times == 1000:
            error2 = loss(X, Y, theta, lamda, size)
            if error1 - error2 < 0.000001:
                break
            error1 = error2
            times = 0
    print(error1)
    return theta


# 主函数
if __name__ == '__main__':
    X11, X12 = create_data(10, 8, 1, 2, 150)
    X21, X22 = create_data(1, 3, 2, 1.5, 150)

    X = []
    X0 = np.ones(300)
    X1 = np.concatenate([X11, X21], axis=0)
    X2 = np.concatenate([X12, X22], axis=0)
    X.append(X0)
    X.append(X1)
    X.append(X2)
    X = np.array(X).T

    Y1 = np.ones(150)
    Y2 = np.zeros(150)
    Y = np.concatenate([Y1, Y2], axis=0)

    theta = np.array([-18.0544146, 2.78338406, 0.79579366])

    theta = np.array(gradient_method(X, Y, theta, 0.0005, 0.001, 300))
    print(theta)
    x = np.linspace(-4, 15, 1000)
    y = [-theta[1] / theta[2] * a - theta[0] / theta[2] for a in x]

    plt.scatter(X11, X12, c='r', marker='+')
    plt.scatter(X21, X22, c='g', marker='+')
    plt.plot(x, y, color='b')
    plt.xlabel("X1")
    plt.ylabel("X2")
    plt.show()