一文彻底理解逻辑回归:从公式推导到代码实现

3,587 阅读14分钟

这是我参与2022首次更文挑战的第4天,活动详情查看:2022首次更文挑战

前言

大家好,这里是 Codeman 。这是本文的第二次修订,在前文的基础之上,这次我增加了很多公式的推导,从数学原理到代码实现,本文提供了一站式服务,希望能帮助读者从根上理解逻辑回归,一劳永逸地解决问题!

多次修订,只源于我精益求精的人生态度。写博客,我是认真的!如果你觉得本文确实对你有帮助,请点个赞支持我一下吧 🌹

正文

逻辑回归在社会和自然科学中应用非常广泛,它其实是一种统计学习方法,因为它的底层方法就是线性回归。它们不同的地方在于:线性回归适用于解决回归问题,而逻辑回归适用于解决分类问题。这是为什么呢?别急,看完本文你就懂了。因此,我们可以说逻辑回归是基于回归的伪回归算法!

image.png

在自然语言处理中,逻辑回归是用于分类的基线监督机器学习算法,它与神经网络也有非常密切的关系,神经网络可以被视为一系列堆叠在一起的逻辑回归分类器。

由于后文涉及到很多数学公式的推导,因此我们先做一个符号约定:

符号含义
𝑚𝑚训练集中样本的数量
𝑛𝑛特征的数量
𝑥𝑥特征/输入变量
𝑦𝑦目标变量/输出变量
(𝑥,𝑦)(𝑥,𝑦)训练集中的样本
(𝑥(𝑖),𝑦(𝑖))(𝑥^{(𝑖)},𝑦^{(𝑖)})𝑖𝑖 个观察样本
𝑥j(i)𝑥_j^{(i)}𝑖𝑖 个观察样本的第 jj 个特征
h学习算法的解决方案或函数,也称为假设(hypothesis)
𝑦^=h(𝑥)\widehat{𝑦}=ℎ(𝑥)预测值

我们先从线性回归开始讲起,一步一步地领略逻辑回归的全貌。

1.线性回归

线性回归是一种使用特征属性的线性组合来预测响应的方法。它的目标是找到一个线性函数,以尽可能准确地描述特征或自变量(xx)与响应值(yy)之间的关系,使得预测值与真实值之间的误差最小化。

image.png

1.1 数学定义

在数学上,线性回归要找的这个线性函数叫回归方程,其定义如下:

h(x(i))=β0+β1x(i)(1.1)h(x^{(i)})=\beta_0+\beta_1x^{(i)}\tag{1.1}

其中,h(x(i))h(x^{(i)})表示第 ii 个样本的预测响应值。b0b_0b1b_1 是回归系数,分别代表回归线的 yy 轴截距和斜率。这种形式通常见于特征只有单个属性的时候,也就是一元线性回归,我们初高中所学的就是这种。

在机器学习中,通常每个样本都有 nn 个特征属性,每个特征 xix_i 都有一个对应的权值 wiw_i,此时我们需要的就是多元线性回归:

h(𝑥)=𝑤0x0+𝑤1𝑥1+𝑤2𝑥2+...+𝑤𝑛𝑥𝑛=wTX(1.2)ℎ(𝑥)=𝑤_0x_0 +𝑤_1𝑥_1 +𝑤_2𝑥_2 +...+𝑤_𝑛𝑥_𝑛=w^TX\tag{1.2}

其中,x0=1x_0=1 没有实义,只是为了方便写成矩阵的形式,w0w_0 则等价于式 (1.1) 中的 β0\beta_0,把 β0\beta_0 融入矩阵中,不仅为了看起来简洁,也是为了方便计算。

损失函数采用平方和损失:

𝑙(𝑥(𝑖))=12(h(𝑥(𝑖))𝑦(𝑖))2(1.3)𝑙(𝑥^{(𝑖)})=\frac{1}{2}(ℎ(𝑥^{(𝑖)})−𝑦^{(𝑖)})^2\tag{1.3}

代价函数定义如下:

J(w)=12𝑖=1𝑚(h(𝑥(𝑖))𝑦(𝑖))2(1.4)J(w)=\frac{1}{2}\sum_{𝑖=1}^𝑚(ℎ(𝑥^{(𝑖)})−𝑦^{(𝑖)})^2\tag{1.4}

损失函数(Loss Function)度量单样本预测的误差,代价函数(Cost Function)度量全部样本的平均误差。常用的代价函数有均方误差、均方根误差、平均绝对误差等。

损失函数的系数 1/2 是为了便于计算,使对平方项求导后的常数系数为 1。

