更多内容关注公众号：数学的旋律

tb店铺搜：FUN STORE玩物社，专业买手挑选送礼好物

引言

    线性模型（linear modal）试图学得一个通过属性的线性组合来进行预测的函数。本文介绍两种经典的线性模型，分别是回归任务中的线性回归（linear regression）与二分类任务中的逻辑回归（logistic regression）。
    如图1，在二维空间中有一些样本点，我们用一条直线对这些点进行拟合，该直线称为最佳拟合直线。线性回归就是根据训练集，寻找对训练样本的最佳拟合直线；逻辑回归则是利用“单位跳跃函数”将线性模型产生的预测值转换为0/1值，从而实现二分类任务。

图1

一、数学预备知识

1.方向导数

如果函数$f(x,y)$在点$P_0(x_0,y_0)$可微分，那么函数在该点沿任一方向$l$的方向导数存在，且有$${\frac{\partial f}{\partial l}|}_{(x_0,y_0)}=f_x(x_0,y_0)\cos \alpha +f_y(x_0,y_0)\cos \beta$$其中$\cos \alpha$，$\cos \beta$是方向$l$的方向余弦。
（α为$l$与x轴正向夹角，β为$l$与y轴正向夹角。方向导数反映的是函数沿任一指定方向的变化率问题）

例1
求函数$z=x{e}^{2y}$在点$P(1,0)$处沿从点$P(1,0)$到点$Q(2,-1)$的方向的方向导数
解：
这里方向$l$即向量$\overrightarrow{PQ}=(1,-1)$的方向，与$l$同向的单位向量为$e_l=(\frac{1}{\sqrt{2}},-\frac{1}{\sqrt{2}})$。
因为函数可微分，且$$Z_x(1,0)=e^{2y}=1$$$$Z_y(1,0)=2x{e}^{2y}=2$$故所求方向导数为$${\frac{\partial Z}{\partial l}|}_{(1,0)}=1\times \frac{1}{\sqrt{2}}+2\times (-\frac{1}{\sqrt{2}})=-\frac{\sqrt{2}}{2}$$

2.梯度

与方向导数有关联的一个概念是函数的梯度。在二元函数的情形，设函数$f(x,y)$在平面区域D内具有一阶连续偏导数，则对于每一点$P_0(x_0,y_0)\in D$，都可定出一个向量$$f_x(x_0,y_0)\mathbf{i}+f_y(x_0,y_0)\mathbf{j}$$这向量称为函数$f(x,y)$在点$P_0(x_0,y_0)$的梯度，记作$\mathbf{grad}f(x_0,y_0)$或$\nabla f(x_0,y_0)$，即$$\mathbf{grad}f(x_0,y_0)=\nabla f(x_0,y_0)=f_x(x_0,y_0)\mathbf{i}+f_y(x_0,y_0)\mathbf{j}$$其中$\nabla =\frac{\partial }{\partial x}\mathbf{i}+\frac{\partial }{\partial y}\mathbf{j}$称为（二维的）向量微分算子或Nabla算子，$\nabla f=\frac{\partial f}{\partial x}\mathbf{i}+\frac{\partial f}{\partial y}\mathbf{j}$。
如果函数$f(x,y)$在点$P_0(x_0,y_0)$可微分，$e_l=(\cos \alpha ,\cos \beta )$是与方向$l$同向的单位向量，则$${\frac{\partial f}{\partial l}|}_{(x_0,y_0)}=f_x(x_0,y_0)\cos \alpha +f_y(x_0,y_0)\cos \beta =\mathbf{grad}f(x_0,y_0)\cdot e_l=|\mathbf{grad}f(x_0,y_0)|\cos \theta$$其中$\theta =<\mathbf{grad}f(x_0,y_0),e_l>$。
这一关系式表明了函数在一点的梯度与函数在这点的方向导数间的关系。特别的：
①$\theta =0$，即方向$e_l$与梯度$\mathbf{grad}f(x_0,y_0)$的方向相同时，函数$f(x,y)$增加最快。此时，函数在这个方向的方向导数达到最大值，这个最大值就是梯度$\mathbf{grad}f(x_0,y_0)$的模，即$${\frac{\partial f}{\partial l}|}_{(x_0,y_0)}=|\mathbf{grad}f(x_0,y_0)|$$②$\theta =\pi$，即方向$e_l$与梯度$\mathbf{grad}f(x_0,y_0)$的方向相反时，函数f$(x,y)$减少最快，函数在这个方向的方向导数达到最小值，即$${\frac{\partial f}{\partial l}|}_{(x_0,y_0)}=-|\mathbf{grad}f(x_0,y_0)|$$③$\theta =\frac{\pi }{2}$，即方向$e_l$与梯度$\mathbf{grad}f(x_0,y_0)$的方向正交时，函数的变化率为零，即$${\frac{\partial f}{\partial l}|}_{(x_0,y_0)}=|\mathbf{grad}f(x_0,y_0)|\cos \theta =0$$例2
设$f(x,y)=\frac{1}{2}(x^2+y^2)$，$P_0(1,1)$，求
①$f(x,y)$在$P_0$处增加最快的方向以及$f(x,y)$沿这个方向的方向导数；
②$f(x,y)$在$P_0$处减少最快的方向以及$f(x,y)$沿这个方向的方向导数；
③$f(x,y)$在$P_0$处的变化率为零的方向；
解：
①$f(x,y)$在$P_0$处沿$\nabla f(1,1)$的方向增加最快$$\nabla f(1,1)=(x\mathbf{i}+y\mathbf{j}){|}_{(1,1)}=\mathbf{i}+\mathbf{j}$$故所求方向可取为$$n=\frac{\nabla f(1,1)}{|\nabla f(1,1)|}=\frac{1}{\sqrt{2}}\mathbf{i}+\frac{1}{\sqrt{2}}\mathbf{j}$$方向导数为$${\frac{\partial f}{\partial n}|}_{(1,1)}=|\nabla f(1,1)|=\sqrt{2}$$②$f(x,y)$在$P_0$处沿$-\nabla f(1,1)$的方向减少最快，这个方向可取为$$n_1=-n=-\frac{1}{\sqrt{2}}\mathbf{i}-\frac{1}{\sqrt{2}}\mathbf{j}$$方向导数为$${\frac{\partial f}{\partial n_1}|}_{(1,1)}=-|\nabla f(1,1)|=-\sqrt{2}$$③$f(x,y)$在$P_0$处沿垂直于$\nabla f(1,1)$的方向变化率为零，这个方向是$$n_2=-\frac{1}{\sqrt{2}}\mathbf{i}+\frac{1}{\sqrt{2}}\mathbf{j}$$或$$n_3=\frac{1}{\sqrt{2}}\mathbf{i}-\frac{1}{\sqrt{2}}\mathbf{j}$$

