多变量线性回归模型

在许多线性回归问题中，训练集的数据特征不止一个，因此我们需要考虑一个包含多个输入变量的模型算法，即多变量线性回归模型。在继续进行之前，我们还是规定一下记号

$n$ ：训练样本的特征数
$\boldsymbol{x^{(i)}}$ ：训练集中索引为 $i$ 的输入数据
$x_j^{(i)}$ ： $\boldsymbol{x^{(i)}}$ 的第 $j$ 个特征量

在单变量线性回归模型中，假设函数的形式为

$h_{\theta}(x)=\theta_0+\theta_1x$

当有 $n$ 个变量时，根据线性函数的定义，我们规定假设函数为

$h_{\theta}(x)=\theta_0+\theta_1x_1+\theta_2x_2+\cdots+\theta_nx_n$

为了让上式看起来更加美观，我们规定 $x_0=1$ ，因此可以改写上式为

$h_{\theta}(\boldsymbol{x})=\boldsymbol{\theta}^{T}\boldsymbol{x}$

其中涉及的向量解释为

$\boldsymbol{x}=\begin{bmatrix}x_0\\x_1\\x_2\\\vdots\\x_n\end{bmatrix}\in \mathbb{R}^{n+1},\quad \boldsymbol{\theta}=\begin{bmatrix}\theta_0\\\theta_1\\\theta_2\\\vdots\\\theta_n\end{bmatrix}\in \mathbb{R}^{n+1}$

类比单变量线性回归模型，多变量线性回归模型的代价函数可以定义如下

$J(\boldsymbol{\theta})=\frac{1}{2m}\sum_{i=1}^{m}\Big(h_{\theta}(\boldsymbol{x^{(i)}})-y^{(i)}\Big)^2$

多元梯度下降法

核心公式

多变量线性回归模型的梯度下降算法可以类比单变量的情形，有

$\theta_{j}:=\theta_j-\alpha\frac{\partial}{\partial \theta_j}J(\boldsymbol{\theta}),\quad(j=0,1,\cdots,n)$

展开代价函数并求偏导，有更简洁的形式

$\theta_j:=\theta_j-\alpha\cdot\dfrac{1}{m}\displaystyle\sum_{i=1}^{m}\Big(h_{\theta}(\boldsymbol{x}^{(i)})-y^{(i)}\Big)x_j^{(i)},\quad(j=0,1,\cdots,n)$

上式是多元梯度下降法的核心迭代公式。特别地，当 $j=0$ 时，由于在上文中已经规定了 $x_0=1$ ，故此迭代公式仍然成立。仍然强调：这些变量的更新是同时进行的。

特征缩放

我们假设在房屋价格预测问题中，存在如下两个变量(特征量)

房屋面积 $x_1$ ：尺寸在0-2000平方英尺之间
卧室数量 $x_2$ ：数量在1-5个之间

如果画出这两个变量对应的参数 $\theta_1,\theta_2$ 和代价函数 $J$ 的等值线图，会发现由于特征量的量级不同，整个图像会十分扁平化。这会导致梯度下降法的效率降低。

为了解决这一问题，可以使用特征缩放(Feature Scaling)，通过归一化使得特征量的量级近似于同一水平。这在直观上表示为对应等值线图像“更圆了”。

[{"url":"/img/mlearning/m3f1.png","alt":"特征缩放","title":""}]

更一般的，我们在特征缩放时，总是想让每一个特征量 $x_i$ 大致处于 $[-1,1]$ 区间内。这一步可以通过归一化操作实现。

有时候我们也会进行均值归一化(Mean Normarlization)。对于特征量 $x_i$ ，如果其取值的数学期望(均值)为 $\mu_i$ ，那么我们就用 $x_i-\mu_i$ 替换 $x_i$ ，确保新特征量的均值大致为0，然后再进行数值的归一化。

不要对 $x_0=1$ 使用均值归一化。

例如，对于一个在 $[100,200]$ 均匀取值的特征量 $x_1$ ，可以进行如下操作

$x^{(i)}_1:=\frac{x^{(i)}_1-150}{200-100}$

