matlab 最小二乘法拟合_基于梯度下降与最小二乘法解线性回归模型—机器学习(6)...

1. 线性回归求解

在前一篇中，了解线性回归原理及模型，并调用开源工具了线性回归的目标函数的求解，并不知道内部的真正原理。接下来祥细对目标函数的求解过程和最小二乘法进行讲解。

1.1 使用梯度下降来确定线性回归的参数

最优化问题在机器学习中有非常重要的地位，很多机器学习算法最后都归结为求解最优化问题。在各种最优化算法中，梯度下降法是最简单、最常见的一种，在深度学习的训练中被广为使用。

1.1.1最优化问题

最优化问题是求解函数极值的问题，包括极大值和极小值。相信对这个问题都不陌生，在初中时就学会了求解二次函数的极值(抛物线的顶点)，高中时学习了幂函数，指数函数，对数函数，三角函数，反三角函数等各种类型的函数，求函数极值的题更是频频出现。这些方法都采用了各种各样的技巧，没有一个统一的方案。

真正的飞跃发生在大学时，微积分为我们求函数的极值提供了一个统一的思路：找函数的导数等于0的点，因为在极值点处，导数必定为0。这样，只要函数的可导的，我们就可以用这个万能的方法解决问题，幸运的是，在实际应用中我们遇到的函数基本上都是可导的。

在机器学习之类的实际应用中，我们一般将最优化问题统一表述为求解函数的极小值问题，即：

其中称为优化变量，称为目标函数。极大值问题可以转换成极小值问题来求解，只需要将目标函数加上负号即可:

一个优化问题的全局极小值   是指对于可行域里所有的，有：

即全局极小值点处的函数值不大于任意一点处的函数值。局部极小值   定义为存在一个     邻域，对于在邻域内：

并且在可行域内的所有，有：

即局部极小值点处的函数值比一个局部返回内所有点的函数值都小。在这里，我们的目标是找到全局极小值。不幸的是，有些函数可能有多个局部极小值点，因此即使找到了导数等于0的所有点，还需要比较这些点处的函数值。

1.1.2 梯度下降思想

对于梯度下降，可以形象地理解为一个人下山的过程。假设现在有一个人在山上，现在他想要走下山，但是他不知道山底在哪个方向，怎么办呢？显然可以想到的是，一定要沿着山高度下降的地方走，不然就不是下山而是上山了。山高度下降的方向有很多，选哪个方向呢？这个人比较有冒险精神，他选择最陡峭的方向，即山高度下降最快的方向。现在确定了方向，就要开始下山了。又有一个问题来了，在下山的过程中，最开始选定的方向并不总是高度下降最快的地方。假设这个人比较聪明，他每次都选定一段距离，每走一段距离之后，就重新确定当前所在位置的高度下降最快的地方。这样，这个人每次下山的方向都可以近似看作是每个距离段内高度下降最快的地方。现在将这个思想引入线性回归，在线性回归中，要找到参数矩阵，使得损失函数最小。如果把损失函数看作是这座山，山底不就是损失函数最小的地方吗，那求解参数矩阵的过程，就是人走到山底的过程 下山的例子中每一段路的距离取名叫学习率(Learning Rate，也称步长表示)，把一次下山走一段距离叫做一次迭代。算法详细过程：

1.确定定参数的初始值，计算损失函数的偏导数
2.将参数代入偏导数计算出梯度。若梯度为 0，结束；否则转到 3
3.用步长乘以梯度，并对参数进行更新
4.重复 2-3

1.1.3 导数与梯度

由于实际应用中一般都是多元函数，因此我们跳过一元函数，直接介绍多元函数的情况。梯度是导数对多元函数的推广，它是多元函数对各个自变量偏导数形成的向量。多元函数的梯度定义为：

其中称为梯度算子，它作用于一个多元函数，得到一个向量

可导函数在某一点处取得极值的必要条件是梯度为0，梯度为0的点称为函数的驻点，这是疑似极值点。需要注意的是，梯度为0只是函数取极值的必要条件而不是充分条件，即梯度为0的点可能不是极值点。

至于是极大值还是极小值，要看二阶导数Hessian矩阵，Hessian矩阵我们将在后面的文章中介绍，这是由函数的二阶偏导数构成的矩阵。这分为下面几种情况：如果Hessian矩阵正定，函数有极小值 如果Hessian矩阵负定，函数有极大值 如果Hessian矩阵不定，则需要进一步讨论

这和一元函数的结果类似，Hessian矩阵可以看做是一元函数的二阶导数对多元函数的推广。一元函数的极值判别法为，假设在某点处导数等于0，则：如果二阶导数大于0，函数有极小值 如果二阶导数小于0，函数有极大值 如果二阶导数等于0，情况不定

