实验目的

熟悉和掌握单变量线性回归算法
熟悉和掌握批处理梯度下降算法
熟悉和掌握多变量线性回归算法

实验要求

采用Python、Matlab等高级语言进行编程，推荐优先选用Python语言
核心模型和算法需自主编程实现，不得直接调用Scikit-learn、PyTorch等成熟框架的第三方实现
代码可读性强：变量、函数、类等命名可读性强，包含必要的注释

实验原理

一元线性回归模型

即 $D=\left\{\left(x_i, y_i\right)\right\}_{i=1}^m$, 其中 $x_i \in R$. 对离散属性, 若属性值间存在 “序” 关系, 可以通过连续化将其转化为连续值。线性回归试图学得
$$
f\left(x_i\right)=w x_i+b \text {, 使得 } f\left(x_i\right) \cong y_i
$$

均方误差是回归任务中最常用的性能度量, 因此我们可试图让均方误差最小化, 即:
$$
\left(w^*, b^*\right)=\underset{(w, b)}{\arg \min } \sum_{i=1}^m\left(f\left(x_i\right)-y_i\right)^2=\underset{(w, b)}{\arg \min } \sum_{i=1}^m\left(w x_i+b-y_i\right)^2
$$

均方误差对应了常用的欧几里得距离。基于均方误差最小化来进行模型求解的方法称为 “最小二乘法”。在线性回归中, 最小二乘法就是试图找到一条直线, 使所有样本到直线上的欧氏距离之和最小。

求解 $w$ 和 $b$ 使 $E_{(w, b)}=\sum_{i=1}^m\left(w x_i+b-y_i\right)^2$ 最小化的过程, 称为线性回归模型的最小二乘 “参数估计”。我们可将 $E_{(w, b)}$ 分别对 $w$ 和 $b$ 求导, 得到:
$$
\begin{gathered}
\frac{\partial E_{(w, b)}}{\partial w}=2\left(w \sum_{i=1}^m x_i^2-\sum_{i=1}^m\left(y_i-b\right) x_i\right), \\
\frac{\partial E_{(w, b)}}{\partial w}=2\left(m b-\sum_{i=1}^m\left(y_i-w x_i\right)\right),
\end{gathered}
$$

然后令上式为 0 , 可以得到 $w$ 和 $b$ 最优解的闭式解:
$$
\begin{gathered}
w=\frac{\sum_{i=1}^m y_i\left(x_i-\bar{x}\right)}{\sum_{i=1}^m x_i^2-\frac{1}{m}\left(\sum_{i=1}^m x\right)^2}, \
b=\frac{1}{m} \sum_{i=1}^m\left(y_i-w x_i\right),
\end{gathered}
$$

其中, $\bar{x}=\frac{1}{m} \sum_{i=1}^m x$ 为 $x$ 的均值。

多元线性回归模型

样本由 $d$ 个属性描述, 此时我们试图学得
$$
f\left(x_i\right)=w^T x_i+b \text {, 使得 } f\left(x_i\right) \cong y_i
$$

类似的, 可利用最小二乘法来对 $w$ 和 $b$ 进行估计。
$$
X=\left(\begin{array}{ccccc}
x_{11} & x_{12} & \cdots & x_{1 d} & 1 \\
x_{21} & x_{22} & \cdots & x_{2 d} & 1 \\
\vdots & \vdots & \ddots & \vdots & 1 \\
x_{m 1} & x_{m 2} & \cdots & x_{m d} & 1
\end{array}\right)=\left(\begin{array}{cc}
x_1^T & 1 \\
x_2^T & 1 \\
\vdots & \vdots \\
x_m^T & 1
\end{array}\right)
$$

再把标记也写成向量形式 $y=\left(y_1 ; y_2 ; \ldots ; y_m\right)$, 有
$$
\widehat{w}^*=\underset{\widehat{w}^*}{\arg \min } \sum_{i=1}^m(y-X \widehat{w})^T(y-X \widehat{w})
$$

