基于Jupyter python实现梯度下降和最小二乘法

一、利用梯度下降算法求解多元线性回归方程

梯度下降法 (Gradient Descent)

接下来我们使用梯度下降法来求解上文中的目标函数 J(θ) ​,步骤如下
1.随机初始化 ​θ
2.沿着负梯度方向迭代，使更新后的 θ​ 令 J(θ) ​ 更小，α 为学习率、步长
3.当下降到无法下降或某个定义的极小值时，则停止下降。
注：梯度下降的最终点并非是全局最小点，可能是一个局部最小点

import numpy as np
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
data=np.genfromtxt('D:\\data.csv',delimiter=',')   #导入数据

data.shape

(10, 3)

data

array([[ 10.,  80., 469.],
       [  8.,   0., 366.],
       [  8., 200., 371.],
       [  5., 200., 208.],
       [  7., 300., 246.],
       [  8., 230., 297.],
       [  7.,  40., 363.],
       [  9.,   0., 436.],
       [  6., 330., 198.],
       [  9., 180., 364.]])

x_data=data[:,:-1]
y_data=data[:,2]
#定义学习率、斜率、截据a
#设方程为y=theta1x1+theta2x2+theta0
lr=0.00001
theta0=0
theta1=0
theta2=0
#定义最大迭代次数，因为梯度下降法是在不断迭代更新k与b
epochs=10000
#定义最小二乘法函数-损失函数（代价函数）
def compute_error(theta0,theta1,theta2,x_data,y_data):
    totalerror=0
    for i in range(0,len(x_data)):#定义一共有多少样本点
        totalerror=totalerror+(y_data[i]-(theta1*x_data[i,0]+theta2*x_data[i,1]+theta0))**2
    return totalerror/float(len(x_data))/2
#梯度下降算法求解参数
def gradient_descent_runner(x_data,y_data,theta0,theta1,theta2,lr,epochs):
    m=len(x_data)
    for i in range(epochs):
        theta0_grad=0
        theta1_grad=0
        theta2_grad=0
        for j in range(0,m):
            theta0_grad-=(1/m)*(-(theta1*x_data[j,0]+theta2*x_data[j,1]+theta2)+y_data[j])
            theta1_grad-=(1/m)*x_data[j,0]*(-(theta1*x_data[j,0]+theta2*x_data[j,1]+theta0)+y_data[j])
            theta2_grad-=(1/m)*x_data[j,1]*(-(theta1*x_data[j,0]+theta2*x_data[j,1]+theta0)+y_data[j])
        theta0=theta0-lr*theta0_grad
        theta1=theta1-lr*theta1_grad
        theta2=theta2-lr*theta2_grad
    return theta0,theta1,theta2
#进行迭代求解
theta0,theta1,theta2=gradient_descent_runner(x_data,y_data,theta0,theta1,theta2,lr,epochs)
print('结果：迭代次数：{0} 学习率：{1}之后 a0={2},a1={3},a2={4},代价函数为{5}'.format(epochs,lr,theta0,theta1,theta2,compute_error(theta0,theta1,theta2,x_data,y_data)))
print("多元线性回归方程为:y=",theta1,"X1+",theta2,"X2+",theta0)
#画图
ax=plt.figure().add_subplot(111,projection='3d')
ax.scatter(x_data[:,0],x_data[:,1],y_data,c='r',marker='o')
x0=x_data[:,0]
x1=x_data[:,1]
#生成网格矩阵
x0,x1=np.meshgrid(x0,x1)
z=theta0+theta1*x0+theta2*x1
#画3d图
ax.plot_surface(x0,x1,z)
ax.set_xlabel('area')
ax.set_ylabel('distance')
ax.set_zlabel("monney")
plt.show()

结果：迭代次数：10000 学习率：1e-05之后 a0=5.3774162274868,a1=45.0533119768975,a2=-0.19626929358281256,代价函数为366.7314528822914
多元线性回归方程为:y= 45.0533119768975 X1+ -0.19626929358281256 X2+ 5.3774162274868

二、最小二乘法法求解多元线性回归方程的python代码实现

最小二乘法概念：

第一，线性回归的模型假设。这是最小二乘方法的优越性前提，否则不能推出最小二乘是最佳（即方差最小）的无偏估计，具体请参考高斯-马尔科夫定理。特别地，当随机噪声服从正态分布时，最小二乘与最大似然等价。
第二，由于是线性回归/拟合，因此可以很容易的求出全局最优的闭式解（close form solution）即：全局最优解，也即我们通常看到的那几个矩阵形式，给了（input data）可以一步到位计算出拟合参数，而不是像梯度下降法或者牛顿法那样一点点地迭代优化调参，最后到达极值点，即：局部最优解，二者得到的参数值可能会相同或特别接近。

