李宏毅机器学习（二）

1、理解偏差(bias)和方差(variance)

首先，我们来理解什么是偏差(bias)和方差(variance)：

机器学习目的是找到一个最好函数 $\hat f$f^​，能够非常准确根据输入的数据得到所期望的结果。为了得到这个最好的函数，我们需要一堆训练数据，找到函数$f^*$f∗，$f^*$f∗是最优函数$\hat f$f^​的估测。衡量两者之间的误差，主要来自两个方面：偏差(bias)和方差(variance)。

从统计学来看，偏差和方差的来源：

• 目的：估计一个变量$x$x的平均值

  • 假设$x$x的平均值是$\mu$μ
  • 假设$x$x的方差是$\sigma^2$σ2

• 估测$\mu$μ的方法：

  • 取样本，样本为$N$N个点：$\{x^1, x^2,…, x^N\}${x1,x2,…,xN}
  • 将$N$N个点求平均值：$m = \frac1N\sum_Nx^n \neq \mu$m=N1​∑N​xn̸​=μ
  • 再取多次样本，得到$m_1,m_2,m_3...$m1​,m2​,m3​...
  • 然后求$m$m的期望值：$E[m] = E[\frac1N \sum_nx^n] = \frac1N\sum_nE[x^n]=\mu$E[m]=E[N1​∑n​xn]=N1​∑n​E[xn]=μ
  • 上述的$m_n$mn​都是无偏差的(unbiased)，但是散落在$\mu$μ周围，存在方差$Var[m] = \frac{\sigma^2}{N}$Var[m]=Nσ2​
  • 方差的值取决于样本大小$N$N，如果$N$N比较大，方差就比较小；反之，方差比较大。

• 估测$\sigma^2$σ2的方法：

  • $\sigma^2$σ2的估测值是$s^2 = \frac1N\sum_n(x^n-m)^2$s2=N1​∑n​(xn−m)2
  • 此时多个不同的$s_n^2$sn2​会散落在$\mu$μ周围，存在偏差：$E[s^2] = \frac{N-1}{N}\sigma^2 \neq \sigma^2$E[s2]=NN−1​σ2̸​=σ2.这个意思应该就是样本数据相比总体会更加紧密聚集在均值周围。【注：从中心极限定理中的样本估计整体可以理解这一问题。】
  • 如果把样本大小$N$N加大，方差估测值$s^2$s2与方差真实值$\sigma^2$σ2的差距会变小。

从机器学习寻找最优函数$\hat f$f^​来看，偏差和方差来源：

• 为了找最优函数$\hat f$f^​，我们通过训练数据找到了一些函数$f^*$f∗，这些函数$f^*$f∗与最优函数在预测数据的时候会产生偏差，这就是误差的来源之一偏差(bias)了。
• 假设找到的所有函数$f^*$f∗的期望值：$E[f^*]=\overline f$E[f∗]=f​，这些函数$f^*$f∗之间又会存在偏离，这些偏离就是方差(variance)。
• 以打靶为例，因为准星没对准靶心而产生的偏离就是偏差，准星对准了，但是因为其它因素导致偏离就是方差。讲义中下面图就是一个很好的例子：
  

2、偏差(bias)、方差(variance)与Model的关系

2.1 重温模型(Model)和函数(Function)

在上次学习中，没弄清楚Model和function的关系。直到在这次学习中才明白两者的差异：

• Model：我们在面对一个机器学习问题时设计的一个函数形式，例如：$y = b + w.x$y=b+w.x

  其中的$b$b和$w$w都是未知的，我们需要求的也就是它们的值。

• Function：针对设计好的模型，我们使用不同的训练数据，会求到不同的$w$w和$b$b的值，每一组值就对应一个函数(Function)。

所以能够使用同一个Model，得到不同的Function。

2.2 过拟合、欠拟合与偏差(bias)、方差(variance)

首先说明三类函数： 理想的最优函数：$\hat f$f^​，一个Model通过训练数据获取的函数：$f^*$f∗，函数$f^*$f∗的期望： $\overline f$f​。