这样就能大致保证新特征量符合特征缩放的要求。有时候我们不需要严格按照要求进行缩放，只要所有的特征量量级相近，梯度下降法的效率就不会受到大的影响。

学习率 $\alpha$

在梯度下降的过程中， $J(\boldsymbol{\theta})$ 理应是随着迭代次数的变多而不断减小的，直到逼近最小值点或局部最小值点。但是，随着迭代次数到达某一阈值后，代价函数 $J(\boldsymbol{\theta})$ 可能因为局部收敛而下降的越来越慢。此时无休止的进行下去是没有太大意义的。

[{"url":"/img/mlearning/m3f2.png","alt":"代价函数的收敛趋势","title":""}]

我们不仅可以从 $\min J(\boldsymbol{\theta})$ 和迭代次数的图像看出这种收敛趋势，也可以通过收敛判断程序辨认当前算法的收敛情况，比如在每一次迭代后判断

$\Delta J(\boldsymbol{\theta})\leqslant \varepsilon,\quad \varepsilon \text{ is a small number}$

在上一篇学习笔记中也提到——学习率 $\alpha$ 过大时，梯度下降算法可能会因为每一步的迭代步长过大从而不收敛。此时选用较小的学习率是一个较优选择。

可以证明，对于充分小的学习率 $\alpha$ ，代价函数可以随着迭代次数增加而不断下降。
然而，设置过小的 $\alpha$ 可能导致梯度下降算法极其缓慢，不满足实际需求。

因此，在最开始设置一个合理的较小学习率，基于此，根据算法的效率不断调整 $\alpha$ 的值是梯度下降算法中常用的一个优化策略。

创建新特征量

还是以我们最熟悉的房屋价格预测问题为例，我们考虑如下两个特征量

房屋占地的长度 $x_1$
房屋占地的宽度 $x_2$

按照传统的线性回归模型，我们可以构造假设函数为

$h_{\theta}(x_1,x_2)=\theta_0+\theta_1x_1+\theta_2x_2$

我们仔细推敲一下这个问题背后的逻辑——房屋购买者往往看重的是房屋的占地面积，而不是房屋占地的长和宽。事实上，对于相关问题的已有的数据，研究证明：使用 $x_1,x_2$ 作为变量的线性回归模型预测性极差。

因此，我们可以创造一个新的变量 $x=x_1x_2$ 表示房屋占地面积，而修改假设函数为

$h_{\theta}(x)=\theta_0+\theta_1x$

这种优化方法称作创建新特征量。

多项式回归

线性函数在拟合一些分布特殊的训练集数据时往往展现出欠佳的质量。我们可以考虑多项式函数

$h(x)=\theta_0+\theta_1x+\theta_2x^2+\cdots+\theta_nx^n$

这个假设函数是 $n$ 次多项式，因此使用此类假设函数的模型称为多项式回归模型(Polynomial Regression)。下图展现了使用二次函数和三次函数进行拟合的情况。

[{"url":"/img/mlearning/m3f3.png","alt":"多项式回归","title":""}]

以三次多项式回归为例，我们可以换个视角看待。设 $x=x_1$ 表示房屋的一维尺寸，那么 $x^2=x_2$ 可以看做房屋的占地， $x^3=x_3$ 可以看做房屋的体积。把这三个量看做三个特征量，那么多项式回归就又回到了线性回归模型上。

需要注意的是，上例中 $x_1,x_2,x_3$ 由于分别是 $x$ 的不同幂次，因此其量级存在着较大差异。如果使用梯度下降法，必须使用特征缩放保证算法的效率和质量。

使用次数为几的多项式进行回归拟合是一个很好的问题，有两个考虑维度

考虑函数图像和实际问题的符合性

二次多项式具有对称性，因此不会一直单调增加或单调减少
三次多项式具有奇性，整体变化趋势统一

避免过高次数多项式带来的过拟合(Overfitting)

在房屋价格问题中，常识告诉我们——房屋占地越大，价格总是会越大，因此采用诸如二次多项式等是不符合常理的。此时可以采用平方根多项式或三次多项式等。

