线性回归最小二乘法与梯度下降法数学推导及 Python 实现

本文将介绍:

  • 线性回归最小二乘法代数推导及实现
  • 线性回归最小二乘法矩阵推导及实现
  • 线性回归梯度下降法数学推导及实现

* 以下内容,未经允许,禁止复制转载。但可以在注明出处前提下引用。

线性回归是一种十分基础的机器学习算法。但是基础并不意味着不重要,相反,基础往往才是最重要的。很多人入门机器学习,第一个学习的就是线性回归。例如,吴恩达的机器学习课程中,也是以线性回归开始。

对于一个线性回归问题,一般来讲有 2 种解决方法,分别是:最小二乘法和梯度下降法。其中最小二乘法又分为两种求解思路:代数求解和矩阵求解。

接下来,我将梳理线性回归问题的求解过程,以及使用 Python 进行编程实现的完整过程。

最小二乘法代数求解推导

最小二乘法是线性回归问题最常用的解决方法。我们这里使用到的是普通最小二乘法,英文为:Ordinary Least Squares,简称 OLS。最小二乘法中的「二乘」代表平方,最小二乘也就是最小平方。而这里的平方就是指代平方损失函数。首先,假设我们需要求解的是一元线性回归问题,及其对应的线性方程如下:

equation--10-

简单来讲,最小二乘法也就是求解平方损失函数最小值的方法。那么,到底该怎样求解呢?这就需要使用到高等数学中的知识。推导如下。

我们所知,平方损失函数的公式为:

equation--7-

由于最终的目标是求取平方损失函数 $min(f)$ 最小时对应的 $w$。所以,先求 $f$ 的一阶偏导数:

equation--8-

然后,我们令 $\frac{\partial f}{\partial w_{0}}=0$ 以及 $\frac{\partial f}{\partial w_{1}}=0$,解得:

equation--9-

至此,我们就用代数的方式求解出线性方程的参数了。其中,$w_0$ 经常习惯称之为截距项。

最小二乘法代数求解实现

接下来,我们依据上面的推导过程,使用 Python 进行实现。这里,假设进行线性拟合的数据如下:

x = [55, 71, 68, 87, 101, 87, 75, 78, 93, 73]
y = [91, 101, 87, 109, 129, 98, 95, 101, 104, 93]

然后,根据公式实现参数计算函数 ols_algebra(x, y)

def ols_algebra(x, y):
    n = len(x)
    w1 = (n*sum(x*y) - sum(x)*sum(y))/(n*sum(x*x) - sum(x)*sum(x))
    w0 = (sum(x*x)*sum(y) - sum(x)*sum(x*y))/(n*sum(x*x)-sum(x)*sum(x))
    return w0, w1

最后,就可以得到拟合后的参数值。

import numpy as np

ols_algebra(np.array(x), np.array(y))

>>> (44.25604341391219, 0.7175629008386778)

那么,所对应的线性方程为:$y = 0.718 * x + 44.256$。

最小二乘法矩阵求解推导

除了上面的代数推导过程。通常在实际情况中,更多会使用矩阵的方式来求线性拟合参数。因为在大规模运算中,矩阵的计算效率更高。

矩阵的推导会稍微复杂一些。同样,根据一元线性函数的表达式 $ y(x, w) = w_0 + w_1x$。那么,转换成矩阵形式为:

equation

即:

equation--1-

然后,平方损失函数转换为矩阵形式:

equation--4-

此时,对矩阵求偏导数得到:

equation--5-

当矩阵 $X^TX$ 满秩(不满秩无法通过 OLS 方法求解)时,$(X^TX)^{-1}X^TX=E$,且 $EW=W$。所以,$(X^TX)^{-1}X^TXW=(X^TX)^{-1}X^Ty$。最终得到:

equation--6-

以上,我们完成了线性回归矩阵求解推导。

最小二乘法矩阵求解实现

使用矩阵方法求解时的代码实现过程会更加简单,因为最终参数的求解就只有一个公式。实现 ols_matrix(x, y) 函数如下:

def ols_matrix(x, y):
    w = (x.T * x).I * x.T * y
    return w

同样,这里使用上面提供的示例数据。注意,由于是矩阵运算,所以需要对原数据进行形状变换,并添加截距项系数 1,否则会因为形状问题无法计算。

x = np.array(x).reshape(len(x), 1)
x = np.concatenate((np.ones_like(x), x), axis=1)
x = np.matrix(x)

y = np.array(y).reshape(len(y), 1)
y = np.matrix(y)

ols_matrix(x, y)

>>> matrix([[44.25604341], [ 0.7175629 ]])

可以看到,这里得到的结果与代数求解方法一致。

上面我们使用最小二乘法来求解线性方程系数,得到的结果是精确的解析解。最后,尝试将拟合的直线绘制成图,查看拟合效果。这里,四舍五入保留参数的 3 为小数。

$$y = 44.256 + 0.718 * x$$

import matplotlib.pyplot as plt

plt.scatter(x, y)
plt.plot(np.array([50, 110]), np.array([50, 110]) * 0.718 + 44.256, 'r')

--

梯度下降法求解数学推导

介绍完最小二乘法,下面来说一说梯度下降法的求解过程,关于梯度下降法的原理和过程这里不再解释。梯度下降法作为迭代优化方法,需要提前定义好损失函数。这里选择均方误差 MSE 作为损失函数。

equation--13-

其中,$y_{i}$ 表示真实值,$\hat y_{i}$ 表示近似值,$n$ 则表示值的个数。近似值则通过当前已知计算,那么:

equation--14-

此时,我们求解参数和截距项对应的梯度:

equation--15-

那么,根据梯度就可以执行迭代更新,从而找到最优参数:

equation--16-

梯度下降法求解代码实现

相比于前面的最小二乘法实现过程,梯度下降法的 Python 实现会稍微麻烦一些。不过通过上面的公式,我们也能很快地完成代码编写:

def gradient_descent(x, y, lr, num_iter):
    
    w1 = 0  # 初始参数为 0
    w0 = 0  # 初始参数为 0
    
    for i in range(num_iter):  # 梯度下降迭代
        # 计算近似值
        y_hat = (w1 * x) + w0
        # 计算参数对应梯度
        w1_gradient = -(2/len(x)) * sum(x * (y - y_hat))
        w0_gradient = -(2/len(x)) * sum(y - y_hat)
        # 根据梯度更新参数
        w1 -= lr * w1_gradient
        w0 -= lr * w0_gradient

    return w1, w0

同样,使用上面的示例数据完成参数计算:

w1, w0 = gradient_descent(np.array(x), np.array(y), lr=0.00001, num_iter=100)
w1, w0

>>> (1.264673601681181, 0.017892580143733863)

因为是迭代方法,所以需要设置学习率和最大迭代次数。你可能发现,这里得到的结果与最小二乘法不太一致,这是不同方法带来的普遍现象。同样绘制出拟合直线:

plt.scatter(x, y)
plt.plot(np.array([50, 110]), np.array([50, 110]) * w1 + w0, 'r')

----1-

肉眼可见,最终的拟合情况不如最小二乘法。当然,如果随机初始参数以及调整合适的迭代次数及学习率,优化迭代结果可能会更理想一些。

小结

本文系统地介绍了线性回归问题的两种不同的求解方法,分别是:最小二乘法和梯度下降法。同时,对 3 种求解过程进行了编程实现,希望对机器学习初学者带来启发。如何有相关疑问,欢迎交流探讨。