机器学习梯度下降法举例

本文是参考Andrew Ng网上课程和别人总结的机器学习笔记，作为笔者学习的参考，笔者主要实现了一些程序问题。

持续更新中…

更新…
先以最为简单的例子举例，比如有一个函数为$$y=x^2$$现在要求当x等于多少的时候y最小。
采用梯度下降法，即往梯度下降的方向走,已知梯度为$$\frac{dy}{dx}=2x$$所以梯度下降的方向即$-2x$,所以参数$x$的更迭方式为$$x=x+(-2x)=x-2x$$但通常我们加入一个表示步长的参数$\eta$,所以更迭变为$$x=x-2\eta x$$
于是代码如下：

import numpy as np
import time

s = time.time()

x = 10000*np.random.randn(1)  #产生服从正太分布的一个数(均值=0，方差=1)，扩大10000倍
eta = 0.001
i = 1
while True:
    y = x**2
    x -= eta*2*x
    if y<=0.0001:
        print('x is %.6f\ny is %.6f\nand steps are %d'%(x,y,i))
        break
    i += 1
e = time.time()
print('time is %.4f second'%(e-s))

对于两个参数，也是一样的求法，如$$z=x^2+y+xy$$
参数更迭：$$x=x-\eta \frac{dz}{dx}=x-\eta (2x+y);y=y-\eta \frac{dz}{dy}=y-\eta (1+x)$$

import numpy as np
import time

s = time.time()

x,y = 10*np.random.randn(2)  #产生服从正太分布的一个数(均值=0，方差=1)，扩大10000倍
eta = 0.001
i = 1

while True:
    z = x**2+y+x*y
    x -= eta*(2*x+y)
    y -= eta*(1+x)
    if z<=-10:
        print('x is %.6f\ny is %.6f\nz is %.6f\nand steps are %d'%(x,y,z,i))
        break
    i += 1
e = time.time()
print('time is %.4f second'%(e-s))

1 两个参数梯度下降法

现在有一组数据，plot之后如下图

x = np.array([1,3,4,7,9,10,11,18])
y = np.array([5,12,12,20,32,30,33,60])

现在问：当x等于8的时候，y的取值是多少？

首先看这个图，具有线性，于是我们假设先用一条直线拟合（即一次函数）
采用梯度下降法
我们的hypothesis为（两个参数${\theta }_0和{\theta }_1$）

loss function或说cost function或说目标函数（objective function）为
$$J({\theta }_0,{\theta }_1)=\frac{1}{2m}\sum _{i=1}^m(h_{\theta }(x^{(i)})-y^{(i)}{)}^2$$
有两个参数，使用梯度下降法参数更新为

用python写code

import numpy as np
import matplotlib.pyplot as plt

x = np.array([1,3,4,7,9,10,11,18])
y = np.array([5,12,12,20,32,30,33,60])

def Gradient_descent(x):
    alpha = 0.001; m = len(x)
    theta0,theta1 = np.random.randint(1,6,2) #initialize two parameters with int between [1,6)
    while(True):
        h = theta0 + theta1*x
        J = 1/(2*m)*np.sum(((h-y)**2))
        theta0 = theta0 - alpha*1/m*np.sum((h-y))
        theta1 = theta1 - alpha*1/m*np.sum(((h-y)*x))
        print(J,theta0,theta1)
        if J <=5:
            return theta0, theta1
            break

if __name__ == "__main__":
    theta0,theta1 = Gradient_descent(x)
    #plot data
    x1 = np.linspace(np.min(x),np.max(x)+1)
    y1 = theta0 + theta1*x1
    plt.scatter(x,y,marker='x')
    plt.plot(x1,y1)
    plt.show()

结果如下

43.375 3.00675 2.083625
36.64715182128906 3.012834703125 2.159869203125
31.055311490596658 3.0183128984472654 2.2293845909140626
26.407644730247434 3.02323818189537 2.2927655520599455
22.544711758531772 3.0276594149910023 2.3505539698782667
19.333996925479305 3.03162114306322 2.4032438603746296
16.665362113113904 3.0351639765197063 2.45128560060768
14.447267906771756 3.038324938438401 2.4950897835393393
12.603632886217634 3.0411377814545903 2.5350307323665024
11.071223279127338 3.0436332766557497 2.571449704413933
9.797483416187205 3.045839476956834 2.604657812010998
8.738732550924968 3.0477819572102907 2.634938685352499
7.858666177326327 3.0494840331059296 2.6625509001354555
7.127110425582591 3.050966960734257 2.687730190750377
6.5189867994023265 3.0522501185213633 2.710691467970093
6.013451735119135 3.0533511730925773 2.731630658405858
5.593181460948169 3.0542862304845384 2.7507263814749403
5.243777619977898 3.0550699739999385 2.768141478233133
4.9532732638769845 3.055715789884853 2.784024405157705

第一列是平方误差，可以看到误差在逐渐减小。第二第三列分别是${\theta }_0和{\theta }_1$
注意：code中的alpha是learning rate，这里我取0.001；此处不能取很大，比如取1则结果会越来越差，平方误差会越来越大。

另外，这里设置了J<5时输出，也可以设置<0.01输出，但时间久一点，或者通常的，可以设置迭代的次数，比如对平方误差损失函数迭代计算100次就输出。