在这里我们可能会问：直接求函数的导数/梯度，然后令导数/梯度为0，解方程，问题不就解决了吗？事实上没这么简单，因为多元方程的导数可能很难解。如函数:

分别对x和y求偏导数，并令它们为0，得到下面的方程组：

这个方程非常难以求解.精确的求解不太可能，因此只能求近似解，这称为数值计算。工程上实现时通常采用的是迭代法，它从一个初始点x_{0}开始，反复使用某种规则从x_{k}   移动到下一个点x_{k+1}  ，构造这样一个数列，直到收敛到梯度为0的点处。即有下面的极限成立：

这些规则一般会利用一阶导数信息即梯度；或者二阶导数信息即Hessian矩阵。这样迭代法的核心是得到这样的由上一个点确定下一个点的迭代公式：

只要没有到达梯度为0的点，则函数值会沿着序列递减，最终会收敛到梯度为0的点，这就是梯度下降法。迭代终止的条件是函数的梯度值为0(实际实现时是接近于0)，此时认为已经达到极值点。注意我们找到的是梯度为0的点，这不一定就是极值点。梯度下降法只需要计算函数在某些点处的梯度，实现简单，计算量小。

1.1.4 实现细节问题

• 初始值的设定:对于不带约束条件的优化问题，我们可以将初始值设置为0，或者设置为随机数，对于神经网络的训练，一般设置为随机数，这对算法的收敛至关重要;
• 学习率的设定:学习率设置为多少，也是实现时需要考虑的问题。最简单的，我们可以将学习率设置为一个很小的正数，如0.001。另外，可以采用更复杂的策略，在迭代的过程中动态的调整学习率的值。比如前1万次迭代为0.001，接下来1万次迭代时设置为0.0001。

1.2 使用最小二乘法来确定线性回归的参数

最小二乘法(又称最小平方法)是一种数学优化技术。它通过最小化误差的平方和寻找数据的最佳函数匹配。

1.2.1 最小二乘法的原理

公式形式：

目标函数观测值理论值

我们的目标是得到使目标函数最小化时候的拟合函数的模型。如有m个只有一个特征的样本：

样本采用下面的拟合函数：θθθ样本有一个特征x，对应的拟合函数有两个参数θ和θ需要求出。目标函数为：

当使θθ最小时，求出使θθ最小时的θ和θ，这样拟合函数就得出了。那么，最小二乘法怎么才能使θθ最小呢？

1.2.2 最小二乘法的原理矩阵求解法

最小二乘法的解法有代数法、矩阵法，一般计算机使用矩阵法。假设函数的矩阵表达方式为：

其中，假设函数为的向量,为的向量，里面有个代数法的模型参数。X为维的矩阵。代表样本的个数，代表样本的特征数。损失函数定义为:

其中Y是样本的输出向量，维度为. 在这主要是为了求导后系数为1，方便计算. 根据最小二乘法的原理，我们要对这个损失函数对向量求导取0:

根据这个公式求出 ：

只要给了数据,就可以这个这个公式求出参数向量：

1.2.3 最小二乘法的局限性和适用场景

• 首先，最小二乘法需要计算的逆矩阵，有可能它的逆矩阵不存在，这样就没有办法直接用最小二乘法了，此时梯度下降法仍然可以使用;
• 当样本特征n非常的大的时候，计算的逆矩阵是一个非常耗时的工作(的矩阵求逆)，甚至不可行。此时以梯度下降为代表的迭代法仍然可以使用;
• 如果拟合函数不是线性的，这时无法使用最小二乘法，需要通过一些技巧转化为线性才能使用，此时梯度下降仍然可以用;
• 当样本量m很少，小于特征数n的时候，这时拟合方程是欠定的，常用的优化方法都无法去拟合数据。当样本量m等于特征数n的时候，用方程组求解就可以了。当m大于n时，拟合方程是超定的，也就是我们常用与最小二乘法的场景了。

2 线性回归的梯度下降算法使用

通过案例说明使用梯度下降来确定线性回归的参数

2.1 创建一个遵循分布的线性回归函数和模拟数据

import numpy as npimport matplotlib.pyplot as pltfrom matplotlib import cmfrom mpl_toolkits.mplot3d import Axes3Dplt.style.use('ggplot')np.random.seed(1)X = np.random.randint(0, 100, 100)# Target variable# y = 2x - 50 + ϵ (noise)y = 2*X - 50 + np.random.uniform(-15, 15, 100)plt.scatter(X, y)plt.xlabel("x")plt.ylabel("y")plt.show()

2.2目标是最小化平方误差的平均值MSE，其计算公式如下

• 开始尝试很多值 和 个，并查看哪一个MSE值最低， 尝试一下介于-100和100之间的每个整数值 和 个，并绘制结果图：

