多变量线性回归 Multiple Variable Linear Regression

扩展我们的回归模型以支持多特征
- 扩展数据结构以支持多特征
- 重写预测、代价和梯度的相关程序以支持多特征
- 使用 NumPy 的 np.dot 函数对实现进行向量化处理，以提高速度和简化代码
使用以下工具：
- NumPy，一个流行的科学计算库
- Matplotlib，一个流行的数据绘图库

import copy, math
import numpy as np
import matplotlib.pyplot as plt
plt.style.use('./deeplearning.mplstyle')
np.set_printoptions(precision=2)  # reduced display precision on numpy arrays

符号 Notation

以下是一些记号的总结，更新以适应多特征场景。

机器学习通用符号 General Notation in Machine Learning

符号 (Notation)	描述 (Description)
$x_j$	第 $j$ 个特征（feature）
$n$	特征的总数量（number of features）
$\vec{x}^{(i)}$	第 $i$ 个训练样本的特征集合（features of the $i^{th}$ training example）
$x_j^{(i)}$	第 $i$ 个训练样本中第 $j$ 个特征的值（value of feature $j$ in the $i^{th}$ training example）

通用符号 (General Notation)

记号 (Notation)	描述 (Description)	Python（如适用）
$a$	标量 (scalar), 非加粗
$\mathbf{a}$	向量 (vector), 加粗
$\mathbf{A}$	矩阵 (matrix), 加粗大写

回归相关符号 (Regression Notation)

记号 (Notation)	描述 (Description)	Python（如适用）
$\mathbf{X}$	训练样本矩阵 (training example matrix)	`X_train`
$\mathbf{y}$	训练样本目标值 (training example targets)	`y_train`
$\mathbf{x}^{(i)}, y^{(i)}$	第 $i$ 个训练样本 ( $i_{th}$ Training Example)	`X[i]`, `y[i]`
$m$	训练样本数 (number of training examples)	`m`
$n$	每个样本的特征数 (number of features in each example)	`n`
$\mathbf{w}$	参数：权重 (parameter: weight)	`w`
$b$	参数：偏置 (parameter: bias)	`b`
$f_{\mathbf{w},b}(\mathbf{x}^{(i)})$	模型在第 $i$ 个样本 $\mathbf{x}^{(i)}$ 上的预测值，公式为： $f_{\mathbf{w},b}(\mathbf{x}^{(i)}) = \mathbf{w} \cdot \mathbf{x}^{(i)} + b$	`f_wb`

问题描述 Problem Statement

你将以房价预测作为动机示例。训练数据集包含三个样本，每个样本有四个特征（面积、卧室数量、楼层数量和房龄），如下表所示。注意，与之前的实验不同，面积以平方英尺 (sqft) 为单位，而不是以千平方英尺 (1000 sqft) 为单位。这引发了一个问题，你将在下一次实验中解决！

面积 (平方英尺)	卧室数量	楼层数量	房龄	价格 (千美元)
2104	5	1	45	460
1416	3	2	40	232
852	2	1	35	178

你将使用这些数据构建一个线性回归模型，以预测其他房屋的价格。例如，预测一套面积为 1200 平方英尺、3 间卧室、1 层楼、房龄 40 年的房屋价格。

请运行以下代码单元来创建 X_train 和 y_train 变量。

X_train = np.array([[2104, 5, 1, 45], [1416, 3, 2, 40], [852, 2, 1, 35]])
y_train = np.array([460, 232, 178])

包含示例的矩阵 $X$ Matrix X containing our examples

类似于上表，示例存储在一个 NumPy 矩阵 X_train 中。矩阵的每一行表示一个示例。当你有 $m$ 个训练示例（在我们的例子中 $m$ 为 3），并且有 $n$ 个特征（在我们的例子中 $n$ 为 4）时， $\mathbf{X}$ 是一个维度为 ( $m$ , $n$ ) 的矩阵（ $m$ 行， $n$ 列）。

$\mathbf{X} = \begin{pmatrix} x^{(0)}_0 & x^{(0)}_1 & \cdots & x^{(0)}_{n-1} \\ x^{(1)}_0 & x^{(1)}_1 & \cdots & x^{(1)}_{n-1} \\ \cdots \\ x^{(m-1)}_0 & x^{(m-1)}_1 & \cdots & x^{(m-1)}_{n-1} \end{pmatrix}$

符号说明：

$\mathbf{x}^{(i)}$ 是包含示例 $i$ 的向量。 $\mathbf{x}^{(i)} = (x^{(i)}_0, x^{(i)}_1, \cdots,x^{(i)}_{n-1})$
$x^{(i)}_j$ 是示例 $i$ 中的第 $j$ 个元素。括号中的上标表示示例编号，而下标表示元素编号。

