利用Python编程，分别使用梯度下降法和最小二乘法求解多元函数

分别使用梯度下降法和最小二乘法求解多元函数并进行比较，这里使用jupyter notebook平台进行Python编程

一、题目描述

为求得某个地区的商品店的月营业额与店铺的面积、该店距离车站距离的相关性大小，以店铺面积、距离车站的距离、以及月营业额建立线性回归方程，并求解方程，得到相关系数。

将表中的数据录入到Excel中，编程时会用到这些数据，如下：

二、使用梯度下降法求解多元函数

（一）梯度下降法基本原理

梯度下降法又称最速下降法，是求解无约束最优化问题的一种最常用的方法，在对损失函数最小化时经常使用。梯度下降法是一种迭代算法。选取适当的初值x(0)，不断迭代，更新x的值，进行目标函数的极小化，直到收敛。由于负梯度方向时使函数值下降最快的方向，在迭代的每一步，以负梯度方向更新x的值，从而达到减少函数值的目的。

（二）梯度下降法公式

（三）Python代码实现

1、导入要使用到的库、定义变量并赋值

import numpy as np
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
data=np.genfromtxt('D:/area_distance.csv',delimiter=',')
x_data=data[:,:-1]
y_data=data[:,2]
#定义学习率、斜率、截据
#设方程为y=theta1x1+theta2x2+theta0
lr=0.00001
theta0=0
theta1=0
theta2=0
#定义最大迭代次数，因为梯度下降法是在不断迭代更新k与b
epochs=10000

2、代价函数

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

3、使用梯度下降法求解多元函数系数

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

4、打印系数和方程

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)

5、绘制线性回归方程拟合图

#画图
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("Monthly turnover")
plt.show()

6、完整代码

import numpy as np
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
data=np.genfromtxt('D:/area_distance.csv',delimiter=',')
x_data=data[:,:-1]
y_data=data[:,2]
#定义学习率、斜率、截据
#设方程为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("Monthly turnover")
plt.show()

7、运行结果

三、使用最小二乘法求解多元函数

（一）最小二乘法基本原理

在我们研究两个变量（x,y）之间的相互关系时，通常可以得到一系列成对的数据（x1,y1.x2,y2… xm,ym）；将这些数据描绘在x -y直角坐标系中，若发现这些点在一条直线附近，可以令这条直线方程为（式1-1）
其中：a0、a1 是任意实数
为建立这直线方程就要确定a0和a1，应用《最小二乘法原理》，将实测值Yi与利用计算值Yj（Yj=a0+a1Xi）（式1-1）的离差（Yi-Yj）的平方和 最小为“优化判据”。
令：φ = （式1-2)
把（式1-1）代入（式1-2）中得：
φ = （式1-3)
当 最小时，可用函数 φ 对a0、a1求偏导数，令这两个偏导数等于零。
∑2(a0 + a1Xi - Yi）=0（式1-4)
∑2Xi（a0 +a1Xi - Yi）=0（式1-5)
亦即：
na0 + （∑Xi ) a1 = ∑Yi （式1-6)
（∑Xi ) a0 + （∑Xi^2 ) a1 = ∑（Xi*Yi) （式1-7)
得到的两个关于a0、 a1为未知数的两个方程组，解这两个方程组得出：
a0 = （∑Yi) / n - a1（∑Xi) / n （式1-8)
a1 = [n∑(Xi Yi) - （∑Xi ∑Yi)] / (n∑Xi^2 -∑Xi∑Xi)（式1-9)
这时把a0、a1代入（式1-1）中， 此时的(式1-1）就是我们回归的一元线性方程即：数学模型。
在回归过程中，回归的关联式不可能全部通过每个回归数据点（x1,y1. x2,y2…xm,ym），为了判断关联式的好坏，可借助相关系数“R”，统计量“F”，剩余标准偏差“S”进行判断；“R”越趋近于 1 越好；“F”的绝对值越大越好；“S”越趋近于 0 越好。
R = [∑XiYi - m （∑Xi / m）（∑Yi / m)]/ SQR{[∑Xi2 - m （∑Xi / m)2][∑Yi2 - m （∑Yi / m)2]} （式1-10) *
在（式1-10）中，m为样本容量，即实验次数；Xi、Yi分别为任意一组实验数据X、Y的数值。

（二）Python代码实现

1、自行推导过程编写代码

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import seaborn as sns

%matplotlib inline
data = np.genfromtxt("D:/area_distance.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={:.1f}".format(a1))
print("系数a2={:.1f}".format(a2))
print("截距为={:.1f}".format(b))
print("多元线性回归方程为:y={:.1f}".format(a1),"X1+ {:.1f}".format(a2),"X2+{:.1f}".format(b))
#画出线性回归分析图
data1=pd.read_excel('D:\面积距离车站数据.xlsx')
sns.pairplot(data1, x_vars=['area','distance'], y_vars='Y', height=3, aspect=0.8, kind='reg')  
plt.show() 
#求月销售量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("营业额-营业额平均值的和为:{:.1f}".format(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("预测营业额-预测营业额平均值的和为:{:.1f}".format(y11_y2))
#求月销售额y-他的平均值的平方和
Syy=0#y-y1的值的平方和
for i in range(0,len(Y)):
    Syy=Syy+((Y[i]-y1)*(Y[i]-y1))
print("Syy={:.1f}".format(Syy))
#求y1-y1平均的平方和
Sy1y1=0
for i in range(0,len(sales1)):
    Sy1y1=Sy1y1+((sales1[i]-y2)*(sales1[i]-y2))
print("Sy1y1={:.1f}".format(Sy1y1))
#（y1-y1平均）*（y-y平均）
Syy1=0
for i in range(0,len(sales1)):
    Syy1=Syy1+((Y[i]-y1)*(sales1[i]-y2))
print("Syy1={:.1f}".format(Syy1))
#求y-y1的平方Se
Se=0
for i in range(0,len(sales1)):
    Se=Se+((Y[i]-y2)*(Y[i]-y2))
print("Se={:.1f}".format(Se))
#求R
R=Syy1/((Syy*Sy1y1)**0.5)
R2=R*R
print("R2={:.4f}".format(R2))

运行结果

2、采用sklearn库

import pandas as pd
import seaborn as sns  
import matplotlib.pyplot as plt
from sklearn import linear_model
%matplotlib inline

data=pd.read_excel('D:\面积距离车站数据.xlsx')      #导入数据
X=data[['area','distance']]           #x,y取值
y=data['Y']    
y=data.Y    
model = linear_model.LinearRegression()
model.fit(X,y)
print("系数a1={:.1f}".format(model.coef_[0]))
print("系数a2={:.1f}".format(model.coef_[1]))
print("截距为={:.1f}".format(model.intercept_))
print("多元线性回归方程为:y={:.1f}".format(model.coef_[0]),"X1+ {:.1f}".format(model.coef_[1]),"X2+{:.1f}".format(model.intercept_))
sns.pairplot(data, x_vars=['area','distance'], y_vars='Y', height=3, aspect=0.8, kind='reg')  
plt.show()  

运行结果

（四）将梯度下降法与最小二乘法进行对比

这里给出利用Excel表格的数据分析功能得出的数据：




通过对两种方法的对比分析，不难得出结论：最小二乘法计算量偏大，而且求逆矩阵相当耗费时间，求逆矩阵也会存在数值不稳定的情况，相比之下梯度下降法可以看做是一种更简单的最小二乘法最后一步解方程的方法，梯度下降法虽然有一定的缺点，但是计算量不大，在数据量比较大的时候选择梯度下降法比较好一点。