3.带佩亚诺余项的泰勒公式

若$f(x)$在$x=a$有$n$阶微商，即$f^{(n)}(a)$存在，则当$x\to a$时有$$f(x)=f(a)+f^{\prime }(a)(x-a)+\frac{f^{\prime \prime }(a)}{2!}(x-a{)}^2+\cdots +\frac{f^{(n)}(a)}{n!}(x-a{)}^n+\mathrm{o}((\mathrm{x}-\mathrm{a}{)}^{\mathrm{n}})$$其中多项式$$f(a)+f^{\prime }(a)(x-a)+\cdots +\frac{f^{(n)}(a)}{n!}(x-a{)}^n$$称为$f(x)$在$a$点的$n$阶泰勒多项式。

4.二元函数的泰勒公式

设$f(x,y)$在$P_0(x_0,y_0)$的某领域$\mathrm{O}({\mathrm{P}}_0)$内有直到$n+1$阶连续偏导数，则对$\mathrm{O}({\mathrm{P}}_0)$内任一点$(x_0+\Delta x,y_0+\Delta y)$，存在$\theta \in (0,1)$，使得$$f(x_0+\Delta x,y_0+\Delta y)=\sum _{k=0}^n\frac{1}{k!}{\left(\frac{\partial }{\partial x}\Delta x+\frac{\partial }{\partial y}\Delta y\right)}^{k}f(x_0,y_0)+R_n$$其中$$R_n=\frac{1}{(n+1)!}{\left(\frac{\partial }{\partial x}\Delta x+\frac{\partial }{\partial y}\Delta y\right)}^{n+1}f(x_0+\theta \Delta x,y_0+\theta \Delta y)$$当$n=1$，忽略余项$R_n$，即得函数的一阶泰勒展开式$$f(x_0+\Delta x,y_0+\Delta y)\approx f(x_0,y_0)+{f^{\prime }}_x(x_0,y_0)\Delta x+{f^{\prime }}_y(x_0,y_0)\Delta y$$利用泰勒公式，可以用$n$阶多项式逼近函数$f(x,y)$，$n$越大则误差越小，近似程度越好。

二、线性回归