显示输入数据。

# 数据存储在 NumPy 数组/矩阵中 data is stored in numpy array/matrix
print(f"X Shape: {X_train.shape}, X Type:{type(X_train)})")
print(X_train)
print(f"y Shape: {y_train.shape}, y Type:{type(y_train)})")
print(y_train)

X Shape: (3, 4), X Type:<class 'numpy.ndarray'>)
[[2104    5    1   45]
 [1416    3    2   40]
 [ 852    2    1   35]]
y Shape: (3,), y Type:<class 'numpy.ndarray'>)
[460 232 178]

参数向量 $\mathbf{w}$ 和 $b$

\mathbf{w} 是一个包含 n 个元素的向量。
- 每个元素对应一个特征的参数。
- 在我们的数据集中， $n$ 为 4。
- 理论上，我们将其表示为一个列向量：

$\mathbf{w} = \begin{pmatrix} w_0 \\ w_1 \\ \cdots\\ w_{n-1} \end{pmatrix}$

$b$ 是一个标量参数（scalar parameter）。

b_init = 785.1811367994083
w_init = np.array([ 0.39133535, 18.75376741, -53.36032453, -26.42131618])
print(f"w_init shape: {w_init.shape}, b_init type: {type(b_init)}")

w_init shape: (4,), b_init type: <class 'float'>

线性回归的正规方程 Normal Equation for Linear Regression

求解过程

代价函数
线性回归通过最小化均方误差（Mean Squared Error, MSE）的代价函数来确定参数 $\mathbf{w}$ 和 $b$ 。代价函数为：
$J(\mathbf{w}) = \frac{1}{2m} \sum_{i=1}^m \left( f_{\mathbf{w}, b}(\mathbf{x}^{(i)}) - y^{(i)} \right)^2$
其中：
$f_{\mathbf{w}, b}(\mathbf{x}) = \mathbf{w}^T \mathbf{x} + b$
假设通过将偏置 $b$ 包含在 $\mathbf{w}$ 中，可以将模型简化为：
$f_{\mathbf{w}}(\mathbf{x}) = \mathbf{w}^T \mathbf{x}$
添加一列全为 1 的列向量到矩阵 $\mathbf{X}$ ，代价函数重写为：
$J(\mathbf{w}) = \frac{1}{2m} \| \mathbf{X} \mathbf{w} - \mathbf{y} \|^2$
展开代价函数
代价函数的平方误差项展开为：
$J(\mathbf{w}) = \frac{1}{2m} (\mathbf{X} \mathbf{w} - \mathbf{y})^T (\mathbf{X} \mathbf{w} - \mathbf{y})$
展开括号：
$J(\mathbf{w}) = \frac{1}{2m} \left( \mathbf{w}^T \mathbf{X}^T \mathbf{X} \mathbf{w} - 2 \mathbf{y}^T \mathbf{X} \mathbf{w} + \mathbf{y}^T \mathbf{y} \right)$
对 $\mathbf{w}$ 求导
为了找到使代价函数最小的 $\mathbf{w}$ ，需要对 $J(\mathbf{w})$ 对 $\mathbf{w}$ 求导，并令其等于 0：
$\frac{\partial J(\mathbf{w})}{\partial \mathbf{w}} = \frac{1}{m} \left( \mathbf{X}^T \mathbf{X} \mathbf{w} - \mathbf{X}^T \mathbf{y} \right)$
令导数为 0：
$\mathbf{X}^T \mathbf{X} \mathbf{w} = \mathbf{X}^T \mathbf{y}$
求解 $\mathbf{w}$
假设 $\mathbf{X}^T \mathbf{X}$ 可逆，左乘 $(\mathbf{X}^T \mathbf{X})^{-1}$ ：
$\mathbf{w} = (\mathbf{X}^T \mathbf{X})^{-1} \mathbf{X}^T \mathbf{y}$
条件

可逆性： $\mathbf{X}^T \mathbf{X}$ 必须是满秩矩阵，即没有冗余特征，否则矩阵不可逆。此时需使用正则化方法（如岭回归）。
时间复杂度：矩阵求逆的时间复杂度为 $O(n^3)$ ，当特征数 $n$ 很大时计算代价较高。

多变量模型预测 Model Prediction With Multiple Variables

多变量模型的预测由以下线性模型给出：

$\begin{aligned} (1) & f_{\mathbf{w},b}(\mathbf{x}) = w_0x_0 + w_1x_1 +... + w_{n-1}x_{n-1} + b \end{aligned}$

