CONTENTS

机器学习基础 第二部分:多元线性回归Draft

欢迎来到我的 C 语言机器学习基础系列的第 2 / N 部分!如果你还没看过,请先阅读第一部分——在本文中,我将假定你已熟悉其中的概念。

超越单一变量

我们的第一个机器学习模型——单变量线性回归模型——从技术上讲确实是机器学习,但平心而论,它并不那么令人印象深刻。那么我们费那么大劲的意义何在?是为了给更复杂的模型奠定基础,比如我们今天要探讨的多变量线性模型。

多元线性模型

多元线性模型(MLM)本质上与其单变量版本相同,区别在于它允许多个输入变量共同影响输出结果。在上一部分中,我们仅尝试基于房屋面积来预测房价。我们都知道,阿肯色州一套 5000 的房屋比圣何塞的同面积房屋便宜得多,但我们的模型却会为两者输出相同的价格。

这正是我们意识到需要让模型能够混合处理不同类别的变量的时候——每个变量都可能对预测产生不同的影响。

问题陈述

首先定义我们的问题:

我们获得以下输入变量及其对应的房屋价格 (以千美元计)的数据:

:建筑面积(以 计)

:地块宽度(以英尺计)

:地块深度(以英尺计)

:距海岸距离(以英里计)

给定一组新数据,预测房屋价格。

我们知道这些因素中的每一个都会对价格产生影响。有些(如建筑面积和距海岸距离)影响会非常大,而其他因素的影响则较小。

那么,我们如何将这些关系表达为一个模型呢?让我们尝试通过观察一些图表并进行估算来构建一个模型。

图表:房屋面积与价格关系

Loading...

可以看出,房屋面积与价格之间存在正相关关系。趋势线由 建模得出。

绘图:地块宽度与价格关系

Loading...

此处同样存在微弱的正相关关系。趋势线为

图表:地块深度与价格关系

Loading...

与地块宽度相比,此处呈现出更强的正相关性,表明购房者更看重纵深较长的地块形态。趋势线为

图表:离岸距离与价格关系

Loading...

正如我们所料,离海岸越远,价格越低。我们还观察到,随着离岸距离的增加,这种影响逐渐减弱。趋势线为

模型的初步猜测

趋势线告诉我们每个特征与房价之间的大致关系。我们可以将这些趋势线组合起来,得到一个近似的模型。让我们尝试将所有趋势线相加,看看是否能得到一个合理的模型。

直观上,这意味着什么?我们有一个 650 万美元的基准值(偏差),这代表了一个面积为 0 平方英尺、没有地块、位于海边的房子的理论价值。等等……这听起来不太对。我们或许应该取这些 y 轴截距的平均值,因为每个图都与其他图有重叠。新模型如下:

好吧……现在是 160 万美元。但这也许还算合理,因为它正好在海边。模型的其他部分似乎说得通。它表明:

每当以下项增加 1 个单位时 房价变化
建筑面积 340 美元
地块宽度 3300 美元
地块深度 5840 美元
离海岸距离 -8070 美元

虽然合理,但这些只是基于粗略计算的假设。 要真正解决我们的问题,我们需要找到能根据这 4 个特征预测房价的最优 MLM。用数学术语来说,我们需要找到权重 和偏差 ,使得

能以最小误差预测房价。

代码中的多元线性模型

与第一部分类似,我们先定义模型。由于现在有多个权重,我们需要使用一个 w 数组。

struct mlinear_model {
	int num_weights;
	double *w;
	double b;
};

你可能已经提前想到,如果我们将 w 表示为一个向量,x 也表示为向量,那么模型的输出就可以表示为 的点积,再加上

嗯,这样清晰多了,对吧?不过,这种简洁不会持续太久……

这里我再进一步,将我们的 维向量表示为一个 的矩阵。

你很快就会明白这样做的好处。我们先定义 matrix 结构体。

// matrix.h

// 这样以后我们可以方便地将 `double` 替换为其他浮点类型
typedef double mfloat;

typedef struct {
	// 行主序
	mfloat *buf;
	int rows;
	int cols;
} matrix;

现在,模型看起来是这样的:

// main.c
struct mlinear_model {
	matrix w;
	mfloat b;
};

注意,num_weights 现在已包含在矩阵信息中。

由于我们改用矩阵,我们需要使用矩阵点积而不是向量点积。我在代码中添加了一些注释,以防你需要复习线性代数。

// matrix.h

