而在这个笔记中，我想谈的是 Logistic Regression( 逻辑回归)。

假定我们有如下两个clusters的数据，一堆是红的那些点，一堆是绿色的那堆，它们拥有如下的分布：

我们规定，其中绿色的点代表了某种正逻辑，用“1” 表示；红色代表负逻辑，用“0"表示，称为数据的标定（label）。显而易见，这两个clusters在特征空间(feature space)上是线性可分的（ 不可分的以后会提及）。意味着我们可以用类似如下的一条直线来分开两个clusters，这样再有一个新的未归类的点出现，我们可以根据它落在线的哪一边来给它着色，也就是分类。

既然是叫 Logistic Regression，那么肯定是要从这堆数据中想办法拟合出这条线最合理的位置和走向的，但具体怎么实现呢？ 我们再接着看，假如我们沿着如图的垂线方向来看点的label值的分布。

理想情况下应该会是如下这样的一个阶越分布，即线的一边全是正样本，一边全是负样本，记住绿=1，红=0：

看上去很美好，但头疼的是这玩意在关键的分隔点上不可导，会给我们后面的工作带来一系列麻烦( 但实际上仍然有分类器利用了这种分布，叫 perceptron)。
于是我们祭出了一条祖传的无比优雅的近似曲线——Logistic function:
$$h(x)=\frac{1}{1+e^{-x}}$$

其导数形式如下，有兴趣者自行推导：
$$h'(x)=h(x)(1-h(x))$$

假设这条直线方程是：
$$wx+b=0$$ 那么如下的式子可以看作对任意数据点\(x\)离开这条直线垂直距离：
$$\theta=w^Tx+b$$ 所以这里我们得到一个式子：
$$P(y=1|x)=h(w^Tx+b)$$ 这里的\(P(y=1|x)\)是，当给定一个输入数据x时，x属于正样本区间y=1的后验概率分布，也就是我们这里的classification所需要的映射关系。顺便一提，这种研究\(P(y|x)\)而非\(P(y，x)\)的机器学习方法，被称作***Discriminative model***，而后者被称作***Generative model***，这两者的区别和特点是非常有趣的一个话题，现按下不表。

\(P(y=1|x)\)显而易见是w,b的函数，所以我们写作\(P(y=1|x;w,b)\)。实际训练过程，可以看作对w,b的取值进行优化，使\(P(y=1|x;w,b)\)最符合我们的训练集数据\( \{x_i,y_i \}\)。

于是问题来了，如何描述这种符合程度？

这里引入了损失函数（loss function）的概念。

对所有的\(x_i\)和\(y_i\)，我们希望以下这个总概率最大化(参考高中数学知识，独立随机事件的总概率)：

$$P(Y|X)=\Pi_i P(y=y_i|x_i)$$ 即：
$$w,b=argmax_{w,b}\Pi_i P(y=y_i|x_i;w,b)$$

\(P(y=y_i|x_i;w,b)\)实际上可以写成如下表示：
$$P(y=y_i|x_i;w,b)=h(w^Tx_i+b)^y_i(1-h(w^Tx_i+b))^{1-y_i}$$ 请自己体会。
但处于某种原因，直接对\(P(Y|X)\)是很困难的，于是我们使用了如下的trick（认真上过高数课的同学并不会陌生），转而优化\(L(w,b)\)：
$$L(w,b)=-\sum_i log(h(w^Tx_i+b)^{y_i}(1-h(w^Tx_i+b))^{1-y_i})$$ 但同时优化两个参数w,b是㯊麻烦的（求两遍导，更新两个参量），所以又使用了如下的 bias trick:
$$w'=[w,b], x'=[x,1]$$ $$w^Tx+b=w'^Tx'$$ 简化后有没有神清气爽？
$$L(w)=-\sum_i log(h(w^Tx_i)^{y_i}(1-h(w^Tx_i))^{1-y_i})$$

拿到了\(L(w)\)的表达式，下面该如何找到正确的w，使得\(L(w)\)最小呢（注意，求最大变成了求最小）？

祭出凸优化中神挡杀神，佛挡杀佛的梯度下降法（gradient descent）。



图片来源
思路简单说就是，虽然我不知道w取坐标轴上什么位置好，但我可以一小步一小步的挪过去，每挪一步都向着当前所在位置坡度最陡的地方去。于是我们就这样慢慢的圆润的滚道谷底了。

可以证明，“坡度最陡”的方向，就是该点导数方向的反方向，这本书的这个章节有极好的详述。于是我们要对\(L(w)=-\sum_i log(h(w^Tx_i)^{y_i}(1-h(w^Tx_i))^{1-y_i})\)求导：
$$\frac{\partial L(w)}{\partial w}=\sum_i (h(w^Tx_i)-y_i)x_i$$ 得益于sigmoid函数良好的求导特性和log函数的特性，这个导数表达非常简洁。