令 $E_{\widehat{w}}=(y-X \widehat{w})^T(y-X \widehat{w})$, 对 $\widehat{w}$ 求导得到
$$
\frac{\partial E_{\hat{w}}}{\partial \widehat{w}}=2 X^T(X \widehat{w}-y) \text {. }
$$

令上式为 0 可得 $\widehat{w}$ 的最优解的闭式解, 但由于设计矩阵逆的计算, 比单变量情形要复杂一些。下面我们做一个简单的讨论。
当 $X^T X$ 为满秩矩阵或正定矩阵时, 令上式为 0 , 可得
$$
\widehat{w}=\left(X^T X\right)^{-1} X^T y
$$

其中 $\left(X^T X\right)^{-1}$ 是矩阵 $\left(X^T X\right)$ 的逆矩阵, 令 $\hat{x}_i=\left(x_i ; 1\right)$. 则最终学得的多元线性回归模型为
$$
f\left(\hat{x}_i\right)=\hat{x}_i^T\left(X^T X\right)^{-1} X^T y
$$

然而, 现实任务中 $X^T X$ 往往不是满秩矩阵, 例如再许多任务中我们会遇到大量的变量,其数目甚至超过样例数, 导致 $X$ 的列数多于行数, $X^T X$ 显然不满秩。此时可解出多个 $\widehat{w}$, 它们都能使均方误差最小化。常见的方法还有引入正则项。

实验内容

单变量线性回归

采用数据集 “data/regress_data1.csv”进行单变量线性回归实验
借助 matplotlib 画出原始数据分布的散点图（x=”人口”, $y=$ “收益”）
以最小平方误差为目标函数, 构造模型的损失（误差）计算函数:
$$
J(w)=\frac{1}{2 m} \sum_{i=1}^m\left(h\left(x^{(i)}\right)-y^{(i)}\right)^2
$$

其中, $h(x)=w^T X=w_0 x_0+w_1 x_1+w_2 x_2+\ldots+w_n x_n$

实现批量梯度下降算法（Batch Gradient Decent）用于优化线性回归模型:
$$
w_j=w_j-\alpha \frac{1}{m} \sum_{i=1}^m \frac{\partial}{\partial w_j} J(w)
$$
其中, $w_j$ 为参数向量, $m$ 样本量。
提示:
批量梯度下降算法可以参考第 5 章神经网络的优化算法（P104）
线性回归模型的偏置 $\mathrm{b}$ 可吸收进参数向量 $\mathrm{w}$, 从而采用向量形式 (P55)
采用上述批量梯度下降法, 优化单变量线性回归模型。其中, 迭代轮数 epoch 设定为 1000 轮, 学习率设定为 0.01 , 参数初始化为 0 。

代码输出:

优化结束时的损失值和模型参数
将模型拟合的直线和原始数据散点图画到同一张图中

多变量线性回归

采用数据集 “data/regress_data2.csv”进行多变量线性回归实验, 通过房子的大小和房间数量两个变量回归房子的价格。
对数据进行特征归一化： $z_i=\frac{x_i-\mu}{\sigma}$
采用上述批量梯度下降法, 优化多变量线性回归模型。其中, 迭代轮数 epoch 设定为 1000 轮, 学习率设定为 0.01 , 参数初始化为 0 。

代码输出:

优化结束时的损失值和模型参数
画图输出训练误差（损失）随着迭代轮数 epoch 的变化曲线

实验代码和结果

import pandas as pd
import matplotlib.pyplot as plt
import numpy as np


def compute_loss(X, Y, theta): #损失函数的计算
    m = Y.size
    loss = 0
    for i in range(0, m):
        x = X[i, 1]
        y = Y[i]
        loss += (y - (theta[1] * x + theta[0])) ** 2
    loss = loss / (2.0 * m)
    return loss