// 获取矩阵中的元素
mfloat
matrix_get(matrix m, int row, int col)
{
    return m.buf[row * m.cols + col];
}
// 设置矩阵中的元素
void
matrix_set(matrix m, int row, int col, mfloat val)
{
    m.buf[row * m.cols + col] = val;
}

// out = m1 点乘 m2
void
matrix_dot(matrix out, const matrix m1, const matrix m2)
{
	// 遍历第一个矩阵的第 i 行
    for (int row = 0; row < m1.rows; row++) {
		// 遍历第二个矩阵的第 j 列
        for (int col = 0; col < m2.cols; col++) {
			// 执行向量点积,并将结果放入 out[i][j]
            mfloat sum = 0.0;
			// m1.cols == m2.rows,因此 k 遍历所有对应元素
            for (int k = 0; k < m1.cols; k++) {
                mfloat x1 = matrix_get(m1, row, k);
                mfloat x2 = matrix_get(m2, k, col);
                sum += x1 * x2;
            }
            matrix_set(out, row, col, sum);
        }
    }
}

首先,我们注意矩阵 点积的以下性质:

  1. 仅当 X.cols == W.rows 时才能计算。
  2. 结果矩阵的维度为 (X.rows, W.cols)。

在本例中,

因此, 需要是一个行向量,而 需要是一个列向量,这样才能复现向量点积的行为。

你可以将矩阵点积理解为:取出第二个矩阵的每一列,将其转置为行,然后与第一个矩阵的对应行进行点积。

查看代码,验证上述乘法是否确实会返回所声称的结果。

现在,我们有了编写预测函数的工具。

// mlr.c

// x 是一个行向量,即 (1 x n) 矩阵
mfloat predict(struct mlinear_model model, const matrix x) {
	// (1 x n) . (n x 1) => (1 x 1) 矩阵,即一个数字
	mfloat result[1][1] = {0.0};
	matrix tmp = {.buf = result, .rows = 1, .cols = 1};
	// 将 tmp 设为点积结果
	matrix_dot(tmp, x, model.w);
	return tmp.buf[0] + model.b;
}

## 优化模型

我们有两个数组:输入数据,格式为 $(m \times 4)$ 矩阵。

$$
X = \begin{bmatrix} x_1^1 & x_2^1 & x_3^1 & x_4^1 \\\\ x_1^2 & x_2^2 & x_3^2 & x_4^2 \\\\ \vdots & \vdots & \vdots & \vdots \\\\ x_1^m & x_2^m & x_3^m & x_4^m \end{bmatrix}
$$

每一行代表一个样本,每列中的值是上述描述的特征。

以及对应的房屋价格 $Y$,单位为千美元。

$$
Y = \begin{bmatrix} y^1 \\\\ y^2 \\\\ y^3 \\\\ y^4 \\\\ \vdots \\\\ y^m \end{bmatrix}
$$
我们得到了 100 个数据样本,因此 $m = 100$。我将把它们存储在 `data.h` 中。