我们的目标是要找到一组 𝑤(𝑤0,𝑤1,𝑤2,...,𝑤𝑛)𝑤(𝑤_0,𝑤_1,𝑤_2,...,𝑤_𝑛),使得代价函数 J(w)J(w) 最小,即最小化 𝐽(𝑤)𝑤\frac{\partial {𝐽(𝑤)}}{\partial 𝑤}

下面我们将详细描述用最小二乘法求 𝑤𝑤 的推导过程。

1.2 最小二乘法

为了方便叙述,我们将J(w)J(w)用矩阵形式表达,即:

J(𝑤)=12(𝑋𝑤𝑌)2=12(𝑋𝑤𝑌)𝑇(𝑋𝑤𝑌)(1.5)J(𝑤)=\frac{1}{2}(𝑋𝑤−𝑌)^2=\frac{1}{2}(𝑋𝑤−𝑌)^𝑇(𝑋𝑤−𝑌)\tag{1.5}

其中,𝑋𝑋𝑚𝑚𝑛+1𝑛+1 列的矩阵(第一列全为 11,即式 (1.2) 中的 x0x_0),𝑤𝑤𝑛+1𝑛+111 列的矩阵(包含了 𝑤0𝑤_0 ),𝑌𝑌𝑚𝑚11 列的矩阵。

X=[1x1(1)x2(1)x3(1)xn(1)1x1(2)x2(2)x3(2)xn(2)1x1(m)x2(m)x3(m)xn(m)],w=[w0w2wn],Y=[y(1)y(2)y(m)]X=\left[\begin{array}{cccccc} 1 & x_{1}^{(1)} & x_{2}^{(1)} & x_{3}^{(1)} & \ldots & x_{n}^{(1)} \\ 1 & x_{1}^{(2)} & x_{2}^{(2)} & x_{3}^{(2)} & \ldots & x_{n}^{(2)} \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ 1 & x_{1}^{(m)} & x_{2}^{(m)} & x_{3}^{(m)} & \ldots & x_{n}^{(m)} \end{array}\right],\quad w=\left[\begin{array}{c}w_0\\w_2\\\vdots \\w_n \end{array}\right],\quad Y=\left[\begin{array}{c}y^{(1)} \\y^{(2)} \\\vdots \\y^{(m)}\end{array}\right]

对式 (1.5) 求导,可得:

𝐽(𝑤)𝑤=12𝑤(wTXTXwYTXwwTXTY+YTY)(1.6)\frac{\partial {𝐽(𝑤)}}{\partial 𝑤}=\frac{1}{2}\frac{\partial}{\partial 𝑤}(w^TX^TXw-Y^TXw-w^TX^TY+Y^TY)\tag{1.6}

YTXw=(wTXTY)TY^TXw=(w^TX^TY)^T𝑑𝑋𝑇𝑋𝑑𝑋=2𝑋\frac{𝑑𝑋^𝑇𝑋}{𝑑𝑋}=2𝑋,所以

𝐽(𝑤)𝑤=12𝑤(wTXTXw2wTXTY+YTY)=12(2XTXw2XTY+0)=XTXwXTY(1.7)\begin{aligned} \frac{\partial {𝐽(𝑤)}}{\partial 𝑤} &=\frac{1}{2}\frac{\partial}{\partial 𝑤}(w^TX^TXw-2w^TX^TY+Y^TY)\\ &=\frac{1}{2}(2X^TXw-2X^TY+0)\\ &=X^TXw-X^TY \end{aligned}\tag{1.7}

𝐽(𝑤)𝑤=0\frac{\partial {𝐽(𝑤)}}{\partial 𝑤}=0,则:

𝑤=(𝑋𝑇𝑋)1𝑋𝑇𝑌(1.8)𝑤=(𝑋^𝑇𝑋)^{−1}𝑋^𝑇𝑌\tag{1.8}

由上式可知,最小二乘法需要计算(𝑋𝑇𝑋)1(𝑋^𝑇𝑋)^{−1},但是矩阵求逆的时间复杂度为 𝑂(𝑛3)𝑂(𝑛3),因此当特征数量 𝑛𝑛 较大时,其运算代价非常大。所以这种方法只适用于特征数量较少的线性模型,不适用于其他模型。

现代机器学习中常用的参数更新方法是梯度下降法。

1.3 梯度下降法

根据梯度下降的每一步中用到的样本数量,可以将梯度下降法分为以下三类:

类别样本数量公式
批量梯度下降(Batch Gradient Descent,BGD)所有wj:=wjα1mi=1m((h(x(i))y(i))xj(i))w_{j}:=w_{j}-\alpha \frac{1}{m} \sum_{i=1}^{m}\left(\left(h\left(x^{(i)}\right)-y^{(i)}\right) \cdot x_{j}^{(i)}\right)
随机梯度下降(Stochastic Gradient Descent,SGD)一个wj:=wjα(h(x(i))y(i))xj(i)w_{j}:=w_{j}-\alpha\left(h\left(x^{(i)}\right)-y^{(i)}\right) \cdot x_{j}^{(i)}
小批量梯度下降(Mini-Batch Gradient Descent,MBGD)部分wj:=wjα1bk=ii+b1(h(x(k))y(k))xj(k)w_{j}:=w_{j}-\alpha \frac{1}{b} \sum_{k=i}^{i+b-1}\left(h\left(x^{(k)}\right)-y^{(k)}\right) x_{j}^{(k)}

其中,α\alpha 称为学习率。根据公式,不难看出,BGD 和 SGD 其实是 MBGD 的 bb 取值为 mm11 时的特殊情况。我们以 SGD 为例,推导一遍参数的更新过程。

wjJ(w)=wj12(h(x(i))y(i))2=212(h(x(i))y(i))wj(h(x(i))y(i))=(h(x(i))y(i))wj(i=0n(wixi(i)y(i)))=(h(x(i))y(i))xj(i)(1.9)\begin{aligned} \frac{\partial}{\partial w_{j}} J(w) &=\frac{\partial}{\partial w_{j}} \frac{1}{2}\left(h\left(x^{(i)}\right)-y^{(i)}\right)^{2}\\ &=2 \cdot \frac{1}{2}\left(h\left(x^{(i)}\right)-y^{(i)}\right) \cdot \frac{\partial}{\partial w_{j}}\left(h\left(x^{(i)}\right)-y^{(i)}\right) \\ &=\left(h\left(x^{(i)}\right)-y^{(i)}\right) \cdot \frac{\partial}{\partial w_{j}}\left(\sum_{i=0}^{n}\left(w_{i} x_{i}^{(i)}-y^{(i)}\right)\right) \\ &=\left(h\left(x^{(i)}\right)-y^{(i)}\right) x_{j}^{(i)} \end{aligned}\tag{1.9}

wj=wjαwjJ(w)w_j=w_j-\alpha\frac{\partial}{\partial w_{j}} J(w),所以

wj:=wjα(h(x(i))y(i))xj(i)(1.10)w_{j}:=w_{j}-\alpha\left(h\left(x^{(i)}\right)-y^{(i)}\right) \cdot x_{j}^{(i)}\tag{1.10}

1.4 回归的评价指标

监督学习主要分回归和分类两种,二者的评价指标是截然不同的。我们先在此介绍回归的评价指标。

首先是误差要越小越好,主要是下面几种:

指标公式
均方误差(Mean Square Error,MSE)MSE=1mi=1m(y(i)y^(i))2MSE=\frac{1}{m} \sum_{i=1}^{m}\left(y^{(i)}-\widehat{y}^{(i)}\right)^{2}
均方根误差(Root Mean Square Error,RMSE)RMSE(y,y^)=1mi=1m(y(i)y^(i))2RMSE(y, \widehat{y})=\sqrt{\frac{1}{m} \sum_{i=1}^{m}\left(y^{(i)}-\widehat{y}^{(i)}\right)^{2}}
平均绝对误差(Mean Absolute Error,MAE)MAE(y,y^)=1mi=1ny(i)y^(i) MAE(y, \widehat{y})=\frac{1}{m} \sum_{i=1}^{n}\|y^{(i)}-\widehat{y}^{(i)}\|

其次是 R2R^2 要越大越好,它越接近于 1,说明模型拟合得越好。计算公式如下:

R2(y,y^)=1i=0m(y(i)y^(i))2i=0m(y(i)yˉ)2=1SSESST=SSRSST(1.11)R^{2}(y, \widehat{y})=1-\frac{\sum_{i=0}^{m}\left(y^{(i)}-\widehat{y}^{(i)}\right)^{2}}{\sum_{i=0}^{m}\left(y^{(i)}-\bar{y}\right)^{2}}=1-\frac{\mathrm{SSE}}{\mathrm{SST}}=\frac{\mathrm{SSR}}{\mathrm{SST}}\tag{1.11}

其中,SSR=i=0m(y^(i)yˉ)2SSR=\sum_{i=0}^{m}\left(\widehat{y}^{(i)}-\bar{y}\right)^{2},SSE=i=0m(y(i)y^(i))2SSE=\sum_{i=0}^{m}\left(y^{(i)}-\widehat{y}^{(i)}\right)^{2},SST=i=0m(y(i)yˉ)2SST=\sum_{i=0}^{m}\left(y^{(i)}-\bar{y}\right)^{2}

SST 是 Sum of Squares Total 的缩写,含义是总平方和。它是变量 yy 与其平均值 yˉ\bar{y} 之差的平方和。我们可以将其视为对于变量 yy 在其平均值 yˉ\bar{y} 周围的分散程度的一种度量,很像描述性统计中的方差。