def gradient_descent(X, Y, theta, lr, epoch): #BGD实现 
    m = Y.size
    t = theta.size
    J_history = np.zeros(epoch)

    for iter in range(epoch):
        for j in range(t):
            theta[j] = theta[j] - lr * np.sum((np.dot(X, theta) - Y) * X[:, j]) / m

        J_history[iter] = compute_loss(X, Y, theta)
    return theta, J_history


path = 'regress_data1.csv'
data = pd.read_csv(path)
X = data[['人口']]
Y = data[['收益']].values.ravel()  # Flatten Y to 1D array
m = X.size
X = np.c_[np.ones(m), X]  # Add a column of ones for the bias term
theta = np.zeros(2)
epoch = 1000
lr = 0.01
theta, J_history = gradient_descent(X, Y, theta, lr, epoch)
ine1,=plt.plot(X[:,1],np.dot(X,theta),color='red')
plt.scatter(data[['人口']],data[['收益']]);plt.xlabel("population");plt.ylabel("profit");
#plt.plot(range(1000),J_history)
plt.show()
print(theta,J_history[999])

$$
\begin{aligned}
& b=-3.25088222 \\
& w=1.12836314
\end{aligned}
$$

训练 1000 次后的损失为 4.514833339953508

我们只需对前面的代码做简单的修改即可用于多变量线性回归，将计算损失函数的改为

1	loss += (y - (theta[1] * x1 + theta[0]+theta[2]x2)) * 2

对数据进行归一化

1
2
3

X=(X-np.average(X))/np.std(X)
Y=(Y-np.average(Y))/np.std(Y)
Z=(Z-np.average(Z))/np.std(Z)

$$
\begin{aligned}
& \mathrm{b}=-8.21860310 \mathrm{e}-17 \\
& \mathrm{w} 1=8.79149410 \mathrm{e}-01 \\
& \mathrm{w} 2=4.75747766 \mathrm{e}-02
\end{aligned}
$$

迭代 1000 次的损失为 0.1335413413355335

小结或讨论

线性回归是基于假设目标变量与特征变量之间存在线性关系，然后通过拟合最佳的线性函数来预测目标变量。

线性回归建立了输入特征（自变量）与输出变量（因变量）之间的线性关系，通过拟合最佳的线性模型来预测输出变量的值。

线性回归的主要目标是通过学习样本数据集中的特征和目标变量之间的线性关系，来进行预测和解释。线性回归可以用于预测连续数值型的目标变量，也可以用于分析特征对目标变量的影响程度和方向。线性回归的目标是找到最佳的权重和截距，使得模型的预测值与实际值之间的差异最小化。

梯度下降法：①梯度下降法是一种迭代优化算法，通过不断调整模型参数来最小化损失函数。在线性回归中，梯度下降法通过计算损失函数对模型参数的梯度，并根据梯度的方向和大小更新参数值，直到达到最小化损失的目标。②梯度下降法使用学习率来控制每次迭代中参数的更新步长。学习率需要手动选择，并且影响算法的收敛速度和稳定性。③梯度下降法可能陷入局部最优解而不是全局最优解，特别是在非凸的损失函数中。这可以通过选择合适的学习率和初始化参数来缓解。④梯度下降法的收敛速度取决于学习率的选择和损失函数的形状。较小的学习率可能导致收敛速度较慢，而较大的学习率可能导致无法收敛或发散。⑤梯度下降法相对于最小二乘法对异常值具有一定的鲁棒性，但仍然可能受到大的异常值的影响。⑥梯度下降法可以主要分为批量梯度下降法、随机梯度下降法和mini-batch梯度下降法，这三种梯度下降法的选择取决于数据集的规模、计算资源的限制以及优化效果的要求。⑦梯度下降法主要适用于目标变量和特征变量之间存在近似线性关系，数据集规模较大，无法直接计算矩阵的逆矩阵，数据集中存在噪声或异常值等情况。