```c
// data.h

#define NUM_FEATURES 4
#define NUM_SAMPLES 100

static mfloat X_data[NUM_SAMPLES][NUM_FEATURES] = { /* 数据省略 */ };

static mfloat Y_data[NUM_SAMPLES] = { /* 数据省略 */ };

static matrix X = {.buf = (mfloat *)X_data,
                   .rows = NUM_SAMPLES,
                   .cols = NUM_FEATURES};

static matrix Y = {
    .buf = (mfloat *)Y_data, .rows = NUM_SAMPLES, .cols = 1};

现在,是优化的时候了!

这意味着我们需要找到参数 的模型,使得所有数据样本的误差最小化。但我们的误差是什么?实际上,我们可以沿用第 1 部分中的定义,因为我们仍然是在比较两个数字

这是平方差之和均值,或者说预期 值与实际 值之间的平均平方差。这个值越接近 ,我们的模型就越好。让我们用 来重写

这里, 指的是 的第 行。 是第 个样本与权重的点积。加上偏置 后,我们得到 ,即预测值。

在代码中:

// 返回一个表示矩阵 m 第 i 行的单行矩阵
matrix matrix_get_row(matrix m, int i) {
  return (matrix){
      .buf = m.buf + i * m.cols, .rows = 1, .cols = m.cols};
}

mfloat cost(struct mlinear_model *model, matrix X,
                    matrix Y) {
  mfloat cost = 0.0;
  for (int i = 0; i < X.rows; i++) {
    matrix x_i = matrix_get_row(X, i);
    mfloat f_wb = predict(model, x_i);
    mfloat diff = matrix_get(Y, 0, i) - f_wb;
    cost += diff * diff;
  }
  return cost / (2.0 * X.rows);
}

猜测的效果如何?

让我们回到之前的猜测模型,看看它的表现。

int main() {
  matrix W = matrix_new(X.cols, 1);
  struct mlinear_model model = {.W = W, .b = 0.0};
  printf("零模型的成本: %f\n", cost(&model, X, Y));
  matrix_set(W, 0, 0, 0.34);
  matrix_set(W, 1, 0, 3.3);
  matrix_set(W, 2, 0, 5.84);
  matrix_set(W, 3, 0, -8.07);
  model.b = 1634.5;
  printf("猜测模型的成本: %f\n",
         cost(&model, X, Y));
}

输出:

零模型的成本:        3340683.781483
猜测模型的成本: 2641267.466911

谁能想到呢!它比瞎猜稍微好一点。回到梯度下降。

让我们引入微积分。我们需要对 关于 求导,以便知道应该朝哪个方向,以多大的幅度移动权重和偏置来减少误差。

很简单,对吧?但对 就没那么简单了。由于 是一个矩阵, 中的每个权重对输出都有单独的影响。这意味着我们需要为每个权重 计算单独的导数。

那个 项来自于对 求导。对于一个给定的权重 和第 行,会有一个形如 的项。由于 是关于 的偏导数中的一个常数系数,根据链式法则,它会被提取出来。

既然我们已经用矩阵来思考了,让我们写出 关于整个矩阵 的导数

简化得很漂亮!你可能开始明白为什么我们在机器学习中喜欢矩阵了。这是思考批量计算的一种好方法。

现在,来看代码

struct grad {
	matrix dJ_dW;
	mfloat dJ_db;
};

/* 计算梯度并将结果写入 out。 */
void compute_gradient(struct grad *out,
                      struct mlinear_model *model, const matrix X,
                      const matrix Y) {
  int m = X.rows;  // 样本数量
  int n = X.cols;  // 特征数量

  // 使用 tmp 存储 X 的每一行
  matrix tmp = matrix_new(1, n);
  for (int i = 0; i < m; i++) {
    // tmp = X^(i)
    matrix curr_row = matrix_get_row(X, i);
    // y_hat = (X^(i) dot W) + b
    mfloat y_hat = predict(model, curr_row);
    // yi = y^(i)
    mfloat yi = matrix_get(Y, 0, i);
    // 括号内的项
    mfloat err = y_hat - yi;

    /*
     * 对于 dJ_dW,我们需要将误差
     * 乘以当前行,并将其加到累加和中
     */

    // tmp = X^(i) * (y_hat^(i) - y^(i))
    matrix_scalar_mul(tmp, curr_row, err);
    // dJ_dW += tmp
    matrix_ip_T_add(out->dJ_dW, tmp);

    // dJ_db += (y_hat^(i) - y^(i))
    out->dJ_db += err;
  }

  /*
   * 我将把 2/m 替换为 1/m,因为 2
   * 可以在下一步中合并到 alpha 里。
   */

  // dJ/db = (dJ/db) / m
  out->dJ_db /= m;
  // dJ/dW = (dJ/dW) / m
  matrix_scalar_ip_mul(out->dJ_dW, 1.0 / m);
  matrix_del(tmp);
}

注意:矩阵函数中的 ip 表示结果赋值给第一个参数。否则,结果赋值给作为第一个参数传入的缓冲区。

新的 matrix_* 函数相对简单,为了缩短文章,我将省略这里的代码。实际上,这些函数是由 #define 宏生成的,所以仓库里甚至没有完整的源代码。请期待第二部分,其中会详细介绍!

我鼓励你仔细阅读代码,并验证它是否与数学公式匹配。

现在我们有了梯度,实现梯度下降就很简单了。回顾一下第一部分的算法:

重复直到收敛:

  1. 设置初始权重和偏置值:
  2. 沿着与导数相反的方向移动变量, 步长为
  3. 转到步骤 2。

代码实现:

void gradient_descent(struct mlinear_model *model, const matrix X,
                      const matrix Y, const int num_iterations,
                      const mfloat alpha) {
  // 可重用的梯度缓冲区
  int n = X.cols, m = X.rows;
  matrix dJ_dW = matrix_new(n, 1);
  struct grad tmp_grad = {.dJ_dW = dJ_dW, .dJ_db = 0.0};

  for (int i = 0; i < num_iterations; i++) {
    // 记录进度
    if (i % (num_iterations >> 4) == 0) {
      printf("\t第 %d 次迭代的成本: %f\n", i,
             compute_cost(model, X, Y));
    }
    // tmp_grad = 模型当前的梯度
    compute_gradient(&tmp_grad, model, X, Y);
    // dJ/dW *= -alpha
    matrix_scalar_ip_mul(tmp_grad.dJ_dW, -alpha);
    // W += dJ/dW
    matrix_ip_add(model->W, tmp_grad.dJ_dW);
    // b += -alpha * dJ/db
    model->b += -alpha * tmp_grad.dJ_db;
  }
  matrix_del(dJ_dW);
}

就这样!现在我们在数据上运行程序,看看成本如何改善:

int main() {
  // 超参数
  const int num_iterations = 1e7;
  const mfloat alpha = 1e-8;

  int n = X.cols, m = X.rows;
  matrix W = matrix_new(n, 1);
  struct mlinear_model model = {.W = W, .b = 0.0};

  printf("初始成本: %f\n", compute_cost(&model, X, Y));
  gradient_descent(&model, X, Y, num_iterations, alpha);
  printf("最终成本: %f\n", compute_cost(&model, X, Y));
  printf("模型参数:\n");
  matrix_print(model.W);
  printf(" b=%f\n", model.b);
}

输出:

初始成本: 3340683.781483
        第 0 次迭代的成本: 3340683.781483
        第 625000 次迭代的成本: 161369.409722
        第 1250000 次迭代的成本: 161332.630253
        第 1875000 次迭代的成本: 161315.326515
        第 2500000 次迭代的成本: 161298.041198
        第 3125000 次迭代的成本: 161280.768143
        第 3750000 次迭代的成本: 161263.507342
        第 4375000 次迭代的成本: 161246.258784
        第 5000000 次迭代的成本: 161229.022462
        第 5625000 次迭代的成本: 161211.798366
        第 6250000 次迭代的成本: 161194.586488
        第 6875000 次迭代的成本: 161177.386819
        第 7500000 次迭代的成本: 161160.199351
        第 8125000 次迭代的成本: 161143.024075
        第 8750000 次迭代的成本: 161125.860983
        第 9375000 次迭代的成本: 161108.710065
最终成本: 161091.571313
模型参数:
W=[ 2.1185e-01 ]
[ 5.8627e+00 ]
[ 4.5904e+00 ]
[ -1.1702e+01 ]

 b=5.310395

太棒了!我们可以看到每次迭代成本都在下降。最终的模型看起来像

优化优化器

你可能已经注意到,即使梯度下降进行到第一千万次迭代,成本仍在下降。为什么找到最优参数如此困难,尤其是在我们只有4个权重需要优化的情况下?

很大程度上是由于特征分布范围的差异。看看这个箱线图:

Loading...

我们可以看到,房屋面积的权重需要移动的幅度远大于其他变量的权重。然而,当我们更新每个权重时,步长 是相同的。这意味着在其他3个特征权重收敛后,我们将不得不等待很长时间,房屋面积的权重才能收敛。

一种解决方案是选择与分布范围成比例的 值,并用对应的 移动每个权重。但这意味着需要存储另一个 alpha 值数组,对于大型模型来说成本很高。

更好的解决方案是实际修改我们的输入数据,使它们具有相似的分布范围。我们可以通过将它们拟合到正态分布来实现。首先,我们计算每个特征的均值。

然后,我们计算每个特征的分布范围,即标准差。

最后,我们将其归一化

将每列减去其均值可以将数据中心化,使得分布在两侧的扩散相似。除以 可以使每列的分布范围大致相同。让我们用矩阵运算来表示这个过程:

// 对输入数据 `X` 进行 Z 分数归一化
void z_score_normalize(matrix X) {
  int n = X.cols, m = X.rows;
  // 用于存储 mu 的缓冲区
  matrix mean = matrix_new(1, n);
  // 用于存储 sigma 的缓冲区
  matrix stdev = matrix_new(1, n);

  // 计算 mu
  for (int i = 0; i < NUM_SAMPLES; i++) {
    matrix_ip_add(mean, matrix_get_row(X, i));
  }
  matrix_scalar_ip_mul(mean, 1.0 / NUM_SAMPLES);

  // 计算 sigma
  matrix buf = matrix_new(1, n);
  for (int i = 0; i < NUM_SAMPLES; i++) {
    matrix row = matrix_get_row(X, i);
    matrix_sub(buf, mean, row);
    matrix_ip_square(buf);
    matrix_ip_add(stdev, buf);
  }
  matrix_ip_sqrt(stdev);

  // 计算 Z
  for (int i = 0; i < NUM_SAMPLES; i++) {
    matrix row = matrix_get_row(X, i);
    matrix_ip_sub(row, mean);
    matrix_ip_div(row, stdev);
  }
}

归一化后的数据看起来如何?

Loading...

我们可以看到,现在所有特征的分布范围都在同一尺度上。这如何改善我们的性能呢?

int main() {
  int n = X.cols, m = X.rows;
  z_score_normalize(X);
  matrix W = matrix_new(n, 1);
  struct mlinear_model model = {.W = W, .b = 0.0};
  printf("初始成本: %f\n", compute_cost(&model, X, Y));
  const int num_iterations = 1e7;
  const mfloat alpha = 1e-8;
  gradient_descent(&model, X, Y, num_iterations, alpha);
  printf("最终成本: %f\n", compute_cost(&model, X, Y));

  printf("模型参数:\nW=");
  matrix_print(model.W);
  printf(" b=%f\n", model.b);
}

输出:

初始成本: 3340683.781483
        第 0 次迭代的成本: 3340683.781483
        第 625000 次迭代的成本: 3305547.933747
        第 1250000 次迭代的成本: 3270852.031317
        第 1875000 次迭代的成本: 3236590.554988
        第 2500000 次迭代的成本: 3202758.054240
        第 3125000 次迭代的成本: 3169349.146943
        第 3750000 次迭代的成本: 3136358.518493
        第 4375000 次迭代的成本: 3103780.920969
        第 5000000 次迭代的成本: 3071611.172296
        第 5625000 次迭代的成本: 3039844.155415
        第 6250000 次迭代的成本: 3008474.817473
        第 6875000 次迭代的成本: 2977498.169013
        第 7500000 次迭代的成本: 2946909.283181
        第 8125000 次迭代的成本: 2916703.294939
        第 8750000 次迭代的成本: 2886875.400289
        第 9375000 次迭代的成本: 2857420.855511
最终成本: 2828334.976402
模型参数:
W=[ 8.3053e+00 ]
[ 2.4385e+00 ]
[ 6.0186e+00 ]
[ -2.2844e+00 ]

 b=227.136534

它甚至更慢了!?但是等等——看看当我们将 改为 时会发生什么,这原本会导致我们的旧模型发散。

初始成本: 3340683.781483
        第 0 次迭代的成本: 3340683.781483
        第 625000 次迭代的成本: 136947.544865
        第 1250000 次迭代的成本: 136947.544865
        第 1875000 次迭代的成本: 136947.544865
        第 2500000 次迭代的成本: 136947.544865
        第 3125000 次迭代的成本: 136947.544865
        第 3750000 次迭代的成本: 136947.544865
        第 4375000 次迭代的成本: 136947.544865
        第 5000000 次迭代的成本: 136947.544865
        第 5625000 次迭代的成本: 136947.544865
        第 6250000 次迭代的成本: 136947.544865
        第 6875000 次迭代的成本: 136947.544865
        第 7500000 次迭代的成本: 136947.544865
        第 8125000 次迭代的成本: 136947.544865
        第 8750000 次迭代的成本: 136947.544865
        第 9375000 次迭代的成本: 136947.544865
最终成本: 136947.544865
模型参数:
W=[ 8.1088e+03 ]
[ 8.1764e+02 ]
[ 6.6840e+02 ]
[ -3.6835e+03 ]

 b=2364.131545

哇!在第一次日志记录内,它就收敛到了最小值!实际上,使用 ,我们只需要大约 次迭代就能收敛,而不是 次!因此,我们可以看到,使用 Z 分数归一化可以将梯度下降的速度提高几个数量级,尤其是在特征具有不同分布范围的情况下。

如果你注意到新的 与原始值相差甚远, 那是因为我们从根本上改变了输入数据。我们现在是基于特征样本的偏差大小进行建模,而不是其绝对值。

总结

以上就是多元线性回归的全部内容。这无疑是本系列迄今为止最难的一篇文章,所以不必畏惧未来!所有代码都可以在这里找到。我鼓励你运行并修改其行为。

如果你有任何问题或意见,欢迎留言或给我发邮件。

✦ 本文的构思、研究、撰写和编辑均未使用大语言模型。