SSR 是 Sum of Squares Regression 的缩写,含义是回归的平方和。它是预测值 y^\widehat{y} 与平均值 yˉ\bar{y} 之差的平方和。我们可以将其视为描述回归线与数据的拟合程度的一种度量。

SSE 是 Sum of Squares Error 的缩写,含义是误差的平方和。它是预测值 y^\widehat{y} 与观测值 yy 之差的平方和。

image.png

从图中不难看出,三者的关系是:SST = SSR + SSE。如果 SSR 的值等于 SST,这意味着我们的回归模型是完美的。

2.逻辑回归

我们前面说过,逻辑回归和线性回归不同的地方在于:线性回归适用于解决回归问题,而逻辑回归适用于解决分类问题。本节我们就讲讲造成这种差异的原因。

2.1 Sigmoid函数

我们知道逻辑回归的目标是训练一个分类器,该分类器可以对输入数据的类别做出决策。以二元逻辑回归为例,对于一个输入 x(x1,x2,...,xn)x(x_1,x_2,...,x_n),我们希望分类器能够输出 yy 是 1(是某类的成员)或 0(不是某类的成员)。也就是说,我们想知道这个输入 xx 是该类成员的概率 P(y=1x)P(y = 1|x)

我们知道逻辑回归的底层就是线性回归,但是线性回归解决的是回归问题,那我们如何将一个回归问题转化为一个分类问题呢?为了方便讨论,我们再对式 (1.2) 做一点形式变换,即令

z=h(𝑥)=𝑤0+𝑤1𝑥1+𝑤2𝑥2+...+𝑤𝑛𝑥𝑛=wTX(2.1)z=ℎ(𝑥)=𝑤_0 +𝑤_1𝑥_1 +𝑤_2𝑥_2 +...+𝑤_𝑛𝑥_𝑛=w^TX\tag{2.1}

有些地方写的是 z=𝑤𝑇𝑥+𝑏z=𝑤^𝑇𝑥+𝑏,其实是一样的,因为 bb 可以融入到 w0w_0 中。

我们观察式 (2.1),不难发现,zz 的输出范围没有任何限制,即 (,+)(-∞, +∞)。而作为一个分类器,我们需要输出的是位于 0 和 1 之间的合法概率值。而为了完成这一步转变,我们就需要 Sigmoid函数 σ(z)σ(z)。Sigmoid 函数(命名是因为它的图像呈𝑆形)也称为逻辑函数,逻辑回归的名称也由此而来。

image.png

Sigmoid 有许多优点,它能将任意实数映射到 [0,1][0,1] 范围内,而且任意阶可导。它在 0 附近几乎是线性的,但在两端趋于平缓,它能将异常值压向 0 或 1。

σ(z)=11+ez=11+exp(z)(2.2)\sigma(z)=\frac{1}{1+e^{-z}}=\frac{1}{1+\exp (-z)}\tag{2.2}

我们将式 (2.1) 的值代入上式,就能得到一个介于 0 和 1 之间的概率。对于一个二元逻辑回归,我们只需要确保 P(y=1)+P(y=0)=1P(y = 1)+P(y = 0)=1 即可。我们可以这样做:

P(y=1)=σ(z)=11+exp(z)P(y=0)=1σ(z)=111+exp(z)=exp(z)1+exp(z)(2.3)\begin{aligned} P(y=1) &=\sigma(z)=\frac{1}{1+\exp (-z)} \\ P(y=0) &=1-\sigma(z)=1-\frac{1}{1+\exp (-z)}=\frac{\exp (-z)}{1+\exp (-z)} \end{aligned}\tag{2.3}

根据 Sigmoid 函数的性质:

1σ(x)=σ(x)(2.4)1−σ(x)=σ(−x)\tag{2.4}

我们也可以将 P(y=0)P(y = 0) 表示为 σ(z)σ(−z)

好了,现在我们已经有了概率值,那么概率值大于多少我们认为它是属于 1 类呢?通常,如果概率 P(y=1x)P(y = 1|x) 大于 0.50.5,我们认为是,否则就不是。这个 0.5 被称为决策边界