1.最小二乘法

    给定训练集$D=\{(x_1,y_1),(x_2,y_2),\cdots ,(x_N,y_N)\}$，其中$x_i=(x_i^{(1)},x_i^{(2)},\cdots ,x_i^{(n)}{)}^T$，$x_i^{(j)}$是第$i$个实例的第$j$个特征。线性回归试图学得$$f(x_i)=w^T{x}_i+b，使得f(x_i)\approx y_i$$其中$w=(w_1,w_2,...,w_n{)}^T$。
    如何确定$w$和$b$呢？一个常用的方法是极小化其平方损失函数，即使均方误差最小化，得$$(w^{\ast },b^{\ast })=arg\underset{(w,b)}{\min }\sum _{i=1}^N(f(x_i)-y_i{)}^2=arg\underset{(w,b)}{\min }\sum _{i=1}^N(y_i-w^T{x}_i-b{)}^2$$    基于均方误差最小化来进行模型求解的方法称为“最小二乘法”。均方误差有非常好的几何意义，它对应了欧几里得距离，在线性回归中，最小二乘法就是试图找到一条直线，使得所有样本到直线上的欧几里得距离之和最小。

2.求解w和b

    为了讨论方便，我们把$w$和$b$吸收入向量形式$\hat{w}=(w^T;b{)}^T$。把训练集$D$表示为一个$N\times (n+1)$大小的矩阵$X$，其中每行对应于一个实例$x_i$，该行前$n$个元素对应于$x_i$的$n$个属性，最后一个元素恒置为1。即

再把类别也写成向量形式$y=(y_1,y_2,\cdots ,y_N{)}^T$，则有$${\hat{w}}^{\ast }=arg\underset{\hat{w}}{\min }(y-X\hat{w}{)}^T(y-X\hat{w})$$令$E=(y-X\hat{w}{)}^T(y-X\hat{w})$，对$\hat{w}$求导得到$$E_{\hat{w}}=2X^T(X\hat{w}-y)$$令上式为零可求得$\hat{w}$的最优解。当$X^{T}X$为满秩矩阵时，可解得$$\hat{w}=(X^{T}X{)}^{-1}{X}^{T}y(1)$$从而求解出$w$和$b$，得到线性回归模型：$$y=w^{T}x+b$$

3.岭回归（ridge regression）

    然而，现实任务中$X^{T}X$往往不是满秩矩阵，无法计算$(X^{T}X{)}^{-1}$。为了解决这个问题，统计学家引入岭回归。简单说来，岭回归就是在矩阵$X^{T}X$上加入一个$\lambda I$，使得能对$X^{T}X+\lambda I$求逆。其中$I$是$N\times N$的单位矩阵。则$$\hat{w}=(X^{T}X+\lambda I{)}^{-1}{X}^{T}y$$    这里通过引入$\lambda$来限制所有$w$之和，通过引入该惩罚项，能够减少不重要的参数，这个技术在统计学中也叫做缩减（shrinkage）。可以通过交叉验证法，选取不同的$\lambda$进行测试，最终得到一个使预测误差最小化的λ。

4.代码实现（python）

以下代码来自Peter Harrington《Machine Learing in Action》。
代码如下（保存为regression.py）：

# -- coding: utf-8 --
from numpy import *

def loadDataSet(fileName):
    dataMat = []; labelMat = []
    fr = open(fileName)
    for line in fr.readlines():
        lineArr = line.strip().split()
        dataMat.append([float(lineArr[0]), 1.0]) # dataMat为吸收后的X，即将特征值最后一个元素恒置为1
        labelMat.append(float(lineArr[1]))       # labelMat为对应的标记类别
    return dataMat,labelMat

def standRegres(xArr,yArr):
    # xArr为特征值，yArr为对应的标记类别，该函数用于计算w和b
    xMat = mat(xArr); yMat = mat(yArr).T
    xTx = xMat.T*xMat               # 计算xTx
    if linalg.det(xTx) == 0.0:
        print "This matrix is singular, cannot do inverse"
        return
    ws = xTx.I * (xMat.T*yMat)      # 若xTx不为零，则根据式(1)计算w^
    return ws

数据集格式如下（保存为ex0.txt）：

0.067732	3.176513
0.427810	3.816464
0.995731	4.550095
0.738336	4.256571
0.981083	4.560815
0.526171	3.929515

运行命令如下：

最后求出的$ws$即$\hat{w}$。

三、逻辑回归

逻辑回归是一种分类模型，本文主要介绍二分类的逻辑回归。