在实际程序中，我们只需要算出每个数据点对应的\((h(w^Tx_i)-y_i)x_i\)，然后求和，就能得到当前w的导数表达。

但是这样真的好吗？

首先训练集的数据量一般都非常大，达到成千上万的地步，对如此大的数组频繁反复求和从计算的角度讲是非常不利的。其次，更重要的一点，我们无法保证整个数据集的数据都是有效数据，其中必然包含噪声，而把这些噪声算入我们每次计算的导数中，是非常危险的。

于是，这里我们使用另一个思路（所以我觉得机器学习里边的方法论实在是占据太多地盘了。。），每次求导时从整个数据集中随机抽取一小部分数据，用这一小部分数据作为样本来计算导数

这就是在各种机器学习算法中一路通杀的大名鼎鼎的 stochastic gradient descent，简称SGD，这种形式也叫 mini batch 法。无论是简单如Logistic regression，还是复杂如深度学习网络，优化思路基本上都是用SGD，或者说没有脱离SGD的模子。

光说不练从来不是科创的风格，下面讲一讲 Logistic Regression 的程序实现和实验。

在实验中我使用 Python，在 scientific computing 上，Python一直是主流，机器学习也不能免俗，而 numpy 提供的便捷的向量化运算，给诸多机器学习的算法实施提供了便利。另外强烈推荐ML新手装 sklearn 这个package。

导入库：

import numpy as np
from numpy import linalg as LA
import matplotlib.pyplot as plt
import sklearn.datasets
import sklearn.cross_validation

以下是计算梯度的函数，使用了numpy的向量化计算特性：

def calculate_gradient(w,x_batch,y_batch):
    sigmoid=1/(1+np.dot(x_batch,np.transpose(w)))
    dL=np.dot(sigmoid-y_batch,x_batch)/y_batch.size
    return dL

计算Loss function，注意，数据溢出是Logistic Regression程序实现的一个主要问题，因为exp函数的输出，用float 表示时，实际上输入被限制在[-750,750]这个区间内，不做处理的话基本上肯定会上溢。这个问题的解决同样在这本书中有讲：

def calculate_loss(w,x_all,y_all):
    ### Avoid Overflow! ###
    pos_index=np.where(y_all==1)    
    neg_index=np.where(y_all==0)
    Loss=np.sum(-np.log(1+np.exp(-np.dot(x_all[pos_index,:],np.transpose(w)))))+np.sum(-np.log(1+np.exp(-np.dot(x_all[pos_index,:],np.transpose(w)))))
    return Loss   

训练过程主循环，注意对learning rate采用了annealing，在SGD过程中分段一点点减小步长，不然很容易最后变成在global minimum周围反复徘徊难以收敛：

def train(x_train,y_train,alpha,batch_sz,loss_thresh,Max_iter,w0):
    ### bias trick ###
    w=w0
    data_sz=y_train.size
    x_train_b=np.concatenate((x_train,np.ones((data_sz,1))),axis=1)
    Loss_old=0
    Loss=[]
    stepCnt=0
    ### Run SGD ###
    for iter in range(1,Max_iter):
        ### sample a mini batch ###
        batch=np.arange(data_sz)
        np.random.shuffle(batch)
        x_batch=x_train_b[batch[:batch_sz],:]
        y_batch=y_train[batch[:batch_sz]]
        ### update weight ###
        dL=calculate_gradient(w,x_batch,y_batch)
        w-=alpha*dL
        ### record loss changes ###
        Loss.append(calculate_loss(w,x_train_b,y_train))
        ### learning rate annealing ###
        stepCnt+=1
        if stepCnt==10:
            stepCnt=0
            alpha*=0.8

        ### Check if converge ###
        if abs(Loss[-1]-Loss_old)<loss_thresh: break loss_old="Loss[-1]" return w,loss < code></loss_thresh:>

使用了sklearn的数据生成函数，之前的图表数据皆来源于此：

def make_data():
    centers = [(-10, -10),(10, 10)]
    x, y = sklearn.datasets.make_blobs(n_samples=2000, n_features=2, cluster_std=5.0,
                  centers=centers, shuffle=False, random_state=100)
    x_train, x_test, y_train, y_test = sklearn.cross_validation.train_test_split(x, y, test_size=.4)
    return x_train, x_test, y_train, y_test

主函数，w被初始化为全零向量，mini batch size 是50 :

def main():
    alpha=0.5
    batch_sz=50
    Max_iter=2000
    loss_thresh=1e-5
    w0=[0,0,0]
    x_train, x_test, y_train, y_test = make_data()
    w,Loss = train(x_train,y_train,alpha,batch_sz,loss_thresh,Max_iter,w0)
    plt.plot(Loss)

以下是记录的学习过程中 Loss function 的收敛过程。

得到的w的值为[ 123.01618818 125.42445694 11.78221087]。

画成直线可以看出，跟理想情况非常近似