或者用向量表示为：
$\begin{aligned} (2) f_{\mathbf{w},b}(\mathbf{x}) = \mathbf{w} \cdot \mathbf{x} + b \end{aligned}$

其中， $\cdot$ 表示向量的“点积”（dot product）。

为了演示点积的计算，我们将使用公式 (1) 和公式 (2) 来实现预测。

单个预测的逐元素计算 Single Prediction element by element

我们之前的预测是将一个特征值与对应的参数相乘，然后加上偏置参数 $b$ （bias parameter）。将之前的预测方法直接扩展到多特征的情况，可以通过对公式 (1) 中的每个元素进行循环计算来实现：将每个特征值与对应的参数相乘，最后再加上偏置参数。

def predict_single_loop(x, w, b): 
    """
    使用线性回归进行单个预测

    参数:
      x (ndarray): 形状为 (n,) 的数组，表示具有多个特征的示例
      w (ndarray): 形状为 (n,) 的数组，表示模型参数    
      b (scalar): 标量，表示模型的偏置参数     
      
    返回:
      p (scalar): 标量，表示预测值
    """
    n = x.shape[0]  # 获取特征数量
    p = 0  # 初始化预测值为 0
    for i in range(n):  # 遍历每个特征
        p_i = x[i] * w[i]  # 计算当前特征与对应参数的乘积
        p = p + p_i        # 将乘积累加到预测值
    p = p + b              # 最后加上偏置参数
    return p  # 返回最终预测值

# 从训练数据中获取一行
x_vec = X_train[0,:]
print(f"x_vec 的形状: {x_vec.shape}, x_vec 的值: {x_vec}")

# 进行一次预测
f_wb = predict_single_loop(x_vec, w_init, b_init)
print(f"f_wb 的形状: {f_wb.shape}, 预测值: {f_wb}")

x_vec 的形状: (4,), x_vec 的值: [2104    5    1   45]
f_wb 的形状: (), 预测值: 459.9999976194083

注意 x_vec 的形状。它是一个包含 4 个元素的一维 NumPy 向量，形状为 (4,)。结果 f_wb 是一个标量。

单个预测，向量化 Single Prediction, vector

注意，上述公式 (1) 可以通过公式 (2) 的点积运算实现。我们可以利用向量化操作来加速预测过程。

回顾 Python/NumPy 实验内容，NumPy 的 np.dot() 函数（文档链接）可用于执行向量点积运算。

def predict(x, w, b): 
    """
    使用线性回归进行单个预测

    参数:
      x (ndarray): 形状为 (n,) 的数组，表示具有多个特征的示例
      w (ndarray): 形状为 (n,) 的数组，表示模型参数   
      b (scalar): 标量，表示模型的偏置参数 
      
    返回:
      p (scalar): 标量，表示预测值
    """
    p = np.dot(x, w) + b  # 计算向量点积并加上偏置参数
    return p

# 从训练数据中获取一行
x_vec = X_train[0,:]
print(f"x_vec 的形状: {x_vec.shape}, x_vec 的值: {x_vec}")

# 进行一次预测
f_wb = predict(x_vec, w_init, b_init)
print(f"f_wb 的形状: {f_wb.shape}, 预测值: {f_wb}")

x_vec 的形状: (4,), x_vec 的值: [2104    5    1   45]
f_wb 的形状: (), 预测值: 459.9999976194083

结果和形状与之前使用循环的版本相同。接下来，这些操作将使用 np.dot 来实现。现在预测仅需一条语句即可完成。大多数程序将直接实现预测，而不是调用单独的预测函数。

多变量情况下的成本计算 Compute Cost With Multiple Variables

多变量情况下的成本函数 $J(\mathbf{w},b)$ 的公式为：
$\begin{aligned} (3) & J(\mathbf{w},b) = \frac{1}{2m} \sum\limits_{i = 0}^{m-1} (f_{\mathbf{w},b}(\mathbf{x}^{(i)}) - y^{(i)})^2 \end{aligned}$

其中：
$\begin{aligned} (4) & f_{\mathbf{w},b}(\mathbf{x}^{(i)}) = \mathbf{w} \cdot \mathbf{x}^{(i)} + b \end{aligned}$

与之前的实验不同， $\mathbf{w}$ 和 $\mathbf{x}^{(i)}$ 是支持多特征的向量，而非标量。

以下是公式 (3) 和 (4) 的实现。注意，这里使用了本课程的标准模式，即通过一个 for 循环遍历所有 $m$ 个样本。