正规方程

正规方程法

前面花了大篇笔墨介绍了梯度下降法，其核心在于不断地迭代。当我们已知代价函数 $J(\boldsymbol{\theta})$ 的形式时，使用代数方法暴力的求解函数的最小值也是一妙计，这就引出了正规方程法。正规方程法的正确性基于多元函数极值的必要条件，叙述如下。

对于一个多元函数 $J(\boldsymbol{\theta})=J(\theta_0,\theta_1,\cdots,\theta_n)$ ，其极值点处对于每一个变元的偏导数为0，即

$\frac{\partial}{\partial \theta_i}J(\boldsymbol{\theta}),\quad (i=0,1,\cdots,n)$

当训练集数据数量为 $m$ ，特征量数量为 $n$ 时，每一个训练数据可写成一个 $(n+1)\times 1$ 的列向量。令

$\boldsymbol{X}=(\boldsymbol{x}^{(1)},\boldsymbol{x}^{(2)},\cdots,\boldsymbol{x}^{(m)})^T=\begin{bmatrix}x_0^{(1)} & x_1^{(1)} & \cdots & x_n^{(1)} \\x_0^{(2)} & x_1^{(2)} & \cdots & x_n^{(2)} \\\vdots&\vdots&&\vdots\\x_0^{(m)} & x_1^{(m)} & \cdots & x_n^{(m)}\end{bmatrix}_{m\times(n+1)}$

再令输出向量

$\boldsymbol{y}=\begin{bmatrix}y^{(1)}\\y^{(2)}\\\vdots\\y^{(m)}\end{bmatrix}\in \mathbb{R}^{m}$

线性代数的最小二乘法知识告诉我们，使代价函数最小的参数 $\boldsymbol{\theta}$ 满足以下方程(若存在逆)

$\boldsymbol{\theta}=(\boldsymbol{X}^T\boldsymbol{X})^{-1}\boldsymbol{X}^T\boldsymbol{y}$

该方程称作正规方程(Normal Equation)，解该方程得到参数解的方法称作正规方程法。正规方程通常是由下述方程导出的

$\boldsymbol{X}^T\boldsymbol{X}\boldsymbol{\theta}=\boldsymbol{X}^T\boldsymbol{y}$

这个方程描述了在 $\boldsymbol{X}\boldsymbol{\theta}=\boldsymbol{y}$ 中，使得残差 $\|\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\theta}\|$ 最小的解 $\boldsymbol{\theta}$ 。然而，有时候 $\boldsymbol{X}^T\boldsymbol{X}$ 是不可逆的，这就导致我们无法直接使用 $\boldsymbol{\theta}=(\boldsymbol{X}^T\boldsymbol{X})^{-1}\boldsymbol{X}^T\boldsymbol{y}$ 求解。导致这种情况发生主要有两种可能

存在冗余特征量，他们之间表现出线性相关性(Linearly Dependent)
- 举例： $x_1$ 表示物体的长度(米)， $x_2$ 表示物体的长度(英尺)
拥有太多特征量，如 $m\leqslant n$

解决这一问题的关键手段是删除部分不必要的特征量，或使用正则化(Regularization)。

方法对比

梯度下降法	正规方程法
需要选择 $\alpha$	不需要选择 $\alpha$
需要经过多次迭代	不需要迭代
计算简单	计算 $\boldsymbol{X}^T\boldsymbol{X}$ 的逆(如果存在)，复杂度为 $O(n^3)$
即使 $n$ 很大时仍然有效	$n$ 很大时计算过程很慢

总结来说，当 $n$ 较小时可以选择正规方程法，当 $n$ 较大时可以选择梯度下降法。

代码实现

Python的numpy库可以直接实现矩阵运算，用法如下

import numpy as np

# 双变量线性回归
X = np.array([[1,1],[0,1],[-1,1],[2,1]])
Y = np.array([[-4],[9],[22],[-17]])

# 解正规方程(逆存在)
theta = np.linalg.inv(X.T@X)@X.T@Y
print(theta)