• **欠拟合：**根据上一节学习中可知，Model设计得太简单，会造成欠拟合，此时Model误差（error）比较大，偏差（$\overline f$f​与 $\hat f$f^​）也比较大，但是Model的各个函数（$f^*$f∗）散布比较紧密，即方差（variance）比较小。
  • 原因很简单，简单Model参数比较少，受不同样本数据的影响也就比较小。
  • 以一个极端例子为例，设计Model为：$f(x) = c$f(x)=c，($c$c是常数)。由此可见，不管训练数据是什么，Model的函数是不变的，方差为0，但是离最优的function差距很大。
• **过拟合：**Model设计太复杂，造成过拟合。此时误差（error）也很大，但是偏差($\overline f$f​与$\hat f$f^​)比较小，此时Model的各个函数($f^*$f∗)散布比较松散，方差比较大。
  • 原因类似，复杂Model参数比较多，受不同样本数据的影响也就比较大。
• 方差(variance)与Model复杂程度关系图：
  
• 偏差((bias)与Model复杂程度关系图：
  
  总结，老师在课程中用一张图说明了它们之间的关系：
  

最优的函数，应该就是Error最小的那点。

由此可见，我们在训练模型的时候，应该了解到Error的来源是Bias还是Variance。

2.3 诊断Error以及优化Model

• 在训练的时候，如果训练数据误差比较大，此时就是偏差比较大了，处于欠拟合，需要重新设计你的模型：
  • 输入更多的特征；
  • 设计一个更复杂的模型。
• 如果是方差比较大，训练Model时体现在训练数据偏差小，但是在测试数据上偏差比较大，方差很大。解决办法是：
  • 非常有效的方法：收集更多的数据；(可以创造新数据，调角度、上下、左右颠倒等)
  • 正则化。

3 、回顾：梯度下降

设计一个损失函数$L$L，其中输入是Model产生的函数$f^*$f∗（或者参数$w$w和$b$b），通过梯度下降，能够在损失函数上找到一个最优点，使得损失最小，从而找到最佳函数。但是，使用梯度下降时需要注意的几点：

3.1 调整学习率

需要小心设置学习率$\eta$η，如果学习率太大，可能造成错过全局最优点在之间来回摆动；如果太小，梯度调整速度太慢。因此，需要画图体现Loss的变化：

因此，动态调整learning rate是比较重要的：

• 一个简单且受欢迎的想法是：每几轮训练之后，通过一些因素减少学习率。

  • 训练开始时，距离目标点比较远，所以用比较大的学习率。
  • 在几轮训练之后，接近目标点，所以减少学习率。
  • 例子：学习率$\eta$η随更新次数$t$t的变化：$\eta^t=\eta/\sqrt{t+1}$ηt=η/t+1​

• 学习率不应该只设置一个。

  • 给不同参数不同的学习率：

  • 参数$w$w更新数学推导如下：

    首先，定义$t$t为更新次数，$\eta^t =\frac{\eta}{\sqrt{t+1}}$ηt=t+1​η​是第$t$t次的学习率，$g^t = \frac{dL(\theta^t)}{dw}$gt=dwdL(θt)​是第$t$t次时损失函数对参数$w$w的偏微分，也就是梯度，$\sigma^t$σt是前$t$t项$g^t$gt的均方根。

    1. $w^1=w^0-\frac{\eta^0}{\sigma^0}g^0$w1=w0−σ0η0​g0，$\sigma^0 = \sqrt{(g^0)^2}$σ0=(g0)2​

    2. $w^2 = w^1-\frac{\eta^1}{\sigma^1}g^1$w2=w1−σ1η1​g1，$\sigma^1 = \sqrt{\frac12[(g^0)^2+(g^1)^2]}$σ1=21​[(g0)2+(g1)2]​

    3. $w^3 = w^2-\frac{\eta^2}{\sigma^2}g^2$w3=w2−σ2η2​g2，$\sigma^2 = \sqrt{\frac12[(g^0)^2+(g^1)^2+(g^2)^2]}$σ2=21​[(g0)2+(g1)2+(g2)2]​

       …

    4. $w^{t+1} = w^t-\frac{\eta^t}{\sigma^t}g^t$wt+1=wt−σtηt​gt，$\sigma^t = \sqrt{\frac{1}{t+1}\sum_{i=0}^t(g^t)^2}$σt=t+11​∑i=0t​(gt)2​

  • 最后，推导出式子：$w^{t+1}=w^t-\frac{\eta}{\sqrt{\sum^t_{t=0}(g^i)^2}}g^t$wt+1=wt−∑t=0t​(gi)2​η​gt

3.2 随机梯度下降

损失函数的形式：$L = \sum_n(\hat y^n-(b+\sum w_ix_i^n))^2$L=∑n​(y^​n−(b+∑wi​xin​))2。损失函数确实应该考虑所有训练样本数据的误差和。下面说明梯度下降和随机梯度下降的区别：

• 梯度下降：$\theta^i = \theta^{i-1}-\eta \nabla L(\theta^{i-1})$θi=θi−1−η∇L(θi−1)，训练完所有样本之后，再更新参数；
• 随机梯度下降：随机或按照顺序取一个样本$x^n$xn，只考虑这个样本的损失：
  $L^n = (\hat y^n-(b+\sum w_ix_i^n))^2$Ln=(y^​n−(b+∑wi​xin​))2 ，$\theta^i = \theta^{i-1}-\eta \nabla L(\theta^{i-1})$θi=θi−1−η∇L(θi−1)
  一个样本训练之后，就更新依次参数。

3.3 数据归一化(Feature Scaling)

• **描述：**为了说明这个知识点，按照课程的例子，假设有一个Model：$y = b+w_1x_1+w_2x_2$y=b+w1​x1​+w2​x2​。如果其中的两个特征值$x_1$x1​和$x_2$x2​的取值范围相差比较大，那就把其中一个数大小缩放到与另一个相同。这样做的原因是什么呢？看下图：
  
• 为了说明这个知识点，按照课程的例子，假设有一个Model：$y = b+w_1x_1+w_2x_2$y=b+w1​x1​+w2​x2​。如果其中的两个特征值$x_1$x1​和$x_2$x2​的取值范围相差比较大，那就把其中一个数大小缩放到与另一个相同。这样做的原因是什么呢？看下图：

可以看到如果$x_1$x1​和$x_2$x2​取值范围相差比较大，$w_1$w1​和$w_2$w2​对损失函数的影响不一样，会导致梯度下降时参数更新效率变慢，在更新时还需要考虑学习率的改变。反之，两者规模如果差不多，情况会好很多。

怎么做？

Model的输入特征值是$x^1x^2x^3…x^r…x^R$x1x2x3…xr…xR，每个特征值有一堆样本数据$x^r_i$xir​。

$m_i$mi​表示全部特征值第$i$i组数据的平均值；$\sigma_i$σi​表示全部特征值第$i$i组数据的标准差。

数据归一化的方式：$x^r_i = \frac{x^r_i-m_i}{\sigma_i}$xir​=σi​xir​−mi​​。

4、数学理论

4.1 梯度下降说明

梯度下降的原理其实就是在某个点附近找到一个相对最小的值，如下图所示：

这个的理论基础就是我们所学的泰勒公式了。

4.2 泰勒展开式

• 单变量泰勒展开式

  定义：函数$h(x)$h(x)在$x_0$x0​处的展开式如下所示
  $h(x)= \sum_ {k = 0}^{\infty} \frac {h ^ k（x_0）} {k！}（x-x_0）^ k = h （x_0）+ {h}'（x_0）（x-x_0）+ \frac {h''（x_0）} {2！}（x-x_0）^ 2+...$h(x)=k=0∑∞​k！hk（x0​）​（x−x0​）k=h（x0​）+h′（x0​）（x−x0​）+2！h′′（x0​）​（x−x0​）2+...
  当$x$x趋近于$x_0$x0​的时候，展开式形式换成：$h（x）≈ h（x_0）+ {h}'（x_0）（x-x_0）$h（x）≈h（x0​）+h′（x0​）（x−x0​）

• 多变量泰勒展开式

  课件中给出了展开式形式：
  

• 将泰勒展开式与损失函数相联系
  

图中损失函数随机取一点$(a, b)$(a,b)

• 损失函数的多变量泰勒展开形式：$L(\theta)≈L(a,b)=\frac{dL(a,b)}{d\theta_1}(\theta_1-a)+\frac{dL(a,b)}{d\theta_2}(\theta_2-b)$L(θ)≈L(a,b)=dθ1​dL(a,b)​(θ1​−a)+dθ2​dL(a,b)​(θ2​−b)

• 为了简化，令$s=L(a,b)$s=L(a,b)，$u=\frac{dL(a,b)}{d\theta_1}$u=dθ1​dL(a,b)​，$v=\frac{dL(a,b)}{d\theta_2}$v=dθ2​dL(a,b)​

• 化简后，形式是：$L(\theta)≈s+u(\theta_1-a)+v(\theta_2-b)$L(θ)≈s+u(θ1​−a)+v(θ2​−b)

• 此时，需要找的是小红圈内最小的损失函数，因此需要满足条件：

  $(\theta_1-a)^2$(θ1​−a)2+$(\theta_2-b)^2 \leq d^2$(θ2​−b)2≤d2
  

• 如果要求最小的损失函数，也就是下一个点$(\theta_1,\theta_2)$(θ1​,θ2​)，就变成了求两个向量$(u, v)$(u,v)和$(\theta_1-a,\theta_2-b)$(θ1​−a,θ2​−b)的内积最小值，如图所示，向量$(u, v)$(u,v)已知，$(\theta_1-a,\theta_2-b)$(θ1​−a,θ2​−b)在反方向时，内积最小：
  

• 整理之后，求最小的损失函数的式子如下：
  

  这就是梯度下降了。

• **注意：**只有红色圈圈的半径无穷小，式子才能保证成立。即学习率$\eta$η需要足够小，才能满足式子的条件。

5、梯度下降的限制

陷入局部最优，在局部最优点处，参数就停止更新了，此时虽然梯度为0，但这只是极值而不是最小值。如下图，参数处在一个比较高的Loss的时候，但是梯度已经很平缓了，会造成已经到了最优值的地方的假象。

6、学习总结

总结主要包含对此次学习过程中，对一些问题的理解。

6.1 batch、mini-batch与SGD

• Batch：对整个训练集进行梯度下降时，整个数据集就是一个batch。这时候就是batch梯度下降。
• SGD：当每次只对一个样本进行梯度下降的时候，是 随机梯度下降（SGD）。
• mini-batch：当每次处理样本的个数在上面二者之间，就是 mini batch 梯度下降。
• 优缺点分析：当数据集很大时，训练算法是非常慢的，和 batch 梯度下降相比，使用 mini batch 梯度下降更新参数更快，有利于更鲁棒地收敛，避免局部最优。和 stochastic 梯度下降相比，使用 mini batch 梯度下降的计算效率更高，可以帮助快速训练模型。
• 代码Mini-Batch SGD代码：(参考代码)

import numpy as np
import os

def get_data(obj_path_name):
    pro_path = os.path.abspath('.')
    data_path = str(pro_path + obj_path_name)
    print data_path
    data = np.loadtxt(data_path)
    x = np.asarray(np.column_stack((np.ones((data.shape[0], 1)), data[:, 0:-1])), dtype='double')
    y = data[:, -1]
    y = np.reshape(y, (data.shape[0], 1))
    print x.shape, y.shape
    return x, y


def inverse_hessian_mat(x):
    # 计算hessian矩阵，并求其逆矩阵
    x_hessian = np.dot(x.T, x)
    inverse_hessian = np.linalg.inv(x_hessian)
    return inverse_hessian


def get_e(x, theta, y):
    y_pre = np.dot(x, theta)
    e = y_pre - y

    return e


def compute_grad(x, e, stand_flag=0):
    # batch法必须做梯度归一化
    print e.T.shape, x.shape
    grad = np.dot(e.T, x)
    # grad = np.dot(e.T, x)/x.shape[0]
    if stand_flag == 1:
        grad = grad/(np.dot(grad, grad.T)**0.5)
    return np.reshape(grad, (x.shape[1], 1))


def get_cost(e):
    # print e.shape
    # 计算当前theta组合点下的各个样本的预测值 y_pre
    cost = np.dot(e.T, e) / 2
    # cost = np.dot(e.T, e) / (2*e.shape[0])
    return cost


def get_cost_l12(e, theta, m, l=1, lmd=10e10):
    # print e.shape
    if l == 1:
        cost = (np.dot(e.T, e) + lmd*sum(abs(theta))) / (2*m)
    elif l == 2:
        cost = (np.dot(e.T, e) + lmd*np.dot(theta.T*theta)) / (2*m)
    else:
        cost = (np.dot(e.T, e)) / (2*m)
    return cost


def update_batch_theta(theta, grad, alpha):
    theta = theta - alpha*grad
    return theta


def update_batch_theta_l12(theta, grad, alpha, m, l=1, lmd=10e10):
    if l == 1:
        theta = theta - alpha * (grad + (lmd/m)*theta)
    elif l == 2:
        theta = theta - alpha * (grad + (lmd/m))
    else:
        theta = theta - alpha * grad
    return theta


def iter_batch(x, theta, y, out_n, out_e_reduce_rate, alpha):
    step = 1
    while step < out_n:
        """计算初始的损失值"""
        if step == 1:
            e = get_e(x, theta, y)
            cost_0 = get_cost(e)
        """计算当前theta组合下的grad值"""
        grad = compute_grad(x, e, stand_flag=1)
        """依据grad更新theta"""
        theta = update_batch_theta(theta, grad, alpha)
        """计算新的损失值"""
        e = get_e(x, theta, y)
        cost_1 = get_cost(e)

        e_reduce_rate = abs(cost_1 - cost_0)/cost_0
        print 'Step: %-6s, cost: %s' % (step, cost_1[0, 0])
        if e_reduce_rate < out_e_reduce_rate:
            flag = 'Finish!\n'
            print flag
            break
        else:
            cost_0 = cost_1
            step += 1

    return theta


def iter_random_batch(x, theta, y, out_n, out_e_reduce_rate, alpha):

    step = 0
    while step < out_n:
        step += 1
        for i in range(x.shape[0]):
            x_i = np.reshape(x[i, ], (1, x.shape[1]))
            y_i = np.reshape(y[i, ], (1, 1))
            """计算初始的损失值"""
            e_0 = get_e(x, theta, y)
            cost_0 = get_cost(e_0)

            """用一个样本，计算当前theta组合下的grad值"""
            e_i = get_e(x_i, theta, y_i)
            grad = compute_grad(x_i, e_i, stand_flag=1)
            """依据grad更新theta"""
            theta = update_batch_theta(theta, grad, alpha)
            """计算新的损失值"""
            e_1 = get_e(x, theta, y)
            cost_1 = get_cost(e_1)

            e_reduce_rate = abs(cost_1 - cost_0)/cost_0
            if e_reduce_rate < out_e_reduce_rate:
                flag = 'Finish!\n'
                print flag
                step = out_n + 1
                break
        print 'Step: %-6s, cost: %s' % (step, cost_1[0,0])

    return theta


def iter_mini_batch(x, theta, y, out_n, out_e_reduce_rate, alpha, batch):
    batch_n = x.shape[0]//batch
    batch_left_n = x.shape[0] % batch

    step = 0
    while step < out_n:
        step += 1
        for i in range(batch_n+1):
            """计算初始的损失值"""
            e_0 = get_e(x, theta, y)
            cost_0 = get_cost(e_0)

            """选取更新梯度用的batch个样本"""
            if i <= batch_n-1:
                start = i*batch
                x_i = np.reshape(x[start:start+batch, ], (batch, x.shape[1]))
                y_i = np.reshape(y[start:start+batch, ], (batch, 1))
            else:
                if batch_left_n != 0:
                    x_i = np.reshape(x[-1*batch_left_n:, ], (batch_left_n, x.shape[1]))
                    y_i = np.reshape(y[-1*batch_left_n:, ], (batch_left_n, 1))

            """用batch个样本，计算当前theta组合下的grad值"""
            e_i = get_e(x_i, theta, y_i)
            grad = compute_grad(x_i, e_i, stand_flag=1)
            """依据grad更新theta"""
            theta = update_batch_theta(theta, grad, alpha)
            """计算新的损失值"""
            e_1 = get_e(x, theta, y)
            cost_1 = get_cost(e_1)

            e_reduce_rate = abs(cost_1 - cost_0)/cost_0
            if e_reduce_rate < out_e_reduce_rate:
                flag = 'Finish!\n'
                print flag
                step = out_n
                break
        print 'Step: %-6s, cost: %s' % (step, cost_1[0, 0])

    return theta


def update_newton_theta(x_hessian_inv, theta, grad, alpha):
    theta = theta - alpha*np.dot(x_hessian_inv, grad)
    # print 'New Theta --> ', theta1
    return theta


def iter_newton(x, theta, y, out_n, out_e_reduce_rate, alpha):
    e = get_e(x, theta, y)
    cost_0 = get_cost(e)
    x_hessian_inv = inverse_hessian_mat(x)

    step = 1
    while step < out_n:
        grad = compute_grad(x, e)
        theta = update_newton_theta(x_hessian_inv, theta, grad, alpha)
        e = get_e(x, theta, y)
        cost_1 = get_cost(e)

        print 'Step: %-6s, cost: %s' % (step, cost_1[0,0])
        e_reduce_rate = abs(cost_1 - cost_0)/cost_0
        if e_reduce_rate < out_e_reduce_rate:
            flag = 'Finish!\n'
            print flag
            break
        else:
            cost_0 = cost_1
            step += 1

    return theta


def iter_batch_linesearch(x, theta_0, y, out_n, out_e_reduce_rate, alpha):
    e_0 = get_e(x, theta_0, y)
    cost_0 = get_cost(e_0)

    step = 1
    while step < out_n:
        grad = compute_grad(x, e_0, stand_flag=1)
        theta_1, cost_1, e_1, count = line_search(x, y, alpha, theta_0, grad)
        e_reduce_rate = abs(cost_1 - cost_0)/cost_0
        print 'Step: %-6s, count: %-4s, cost: %s' % (step, count, cost_1[0, 0])
        if e_reduce_rate < out_e_reduce_rate:
            flag = 'Finish!\n'
            print flag
            break
        else:
            cost_0 = cost_1
            theta_0 = theta_1
            e_0 = e_1
            step += 1

    return theta_1


def line_search(x, y, alpha, theta_0, grad):
    # 不更新梯度，在当前梯度方向上，迭代100次寻找最佳的下一个theta组合点，
    e_0 = get_e(x, theta_0, y)
    cost_0 = get_cost(e_0)
    max_iter, count, a, b = 1000, 0, 0.8, 0.5

    while count < max_iter:
        # 随count增大，alpha减小，0.8*0.8*0.8*0.8*...
        alpha_count = pow(a, count) * alpha
        theta_1 = theta_0 - alpha_count*grad
        e_1 = get_e(x, theta_1, y)
        cost_1 = get_cost(e_1)
        # 当前theta组合下的梯度为grad, grad == tan(w) == 切线y变化量／theta变化量， 推出 切线y变化量 = grad * theta变化量
        grad_dy_chage = abs(np.dot(grad.T, theta_1-theta_0))
        # 实际y变化量为cost0-cost1, 不能加绝对值，如果加了就不收敛了。只有cost_y_chage>0且大于b*grad_dy_chage时才符合条件
        cost_y_change = cost_0-cost_1
        # 当前梯度方向上，实际y变化量 占 切线y变化量的比率为b，b越大越好，至少大于0.5才行，实际损失减小量要占 dy的一半以上。
        if cost_y_change > b*grad_dy_chage:
            break
        else:
            count += 1
    return theta_1, cost_1, e_1, count


def iter_bfgs(x, theta_0, y, out_n, out_e_reduce_rate, dfp):
    e_0 = get_e(x, theta_0, y)
    cost_0 = get_cost(e_0)
    grad_0 = np.reshape(compute_grad(x, e_0), (x.shape[1], 1))
    hk = np.eye(x.shape[1])

    step = 1
    while step < out_n:
        dk = -1*np.dot(hk, grad_0)
        theta_1, cost_1, e_1, l_count = armijo_line_search(x, y, theta_0, grad_0, dk, cost_0)
        grad_1 = compute_grad(x, e_1)
        # print theta_1
        # print grad_1
        yk = grad_1 - grad_0
        sk = theta_1 - theta_0

        condition = np.dot(sk.T, yk)
        # 更新以此theta就更新一次hk
        if dfp == 1:
            if condition > 0:
                hk = hk - X2(X3(hk, yk, yk.T), hk)/X3(yk.T, hk, yk)+(X2(sk, sk.T))/condition
        else:
            if condition > 0:
                hk = hk + (1+X3(yk.T, hk, yk)/condition)*(X2(sk, sk.T)/condition)-(X3(sk, yk.T, hk)+X3(hk, yk, sk.T))/condition
        e_reduce_rate = abs(cost_1 - cost_0) / cost_0
        print 'Step: %-6s, l_count: %-6s, cost: %s' % (step, l_count, cost_1[0,0])
        if e_reduce_rate < out_e_reduce_rate:
            flag = 'Finish!\n'
            print flag
            break
        cost_0, grad_0, theta_0 = cost_1, grad_1, theta_1
        step += 1

    return theta_1


def iter_lbfgs(x, theta_0, y, out_n, out_e_reduce_rate):
    limit_n = 5
    ss, yy = [], []

    step = 1
    while step < out_n:

        if step == 1:
            e_0 = get_e(x, theta_0, y)
            cost_0 = get_cost(e_0)
            grad_0 = compute_grad(x, e_0)
            dk = -1*grad_0

        theta_1, cost_1, e_1, l_count = armijo_line_search(x, y, theta_0, grad_0, dk, cost_0)
        grad_1 = compute_grad(x, e_1)

        if len(ss) > limit_n and len(yy) > limit_n:
            del ss[0]
            del yy[0]

        yk = grad_1 - grad_0
        sk = theta_1 - theta_0
        ss.append(sk)
        yy.append(yk)

        qk = grad_1
        k = len(ss)
        condition = X2(yk.T, sk)
        alpha = []

        for i in range(k):
            # t 4->0 倒着计算，倒着存alpha
            t = k-i-1
            pt = 1 / X2(yy[t].T, ss[t])
            alpha.append(pt*X2(ss[t].T, qk))

            qk = qk - alpha[i] * yy[t]

        z = qk

        for i in range(k):
            t = k - i - 1
            # i 0->4 正着计算，正着存beta
            pi = 1 / (X2(yy[i].T, ss[i])[0, 0])
            # pi数
            beta = pi*X2(yy[i].T, z)
            # beta[i]数
            z = z+ss[i]*(alpha[t] - beta)

        if condition > 0:
            dk = -z

        e_reduce_rate = abs(cost_1 - cost_0) / cost_0
        print 'Step: %-6s, l_count: %-6s, cost: %s' % (step, l_count, cost_1[0, 0])
        if e_reduce_rate < out_e_reduce_rate:
            flag = 'Finish!\n'
            print flag
            break
        cost_0, grad_0, theta_0 = cost_1, grad_1, theta_1
        step += 1

    return theta_1


def armijo_line_search(x, y, theta_0, grad_0, dk_0, cost_0):
    # print theta_0.shape, grad_0.shape, dk_0.shape, cost_0.shape
    # 不更新梯度，在当前梯度方向上，迭代100次寻找最佳的下一个theta组合点，
    max_iter, count, countk, a, b = 100, 0, 0, 0.55, 0.4

    while count < max_iter:
        """
        batch方法使用梯度方向grad更新theta，newton等使用牛顿方向dk来更新theta
        newton:hessian逆*grad; dfp:当前步dfpHk; bfgs:当前步bfgsHk
        当前梯度方向上，实际y变化量 占 切线y变化量的比率为b，b越大越好，至少大于0.5才行，实际损失减小量要占 dy的一半以上。
        """

        """更新theta"""
        theta_1 = theta_0 + pow(a, count)*dk_0
        """计算损失"""

        e_1 = get_e(x, theta_1, y)
        cost_1 = get_cost(e_1)
        cost_y_change = cost_1 - cost_0
        dy_change = b * pow(a, count) * np.dot(grad_0.T, dk_0)

        # if cost_1 < cost_0 + b * pow(a, count) * np.dot(grad_0.T, dk_0):
        if cost_y_change < dy_change:
            """如果一直不满足条件，那么一直没有将count赋值给countk，countk仍为0"""
            countk = count
            break
        count += 1

    theta_1 = theta_0 + pow(a, countk) * dk_0
    e_1 = get_e(x, theta_1, y)
    cost_1 = get_cost(e_1)
    return theta_1, cost_1, e_1, count


if __name__ == '__main__':

    x1, y1 = get_data('/optdata/line_sample_data.txt')
    # x1, y1 = get_data('/optdata/airfoil_self_noise.txt')
    theta1 = np.zeros(x1.shape[1]).reshape(x1.shape[1], 1)

    res_theta = iter_batch(x1, theta1, y1, out_n=1e5, out_e_reduce_rate=1e-6, alpha=0.02)
    # res_theta = iter_random_batch(x1, theta1, y1, out_n=1e5, out_e_reduce_rate=1e-6, alpha=0.02)
    # res_theta = iter_mini_batch(x1, theta1, y1, out_n=1e5, out_e_reduce_rate=1e-6, alpha=0.01, batch=10)
    # res_theta = iter_batch_linesearch(x1, theta1, y1, out_n=1e5, out_e_reduce_rate=1e-6, alpha=1)
    # res_theta = iter_newton(x1, theta1, y1, out_n=1e5, out_e_reduce_rate=1e-6, alpha=1)

    # res_theta = iter_bfgs(x1, theta1, y1, out_n=100, out_e_reduce_rate=1e-6, dfp=1)
    # res_theta = iter_bfgs(x1, theta1, y1, out_n=100, out_e_reduce_rate=1e-6, dfp=0)
    # res_theta = iter_lbfgs(x1, theta1, y1, out_n=100, out_e_reduce_rate=1e-6)
    print 'Res_theta:', np.reshape(res_theta, (1, res_theta.shape[0]))[0]

    # def iter_bfgs(x, theta_0, y, out_n, out_e_reduce_rate):

6.2 交叉验证

交叉验证的基本思想是把在某种意义下将原始数据(dataset)进行分组，一部分做为训练集(train set)，另一部分做为验证集(validation set or test set)，首先用训练集对分类器进行训练，再利用验证集来测试训练得到的模型(model)，以此来做为评价分类器的性能指标。代码如下所示：

from sklearn.cross_validation import cross_val_score # K折交叉验证模块
from sklearn.datasets import load_iris 	# iris数据集
from sklearn.model_selection import train_test_split # 分割数据模块
from sklearn.neighbors import KNeighborsClassifier # K最近邻(kNN，k-NearestNeighbor)分类算法
 
#加载iris数据集
iris = load_iris()
X = iris.data
y = iris.target
 
#分割数据并
X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=4)
 