decision(x)={1 if P(y=1x)>0.50 otherwise (2.5)\operatorname{decision}(x)= \begin{cases}1 & \text { if } P(y=1 \mid x)>0.5 \\ 0 & \text { otherwise }\end{cases}\tag{2.5}

总结:逻辑回归的总体思路就是,先用逻辑函数把线性回归的结果 (-∞,∞) 映射到 (0,1),再通过决策边界建立与分类的概率联系。

逻辑函数还有一个很好的特性就是,其导函数可以转化成本身的一个表达式,推导过程如下:

σ(z)=(11+ez)=ez(1+ez)2=1+ez1(1+ez)2=1(1+ez)(11(1+ez))=σ(z)(1σ(z))(2.6)\begin{aligned} \sigma^{\prime}(z)&=\left(\frac{1}{1+e^{-z}}\right)^{\prime} \\ &=\frac{e^{-z}}{\left(1+e^{-z}\right)^{2}} \\ &=\frac{1+e^{-z}-1}{\left(1+e^{-z}\right)^{2}} \\ &=\frac{1}{\left(1+e^{-z}\right)}\left(1-\frac{1}{\left(1+e^{-z}\right)}\right) \\ &=\sigma(z)(1-\sigma(z)) \end{aligned}\tag{2.6}

在二分类模型中,事件发生与不发生的概率之比 𝑝1𝑝\frac{𝑝}{1−𝑝} 称为事件的几率(odds)。令σ(z)=p\sigma(z)=p,解得:

z=log(p1p)(2.7)z=log(\frac{p}{1-p})\tag{2.7}

也就是说,线性回归的结果(即 zz )等于对数几率

2.2 代价函数

下面我们讲讲模型的参数如何更新。我们使用极大似然估计法来求解。极大似然估计法利用已知样本结果,反推最有可能导致这样结果的原因。我们的目标就是找到一组参数,使得在这组参数下,我们的样本的似然度(概率)最大。

对于一个二分类模型,已知 P(y=1)=σ(z),P(y=0)=1σ(z)P(y=1)=\sigma(z),P(y=0)=1-\sigma(z), 则似然函数为:

L(w)=i=1mP(y(i)x(i);w)=i=1m(σ(x(i)))y(i)(1σ(x(i)))1y(i)(2.8)L(w)=\prod_{i=1}^{m} P\left(y^{(i)} \mid x^{(i)} ; w\right)=\prod_{i=1}^{m}\left(\sigma\left(x^{(i)}\right)\right)^{y^{(i)}}\left(1-\sigma\left(x^{(i)}\right)\right)^{1-y^{(i)}}\tag{2.8}

等式两边同时取对数:

l(w)=logL(w)=i=1m(y(i)log(σ(x(i)))+(1y(i))log(1σ(x(i))))(2.9)l(w)=\log L(w)=\sum_{i=1}^{m}\left(y^{(i)} \log \left(\sigma\left(x^{(i)}\right)\right)+\left(1-y^{(i)}\right) \log \left(1-\sigma\left(x^{(i)}\right)\right)\right)\tag{2.9}

则代价函数为:

J(w)=1ml(w)(2.10)J(w)=-\frac{1}{m} l(w)\tag{2.10}

其实,式(2.9)前面再加一个负号,就是我们常用的损失函数——交叉熵损失(cross-entropy loss)

代价函数之所以要加负号,是因为机器学习的目标是最小化损失函数,而极大似然估计法的目标是最大化似然函数。那么加个负号,正好使二者等价。

2.3 梯度下降法求解

求解逻辑回归的方法有非常多,我们这里主要了解一下梯度下降法。

将式 (2.2) 代入式 (2.10),得:

y(i)log(σ(x(i)))+(1y(i))log(1σ(x(i)))=y(i)log(11+ewTx(i))+(1y(i))log(111+ewTx(i))=y(i)log(1+ewTx(i))(1y(i))log(1+ewTx(i))J(w)=1mi=1m(y(i)log(1+ewTx(i))(1y(i))log(1+ewTx(i)))(2.11)\begin{aligned} \because & y^{(i)} \log \left(\sigma\left(x^{(i)}\right)\right)+\left(1-y^{(i)}\right) \log \left(1-\sigma\left(x^{(i)}\right)\right)\\ &=y^{(i)} \log \left(\frac{1}{1+e^{-w^{T} x^{(i)}}}\right)+\left(1-y^{(i)}\right) \log \left(1-\frac{1}{1+e^{-w^{T} x^{(i)}}}\right)\\ &=-y^{(i)} \log \left(1+e^{-w^{T} x^{(i)}}\right)-\left(1-y^{(i)}\right) \log \left(1+e^{w^{T} x^{(i)}}\right)\\\therefore J(w)&=-\frac{1}{m}\sum_{i=1}^{m}\left(-y^{(i)} \log \left(1+e^{-w^{T} x^{(i)}}\right)-\left(1-y^{(i)}\right) \log \left(1+e^{w^{T} x^{(i)}}\right)\right) \end{aligned}\tag{2.11}

J(w)J(w) 求偏导,得:

wjJ(w)=wj(1mi=1m(y(i)log(1+ewTx(i))(1y(i))log(1+ewTx(i))))=1mi=1m(y(i)xj(i)ewTx(i)1+ewTx(i)(1y(i))xj(i)ewTx(i)1+ewTx(i))=1mi=1m(y(i)σ(x(i)))xj(i)=1mi=1m(σ(x(i))y(i))xj(i)(2.12)\begin{aligned} \frac{\partial}{\partial w_{j}} J(w)&=\frac{\partial}{\partial w_{j}}\left(-\frac{1}{m} \sum_{i=1}^{m}\left(-y^{(i)} \log \left(1+e^{-w^{T} x^{(i)}}\right)-\left(1-y^{(i)}\right) \log \left(1+e^{w^{T} x^{(i)}}\right)\right)\right) \\ &=-\frac{1}{m} \sum_{i=1}^{m}\left(-y^{(i)} \frac{-x_{j}^{(i)} e^{-w^{T} x^{(i)}}}{1+e^{-w^{T} x^{(i)}}}-\left(1-y^{(i)}\right) \frac{x_{j}^{(i)} e^{w^{T} x^{(i)}}}{1+e^{w^{T} x^{(i)}}}\right) \\ &=-\frac{1}{m} \sum_{i=1}^{m}\left(y^{(i)}-\sigma\left(x^{(i)}\right)\right) x_{j}^{(i)} \\ &=\frac{1}{m} \sum_{i=1}^{m}\left(\sigma\left(x^{(i)}\right)-y^{(i)}\right) x_{j}^{(i)} \end{aligned}\tag{2.12}

所以:

wj:=wjwjJ(w)=wjα1mi=1m(σ(x(i))y(i))xj(i)(2.13)w_{j}:=w_{j}-\frac{\partial}{\partial w_{j}} J(w)=w_{j}-\alpha \frac{1}{m} \sum_{i=1}^{m}\left(\sigma\left(x^{(i)}\right)-y^{(i)}\right) x_{j}^{(i)}\tag{2.13}

2.4 逻辑回归的分类

逻辑回归对特征变量(x)和分类响应变量(y)之间的关系进行建模,在给定一组预测变量的情况下,它能给出落入特定类别响应水平的概率。也就是说,你给它一组数据(特征),它告诉你这组数据属于某一类别的概率。根据分类响应变量(y)的性质,我们可以将逻辑回归分为三类:

  • 二元逻辑回归(Binary Logistic Regression) 当分类结果只有两种可能的时候,我们就称为二元逻辑回归。例如,考试通过或未通过,回答是或否,血压高或低。
  • 名义逻辑回归(Nominal Logistic Regression) 当存在三个或更多类别且类别之间没有自然排序时,我们就称为名义逻辑回归。例如,企业的部门有策划、销售、人力资源等,颜色有黑色、红色、蓝色、橙色等。
  • 序数逻辑回归(Ordinal Logistic Regression) 当存在三个或更多类别且类别之间有自然排序时,我们就称为序数逻辑回归。例如,评价有好、中、差,身材有偏胖、中等、偏瘦。注意,类别的排名不一定意味着它们之间的间隔相等。

2.5 Softmax Regression

我们以上讨论都是以二元分类为前提展开的,但有时可能有两个以上的类,例如情绪分类有正面、负面或中性三种,手写数字识别有 0-9 总共十种分类。

在这种情况下,我们需要多元逻辑回归,也称为 Softmax Regression。在多元逻辑回归中,我们假设总共有 KK 个类,每个输入 xx 只能属于一个类。也就是说,每个输入 xx 的输出 yy 将是一个长度为 KK 的向量。如果类 cc 是正确的类,则设置yc=1\mathbf{y}_c=1,其他元素设置为 0,即yc=1yj=0jc\mathbf{y}_{c}=1,\mathbf{y}_{j}=0\quad\forall j \neq c。这种编码方式称为 one-hot 编码。此时,分类器的工作是产生一个估计向量。对于每个类 kk,产生一个概率估计 P(yk=1x)P(y_k = 1|x)

多元逻辑分类器采用 sigmoid 的泛化函数 softmax 来计算 P(yk=1x)P(y_k = 1|x)。它能将 KK 个任意值的向量 z=[z1,z2,,zK]z = [z_1,z_2,\dots,z_K] 映射到 (0,1) 范围内的概率分布上,并且所有值的总和为 1 。其定义如下:

softmax(zi)=exp(zi)j=1Kexp(zj)1iK(2.14)\operatorname{softmax}\left(\mathbf{z}_{i}\right)=\frac{\exp \left(\mathbf{z}_{i}\right)}{\sum_{j=1}^{K} \exp \left(\mathbf{z}_{j}\right)} 1 \leq i \leq K\tag{2.14}

因此,输入向量 z=[z1,z2,,zK]z = [z_1,z_2,\dots,z_K] 的 softmax输出如下:

softmax(z)=[exp(z1)i=1Kexp(zi),exp(z2)i=1Kexp(zi),,exp(zK)i=1Kexp(zi)](2.15)\operatorname{softmax}(\mathbf{z})=\left[\frac{\exp \left(\mathbf{z}_{1}\right)}{\sum_{i=1}^{K} \exp \left(\mathbf{z}_{i}\right)}, \frac{\exp \left(\mathbf{z}_{2}\right)}{\sum_{i=1}^{K} \exp \left(\mathbf{z}_{i}\right)}, \ldots, \frac{\exp \left(\mathbf{z}_{K}\right)}{\sum_{i=1}^{K} \exp \left(\mathbf{z}_{i}\right)}\right]\tag{2.15}

举个例子,若 z=[0.6,1.1,1.5,1.2,3.2,1.1]\mathbf{z}=[0.6,1.1,-1.5,1.2,3.2,-1.1],则

softmax(z)=[0.055,0.090,0.006,0.099,0.74,0.010]\operatorname{softmax}(\mathbf{z})=[0.055,0.090,0.006,0.099,0.74,0.010]

与 sigmoid 一样,softmax 也具有将异常值压向 0 或 1 的特性。因此,如果其中一个输入远大于其他输入,它将倾向于将其概率推向 1,并抑制较小输入的概率。

由于现在有 KK 个类别,而且每一个类别都有一个单独的权重向量 wkw_k,则每个输出类别的概率估计 yk^\widehat{y_k} 应该这样计算:

P(yk=1x)=exp(wkTx)j=1Kexp(wjTx)(2.16)P\left(\mathbf{y}_{k}=1 \mid \mathbf{x}\right)=\frac{\exp \left(\mathbf{w}_{k}^T \cdot \mathbf{x}\right)}{\sum_{j=1}^{K} \exp \left(\mathbf{w}_{j}^T \cdot \mathbf{x}\right)}\tag{2.16}

上式看起来我们将分别每个类别计算输出。但是,我们可以将 KK 个权重向量表示为权重矩阵 WWWW 的第 kk 列对应于权重向量 wkw_k。因此,WW 的形状为 [(n+1)×K][(n+1) \times K],其中,KK 是输出类的数量,nn 是输入特征的数量,加 11 对应的就是偏置。这样一来,我们还是可以通过一个优雅的公式计算 y^\widehat{y}

y^=softmax(WTX)(2.17)\widehat{y}=\operatorname{softmax}{(W^TX)}\tag{2.17}

3. 代码实现(Pytorch)

class LogisticRegression(torch.nn.Module):

    def __init__(self, num_features):
        super(LogisticRegression, self).__init__() #调用父类的构造函数
        self.linear = torch.nn.Linear(num_features, 1)
        
        self.linear.weight.detach().zero_() #权值初始化为0
        self.linear.bias.detach().zero_() #偏置初始化为0
        
    def forward(self, x):
        logits = self.linear(x)
        probas = torch.sigmoid(logits)
        return probas

在 Pytorch 中可以通过继承torch.nn.Module类来实现自定义模型,在__init__中定义每一层的构造,在forward中定义每一层的连接关系,是实现模型功能的核心,且必须重写,否则模型将由于无法找到各层的连接关系而无法执行。

torch.nn.Linear对传入的数据做线性变换,weightbias 是它的两个变量,分别代表学习的权值和偏置。

4. 鸢尾花分类

我们以鸢尾花数据的分类为例,做一个简单地用 Logistic 回归进行分类的任务。

  • 数据集获取与划分
import matplotlib.pyplot as plt
import numpy as np
from io import BytesIO

import torch
import torch.nn.functional as F

ds = np.lib.DataSource()
fp = ds.open('http://archive.ics.uci.edu/ml/machine-learning-databases/iris/iris.data')

x = np.genfromtxt(BytesIO(fp.read().encode()), delimiter=',', usecols=range(2), max_rows=100)
y = np.zeros(100)
y[50:] = 1

np.random.seed(1)
idx = np.arange(y.shape[0])
np.random.shuffle(idx)
X_test, y_test = x[idx[:25]], y[idx[:25]]
X_train, y_train = x[idx[25:]], y[idx[25:]]
mu, std = np.mean(X_train, axis=0), np.std(X_train, axis=0)
X_train, X_test = (X_train - mu) / std, (X_test - mu) / std

fig, ax = plt.subplots(1, 2, figsize=(7, 2.5))
ax[0].scatter(X_train[y_train == 1, 0], X_train[y_train == 1, 1])
ax[0].scatter(X_train[y_train == 0, 0], X_train[y_train == 0, 1])
ax[1].scatter(X_test[y_test == 1, 0], X_test[y_test == 1, 1])
ax[1].scatter(X_test[y_test == 0, 0], X_test[y_test == 0, 1])
plt.show()

实际的数据共有152行,我们只取前100行。按照下标随机打乱之后来划分训练集和测试集。其中,训练集有75行,测试集有25行。

image.png

  • 模型训练
model = LogisticRegression(num_features=2).to(device)
cost_fn = torch.nn.BCELoss(reduction='sum')
optimizer = torch.optim.SGD(model.parameters(), lr=0.1)

def custom_where(cond, x_1, x_2):
    return (cond * x_1) + ((1-cond) * x_2)
    
def comp_accuracy(label_var, pred_probas):
    pred_labels = custom_where((pred_probas > 0.5).float(), 1, 0).view(-1)
    acc = torch.sum(pred_labels == label_var.view(-1)).float() / label_var.size(0)
    return acc


num_epochs = 10

X_train_tensor = torch.tensor(X_train, dtype=torch.float32, device=device)
y_train_tensor = torch.tensor(y_train, dtype=torch.float32, device=device).view(-1, 1)


for epoch in range(num_epochs):
    
    #### Compute outputs ####
    out = model(X_train_tensor)
    
    #### Compute gradients ####
    cost = cost_fn(out, y_train_tensor)
    optimizer.zero_grad()
    cost.backward()
    
    #### Update weights ####  
    optimizer.step()
    
    #### Logging ####      
    pred_probas = model(X_train_tensor)
    acc = comp_accuracy(y_train_tensor, pred_probas)
    print('Epoch: %03d' % (epoch + 1), end="")
    print(' | Train ACC: %.3f' % acc, end="")
    print(' | Cost: %.3f' % cost_fn(pred_probas, y_train_tensor))


    
print('\nModel parameters:')
print('  Weights: %s' % model.linear.weight)
print('  Bias: %s' % model.linear.bias)

BCELoss计算目标值和预测值之间的二进制交叉熵损失函数,数学公式如下:

Loss=w[ylog(y)+(1y)log(1y)]Loss = -w[ylog(y')+(1-y)log(1-y')]

其中,wwyyyy' 分别表示权值、标签(target)、预测概率值(input probabilities)。reduction取值为sum表明对样本的损失值进行求和。

输出如下:

Epoch: 001 | Train ACC: 0.987 | Cost: 5.581
Epoch: 002 | Train ACC: 0.987 | Cost: 4.882
Epoch: 003 | Train ACC: 1.000 | Cost: 4.381
Epoch: 004 | Train ACC: 1.000 | Cost: 3.998
Epoch: 005 | Train ACC: 1.000 | Cost: 3.693
Epoch: 006 | Train ACC: 1.000 | Cost: 3.443
Epoch: 007 | Train ACC: 1.000 | Cost: 3.232
Epoch: 008 | Train ACC: 1.000 | Cost: 3.052
Epoch: 009 | Train ACC: 1.000 | Cost: 2.896
Epoch: 010 | Train ACC: 1.000 | Cost: 2.758

Model parameters:
  Weights: Parameter containing:
tensor([[ 4.2267, -2.9613]], requires_grad=True)
  Bias: Parameter containing:
tensor([0.0994], requires_grad=True)
  • 模型评估
X_test_tensor = torch.tensor(X_test, dtype=torch.float32, device=device)
y_test_tensor = torch.tensor(y_test, dtype=torch.float32, device=device)

pred_probas = model(X_test_tensor)
test_acc = comp_accuracy(y_test_tensor, pred_probas)

print('Test set accuracy: %.2f%%' % (test_acc*100))

输出如下:

Test set accuracy: 100.00%

image.png

结语

一个非常基础的入门级的机器学习算法,但是要完全讲透还是非常困难,尤其是公式的推导部分,光用文字描述很难讲清楚,但是总体上公式的推导都没问题,只是说第一眼看不是那么容易懂,多看几遍或者自己在纸上写一写还是不难理解的。另外,也请读者发现错误后能在评论区指出,我看到后会及时更正。

全文写下来洋洋洒洒六千字,也花费了好几天的时间,从自己理解到写出来让别人看懂,这之间差别真的很大,可能也是我不善于表达,词不达意,也请大家见谅。

篇幅这么长,公式这么多,看完本文你可能什么也记不住,但是这两点你一定要记住:

1. Logistic回归不适用于解决回归问题,它适用于解决分类问题。千万不要被它的名称迷惑!

2. Logistic回归 = 线性回归 + Sigmoid 函数。当然了,如果是多元回归的话就是 softmax 函数。

例行公事😂:如果你觉得本文确实对你有帮助,请点个赞支持我一下吧 🌹