批量梯度下降&正规方程

---

前言

对于求解线性回归问题，尤其是多变量线性回归，可以使用批量梯度下降，也可以使用正规方程来求解，本篇主要对这两种方法进行介绍，并比较二者的特点以及使用方法。

---

一、批量梯度下降（Batch Gradient Descent）

1.算法介绍

本篇以多变量线性回归为例，对于给定的$m$组包含$n$个特征的变量 $X$（$x_j^{(i)}$表示第$i$行、第$j$个特征值），和$m$个目标值 $y$，拟合的回归方程为：
$$h_{\theta }(x)={\theta }_0+{\theta }_1{x}_1+{\theta }_2{x}_2+\cdots +{\theta }_n{x}_n$$
为简化公式，令$x_0=1$，则矩阵X的维度变为$m\ast (n+1)$，上式可写为矩阵形式：
$$h_{\theta }(x)={\theta }^{T}X$$
构建的代价函数为（cost function）：
$$J({\theta }_0,{\theta }_1,\cdots ,{\theta }_n)=\frac{1}{2m}\sum _{i=1}^m(h_{\theta }(x^{(i)})-y^{(i)}{)}^2$$
对应的批量梯度下降算法公式如下，其中 $\alpha$ 为学习率，表示每次移动的步长大小，注意更新 $n+1$ 个 $\theta$ 值时，需要同步更新，因此才称之为批量梯度下降。
$$\begin{gathered} {\theta }_j:={\theta }_j-\alpha \frac{\partial J({\theta }_0,{\theta }_1,\cdots ,{\theta }_n)}{\partial {\theta }_j} \\ ={\theta }_j-\alpha \frac{\partial }{\partial {\theta }_j}(\frac{1}{2m}\sum _{i=1}^m(h_{\theta }(x^{(i)})-y^{(i)}{)}^2) \\ ={\theta }_j-\alpha \frac{1}{m}\sum _{i=1}^m((h_{\theta }(x^{(i)})-y^{(i)})\cdot x_j^{(i)}),forj=0,1,\cdots ,n. \end{gathered}$$

2.Python代码

import numpy as np
import random
from numpy import genfromtxt

def getData(dataSet):
    m, n = np.shape(dataSet)
    trainData = np.ones((m, n))
    trainData[:,:-1] = dataSet[:,:-1] #第一列到倒数第二列为X，共n-1列，加上最后一列补为1，共为n列
    trainLabel = dataSet[:,-1] #最后一列为y
    return trainData, trainLabel

def batchGradientDescent(x, y, theta, alpha, m, maxIterations):
    xTrains = x.transpose()
    for i in range(0, maxIterations):
        hypothesis = np.dot(x, theta)
        loss = hypothesis - y
        # print loss
        gradient = np.dot(xTrains, loss) / m
        theta = theta - alpha * gradient
    return theta

def predict(x, theta):
    m, n = np.shape(x)
    xTest = np.ones((m, n+1))
    xTest[:, :-1] = x
    yP = np.dot(xTest, theta)
    return yP

dataPath = r"E:\learning\house.csv" #数据集路径，为避免各特征尺度不同，在导入之前可先将数据集归一化
dataSet = genfromtxt(dataPath, delimiter=',')
trainData, trainLabel = getData(dataSet) #从中提取X，y
m, n = np.shape(trainData) #训练数据集，m行n列，其中最后一列全为1
theta = np.ones(n) #初始化theta
alpha = 0.1 #学习率，可以尝试0.01,0.03,0.1,0.3,1,3,10
maxIteration = 5000 #最大迭代次数
theta = batchGradientDescent(trainData, trainLabel, theta, alpha, m, maxIteration)
x = np.array([[3.1, 5.5], [3.3, 5.9], [3.5, 6.3], [3.7, 6.7], [3.9, 7.1]]) #测试数据集，列数为特征值数
print predict(x, theta) #输出测试集的目标预测值

