前言：
本文主要讲解逻辑回归的二分类和多分类问题。在此，笔者将不做过多引入的讲解，而是直接从理论出发。因为在我看来，网上对这方面的讲解已经很多，写得也非常的好。笔者主要是讲解公式转化为代码这一部分。
关键字：逻辑回归、sigmoid、softmax、Python
sigmoid:
p = 1 1 + e − w T x p=\frac{1}{1+e^{-w^Tx}} p=1+e−wTx1
对于二分类问题，我们使用的是sigmoid函数，采用极大似然估计的方法进行估计。首先二分类的似然函数学过概率论的人都知道，他是这样的
f
(
x
)
=
p
y
i
(
1
−
p
)
(
1
−
y
i
)
f(x)=p^{y_{i}}(1 - p)^{(1 - y_{i})}
f(x)=pyi(1−p)(1−yi)
为什么是极大似然估计？因为极大似然能够估计出样本的分布空间。如果我们的样本具有代表性，那么理论上此时的样本的模型参数也是接近总体的样本的模型参数。
损失函数及其梯度计算：
显而易见，我们可以直接将极大似然估计作为损失函数，当燃了，损失函数一般我们是越小越好。而极大似然估计是越大越好，so？我们要对极大似然取反，最终损失函数应当如下：
L
(
w
)
=
min
w
−
1
m
∏
i
=
1
n
p
y
i
(
1
−
p
)
1
−
y
i
L(w)=\min\limits_{w}-\frac{1}{m}{\prod\limits_{i=1}^n{p^{y_i}(1-p)^{1-y_i}}}
L(w)=wmin−m1i=1∏npyi(1−p)1−yi
为了便于运算，一般情况下我们会对其进行取log对数，然后对其化简最终得到（此处笔者不作化简过程，因为LaTeX公式太难写，笔者都是从word文档写好再转化过来的，读者可参照其他博主的相关文章）：
l
o
g
L
(
w
)
=
min
w
−
1
n
∑
i
=
1
n
[
y
i
log
(
p
(
x
i
)
)
+
(
1
−
y
i
)
l
o
g
(
1
−
p
(
x
i
)
)
logL(w) = {\min\limits_{w}{- \frac{1}{n}{\sum\limits_{i = 1}^{n}\left\lbrack y_{i}{\log(p({x_i}))} + (1 - {y_i})log(1 -p({x_i})) \right.}}}
logL(w)=wmin−n1i=1∑n[yilog(p(xi))+(1−yi)log(1−p(xi))
然后有了损失函数，我们就要利用梯度下降算法对其求解，因此，要对其求导求梯度。我们直接对w求导(求导过程请参照其他文章)
∂
l
o
g
L
(
w
)
∂
w
=
∂
∂
w
(
−
1
n
∑
i
=
1
n
[
y
i
log
(
1
1
+
e
−
w
T
x
i
)
+
(
1
−
y
i
)
log
(
1
−
1
1
+
e
−
w
T
x
i
)
]
)
=
1
n
∑
i
=
1
n
x
i
(
1
1
+
e
−
w
T
x
i
−
y
i
)
\begin{equation} \begin{aligned} \frac{\partial logL(w)}{\partial w} & =\frac{\partial}{\partial w}({{{- \frac{1}{n}{\sum\limits_{i = 1}^{n}\left\lbrack {y_{i}{\log\left( \frac{1}{1 + e^{- w^Tx_i}} \right)} + \left( {1 - y_i} \right){\log\left( {1 - \frac{1}{1+e^{-w^Tx_i}}} \right)}} \right]}}}})\\ &=\frac{1}{n}\sum\limits_{i=1}^nx_i(\frac{1}{1+e^{-w^T{x_i}}}-y_i) \end{aligned} \end{equation}
∂w∂logL(w)=∂w∂(−n1i=1∑n[yilog(1+e−wTxi1)+(1−yi)log(1−1+e−wTxi1)])=n1i=1∑nxi(1+e−wTxi1−yi)
因为在程序实现当中，为了程序的简介和提高运算的速度，我们往往是采用矩阵的形式进行预算。因此，要将上面的公式转化为矩阵相乘。当燃了，眼尖一眼就可以看出来。在转化之前，我们先介绍一下各个变量的矩阵表达
x
=
(
x
1
T
⋮
x
n
T
)
=
(
x
11
⋯
x
1
p
⋮
⋱
⋮
x
n
1
⋯
x
n
p
)
n
∗
p
;
y
=
(
y
1
⋮
y
n
)
n
∗
1
;
p
=
(
p
1
⋮
p
n
)
n
∗
1
x = \begin{pmatrix} {x_{1}}^T\\ \vdots \\ {x_{n}}^T\\ \end{pmatrix} = \begin{pmatrix} x_{11} & \cdots & x_{1p} \\ \vdots & \ddots & \vdots \\ x_{n1} & \cdots & x_{np} \\ \end{pmatrix}_{n*p};y = \begin{pmatrix} y_{1} \\ \vdots \\ y_{n} \\ \end{pmatrix}_{n*1};p = \begin{pmatrix} p_{1} \\ \vdots \\ p_{n} \\ \end{pmatrix}_{n*1}
x=
x1T⋮xnT
=
x11⋮xn1⋯⋱⋯x1p⋮xnp
n∗p;y=
y1⋮yn
n∗1;p=
p1⋮pn
n∗1
n*p代表有n个数据，p代表每个数据有p个维度。里面的单个分量(如
x
1
,
w
1
x_1,w_1
x1,w1)都是以列向量的形式存在.
下面我们对梯度简单转化一下（暂时忽略1/n这个标量）
原式
=
(
x
1
x
2
⋯
x
n
)
(
1
1
+
e
−
w
T
x
1
−
y
1
1
1
+
e
−
w
T
x
2
−
y
2
⋮
1
1
+
e
−
w
T
x
n
−
y
n
)
=
x
T
(
p
−
y
)
\begin{equation} \begin{aligned} 原式 &= \left( {\begin{matrix} x_{1} & x_{2} & \cdots \\ \end{matrix}~~~~~~x_{n}~} \right)\begin{pmatrix} \begin{matrix} {\frac{1}{1 + e^{- w^{T}x_{1}}} - y_{1}} \\ {\frac{1}{1 + e^{- w^{T}x_{2}}} - y_{2}} \\ \vdots \\ \end{matrix} \\ {\frac{1}{1 + e^{- w^{T}x_{n}}} - y_{n}} \\ \end{pmatrix}\\ &=x^T(p-y) \end{aligned} \end{equation}
原式=(x1x2⋯ xn )
1+e−wTx11−y11+e−wTx21−y2⋮1+e−wTxn1−yn
=xT(p−y)
因此
∂
l
o
g
L
(
w
)
∂
w
=
1
n
x
T
(
p
−
y
)
\frac{\partial logL(w)}{\partial w}= \frac{1}{n}x^T(p-y)
∂w∂logL(w)=n1xT(p−y)
代码实现：
注意了，这条直线实际上并不是一条，而是两条，一条是实际上区分数据的，另一条是训练出来的，可以看到基本重合，因此，可以说明训练效果挺好的。
import numpy as np
from scipy import stats
import matplotlib.pyplot as plt
class Logistic_Regression():
def __init__(self,x,y):
'''
:param x: 样本值
:param y: 标签值
'''
#在x的中增加一个全为1的维度。因为后续要将偏置项b写入w之中
self.x=np.insert(x,x.shape[1],1,1)
self.y=y
#m，n为x的维度
self.m,self.n=self.x.shape
#生成一个与（n,1）维的全为1的数组
self.w=np.zeros_like(x,shape=(self.n,1))
def sigmoid(self,x):
'''
概率生成函数
'''
return 1/(np.exp(-x)+1)
def train(self,epoch):
'''
:param epoch: 训练轮次
:return:
'''
for i in range(epoch):
#将x*w后的值代入sigmoid函数生成概率
p=self.sigmoid(self.x@self.w)
#用交叉熵损失计算损失
#1e-5为了防止下溢出
loss=(-1/self.m*(self.y*np.log(p+1e-5)+(1-self.y)*np.log(1-p+1e-5))).sum()
#梯度下降并更新梯度
self.w-=0.0001*self.x.T@(p-self.y)
#设定停止训练的的loss最小
if loss<0.01:
break
def predict(self,x):
#给x增加一个全为1的维度
x=np.insert(x,x.shape[1],1,1)
#计算概率
p=self.sigmoid(x@self.w)
#大于0.5则认为是正类，反之则为负类
result=np.array(p>0.5,dtype=int)
return result
#绘图函数
def plot_figure(x1,x2,norm_x2,y,w):
'''
:param x1: x的第一个维度
:param x2: x的第二个维度
:param norm_x2: 加入噪声之后的x2
:param y: 标签值
:return:
'''
map_color={0:"r",1:"g"}
y=[ map_color[i] for i in y.squeeze()]
plt.scatter(x1,norm_x2,c=y)
w=w.squeeze()
x3=-(w[0]*x1.squeeze()+w[2])/w[1]
plt.plot(x1.squeeze(),x3)
plt.plot(x1,x2)
plt.show()
def main():
x1=np.arange(0,100,1).reshape(-1,1) #生成x的第一个维度
x2=np.arange(0,100,1).reshape(-1,1)#生成x的第二个维度
norm_x2=x2+stats.norm.rvs(1,20,x2.shape)#给第二个维度加上噪声
y=np.array(x2>norm_x2,dtype=int)#对大于norm的取1，反之取0
x = np.concatenate([x1, norm_x2], axis=1)#联立x1,x2 (100,2)
model=Logistic_Regression(x,y) #初始化模型并将数据传入
model.train(100000) #训练
print(model.w) #打印模型参数，最后一个维度为偏置项b
result=model.predict(x) #预测
plot_figure(x1,x2,norm_x2,result,model.w) #将预测结果画出
if __name__ == '__main__':
main()
softmax:
对于多分类问题，原始的sigmoid函数已经不再适用，因此，要引入softmax函数。在引入之前，我们先对w矩阵做一下定义
w
=
(
w
1
T
⋮
w
n
T
)
=
(
w
11
w
12
⋯
x
1
p
w
21
w
22
⋯
x
2
p
⋮
⋮
⋱
⋮
w
k
1
w
k
2
⋯
w
k
p
)
k
∗
p
w = \begin{pmatrix} {w_{1}}^T\\ \vdots \\ {w_{n}}^T \\ \end{pmatrix} = \begin{pmatrix} w_{11} & w_{12}&\cdots & x_{1p} \\ w_{21} & w_{22} & \cdots & x_{2p}\\ \vdots & \vdots&\ddots & \vdots \\ w_{k1} & w_{k2}&\cdots & w_{kp} \\ \end{pmatrix}_{k*p}
w=
w1T⋮wnT
=
w11w21⋮wk1w12w22⋮wk2⋯⋯⋱⋯x1px2p⋮wkp
k∗p
该矩阵第i列代表的是第i类的参数。k代表是一列有多少个元素，实际上对应的就是x的维度，而p则代表是分类的类别数。
而在此之前，我们也需要对标签值y的矩阵进行重定义，它将不再是一个列向量。而是一个矩阵。
y
=
(
y
1
T
y
2
T
⋮
y
n
T
)
=
(
y
11
y
12
⋯
y
1
p
y
21
y
22
⋯
y
2
p
⋮
⋮
⋱
⋮
y
n
1
y
n
2
⋯
y
n
p
)
n
∗
p
y=\begin{pmatrix} y_1^T \\ y_2^T \\ \vdots \\ y_n ^T \end{pmatrix} =\begin{pmatrix} y_{11} & y_{12} & \cdots & y_{1p}\\ y_{21} & y_{22} & \cdots & y_{2p}\\ \vdots & \vdots & \ddots & \vdots\\ y_{n1} & y_{n2} & \cdots & y_{np}\\ \end{pmatrix}_{n*p}
y=
y1Ty2T⋮ynT
=
y11y21⋮yn1y12y22⋮yn2⋯⋯⋱⋯y1py2p⋮ynp
n∗p
其中n代表样本数量，p代表类别数量。
因此，我们定义softmax函数：
g
(
x
)
=
1
∑
i
=
1
k
e
w
T
x
(
e
p
1
e
p
2
⋮
e
p
n
)
g(x)=\frac{1}{\sum\limits_{i=1}^ke^{w^Tx}} \begin{pmatrix} e_{p_1} \\ e_{p_2} \\ \vdots \\ \\ e_{p_n} \end{pmatrix}
g(x)=i=1∑kewTx1
ep1ep2⋮epn
其中，
e
p
i
e_{p_i}
epi代表第i类的概率。所以某一类的概率就可以写成如下：
p
(
y
=
y
j
∣
x
i
;
w
j
)
=
e
w
T
x
i
∑
i
=
1
k
e
w
T
x
i
p(y=y_j|x_i;w_j)=\frac{e^{w^Tx_i}}{\sum\limits_{i=1}^ke^{w^Tx_i}}
p(y=yj∣xi;wj)=i=1∑kewTxiewTxi
以上代表第i个数据，属于第j类的概率。
损失函数：
softmax的损失函数如下：
l
o
g
L
(
w
)
=
−
1
n
[
∑
i
=
1
n
∑
j
=
1
p
y
i
j
l
o
g
(
e
w
j
T
x
i
∑
l
=
1
k
w
l
T
x
i
)
]
logL(w)=-\frac{1}{n}[\sum\limits_{i=1}^n\sum\limits_{j=1}^p{y_{ij} {log(\frac{e^{w_j^Tx_i}}{\sum\limits_{l=1}^kw_l^Tx_i})}}]
logL(w)=−n1[i=1∑nj=1∑pyijlog(l=1∑kwlTxiewjTxi)]
其中
1
{
y
i
=
j
}
1\left\{y_{i}=j\right\}
1{yi=j}是示性函数，为真是取1，反之取0。
有了损失函数，我们就可以依照上面的方法对其进行求导然后利用梯度下降求解啦！
∂
l
o
g
L
(
w
)
∂
w
=
(
∂
l
o
g
L
(
w
)
∂
w
j
1
∂
l
o
g
L
(
w
)
∂
w
j
2
⋯
∂
l
o
g
L
(
w
)
∂
w
j
p
)
\begin{equation} \begin{aligned} \frac{\partial logL(w)}{\partial w}= \begin{pmatrix} \frac{\partial logL(w)}{\partial w_{j1}} & \frac{\partial logL(w)}{\partial w_{j2}} & \cdots \frac{\partial logL(w)}{\partial w_{jp}} \end{pmatrix} \end{aligned} \end{equation}
∂w∂logL(w)=(∂wj1∂logL(w)∂wj2∂logL(w)⋯∂wjp∂logL(w))
这里的
w
j
p
w_{jp}
wjp代表的是对第p个分类的w求偏导，下面我们单独提取出一个进行求导（
w
j
w_{j}
wj表示某一类）
∂
l
o
g
L
(
w
)
∂
w
j
=
∂
∂
w
j
(
−
1
n
[
∑
i
=
1
n
∑
j
=
1
p
y
i
j
l
o
g
(
e
w
j
T
x
i
∑
l
=
1
k
w
l
T
x
i
)
]
)
=
−
1
n
[
∑
i
=
1
n
x
i
(
y
i
j
−
p
(
y
i
=
j
∣
x
i
;
w
j
)
)
]
\begin{equation} \begin{aligned} \frac{\partial logL(w)}{\partial w_{j}}&= \frac{\partial }{\partial w_{j}}\left( -\frac{1}{n}[\sum\limits_{i=1}^n\sum\limits_{j=1}^p{y_{ij} {log(\frac{e^{w_j^Tx_i}}{\sum\limits_{l=1}^kw_l^Tx_i})}}]\right )\\ &=-\frac{1}{n}[\sum\limits_{i=1}^n{x_i}\left(y_{ij} -p(y_i=j|x_i;w_j)\right)] \end{aligned} \end{equation}
∂wj∂logL(w)=∂wj∂
−n1[i=1∑nj=1∑pyijlog(l=1∑kwlTxiewjTxi)]
=−n1[i=1∑nxi(yij−p(yi=j∣xi;wj))]
那么接下来就把它化为矩阵表达就可以完美收工啦(以下转化暂时忽略掉-1/n这个标量)
∂
l
o
g
L
(
w
)
∂
w
j
=
(
x
1
x
2
⋯
x
n
)
(
y
1
j
−
p
1
j
y
2
j
−
p
2
j
⋮
y
n
j
−
p
n
j
)
=
x
T
∗
(
y
j
−
p
j
)
\begin{equation} \begin{aligned} \frac{\partial logL(w)}{\partial w_{j}}&= \begin{pmatrix} x_1 & x_2 & \cdots & x_n \end{pmatrix} \begin{pmatrix} y_{1j}-p_{1j} \\ y_{2j}-p_{2j} \\ \vdots \\ y_{nj}-p_{nj} \\ \end{pmatrix} \\ &=x^T*(y_j-p_j) \end{aligned} \end{equation}
∂wj∂logL(w)=(x1x2⋯xn)
y1j−p1jy2j−p2j⋮ynj−pnj
=xT∗(yj−pj)
因此：
∂
l
o
g
L
(
w
)
∂
w
=
(
∂
l
o
g
L
(
w
)
∂
w
j
1
∂
l
o
g
L
(
w
)
∂
w
j
2
⋯
∂
l
o
g
L
(
w
)
∂
w
j
p
)
=
(
x
T
∗
(
y
j
1
−
p
j
1
)
x
T
∗
(
y
j
2
−
p
j
2
)
⋯
x
T
∗
(
y
j
p
−
p
j
p
)
)
=
x
T
(
y
j
1
−
p
j
1
y
j
2
−
p
j
2
⋯
y
j
p
−
p
j
p
)
=
−
1
n
x
T
(
y
−
p
)
\begin{equation} \begin{aligned} \frac{\partial logL(w)}{\partial w}&= \begin{pmatrix} \frac{\partial logL(w)}{\partial w_{j1}} & \frac{\partial logL(w)}{\partial w_{j2}} & \cdots \frac{\partial logL(w)}{\partial w_{jp}} \end{pmatrix}\\ &=\begin{pmatrix} x^T*(y_{j1}-p_{j1}) & x^T*(y_{j2}-p_{j2}) & \cdots & x^T*(y_{jp}-p_{jp}) & \end{pmatrix}\\ &=x^T \begin{pmatrix} y_{j1}-p_{j1} & y_{j2}-p_{j2} & \cdots & y_{jp}-p_{jp} \end{pmatrix}\\ &=-\frac{1}{n}x^T\left( y-p\right) \end{aligned} \end{equation}
∂w∂logL(w)=(∂wj1∂logL(w)∂wj2∂logL(w)⋯∂wjp∂logL(w))=(xT∗(yj1−pj1)xT∗(yj2−pj2)⋯xT∗(yjp−pjp))=xT(yj1−pj1yj2−pj2⋯yjp−pjp)=−n1xT(y−p)
此处直接将上文忽略的-1/n写回。
代码实现：
import numpy as np
from scipy import stats
import matplotlib.pyplot as plt
class Logistic_Regression():
def __init__(self,x,y):
'''
:param x: 样本值
:param y: 标签值
'''
#在x的中增加一个全为1的维度。因为后续要将偏置项b写入w之中
self.x=np.insert(x,x.shape[1],1,1)
self.y=np.array(y.squeeze(),dtype=int).reshape(-1)
#计算y的类别数
self.class_num=len(set(self.y))
#m，n为x的维度
self.m,self.n=self.x.shape
#构造one_hot 编码的标签值
self.onehot_y=np.zeros_like(y,shape=(y.shape[0],self.class_num))
self.onehot_y[np.arange(self.y.shape[0]),self.y]=1
#生成一个(维度，类别) 也就是 (n,self.class_num)
self.w=np.zeros_like(x,shape=(self.n,self.class_num))
def softmax(self,x):
'''
概率生成函数
x:(400,4)
'''
exp=np.exp(x)
ex_sum=np.sum(exp,axis=1).reshape(-1,1)
result=exp/ex_sum
return result
def train(self,epoch):
'''
:param epoch: 训练轮次
:return:
'''
for i in range(epoch):
#将x*w后的值代入sigmoid函数生成概率
p=self.softmax(self.x@self.w)#(400,4)
#用交叉熵损失计算损失
loss=-(1/self.m)*np.sum((p*self.onehot_y))
print("函数损失:",loss)
#梯度下降并更新梯度
GD=(-1/self.m)*self.x.T@(self.onehot_y-p)
self.w-=0.001*GD
def predict(self,x):
#给x增加一个全为1的维度
x=np.insert(x,x.shape[1],1,1)
#计算概率
p=self.softmax(x@self.w)
#按行求最大值并返回其对应的index
result=np.argmax(p,axis=1)
return result
#绘图函数
def plot_figure(x,y):
map_color={0:"r",1:"g",2:"b",3:"y"}
color=[map_color[i] for i in y]
plt.scatter(x[:,0],x[:,1],c=color)
plt.show()
def main():
#生成第一个团和对应的标签值
class1_x1=np.concatenate([stats.norm.rvs(2,2,(100,1)),stats.norm.rvs(1,2,(100,1))],axis=1)
y1=np.zeros_like(class1_x1,shape=(class1_x1.shape[0],1))
# 生成第二个团和对应的标签值
class1_x2=np.concatenate([stats.norm.rvs(20,2,(100,1)),stats.norm.rvs(1,2,(100,1))],axis=1)
y2 = np.ones_like(class1_x2, shape=(class1_x2.shape[0], 1))
# 生成第三个团和对应的标签值
class1_x3=np.concatenate([stats.norm.rvs(-5,2,(100,1)),stats.norm.rvs(10,2,(100,1))],axis=1)
y3 = np.ones_like(class1_x3, shape=(class1_x3.shape[0], 1))*2
# 生成第四个团和对应的标签值
class1_x4 = np.concatenate([stats.norm.rvs(20, 2, (100, 1)), stats.norm.rvs(10, 2, (100, 1))], axis=1)
y4 = np.ones_like(class1_x4, shape=(class1_x4.shape[0], 1)) * 3
#将所有的x合在一起作为训练数据，将所有的y合在一起作为训练数据
y=np.concatenate([y1,y2,y3,y4],axis=0)
x=np.concatenate([class1_x1,class1_x2,class1_x3,class1_x4],axis=0)
model=Logistic_Regression(x,y) #初始化模型并将数据传入
model.train(100000) #训练
print(model.w) #打印模型参数，最后一个维度为偏置项b
result=model.predict(x) #预测
plot_figure(x,result) #将预测结果画出
if __name__ == '__main__':
main()
结束语：
好啦，以上就是这篇文章的全部内容了。笔者也是在学习中，对于公式如何进行矩阵化的表达是很多大佬都忽略掉的问题。当燃也有可能是我太笨了，才会产生这种问题，亦或者是对线性代数缺乏感性的认识。也许对于很多人而言，可窥一眼而知全貌。但不管怎么说，我都希望能够帮助到与我有着同样的小伙伴。另外，本文公式中有一些不严谨的地方，有任何错误还望能够不吝赐教。就这样把，阿里嘎多（用LateX敲公式真累啊，什么时候优快云能直接上传word文档）。
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