def compute_cost(X, y, w, b): 
    """
    计算成本函数
    参数:
      X (ndarray (m,n)): 数据，包含 m 个样本，每个样本有 n 个特征
      y (ndarray (m,)) : 目标值
      w (ndarray (n,)) : 模型参数  
      b (scalar)       : 模型的偏置参数
      
    返回:
      cost (scalar): 成本
    """
    m = X.shape[0]  # 获取样本数量 m
    cost = 0.0
    for i in range(m):                                
        f_wb_i = np.dot(X[i], w) + b           # (n,) 和 (n,) 的点积结果是标量 (参见 np.dot)
        cost = cost + (f_wb_i - y[i])**2       # 累加预测值与目标值差值的平方，结果是标量
    cost = cost / (2 * m)                      # 计算平均成本，结果是标量    
    return cost

# 使用我们预先选择的最佳参数计算并显示成本
cost = compute_cost(X_train, y_train, w_init, b_init)
print(f'在最佳 w 下的误差成本: {cost}')

在最佳 w 下的误差成本: 1.5578904045996674e-12

预期结果：在最佳 $w$ 下的成本：1.5578904045996674e-12

多变量的梯度下降法

多变量情况下的梯度下降公式为：

$(5) \begin{align*} \text{repeat}&\text{ until convergence:} \; \lbrace \newline & w_j = w_j - \alpha \frac{\partial J(\mathbf{w},b)}{\partial w_j}\; & \text{for j = 0..n-1}\newline & b\ \ = b - \alpha \frac{\partial J(\mathbf{w},b)}{\partial b} \newline \rbrace \end{align*}$

其中， $n$ 是特征数量，参数 $w_j$ 和 $b$ 同时更新，并且：

$(6) \begin{align} \frac{\partial J(\mathbf{w},b)}{\partial w_j} &= \frac{1}{m} \sum\limits_{i = 0}^{m-1} (f_{\mathbf{w},b}(\mathbf{x}^{(i)}) - y^{(i)})x_{j}^{(i)} \end{align}$

$(7) \begin{align} \frac{\partial J(\mathbf{w},b)}{\partial b} &= \frac{1}{m} \sum\limits_{i = 0}^{m-1} (f_{\mathbf{w},b}(\mathbf{x}^{(i)}) - y^{(i)}) \end{align}$

$m$ 是数据集中训练样本的数量。
$f_{\mathbf{w},b}(\mathbf{x}^{(i)})$ 是模型的预测值， $y^{(i)}$ 是目标值。

多变量情况下计算梯度 Compute Gradient with Multiple Variables

以下是计算公式 (6) 和 (7) 的实现方法。这种实现方式有很多种，在本版本中：

外层循环：遍历所有 m 个样本。
- 对于每个样本，可以直接计算并累加 $\frac{\partial J(\mathbf{w},b)}{\partial b}$ 。
- 内层循环：遍历所有 n 个特征：
  - 为每个 $w_j$ 计算 $\frac{\partial J(\mathbf{w},b)}{\partial w_j}$ 。

def compute_gradient(X, y, w, b): 
    """
    计算线性回归的梯度
    参数:
      X (ndarray (m,n)): 数据，包含 m 个样本，每个样本有 n 个特征 (Data, m examples with n features) 
      y (ndarray (m,)) : 目标值 (target values)
      w (ndarray (n,)) : 模型参数 (model parameters)
      b (scalar)       : 模型的偏置参数 (model parameter)
      
    返回:
      dj_dw (ndarray (n,)): 成本函数对参数 w 的梯度 (The gradient of the cost w.r.t. the parameters w.)
      dj_db (scalar):       成本函数对参数 b 的梯度 (The gradient of the cost w.r.t. the parameter b.)
    """
    m, n = X.shape           # 获取样本数 (m) 和特征数 (n)
    dj_dw = np.zeros((n,))   # 初始化对 w 的梯度为零向量
    dj_db = 0.               # 初始化对 b 的梯度为零标量

    for i in range(m):                             
        err = (np.dot(X[i], w) + b) - y[i]         # 计算预测值与目标值的误差
        for j in range(n):                         
            dj_dw[j] = dj_dw[j] + err * X[i, j]    # 累加每个特征的梯度
        dj_db = dj_db + err                        # 累加对偏置参数的梯度
    dj_dw = dj_dw / m                              # 计算平均梯度
    dj_db = dj_db / m                              # 计算平均梯度
        
    return dj_db, dj_dw

# 计算并显示梯度
tmp_dj_db, tmp_dj_dw = compute_gradient(X_train, y_train, w_init, b_init)
print(f'初始 w 和 b 下的 dj_db: {tmp_dj_db}')
print(f'初始 w 和 b 下的 dj_dw: \n {tmp_dj_dw}')

