利用梯度下降法解决回归问题

解决回归问题时我们最先想到的是最小二乘法，但是最小二乘法有一定的弊端：根据损失函数求得其最小值时对应w和b对于小样本没有问题，但是数据量大耗时甚至不可行。由此引出我们接下来要用到的方法：梯度下降法。

1.梯度下降法基本原理

对于一组数据：$$X=(x_1,x_2,\dots ,x_i,\dots ,x_m),Y=(y_1,y_2,\dots ,y_i,\dots ,y_m)$$假设最终找出的预测函数为$y=h(x)$，则我们可以通过如下函数来判断此预测函数是否合适：$$J=\frac{1}{2m}\sum _{i=1}^m(h(x_i)-y_i{)}^2$$一般称其为代价函数。我们需要做的就是找到函数$h(x)$中的参数使得代价函数尽量小。
梯度下降法相当于“无固定路线下坡”。假设我们处于山坡上某点，想要到达谷底。每次我们要走一步下山，那么最可靠的方法便是环顾四周，寻找能一步到达的最低点，持续该过程，最后得到快速下坡路线。
对于函数而言，便是求得该函数对所有参数（变量）的偏导，也就是求梯度。梯度方向是函数上升最快的方向，那么梯度的反方向便是下降最快的方向。
每次更新这些参数，直到到达最低点为止，注意这些参数必须在每一轮一起更新。我们通过下面线性回归的例子来进一步体会。

2.线性回归的例子

我们拿如下样本作为例子：
X = [1,2,3,4,5,6,7,8,9,10]
Y = [1,3,4,5.4,4.8,6.3,8.4,10.1,10.9,11.4]
这组数据的散点图：

首先定义学习率，学习率可以理解为步长，就是下山途中每一步的距离，不能过大也不能过小。过大会导致一步迈过最低点，而过小会导致到达最低点的时间过长，降低效率。

#定义学习率：
alpha = 0.05

对于线性回归，我们的预测函数定义为一次函数即可：

#定义预测函数：
def h(a,x,b):
    return a*x+b

定义代价函数：

#定义代价函数：
def J(a,b,X,Y,m):
    sum = np.dot(h(a,X,b)-Y,h(a,X,b)-Y)
    return sum/(2*m)

接下来是求梯度。此时的代价函数为$J=\frac{1}{2m}{\sum }_{i=1}^m(a{x}_i+b-y_i{)}^2$,则$$\frac{\partial J}{\partial a}=\frac{1}{m}\sum _{i=1}^m(a{x}_i+b-y_i)x_i,\frac{\partial J}{\partial b}=\frac{1}{m}\sum _{i=1}^m(a{x}_i+b-y_i)$$

#定义梯度函数：
def J1(a,b,X,Y,m):
    sum1 = 0
    for j in range(0,m):
        sum1 = sum1+h(a,X[j],b)-Y[j]
    return sum1/m

def J2(a,b,X,Y,m):
    sum2 = np.dot(h(a,X,b)-Y,X)
    return sum2/m

接下来初始化相关参数：

#初始化a,b,m:
a = Y[0]/X[0]
b = 0
m = 10

随后便可进行迭代求解：

#计算代价函数,进行迭代：
while np.abs(J1(a,b,X,Y,m))>=1e-5 or np.abs(J2(a,b,X,Y,m))>=1e-5:
    a = a-alpha*J2(a,b,X,Y,m)
    b = b-alpha*J1(a,b,X,Y,m)

当梯度的两部分均小于1e-5时，可以知道已经进入基本水平阶段，即认为找到最低点，停止循环。
最终输出结果如图：

a,b以及代价函数的值为：1.1509,0.2000,0.1771
完整代码如下：

import matplotlib.pyplot as plt
import numpy as np

#定义学习率：
alpha = 0.05

#定义预测函数：
def h(a,x,b):
    return a*x+b

#定义代价函数：
def J(a,b,X,Y,m):
    sum = np.dot(h(a,X,b)-Y,h(a,X,b)-Y)
    return sum/(2*m)

#定义梯度函数：
def J1(a,b,X,Y,m):
    sum1 = 0
    for j in range(0,m):
        sum1 = sum1+h(a,X[j],b)-Y[j]
    return sum1/m

def J2(a,b,X,Y,m):
    sum2 = np.dot(h(a,X,b)-Y,X)
    return sum2/m

#给出数据集：
X = np.array([1,2,3,4,5,6,7,8,9,10])
Y = np.array([1,3,4,5.4,4.8,6.3,8.4,10.1,10.9,11.4])
plt.scatter(X,Y)#绘制散点图

#初始化a,b,m:
a = Y[0]/X[0]
b = 0
m = 10
#计算代价函数,进行迭代：
while np.abs(J1(a,b,X,Y,m))>=1e-5 or np.abs(J2(a,b,X,Y,m))>=1e-5:
    a = a-alpha*J2(a,b,X,Y,m)
    b = b-alpha*J1(a,b,X,Y,m)

print(a,b,J(a,b,X,Y,m))
x = np.arange(0.5,10.5,0.01)
y = a*x+b
plt.plot(x,y,c='red')
plt.legend(['predict line'])
plt.show()

3.非线性回归的例子

仿照以上例子，我们可以进行对于非线性数据的预测。
例如对于数据：
X = [1,2,3,4]
Y = [3.2,6.5,14.5,20.6]
散点图如图：

此时若选择二次函数为预测函数：$h(x)=a{x}^2+bx+c$,其有三个参数a,b,c，则在写梯度函数时需要三个函数，代码如下：

import matplotlib.pyplot as plt
import numpy as np

#定义学习率：
alpha = 0.01

#定义预测函数：
def h(a,b,c,x):
	return a*x*x+b*x+c

#定义代价函数：
def J(a,b,c,X,Y,m):
	sum = np.dot(h(a,b,c,X)-Y,h(a,b,c,X)-Y)
	return sum/(2*m)

#定义梯度函数：
def J_a(a,b,c,X,Y,m):
	sum1 = np.dot(h(a,b,c,X)-Y,X**2)
	return sum1/m

def J_b(a,b,c,X,Y,m):
	sum2 = np.dot(h(a,b,c,X)-Y,X)
	return sum2/m

def J_c(a,b,c,X,Y,m):
	sum3 = 0
	for t in range(0,m):
		sum3 = sum3+h(a,b,c,X[t])-Y[t]
	return sum3/m

#给出数据集：
X = np.array([1,2,3,4])
Y = np.array([3.2,6.5,14.5,20.6])
plt.scatter(X,Y)#绘制散点图

#初始化a,b,m:
a = 1
b = 0
c = 0
m = 4
#计算代价函数,进行迭代：
while np.abs(J_a(a,b,c,X,Y,m))>=1e-5 or np.abs(J_b(a,b,c,X,Y,m))>=1e-5 or np.abs(J_c(a,b,c,X,Y,m))>=1e-5:
	a = a-alpha*J_a(a,b,c,X,Y,m)				    
	b = b-alpha*J_b(a,b,c,X,Y,m)
	c = c-alpha*J_c(a,b,c,X,Y,m)
																				 

print(a,b,c,J(a,b,c,X,Y,m))
x = np.arange(0.5,4.5,0.01)
y = a*x*x+b*x+c
plt.plot(x,y,c='red')
plt.legend(['predicted'])
plt.show()

结果如图：

求得的a,b,c及代价函数值为：0.7001,2.5195，-0.3494,0.2723