线性回归（梯度下降法，正规方程法）手写推导-CFANZ编程社区

本文主要从数学原理入手，探讨如何通过已知数据，使电脑学习出线性回归方程。在上来开始详细介绍之前，我们先来看一个例子。

这是一个关于银行贷款额度的数据统计表。可以看出，在本例中银行贷款额度与贷款者的工资、年龄有关，换做计算机语言来讲，就是我们的预测函数受到两个变量的控制。如何才能得出本例的线性回归方程呢，首先从数学方面出发，通过从预测函数的函数式着手，对变量进行迭代，以下是详细步骤。

一、损失函数的数学推导：

（1）梯度下降法求线性回归方程

（2）正规方程法求线性回归方程

二、源码展示

一、损失函数的数学推导：

以下为损失函数的数学推导，不是必要内容，可以选择性跳过

在开始推导前，我们可以大致估算出我们的函数模型，函数模型的估算因人而异。不同的函数模型会影响最后的拟合效果，倘若函数模型过于复杂或者维度过高，则会出现过拟合的现象，本文不再予以讨论，如果读者感兴趣可以浏览作者有关正则化的博文，里面有详细的讲解。

实际上，引入矩阵的概念可以极大程度上方便我们的计算，如果不会矩阵的操作可以跳过这一点。通过矩阵的概念显然我们可以得出：预测函数=theta的转置乘以矩阵x

这里运用到了概率论的知识。我们得到了最后一行的结果，目的是要让P最大，显然需要让后面的式子越小越好，因而得出损失函数的数学概念：

对于我们是否带入矩阵操作，得到两个不同的表达式，这两个表达式都可以完成对线性回归方程的计算。

（1）梯度下降法求线性回归方程

需注意，在梯度下降算法中，我们需要定义每次迭代的步长alpha，alpha的大小会影响迭代次数和是否能够找到损失函数的最值，通常取小于1的数。梯度下降算法具有一定的泛用性，在参数极多的情况下也能正常运行。

（2）正规方程法求线性回归方程

在正规方程法当中，我们无需去定义步长alpha的大小，也无需迭代，即可求得参数theta。但相反，正规方程需要求解参数x矩阵的转置，即便numpy函数库可以很方便的计算转置矩阵，但在参数极多的情况下，转置矩阵的运算量会达到一个十分恐怖的数量级，因而不具有泛用性。

二、源码展示

在梯度下降算法中，参数的数量仅涉及到对参数求导、参数迭代的过程，原理上类似。作者是通过notebook环境下编写的代码，下面展示源码。

导入必须的python库，我们用到了pandas、numpy和matplotlib库。

import pandas as pd               #pandas库用来从csv表格中读取数据
import numpy as np
import matplotlib.pyplot as plt

定义求损失函数的函数模块，返回损失函数的函数值。

#求损失函数
def Loss(theta0, theta1, theta2, x1, x2, y):
    h = theta0 + theta1 * x1 + theta2 * x2
    diff = h - y                        #此时h和y都是一个5维矩阵，可以直接做减法，得到每一项与 
                                                  原始值的差值
    L = (diff.sum())**2 / (2 * x1.shape[0])    #x1.shape[0]得到的是矩阵的行数，在5维矩阵中即 
                                                      表示有几个元素
    return L

定义求各个参数偏导的函数模块，为梯度下降做准备。

#求损失函数中theta0的偏导
def partial_loss_theta0(theta0, theta1, theta2, x1, x2, y):
    h = theta0 + theta1 * x1 + theta2 * x2
    diff = h - y                                       #求导的过程用上文推导的式子即可
    partial = diff.sum() / x1.shape[0]
    return partial

#求损失函数中theta1的偏导
def partial_loss_theta1(theta0, theta1, theta2, x1, x2, y):
    h = theta0 + theta1 * x1 + theta2 * x2
    diff = (h - y) * x1                               #求导的过程用上文推导的式子即可
    partial = diff.sum() / x1.shape[0]
    return partial

#求损失函数中theta2的偏导
def partial_loss_theta2(theta0, theta1, theta2, x1, x2, y):
    h = theta0 + theta1 * x1 + theta2 * x2
    diff = (h - y) * x2                               #求导的过程用上文推导的式子即可
    partial = diff.sum() / x1.shape[0]
    return partial

梯度下降核心算法如下。

#梯度下降求损失函数最小值点
def gradient_descent(x1, x2, y, alpha=0.1, theta0 = 1, theta1 = 1, theta2 = 1, steps = 5000):
    num = 0                   #迭代次数
    result=[]
    for step in range(steps):
        num += 1
        update0 = partial_loss_theta0(theta0, theta1, theta2, x1, x2, y)
        update1 = partial_loss_theta1(theta0, theta1, theta2, x1, x2, y)
        update2 = partial_loss_theta2(theta0, theta1, theta2, x1, x2, y)
        
        theta0 = theta0 - alpha * update0                  #对三个参数同时进行迭代
        theta1 = theta1 - alpha * update1
        theta2 = theta2 - alpha * update2
        
        e = Loss(theta0, theta1, theta2, x1, x2, y)
        result.append(e)
        if e < 0.001:
            break
    
    return theta0,theta1,theta2,result

此处作者用到了对数据的特征缩放，目的是让数据更集中、更规则，从而大大提高了运算速度。此处附上notebook上的数据标准化后的结果展示：

ranges = pd.read_csv("ranges.csv")
ranges.Salary = (ranges.Salary - ranges.Salary.mean()) / ranges.Salary.std() #此处对原本的数据进行了标准化
ranges.Age = (ranges.Age - ranges.Age.mean()) / ranges.Age.std()        #标准化后的数据会大大减少运算次数
ranges.Range = (ranges.Range - ranges.Range.mean()) / ranges.Range.std()
ranges                                                   #数据标准化后的结果展示

最后定义main函数，并通过matplotlib画出损失函数梯度下降的图像.

if __name__ == "__main__":
    theta0,theta1,theta2,result = gradient_descent(ranges.Salary, ranges.Age, ranges.Range)
    print("\n线性规划方程为：\n")
    print("h(x) = %.3f + %.3f * x1 + %.3f * x2\n"%(theta0,theta1,theta2))
    n = len(result)
    x = range(n)
    plt.plot(x,result,color='r',linewidth=3)
    plt.title("Loss")
    plt.xlabel("num")
    plt.ylabel("y")
    plt.show()

我们可以看到，经过标准化后的数据，经过大概三十次左右的的计算即可基本拟合出较好的预测函数，效率很高。

一、损失函数的数学推导：

（1） 梯度下降法求线性回归方程

（2）正规方程法求线性回归方程

二、源码展示

（1）梯度下降法求线性回归方程