初始 w 和 b 下的 dj_db: -1.6739251122999121e-06
初始 w 和 b 下的 dj_dw: 
 [-2.73e-03 -6.27e-06 -2.22e-06 -6.92e-05]

多变量的梯度下降 Gradient Descent With Multiple Variables

下面的程序实现了上述公式 (5)。

def gradient_descent(X, y, w_in, b_in, cost_function, gradient_function, alpha, num_iters): 
    """
    批量梯度下降法 (Batch Gradient Descent) 用于优化参数。通过使用学习率 alpha 进行 num_iters 次梯度更新来学习参数 theta。
    
    参数:
      X (ndarray (m,n))   : 数据，包含 m 个样本，每个样本有 n 个特征
      y (ndarray (m,))    : 目标值 (target values)
      w_in (ndarray (n,)) : 初始模型参数 (initial model parameters)
      b_in (scalar)       : 初始偏置参数 (initial model parameter)
      cost_function       : 用于计算成本的函数 (function to compute cost)
      gradient_function   : 用于计算梯度的函数 (function to compute the gradient)
      alpha (float)       : 学习率 (Learning rate)
      num_iters (int)     : 梯度下降的迭代次数 (number of iterations to run gradient descent)
      
    返回:
      w (ndarray (n,)) : 更新后的参数值 (Updated values of parameters)
      b (scalar)       : 更新后的偏置值 (Updated value of parameter)
      """
    
    # 用于存储每次迭代的成本 J 和参数 w 的数组，主要用于后续绘图
    J_history = []
    w = copy.deepcopy(w_in)  # 避免在函数中修改全局 w
    b = b_in
    
    for i in range(num_iters):

        # 计算梯度并更新参数
        dj_db, dj_dw = gradient_function(X, y, w, b)

        # 使用 w, b, alpha 和梯度更新参数
        w = w - alpha * dj_dw
        b = b - alpha * dj_db
      
        # 在每次迭代时保存成本 J
        if i < 100000:      # 防止资源耗尽 
            J_history.append(cost_function(X, y, w, b))

        # 每 10% 的迭代或总迭代次数不足 10 次时，打印成本
        if i % math.ceil(num_iters / 10) == 0:
            print(f"Iteration {i:4d}: Cost {J_history[-1]:8.2f}   ")
        
    return w, b, J_history # 返回最终的 w, b 和 J 历史记录用于绘图

# 初始化参数
initial_w = np.zeros_like(w_init)  # 初始化权重为与 w_init 相同形状的零向量
initial_b = 0.                     # 初始化偏置为 0

# 梯度下降的设置
iterations = 1000                  # 迭代次数
alpha = 5.0e-7                     # 学习率

# 运行梯度下降
w_final, b_final, J_hist = gradient_descent(X_train, y_train, initial_w, initial_b,
                                                    compute_cost, compute_gradient, 
                                                    alpha, iterations)

print(f"通过梯度下降找到的 b 和 w：b = {b_final:0.2f}, w = {w_final} ")

m, _ = X_train.shape  # 获取样本数量
for i in range(m):
    prediction = np.dot(X_train[i], w_final) + b_final  # 计算预测值
    print(f"预测值: {prediction:0.2f}, 目标值: {y_train[i]}")

Iteration    0: Cost  2529.46   
Iteration  100: Cost   695.99   
Iteration  200: Cost   694.92   
Iteration  300: Cost   693.86   
Iteration  400: Cost   692.81   
Iteration  500: Cost   691.77   
Iteration  600: Cost   690.73   
Iteration  700: Cost   689.71   
Iteration  800: Cost   688.70   
Iteration  900: Cost   687.69   
通过梯度下降找到的 b 和 w：b = -0.00, w = [ 0.2   0.   -0.01 -0.07] 
预测值: 426.19, 目标值: 460
预测值: 286.17, 目标值: 232
预测值: 171.47, 目标值: 178

# 绘制成本与迭代次数的关系图
fig, (ax1, ax2) = plt.subplots(1, 2, constrained_layout=True, figsize=(12, 4))
ax1.plot(J_hist)
ax2.plot(100 + np.arange(len(J_hist[100:])), J_hist[100:])
ax1.set_title("Cost vs. iteration");  ax2.set_title("Cost vs. iteration (tail)")
ax1.set_ylabel('Cost')             ;  ax2.set_ylabel('Cost') 
ax1.set_xlabel('iteration step')   ;  ax2.set_xlabel('iteration step') 
plt.show()

png

这些结果并不令人满意! 成本仍在下降，而我们的预测并不十分准确。在下一个实验中，我们将探讨如何改进这一点。