1 对结果分析

画出两个参数和损失函数的关系图如下

我们通过上面梯度下降法找到一条路径使得J越来越小，即找到最小值点。
在上面程序的基础上，如果我们画出theta0,theta1和J的函数的三维图，则如下：

可以看到找了一条路径来找到最小值
代码：

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

x = np.array([1,3,4,7,9,10,11,18])
y = np.array([5,12,12,20,32,30,33,60])

x_scale = []
y_scale = []
z_scale = []
def Gradient_descent(x):
    alpha = 0.001; m = len(x)
    theta0,theta1 = np.random.randint(1,6,2) #initialize two parameters with int between [1,6)
    while(True):
        h = theta0 + theta1*x
        J = 1/(2*m)*np.sum(((h-y)**2))
        theta0 = theta0 - alpha*1/m*np.sum((h-y))
        theta1 = theta1 - alpha*1/m*np.sum(((h-y)*x))
        x_scale.append(theta0)
        y_scale.append(theta1)
        z_scale.append(J)
#        print(J,theta0,theta1)
        if J <=5:
            return theta0, theta1
            break

if __name__ == "__main__":
    theta0,theta1 = Gradient_descent(x)
    #plot data
#    x1 = np.linspace(np.min(x),np.max(x)+1)
#    y1 = theta0 + theta1*x1
#    plt.scatter(x,y,marker='x')
#    plt.plot(x1,y1)
#    plt.show()
    fig = plt.figure()
    ax = fig.add_subplot(111,projection='3d')
    ax.scatter(x_scale,y_scale,z_scale)
    plt.show()

2 三个参数梯度下降法

现在将参数调到三个，即
$$h_{\theta }(x)={\theta }_0+{\theta }_1{x}_1+{\theta }_2{x}_2$$
比如我先假定结果是
$$y=2+3x_1+7x_2$$
先随机给出$x_1和x_2的几个值$，求出y作为真实结果

import matplotlib.pyplot as plt
import numpy as np

x1 = [0,1,2.2,3,3.4,5,8,10,13,15]
x2 = [2,3,4,5,2,3,4,5,18,23]

y = 2+3*np.asfarray(x1) + 7*np.asfarray(x2)

fig = plt.figure()
ax = fig.gca(projection = '3d')

ax.plot(x1,x2,y)
plt.show()

y对应的值为

16
26
36.6
46
26.2
38
54
67
167
208

图如下

因此再总结一下数据是：

x1 = [0,1,2.2,3,3.4,5,8,10,13,15]
x2 = [2,3,4,5,2,3,4,5,18,23]
y = [ 16. ,  26. ,  36.6,  46. ,  26.2,  38. ,  54. ,  67. , 167. ,
       208. ]

现在假设参数${\theta }_0,{\theta }_1,{\theta }_2$都是未知的，要通过上面的数据使用梯度下降法得到这三个参数。
使用同上面两个参数一样的更新方法（这里$x_0=1$）

我J的结果小于1的时候输出值，即

结果的最后几行为：$第一列是J，第二列到第四列分别是{\theta }_0,{\theta }_1,{\theta }_2$

即结果近似是
$$h=4+3x_1+7x_2$$
增加数据量并将J小于0.01才输出，则结果如下
$${\theta }_0=2.077,{\theta }_1=2.997,{\theta }_2=6.997$$
已经很接近真实函数
对于N个参数的拟合与上面一样

使用正规方程Normal Equation求解

对于线性回归问题，其实要得到最小值就是求导等于0即可得到参数
推导过程要用到矩阵的导数，在此不说了，可以参考课程的文本。参数求导=0可以推导得到，参数向量为


代码为

import numpy as np
x1 = [0,1,2.2,3,3.4,5,8,10,13,15]
x2 = [2,3,4,5,2,3,4,5,18,23]
x0 = [1]*len(x1)
y = [ 16. ,  26. ,  36.6,  46. ,  26.2,  38. ,  54. ,  67. , 167. ,
       208. ]
y = np.asfarray(y).T

X = np.mat((x0,x1,x2)).T
theta_left = np.dot(np.dot((np.dot(X.T,X)).I,X.T),y)

最后给出直接调用sklearn包的方法。

通过调用包可以很快的完成，因为所有的方法已经被优秀的人写入包中，代码和解释如下：

# import packages that may be used
import numpy as np
from sklearn import linear_model

#make data linear regression from datasets
x1 = [0,1,2.2,3,3.4,5,8,10,13,15]
x2 = [2,3,4,5,2,3,4,5,18,23]

X = np.array([[0,2],[1,3],[2.2,4],[3,5],[3.4,2],[5,3],[6,4],[10,5],[13,18],[15,23]])

# y = 2+3*x1+7*x2
Y = 2+np.dot(X,np.array([3,2]))

#create linear regression object
reg = linear_model.LinearRegression()

#Train the model using the training sets
reg.fit(X,Y)

#Make predictions using the testing set
test_pred = reg.predict(np.array([[4,2]]))

# The coefficients 
print('coefficients and intercept >> ',(reg.coef_, reg.intercept_))

从输出的结果知道
$${\theta }_0=1.999999,{\theta }_1=3,{\theta }_2=7$$