1.Sigmoid函数

    给定训练集$D=\{(x_1,y_1),(x_2,y_2),\cdots ,(x_N,y_N)\}$，其中$x_i=(x_i^{(1)},x_i^{(2)},\cdots ,x_i^{(n)}{)}^T$，$x_i^{(j)}$是第$i$个实例的第$j$个特征，$y_i\in 0,1$为实例的类别。若用线性回归模型$y=w^{T}x+b$对该训练集进行分类，分类图形如图2线A（此处设$x_i=(x_i^{(1)}{)}^T$）所示：

图2

    从上图看出，线性回归模型产生的预测值$z=w^{T}x+b$是实值，于是，我们需要将实值$z$转换为0/1值，这里我们引入Sigmoid函数，其数学形式是$$g(z)=\frac{1}{1+e^{-z}}$$    图3给出了Sigmoid函数的坐标图。当$z$为0时，Sigmoid函数值为0.5，随着$z$的增大，对应的Sigmoid值将逼近于1；而随着$z$的减少，Sigmoid值将逼近于0。

图3

我们将z=wTx+b代入Sigmoid函数，从而将z值转化为一个接近0或1的数，得到
$$y=\frac{1}{1+e^{-(w^{T}x+b)}}$$

2.模型

    逻辑回归模型所做的假设是，将(1)中$y$视为类后验估计$P(Y=1|x)$，得到$$P(Y=1|x)=\frac{1}{1+e^{-(w^{T}x+b)}}=\frac{e^{w^{T}x+b}}{1+e^{w^{T}x+b}}$$    通过$P(Y=0|x)=1-P(Y=1|x)$求得$$P(Y=0|x)=\frac{1}{1+e^{w^{T}x+b}}$$    对于给定的实例$x$，可以求得$P(Y=1|x)$和$P(Y=0|x)$，比较两个条件概率值的大小，将实例$x$分到概率大的那一类，若概率值相等则可任意判别。

3.损失函数

    这里使用对数损失函数$L(Y,P(Y|X))=-logP(Y|X)$。$P(Y|X)$的取值范围为[0,1]，$P(Y|X)$的值越接近1，则损失函数值越小。
    为了讨论方便，我们把w和b吸收入向量形式$\hat{w}=(w^T;b{)}^T$，令$x_i=(x_i^{(1)},x_i^{(2)},\cdots ,x_i^{(n)},1{)}^T$，则$w^{T}x+b$可简写为$\hat{w}x$。设$P(Y=1|x)=\pi (x)$，$P(Y=0|x)=1-\pi (x)$，则最终逻辑回归模型的损失函数为$$L(\hat{w})=\sum _{i=1}^N[-y_i\log \pi (x_i)-(1-y_i)\log (1-\pi (x_i))]$$$$=-\sum _{i=1}^N[y_i\log \pi (x_i)-y_i\log (1-\pi (x_i))+\log (1-\pi (x_i))](2)$$$$=-\sum _{i=1}^N\left[y_i\log \frac{\pi (x_i)}{1-\pi (x_i)}+\log (1-\pi (x_i))\right]=-\sum _{i=1}^N\left[y_i({\hat{w}}^T{x}_i)-log(1+e^{{\hat{w}}^T{x}_i})\right]$$    对损失函数$L(\hat{w})$求极小值，从而可得$\hat{w}$。
    最后求得的直线应该如图2线B，分类为0的样本点使得z<0，此时$P(Y=1|x)<P(Y=0|x)$；分类为1的样本点使得z>0，此时$P(Y=1|x)>P(Y=0|x)$。

4.梯度下降法

    假设$f(x$)是$R_n$上具有一阶连续偏导数的函数，要求解的无约束优化问题是$$\underset{x\in R_n}{\min }f(x)$$$x^{\ast }$表示目标函数$f(x)$的极小点。
    梯度下降算法是一种迭代算法，选取适当的初值$x^{(0)}$，不断迭代，更新$x$的值，进行目标函数的极小化，直到收敛。由于负梯度方向是使函数值下降最快的方向，在迭代的每一步，以负梯度方向更新$x$的值，从而达到减少函数值的目的。
    由于$f(x)$具有一阶连续偏导数，若第$k$次迭代值为$x^{(k)}$，则可将$f(x)$在$x^{(k)}$附近进行一阶泰勒展开：$$f(x)=f(x^{(k)})+{g_k}^T(x-x^{(k)})$$这里，$g_k=g(x^{(k)})=\nabla f(x^{(k)})$为$f(x)$在$x^{(k)}$的梯度。
    求出第k+1次迭代值x(k+1)：$$x^{(k+1)}=x^{(k)}+{\lambda }_k{p}_k$$其中，$p_k$是搜索方向，取负梯度方向$p_k=-\nabla f(x^{(k)})$，${\lambda }_k$是步长，由一维搜索确定，即${\lambda }_k$使得$$f(x^{(k)}+{\lambda }_k{p}_k)=\underset{\lambda \ge 0}{\min }f(x^{(k)}+\lambda p_k)$$算法1（梯度下降法）