#利用线性代数的矩阵模拟最小二乘法求解法求解多元线性回归方程的系数
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import pandas as pd
%matplotlib inline
data = np.genfromtxt("D:/data.csv",delimiter=",")
X1=data[0:10,0]#自变量温度
X2=data[0:10,1]#因变量销售量
Y=data[0:10,2]#自变量温度
#将因变量赋值给矩阵Y1
Y1=np.array([Y]).T
#为自变量系数矩阵X赋值
X11=np.array([X1]).T
X22=np.array([X2]).T
A=np.array([[1],[1],[1],[1],[1],[1],[1],[1],[1],[1]])#创建系数矩阵
B=np.hstack((A,X11))#将矩阵a与矩阵X11合并为矩阵b
X=np.hstack((B,X22))#将矩阵b与矩阵X22合并为矩阵X
#求矩阵X的转置矩阵
X_=X.T
#求矩阵X与他的转置矩阵的X_的乘积
X_X=np.dot(X_,X)
#求矩阵X与他的转置矩阵的X_的乘积的逆矩阵
X_X_=np.linalg.inv(X_X)
#求解系数矩阵W，分别对应截距b、a1、和a2
W=np.dot(np.dot((X_X_),(X_)),Y1)
b=W[0][0]
a1=W[1][0]
a2=W[2][0]
print("系数a1=",a1)
print("系数a2=",a2)
print("截距为=",b)
print("多元线性回归方程为:y=",a1,"X1+",a2,"X2+",b)
#求月销售量Y的和以及平均值y1
sumy=0#因变量的和
y1=0#因变量的平均值
for i in range(0,len(Y)):
    sumy=sumy+Y[i]
y1=sumy/len(Y)
#求月销售额y-他的平均值的和
y_y1=0#y-y1的值的和
for i in range(0,len(Y)):
    y_y1=y_y1+(Y[i]-y1)
print("销售量-销售量平均值的和为:",y_y1)
#求预测值sales1
sales1=[]
for i in range(0,len(Y)):
    sales1.append(a1*X1[i]+a2*X2[i]+b)
#求预测值的平均值y2
y2=0
sumy2=0
for i in range(len(sales1)):
    sumy2=sumy2+sales1[i]
y2=sumy2/len(sales1)
#求预测值-平均值的和y11_y2
y11_y2=0
for i in range(0,len(sales1)):
   y11_y2=y11_y2+(sales1[i]-y2)
print("预测销售值-预测销售平均值的和为:",y11_y2)
#求月销售额y-他的平均值的平方和
Syy=0#y-y1的值的平方和
for i in range(0,len(Y)):
    Syy=Syy+((Y[i]-y1)*(Y[i]-y1))
print("Syy=",Syy)
#求y1-y1平均的平方和
Sy1y1=0
for i in range(0,len(sales1)):
    Sy1y1=Sy1y1+((sales1[i]-y2)*(sales1[i]-y2))
print("Sy1y1=",Sy1y1)
#（y1-y1平均）*（y-y平均）
Syy1=0
for i in range(0,len(sales1)):
    Syy1=Syy1+((Y[i]-y1)*(sales1[i]-y2))
print("Syy1=",Syy1)
#求R
R=Syy1/((Syy*Sy1y1)**0.5)
R2=R*R
print("判定系数R2=",R2)

系数a1= 41.513478256438404
系数a2= -0.34088268566361957
截距为= 65.32391638894777
多元线性回归方程为:y= 41.513478256438404 X1+ -0.34088268566361957 X2+ 65.32391638894777
销售量-销售量平均值的和为: -1.1368683772161603e-13
预测销售值-预测销售平均值的和为: 2.842170943040401e-14
Syy= 76199.6
Sy1y1= 72026.59388000517
Syy1= 72026.59388000521
判定系数R2= 0.9452358526817103

三、分析

  两个方法求出的多元线性回归方程还是有点误差的。最小二乘法计算一个矩阵的逆是相当耗费时间的, 而且求逆也会存在数值不稳定的情况 (比如对希尔伯特矩阵求逆就几乎是不可能的). 因而这样的计算方法有时不值得提倡.相比之下, 梯度下降法虽然有一些弊端, 迭代的次数可能也比较高, 但是相对来说计算量并不是特别大. 而且, 在最小二乘法这个问题上, 收敛性有保证. 故在大数据量的时候, 反而是梯度下降法 (其实应该是其他一些更好的迭代方法) 更加值得被使用。
 当然, 其实梯度下降法还有别的其他用处, 比如其他找极值问题. 另外, 牛顿法也是一种不错的方法, 迭代收敛速度快于梯度下降法, 只是计算代价也比较高.