plt_beta_0, plt_beta_1 = np.meshgrid(beta_0, beta_1)def calculate_mse(beta_0, beta_1):    y_hat = beta_0 + beta_1 * X    error = y_hat - y    sse = np.sum(error ** 2)    return ((1 / (2 * len(X))) * sse)calculate_mse_v = np.vectorize(calculate_mse)mse = calculate_mse_v(plt_beta_0, plt_beta_1).reshape(len(plt_beta_0), len(plt_beta_1))fig = plt.figure(figsize=(10, 5))ax = fig.gca(projection='3d')surf = ax.plot_surface(plt_beta_0,                       plt_beta_1,                       mse,                       cmap=cm.jet,                       linewidth=0,                       antialiased=False)ax.set_xlabel(r'$\hat{\beta_0}$')ax.set_ylabel(r'$\hat{\beta_1}$')ax.set_zlabel(r'MSE')fig.colorbarsurf

对y的影响明显大于

• 获取值的索引 和 最低的mse并绘制结果图

beta_1_idx, beta_0_idx = np.unravel_index(mse.argmin(), mse.shape)beta_0_hat = beta_0[beta_0_idx]beta_1_hat = beta_1[beta_1_idx]print("y = {} + {}x".format(beta_0_hat, beta_1_hat))# Plot a line for our modelplt.scatter(X, y)plt.plot(    [min(X), max(X)],    [min(X) * beta_1_hat + beta_0_hat, max(X) * beta_1_hat + beta_0_hat],    color='blue')plt.xlabel('x')plt.ylabel('y')plt.show()

• 使用梯度下降求解 梯度下降的目的是根据成本函数给出的凹形向下移动梯度以达到最小值，如上面的3D图形所示。梯度下降是一个迭代过程：

1.初始化 和 个 具有随机值并计算MSE；

2.计算梯度，以便我们朝着最小化MSE的方向发展；

3.调整 和个 带有渐变；

4.使用新的权重获取值 ÿ^ 计算MSE；

5.重复步骤2-4。

使用 α确定应该将梯度向下移动多远的术语。一个α太小会花费更长的时间才能收敛到最小值，而 α 太大可能会开始偏离MSE的最小值。使用以下初始化数据：

1. α = 0.0005
2. 100,000次迭代。
3. 初始化 值为1
4. 初始化 值为1

X = np.array([np.ones(len(X)), X])alpha = 0.0005iterations = 100000beta_0_hat = 1beta_1_hat = 1mses = []for i in range(1, iterations+1):    y_hat = beta_0_hat * X[0] + beta_1_hat * X[1]    error = y_hat - y    sse = np.sum(error ** 2)    mse = ((1 / (2 * len(X.T))) * sse)    mses.append(mse)    gradient = np.dot(X, error) / len(X.T)    beta_0_hat = beta_0_hat - (gradient[0] * alpha)    beta_1_hat = beta_1_hat - (gradient[1] * alpha)    if i % 10000 == 0:        print("Iteration {}, MSE={}, β0={}, β1={}".format(            i, round(mse, 3), round(beta_0_hat, 3), round(beta_1_hat, 3)))plt.scatter(X[1], y)plt.plot(    [min(X[1]), max(X[1])],    [min(X[1]) * beta_1_hat + beta_0_hat, max(X[1]) * beta_1_hat + beta_0_hat],    color='blue')plt.xlabel('x')plt.ylabel('y')plt.show()

迭代10000，MSE = 53.152，β0= -38.628，β1= 1.791 迭代20000，MSE = 34.885，β0= -47.416，β1= 1.929 迭代30000，MSE = 33.986，β0= -49.364，β1= 1.959 迭代40000，MSE = 33.942，β0= -49.796，β1= 1.966 迭代50000，MSE = 33.94，β0= -49.892，β1= 1.967 迭代60000，MSE = 33.94，β0= -49.913，β1= 1.968 迭代70000，MSE = 33.94，β0=- 49.918，β1= 1.968 迭代80000，MSE = 33.94，β0= -49.919，β1= 1.968 迭代90000，MSE = 33.94，β0= -49.919，β1= 1.968 迭代100000，MSE = 33.94，β0= -49.919，β1= 1.968

最后 值是-49.919，而我们的最终值  值是1.968

• 在迭代梯度下降步骤时收集MSE值时，可以绘制图形并查看MSE如何随时间衰减：

plt.figure(figsize=(10,5))plt.plot(np.arange(1, iterations+1), mses, color='green')plt.ylabel("MSE")plt.xlabel("Iteration")plt.show()

可以看到，大部分衰减发生在前20,000次迭代中

3 小结

今天总结了机器学习的线性回归目标函数最小化的两种算法和实例:

• 利用梯度下降方法寻找目标函数最小化；
• 利用最小二乘法方法寻找目标函数最小化。