输入：目标函数$f(x)$，梯度函数$g(x)=\nabla f(x)$，计算精度$\epsilon$
输出：$f(x)$的极小点$x^{\ast }$
①取初始值$x^{(0)}\in R^n$，置$k=0$
②计算$f(x^{(k)})$
③计算梯度$g_k=g(x^{(k)})$，当$||g_k||<\epsilon$，停止迭代，令$x^{\ast }=x^{(k)}$；否则，令$p_k=-\nabla f(x^{(k)})$，求$\lambda k$，使$$f(x^{(k)}+{\lambda }_k{p}_k)=\underset{\lambda \ge 0}{\min }f(x^{(k)}+\lambda p_k)$$④置$x^{(k+1)}=x^{(k)}+{\lambda }_k{p}_k$，计算$f(x^{(k+1)})$，当$||f(x(k+1))-f(x(k))||<\epsilon$或$||x(k+1)-x(k)||<\epsilon$时，停止迭代，令$x^{\ast }=x^{(k+1)}$
⑤否则，置$k=k+1$，转③

5.代码实现（python）

以下代码来自Peter Harrington《Machine Learing in Action》。
    我们需要求式(2)中$L(\hat{w})$的最小值。若令$f(\hat{w})=-L(\hat{w})$，则问题可转换为求$f(\hat{w})$的最大值：$$f(\hat{w})=\sum _{i=1}^N\left[y_i({\hat{w}}^T{x}_i)-log(1+e^{{\hat{w}}^T{x}_i})\right]$$    为了求上式的最大值，我们使用梯度上升算法。算法1（梯度下降法）中，由于与梯度相反方向下降最快，所以$p_k=-\nabla f(x(k))$；而在本例的梯度上升算法，令$p_k=\nabla f(x(k))$，即梯度方向上升最快。先求$\nabla f(\hat{w})$：$$\nabla f(\hat{w})=\sum _{i=1}^N\left[y_i{x}_i-\frac{1}{1+e^{-{\hat{w}}^T{x}_i}}{x}_i\right]=\sum _{i=1}^N{x}_i\left[y_i-\frac{1}{1+e^{-{\hat{w}}^T{x}_i}}\right](3)$$    在本代码中，预先设置好步长$\lambda$的值和迭代次数，则$w^{\ast }=\lambda +\nabla f(\hat{w})$，在循环迭代完成后输出训练好的回归系数。
代码如下（保存为logRegres.py）：

# -- coding: utf-8 --
from numpy import *

def loadDataSet():
    dataMat = []; labelMat = []
    fr = open('testSet.txt')
    for line in fr.readlines():
        lineArr = line.strip().split()
        dataMat.append([float(lineArr[0]), float(lineArr[1]), 1.0]) # dataMat为吸收后的X，即将特征值最后一个元素恒置为1
        labelMat.append(int(lineArr[2])) # labelMat为对应的标记类别
    return dataMat,labelMat

def sigmoid(inX):
    return 1.0/(1+exp(-inX))

def gradAscent(dataMatIn, classLabels):
    dataMatrix = mat(dataMatIn)
    labelMat = mat(classLabels).transpose()
    m,n = shape(dataMatrix)
    alpha = 0.001          # 步长
    maxCycles = 500        # 迭代次数
    weights = ones((n,1))  # 将w的各个元素和b初始化为1
    for k in range(maxCycles):
        h = sigmoid(dataMatrix*weights)
        error = (labelMat - h)
        # 根据式(3)可知dataMatrix.transpose()* error为梯度，则可得（新的weights）=（旧的weights）+步长＊梯度
        weights = weights + alpha * dataMatrix.transpose()* error
    return weights

数据集如下（保存为testSet.txt）：

-1.395634	4.662541	1
-1.322371	7.152853	0
0.423363	11.054677	0
0.406704	7.067335	1
0.667394	12.741452	0
-2.460150	6.866805	1

调用gradAscent函数求得的值即为$\hat{w}$。









以上全部内容参考书籍如下：
李航《统计学习方法》
周志华《机器学习》
Peter Harrington《Machine Learing in Action》
《高等数学第六版下册 同济大学》
《数学分析简明教程 第二版上册》