#建立模型
knn = KNeighborsClassifier()
 
#训练模型
knn.fit(X_train, y_train)
 
#将准确率打印出
print(knn.score(X_test, y_test))
# 0.973684210526     基础验证的准确率

#使用K折交叉验证模块
scores = cross_val_score(knn, X, y, cv=5, scoring='accuracy')
 
#将5次的预测准确率打印出
print(scores)

#将5次的预测准确平均率打印出
print(scores.mean())

6.3 数据归一化

如3.3节所示，如果某个特征的方差比其他特征大几个数量级，那么它就会在学习算法中占据主导位置，导致学习器并不能像我们说期望的那样，从其他特征中学习。数据归一化是让不同维度之间的特征在数值上有一定比较性，可以大大提高分类器的准确性。好处如下：

• 提升模型精度
• 提升收敛速度

6.4 回归模型评价指标

参考链接

对于回归模型效果的判断指标经过了几个过程，从SSE到R-square再到Ajusted R-square, 是一个完善的过程：

• SSE(误差平方和)

  计算公式：$SSE=\sum(\hat y-y^*)^2$SSE=∑(y^​−y∗)2，其中的$\hat y$y^​是真实值，$y^*$y∗是预测值。

• R-square(决定系数)

  计算公式：$R^2 = 1-\frac{(\hat y-y^*)^2}{(\hat y-\overline y)^2}$R2=1−(y^​−y​)2(y^​−y∗)2​，其中$\overline y$y​是$y$y的平均值。

  • 数学理解： 分母理解为原始数据的离散程度，分子为预测数据和原始数据的误差，二者相除可以消除原始数据离散程度的影响。

  • 其实“决定系数”是通过数据的变化来表征一个拟合的好坏。

  • 理论上取值范围（-∞，1], 正常取值范围为[0 1] 。实际操作中通常会选择拟合较好的曲线计算R²，因此很少出现$-\infty$−∞.

  • 越接近1，表明方程的变量对y的解释能力越强，这个模型对数据拟合的也较好

    越接近0，表明模型拟合的越差

    经验值：>0.4， 拟合效果好

  • 缺点：数据集的样本越大，R²越大，因此，不同数据集的模型结果比较会有一定的误差。

• Adjusted R-Square (校正决定系数）

  计算公式：$R\_adjusted=1-\frac{(1-R)^2(n-1)}{n-p-1}$R_adjusted=1−n−p−1(1−R)2(n−1)​，其中$n$n为样本数量，$p$p为特征数量。

  • 为了消除了样本数量和特征数量的影响。