3.正则化优化

正则化线性回归函数的代价函数为：
$$J(\theta )=\frac{1}{2m}\sum _{i=1}^m[(h_{\theta }(x^{(i)})-y^{(i)}{)}^2+\lambda \sum _{j=1}^n{\theta }_j^2]$$
求导后得到对应批量梯度下降算法的公式：
$$\begin{gathered} {\theta }_0:={\theta }_0-\alpha \frac{1}{m}\sum _{i=1}^m[(h_{\theta }(x^{(i)})-y^{(i)})x_0^{(i)}] \\ {\theta }_j:={\theta }_j-\alpha [\frac{1}{m}\sum _{i=1}^m((h_{\theta }(x^{(i)})-y^{(i)})x_j^{(i)})+\frac{\lambda }{m}{\theta }_j],forj=1,2,\cdots ,n. \end{gathered}$$

二、正规方程（Normal Equation）

1.算法介绍

对于上述代价函数，将其用矩阵表示出来，即为$J(\theta )=\frac{1}{2}(X\theta -y{)}^2$，其中$X$为$m$行$n$列，$\theta$为$n$行1列，$y$为$m$行1列，可将其展开写成如下形式：
$$J(\theta )=\frac{1}{2}(X\theta -y{)}^T(X\theta -y)=\frac{1}{2}({\theta }^T{X}^T-y^T)(X\theta -y)=\frac{1}{2}({\theta }^T{X}^{T}X\theta -{\theta }^T{X}^{T}y-y^{T}X\theta +y^{T}y)$$
下面求$J(\theta )$对$\theta$的偏导数，得到：
$$\frac{\partial J(\theta )}{\partial \theta }=X^{T}X\theta -X^{T}y$$
令$\frac{\partial J(\theta )}{\partial \theta }=0$，得到$\theta =(X^{T}X{)}^{-1}{X}^{T}y$
需要注意的是，使用该算法时，需要考虑$X^{T}X$是否可逆，对于某些不可逆的情况（①特征之间不独立；②特征数量大于训练集的数量），该方法不可用。

2.Python代码

import numpy as np
import random
from numpy import genfromtxt

def getData(dataSet):
    m, n = np.shape(dataSet)
    trainData = np.ones((m, n))
    trainData[:,:-1] = dataSet[:,:-1] #第一列到倒数第二列为X，共n-1列，加上最后一列补为1，共为n列
    trainLabel = dataSet[:,-1] #最后一列为y
    return trainData, trainLabel

def normalEqu(x, y):
    a = numpy.linalg.pinv(np.dot(x.T,x))
    b = np.dot(a,x.T)
    theta = np.dot(b,y)
    return theta

def predict(x, theta):
    m, n = np.shape(x)
    xTest = np.ones((m, n+1))
    xTest[:, :-1] = x
    yP = np.dot(xTest, theta)
    return yP

dataPath = r"E:\learning\house.csv" #数据集路径，为避免各特征尺度不同，在导入之前可先将数据集归一化
dataSet = genfromtxt(dataPath, delimiter=',')
trainData, trainLabel = getData(dataSet) #从中提取X，y
m, n = np.shape(trainData) #训练数据集，m行n列，其中最后一列全为1
theta = normalEqu(trainData, trainLabel)
x = np.array([[3.1, 5.5], [3.3, 5.9], [3.5, 6.3], [3.7, 6.7], [3.9, 7.1]]) #测试数据集，列数为特征值数
print predict(x, theta) #输出测试集的目标预测值

3.正则化优化

利用正规方程来求解正则化线性回归模型，方法如下所示：
$$\theta ={\left(X^{T}X+\lambda \left[\begin{array}{ccccc}0& & & & \\ & 1& & & \\ & & 1& & \\ & & & \ddots & \\ & & & & 1\end{array}\right]\right)}^{-1}{X}^{T}y$$

三、二者比较