机器学习 - 线性回归

240

(本系列的 Jupyter Notebook 可以到我的码云下载。)

概述

本篇教程讲述线性回归,那么什么是线性回归呢?

首先举个例子,我们去市场买牛肉,一斤牛肉52块钱,两斤牛肉104块钱,三斤牛肉156块钱,以此类推。也是说牛肉的价格随着牛肉斤数的增加而有规律地增加,这种规律可以用下图表示:

牛肉斤数vs价格

可以看到上述规律可以用一条直线来表述,这就是一个线性模型。用 $x$ 表示牛肉斤数,用 $y$ 表示价格,就得到下面的方程:


$$ y = 52x \tag{1} $$

这个方程就叫做回归方程,$52$ 叫做回归系数,求解回归系数的过程叫做回归

线性回归首先假设自变量和因变量是线性关系,然后通过对现有样本进行回归,进而计算出回归系数以确定线性模型,最后使用这个模型对未知样本进行预测。

实际的回归方程要比上面例子中的复杂一些,它的一般形式如下:


$$ f(x_1, x_2, x_3, \dots, x_n ; \theta_0, \theta_1, \dots, \theta_n) = \theta_0 + \theta_1 x_1 + \theta_2 x_2 + \dots + \theta_n x_n \tag{2} $$

通常使用矩阵和向量的方式来简化书写。为了将上式写成矩阵形式,增加 $x_0 \equiv 1$($x_0$ 恒等于 $1$):


$$ f(x_0, x_1, x_2, x_3, \dots, x_n ; \theta_0, \theta_1, \dots, \theta_n) = x_0 \theta_0 + \theta_1 x_1 + \theta_2 x_2 + \dots + \theta_n x_n \tag{3} $$

令 $ \mathbf{x} = \left[ x_0 \ x_1 \ \dots \ x_n \right]^T $,$ \boldsymbol{\theta} = \left[ \theta_0 \ \theta_1 \ \dots \ \theta_n \right]^T $,就得到回归方程的矩阵形式:


$$ f(\mathbf{x}; \boldsymbol{\theta}) = \boldsymbol{\theta}^T \mathbf{x} = \mathbf{x}^T \boldsymbol{\theta} \tag{4} $$

上式中 $ \boldsymbol{\theta}^T \mathbf{x} = \mathbf{x}^T \boldsymbol{\theta} $ 是因为 $ \boldsymbol{\theta}^T \mathbf{x} $ 的结果是标量((1x1)的矩阵),而标量的转置就是它本身。因此


$$ \mathbf{x}^T \boldsymbol{\theta} = (\boldsymbol{\theta}^T \mathbf{x})^T = scalar^T = scalar \tag{5} $$

我们再来看一个例子,比如有如下的数据集:

所有样本点

容易看出,这些数据是线性的,因此可以使用线性回归来拟合。由于每个样本点只有一个自变量和一个因变量,因此可以用如下回归方程来表示模型:


$$ y = \theta_0 + \theta_1 x \tag{6} $$

在上图中画出几条通过样本集的直线:

通过样本点的直线

可以看到通过样本集的直线有无数条,每一条都由 $\theta_0$ 和 $\theta_1$ 确定。我们的任务就是找到那条最符合数据集的直线,也就是求得最优的 $\theta_0$ 和 $\theta_1$。

上面的讨论同样适用于 n 个自变量的情况,也就是说最终要计算的是最优的 $\boldsymbol{\theta}$ 值。如何计算最优的 $\boldsymbol{\theta}$ 呢?本篇介绍两种方法:最小二乘法、最大似然法。先来看最小二乘法。

最小二乘法

最小二乘法是通过最小化模型预测值(估计值,拟合值)与真实值(观察值)的差距(残差)求得最优 $\boldsymbol{\theta}$ 的一种方法。假设用 $y_i$ 表示一个样本 $i$ 的真实值,$f(\mathbf{x}_i;\boldsymbol{\theta})$ 表示该样本的预测值,那么它们的差距 $\mathcal{L}_i$ 可以表示为欧氏距离:


$$ \mathcal{L}_i = \mathcal{L}_i(y_i,f(\mathbf{x}_i;\boldsymbol{\theta})) = (y_i - f(\mathbf{x}_i;\boldsymbol{\theta}))^2 \tag{7} $$

对于整个样本集,则计算所有样本的真实值与预测值差距(残差平方和)的均值:


$$ \mathcal{L} = \frac{1}{m} \sum_{i=1}^m \mathcal{L}_{i} = \frac{1}{m} \sum_{i=1}^m \mathcal{L}_i(y_i,f(\mathbf{x}_i;\boldsymbol{\theta})) = \frac{1}{m} \sum_{i=1}^m (y_i - f(\mathbf{x}_i;\boldsymbol{\theta}))^2 \tag{8} $$

上式中,$m$ 表示样本的数量。在机器学习中,很多算法都是先找到真实值与预测值的差距,然后通过最小化这个差距来调整模型参数。“真实值与预测值的差距” 在机器学习中有一个专用术语叫做损失函数

最小二乘法就是求 $(8)$ 式取最小值时 $\boldsymbol{\theta}$ 的值(即 $\boldsymbol{\theta}$ 的最优值,用 $\widehat{\boldsymbol{\theta}}$ 表示):


$$ \widehat{\boldsymbol{\theta}} = \underset{\boldsymbol{\theta}}{\arg\min}\ \frac{1}{m} \sum_{i=1}^m \mathcal{L}_i(y_i,f(\mathbf{x}_i;\boldsymbol{\theta})) \tag{9} $$

为了解决最小二乘问题,首先需要复习一些基础知识。

基础知识

标量对向量的导数

所谓标量对向量求导,就是用标量对向量中的每一个元素求导,求得的结果为一个向量:


$$ \frac{\partial{y}}{\partial{\mathbf{x}}} = \begin{bmatrix} \dfrac{\partial{y}}{\partial{x_1}} & \dfrac{\partial{y}}{\partial{x_2}} & \cdots & \dfrac{\partial{y}}{\partial{x_n}} \end{bmatrix} \tag{10} $$

当然也可以为下面的向量:


$$ \frac{\partial{y}}{\partial{\mathbf{x}}} = \begin{bmatrix} \dfrac{\partial{y}}{\partial{x_1}} \\ \dfrac{\partial{y}}{\partial{x_2}} \\ \vdots \\ \dfrac{\partial{y}}{\partial{x_n}} \end{bmatrix} \tag{11} $$

在实际应用中,上面两种形式都可以表示标量对向量的导数,不过为了区分它们,将 $(10)$ 式这样的标记方式叫做分子分布,将 $(11)$ 式这样的标记方式叫做分母分布。我们可以把分子分布理解为把分子固定,而将分母转置。分母分布则可以理解为把分母固定,而将分子转置。说白了,就是谁的分布谁做主

本篇接下来使用的都是分母分布的形式,因此我们将不再讨论分子分布了。

向量对向量的导数

上面标量对向量的导数可以很容易扩展到向量对向量的导数:


$$ \frac{\partial{\mathbf{y}}}{\partial{\mathbf{x}}} = \begin{bmatrix} \dfrac{\partial{y_1}}{\partial{x_1}} & \dfrac{\partial{y_2}}{\partial{x_1}} & \cdots & \dfrac{\partial{y_n}}{\partial{x_1}} \\ \dfrac{\partial{y_1}}{\partial{x_2}} & \dfrac{\partial{y_2}}{\partial{x_2}} & \cdots & \dfrac{\partial{y_n}}{\partial{x_2}} \\ \vdots & \vdots & \ddots & \vdots \\ \dfrac{\partial{y_1}}{\partial{x_n}} & \dfrac{\partial{y_2}}{\partial{x_n}} & \cdots & \dfrac{\partial{y_n}}{\partial{x_n}} \end{bmatrix} \tag{12} $$

标量对矩阵的导数

标量对矩阵的导数即为标量对矩阵的每一个元素求导:


$$ \dfrac{\partial f(\mathbf{A})}{\partial \mathbf{A}} = \begin{bmatrix} \dfrac{\partial f(\mathbf{A})}{\partial a_{11}} & \dfrac{\partial f(\mathbf{A})}{\partial a_{12}} & \cdots & \dfrac{\partial f(\mathbf{A})}{\partial a_{1n}} \\ \dfrac{\partial f(\mathbf{A})}{\partial a_{21}} & \dfrac{\partial f(\mathbf{A})}{\partial a_{22}} & \cdots & \dfrac{\partial f(\mathbf{A})}{\partial a_{2n}} \\ \vdots & \vdots & \ddots & \vdots \\ \dfrac{\partial f(\mathbf{A})}{\partial a_{m1}} & \dfrac{\partial f(\mathbf{A})}{\partial a_{m2}} & \cdots & \dfrac{\partial f(\mathbf{A})}{\partial a_{mn}} \end{bmatrix} \tag{13} $$

公式

假设 $\mathbf{A}$ 是一个与 $\mathbf{x}$ 无关的常数矩阵,则有下面两个向量对向量的导数公式:

公式(一)


$$ \frac{\partial{\mathbf{A} \mathbf{x}}}{\partial{\mathbf{x}}} = \mathbf{A}^T \tag{14} $$

证明:

假设 $\mathbf{A}$ 是一个 $m \times n$ 的常数矩阵,$\mathbf{x}$ 是一个 $n$ 维向量(即 $n \times 1$ 的矩阵),则:


\begin{align*} \frac{\partial{\mathbf{A} \mathbf{x}}}{\partial{\mathbf{x}}} &= \frac{\partial}{\partial{\mathbf{x}}} \left( \begin{bmatrix} a_{11} & a_{12} & \cdots & a_{1n} \\ a_{21} & a_{22} & \cdots & a_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{m1} & a_{m2} & \cdots & a_{mn} \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ \vdots \\ x_n \end{bmatrix} \right) \\ &= \frac{\partial}{\partial{\mathbf{x}}} \begin{bmatrix} \sum\limits_{i=1}^n a_{1i} x_{i} \\ \sum\limits_{i=1}^n a_{2i} x_{i} \\ \vdots \\ \sum\limits_{i=1}^n a_{mi} x_{i} \end{bmatrix} \\ &= \begin{bmatrix} \dfrac{\partial{\sum\limits_{i=1}^n a_{1i} x_{i}}}{\partial{x_{1}}} & \dfrac{\partial{\sum\limits_{i=1}^n a_{2i} x_{i}}}{\partial{x_{1}}} & \cdots & \dfrac{\partial{\sum\limits_{i=1}^n a_{mi} x_{i}}}{\partial{x_{1}}} \\ \dfrac{\partial{\sum\limits_{i=1}^n a_{1i} x_{i}}}{\partial{x_{2}}} & \dfrac{\partial{\sum\limits_{i=1}^n a_{2i} x_{i}}}{\partial{x_{2}}} & \cdots & \dfrac{\partial{\sum\limits_{i=1}^n a_{mi} x_{i}}}{\partial{x_{2}}} \\ \vdots & \vdots & \ddots & \vdots \\ \dfrac{\partial{\sum\limits_{i=1}^n a_{1i} x_{i}}}{\partial{x_{n}}} & \dfrac{\partial{\sum\limits_{i=1}^n a_{2i} x_{i}}}{\partial{x_{n}}} & \cdots & \dfrac{\partial{\sum\limits_{i=1}^n a_{mi} x_{i}}}{\partial{x_{n}}} \end{bmatrix} \\ &= \begin{bmatrix} a_{11} & a_{21} & \cdots & a_{m1} \\ a_{12} & a_{22} & \cdots & a_{m2} \\ \vdots & \vdots & \ddots & \vdots \\ a_{1n} & a_{2n} & \cdots & a_{mn} \end{bmatrix} \\ &= \mathbf{A}^T \tag{15} \end{align*}

证毕。

公式(二)


$$ \frac{\partial{ \mathbf{x}^T \mathbf{A} \mathbf{x} }}{ \partial{ \mathbf{x} }} = 2 \mathbf{A} \mathbf{x} \tag{16} $$

证明:

假设 $\mathbf{x}$ 是一个 $n$ 维向量(即 $n \times 1$ 的矩阵),$\mathbf{A}$ 是一个 $n \times n$ 的对称矩阵,则:


\begin{align*} \frac{\partial{ \mathbf{x}^T \mathbf{A} \mathbf{x} }}{ \partial{ \mathbf{x} }} &= \frac{\partial}{\partial{\mathbf{x}}} \left( \begin{bmatrix} x_1 & x_2 & \cdots & x_n \end{bmatrix} \begin{bmatrix} a_{11} & a_{12} & \cdots & a_{1n} \\ a_{21} & a_{22} & \cdots & a_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{n1} & a_{n2} & \cdots & a_{nn} \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ \vdots \\ x_n \end{bmatrix} \right) \\ &= \frac{\partial}{\partial{\mathbf{x}}} \left( \begin{bmatrix} \sum\limits_{i=1}^n a_{i1} x_i & \sum\limits_{i=1}^n a_{i2} x_i & \cdots & \sum\limits_{i=1}^n a_{in} x_i \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ \vdots \\ x_n \end{bmatrix} \right) \\ &= \frac{\partial}{\partial{\mathbf{x}}} \left( \sum\limits_{i=1}^n a_{i1} x_{i} x_1 + \sum\limits_{i=1}^n a_{i2} x_{i} x_2 + \cdots + \sum\limits_{i=1}^n a_{in} x_{i} x_n \right) \\ &= \frac{\partial}{\partial{\mathbf{x}}} \left( \sum\limits_{i=1}^n \sum\limits_{j=1}^n a_{ij} x_i x_j \right) \\ &= \left[ \begin{matrix} \dfrac{\partial{\sum\limits_{i=1}^n \sum\limits_{j=1}^n a_{ij} x_i x_j}}{\partial{x_1}} \\ \dfrac{\partial{\sum\limits_{i=1}^n \sum\limits_{j=1}^n a_{ij} x_i x_j}}{\partial{x_2}} \\ \vdots \\ \dfrac{\partial{\sum\limits_{i=1}^n \sum\limits_{j=1}^n a_{ij} x_i x_j}}{\partial{x_n}} \end{matrix}\ \ \right] & ① \\ &= \begin{bmatrix} 2 \sum\limits_{i=1}^n a_{1i} x_i \\ 2 \sum\limits_{i=1}^n a_{2i} x_i \\ \vdots \\ 2 \sum\limits_{i=1}^n a_{ni} x_i \end{bmatrix} & ② \\ &= 2 \begin{bmatrix} a_{11} & a_{12} & \cdots & a_{1n} \\ a_{21} & a_{22} & \cdots & a_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{n1} & a_{n2} & \cdots & a_{nn} \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ \vdots \\ x_n \end{bmatrix} \\ &= 2 \mathbf{A} \mathbf{x} \tag{17} \end{align*}

$① \Rightarrow ②$ 的化简过程如下:

$①$ 中的每一项都是 $\dfrac{\partial{\sum\limits_{i=1}^n \sum\limits_{j=1}^n a_{ij} x_i x_j}}{\partial{x_k}}\ \ $ 的形式,其中 $1 \le k \le n$,仔细观察我们可以发现累加的每一项只有在 $i$ 和 $j$ 至少有一个为 $k$ 时导数才不为 $0$ ,也就是说只需要考虑三种情况:

  1. $i$ 和 $j$ 都为 $k$,即 $i = j = k$
  2. $i$ 为 $k$ 但 $j$ 不为 $k$,即 $i = k, j \ne k$
  3. $j$ 为 $k$ 但 $i$ 不为 $k$,即 $i \ne k, j = k$

将 $\sum\limits_{i=1}^n \sum\limits_{j=1}^n a_{ij} x_i x_j$ 按这三种情况展开:


\begin{align*} \sum\limits_{i=1}^n \sum\limits_{j=1}^n a_{ij} x_i x_j &= \underbrace{ a_{ij} x_{i} x_{j} }_{i=j=k} + \sum\limits_{\{j \mid 1 \le j \le n, j \ne k, i = k \}} a_{ij} x_i x_j + \sum\limits_{\{i \mid 1 \le i \le n, i \ne k, j = k \}} a_{ij} x_i x_j \\ &= a_{kk} x_k x_k + \underbrace{ \sum\limits_{\{j \mid 1 \le j \le n, j \ne k\}} a_{kj} x_k x_j + \sum\limits_{\{i \mid 1 \le i \le n, i \ne k \}} a_{ik} x_i x_k }_{对应元素相加} \\ &= a_{kk} x_k^2 + \sum\limits_{\{p \mid 1 \le p \le n, p \ne k\}} \left( a_{kp} x_k x_p + a_{pk} x_{p} x_{k} \right) \tag{18} \end{align*}

因为 $\mathbf{A}$ 是对称矩阵,所以 $a_{kp} = a_{pk}$,因此:


\begin{align*} \text{上式} &= a_{kk} x_k^2 + \sum\limits_{\{p \mid 1 \le p \le n, p \ne k \}} 2 a_{kp} x_k x_p \\ &= a_{kk} x_k^2 + 2 \sum\limits_{\{p \mid 1 \le p \le n, p \ne k \}} a_{kp} x_k x_p \tag{19} \end{align*}

现在对 $\dfrac{\partial{\sum\limits_{i=1}^n \sum\limits_{j=1}^n a_{ij} x_i x_j}}{\partial{x_k}}\ \ $ 求导:


\begin{align*} \frac{\partial{\sum\limits_{i=1}^n \sum\limits_{j=1}^n a_{ij} x_i x_j}}{\partial{x_k}} &= \frac{\partial}{\partial{x_k}} \left( a_{kk} x_k^2 + 2 \sum\limits_{\{p \mid 1 \le p \le n, p \ne k \}} a_{kp} x_k x_p \right) \\ &= 2 a_{kk} x_k + 2 \sum\limits_{\{p \mid 1 \le p \le n, p \ne k \}} a_{kp} x_p \\ &= 2\sum\limits_{p=1}^n a_{kp} x_p \tag{20} \end{align*}

这样就将 $①$ 化简成了 $②$ 。

证毕。

公式(三)


$$ \nabla_{\mathbf{A}^T} f(\mathbf{A}) = \left( \nabla_{\mathbf{A}} f(\mathbf{A}) \right)^T \tag{21} $$

证明:


\begin{align*} \nabla_{\mathbf{A}^T} f(\mathbf{A}) &= \begin{bmatrix} \dfrac{\partial f(\mathbf{A})}{\partial a_{11}} & \dfrac{\partial f(\mathbf{A})}{\partial a_{21}} & \cdots & \dfrac{\partial f(\mathbf{A})}{\partial a_{n1}} \\ \dfrac{\partial f(\mathbf{A})}{\partial a_{12}} & \dfrac{\partial f(\mathbf{A})}{\partial a_{22}} & \cdots & \dfrac{\partial f(\mathbf{A})}{\partial a_{n2}} \\ \vdots & \vdots & \ddots & \vdots \\ \dfrac{\partial f(\mathbf{A})}{\partial a_{1m}} & \dfrac{\partial f(\mathbf{A})}{\partial a_{2m}} & \cdots & \dfrac{\partial f(\mathbf{A})}{\partial a_{nm}} \end{bmatrix} \\ &= \left(\nabla_{\mathbf{A}} f(\mathbf{A})\right)^T \tag{22} \end{align*}

证毕。

正定矩阵与负定矩阵

如果对于任意非零向量 $\mathbf{z}$ 都有 $\mathbf{z}^T \mathbf{M} \mathbf{z} \gt 0$,那么 $\mathbf{M}$ 就是正定矩阵。

如果对于任意非零向量 $\mathbf{z}$ 都有 $\mathbf{z}^T \mathbf{M} \mathbf{z} \lt 0$,那么 $\mathbf{M}$ 就是负定矩阵。

海森矩阵

海森矩阵是由多元函数的二阶导数组成的 $n$ 阶方阵,用 $\mathbf{H}$ 表示,它的形态如下:


$$ \mathbf{H} = \begin{bmatrix} \dfrac{\partial^2 f}{\partial{x_1^2}} & \dfrac{\partial^2 f}{\partial{x_1} \partial{x_2}} & \cdots & \dfrac{\partial^2 f}{\partial{x_1} \partial{x_n}} \\ \dfrac{\partial^2 f}{\partial{x_2} \partial{x_1}} & \dfrac{\partial^2 f}{\partial{x_2^2}} & \cdots & \dfrac{\partial^2 f}{\partial{x_2} \partial{x_n}} \\ \vdots & \vdots & \ddots & \vdots \\ \dfrac{\partial^2 f}{\partial{x_n} \partial{x_1}} & \dfrac{\partial^2 f}{\partial{x_n} \partial{x_2}} & \cdots & \dfrac{\partial^2 f}{\partial{x_n^2}} \end{bmatrix} \tag{23} $$

海参矩阵有如下的性质:

  • 如果 $\mathbf{H}$ 为正定矩阵,那么函数取得最小值。
  • 如果 $\mathbf{H}$ 为负定矩阵,那么函数取得最大值。

矩阵的迹

矩阵的迹定义为 $n \times n$ 的方阵 $A$ 主对角线元素的和,用 $\textbf{tr}$ 表示,即:


$$ \textbf{tr}(\mathbf{A}) = \sum\limits_{i=1}^n a_{ii} = a_{11} + a_{22} + \cdots + a_{nn} \tag{24} $$

性质

性质(一)


$$ \textbf{tr}(c) = c \quad (c \in \mathbb{R}) \tag{25} $$

证明:将 $c$ 看成 $1 \times 1$ 的矩阵,立即得证。

性质(二)


$$ \textbf{tr}(\mathbf{A} + \mathbf{B}) = \textbf{tr}(\mathbf{A}) + \textbf{tr}(\mathbf{B}) \tag{26} $$

证明:$ \textbf{tr}(\mathbf{A + B}) = \sum\limits_{i=1}^n (a_{ii} + b_{ii}) = \sum\limits_{i=1}^n a_{ii} + \sum\limits_{i=1}^n b_{ii} = \textbf{tr}(\mathbf{A}) + \textbf{tr}(\mathbf{B}) $

性质(三)


$$ \textbf{tr}(c\mathbf{A}) = c \cdot \textbf{tr}(\mathbf{A}) \quad (c \in \mathbb{R}) \tag{27} $$

证明:$ \textbf{tr}(c\mathbf{A}) = \sum\limits_{i=1}^n c a_{ii} = c \sum\limits_{i=1}^n a_{ii} = c \cdot \textbf{tr}(\mathbf{A}) $

性质(四)


$$ \textbf{tr}(\mathbf{A}) = \textbf{tr}(\mathbf{A}^T) \tag{28} $$

证明:方阵的转置并不改变对角线上的元素。证毕。

性质(五)


$$ \textbf{tr}(\mathbf{AB}) = \textbf{tr}(\mathbf{BA}) \tag{29} $$

证明:$ \textbf{tr}(\mathbf{AB}) = \sum\limits_{i=1}^m(\mathbf{AB})_{ii} = \sum\limits_{i=1}^m \sum\limits_{j=1}^n a_{ij} b_{ji} = \sum\limits_{j=1}^n \sum\limits_{i=1}^m b_{ji} a_{ij} = \sum\limits_{j=1}^n (\mathbf{BA})_{jj} = \textbf{tr} (\mathbf{BA}) $

性质(六)(循环置换法则)


$$ \textbf{tr}(\mathbf{ABC}) = \textbf{tr}(\mathbf{CAB}) = \textbf{tr}(\mathbf{BCA}) \tag{30} $$

证明:根据性质(五)和矩阵乘法法则,立即得证。注意只能循环置换,而不能随意置换:$ \textbf{tr}(\mathbf{ABC}) \ne \textbf{tr}(\mathbf{ACB}) $

迹的导数

公式(一)


$$ \nabla_A \textbf{tr} (\mathbf{AB}) = \mathbf{B}^T \tag{31} $$

证明:

设 $A$ 为 $m \times n$ 的矩阵,$B$ 为 $n \times m$ 的矩阵,则有:


\begin{align*} \nabla_A \textbf{tr} (\mathbf{AB}) &= \nabla_A \sum\limits_{i=1}^m \sum\limits_{j=1}^n a_{ij} b_{ji} \\ &= \begin{bmatrix} \dfrac{\partial \sum\limits_{i=1}^m \sum\limits_{j=1}^n a_{ij} b_{ji}}{\partial a_{11}} & \dfrac{\partial \sum\limits_{i=1}^m \sum\limits_{j=1}^n a_{ij} b_{ji}}{\partial a_{12}} & \cdots & \dfrac{\partial \sum\limits_{i=1}^m \sum\limits_{j=1}^n a_{ij} b_{ji}}{\partial a_{1n}} \\ \dfrac{\partial \sum\limits_{i=1}^m \sum\limits_{j=1}^n a_{ij} b_{ji}}{\partial a_{21}} & \dfrac{\partial \sum\limits_{i=1}^m \sum\limits_{j=1}^n a_{ij} b_{ji}}{\partial a_{22}} & \cdots & \dfrac{\partial \sum\limits_{i=1}^m \sum\limits_{j=1}^n a_{ij} b_{ji}}{\partial a_{2n}} \\ \vdots & \vdots & \ddots & \vdots \\ \dfrac{\partial \sum\limits_{i=1}^m \sum\limits_{j=1}^n a_{ij} b_{ji}}{\partial a_{m1}} & \dfrac{\partial \sum\limits_{i=1}^m \sum\limits_{j=1}^n a_{ij} b_{ji}}{\partial a_{m2}} & \cdots & \dfrac{\partial \sum\limits_{i=1}^m \sum\limits_{j=1}^n a_{ij} b_{ji}}{\partial a_{mn}} \end{bmatrix} \\ &= \begin{bmatrix} b_{11} & b_{21} & \cdots & b_{n1} \\ b_{12} & b_{22} & \cdots & b_{n2} \\ \vdots & \vdots & \ddots & \vdots \\ b_{1m} & b_{2m} & \cdots & b_{nm} \end{bmatrix} \\ &= \mathbf{B}^T \tag{32} \end{align*}

证毕。

公式(二)


$$ \nabla_A \textbf{tr} (\mathbf{B} \mathbf{A}^T) = \mathbf{B} \tag{33} $$

证明:根据性质(四)和公式(一)立即得:


$$ \nabla_A \textbf{tr} ( \mathbf{B} \mathbf{A}^T) = \nabla_A \textbf{tr} (\left( \mathbf{B} \mathbf{A}^T\right)^T) = \nabla_A \textbf{tr} (\mathbf{A} \mathbf{B}^T) = \left(\mathbf{B}^T\right)^T = \mathbf{B} \tag{34} $$

证毕。

公式(三)


$$ \mathrm{d}(\textbf{tr}\mathbf{A}) = \textbf{tr}(\mathrm{d}\mathbf{A}) \tag{35} $$

证明:


\begin{align*} & \textbf{tr} \mathbf{A} = a_{11} + a_{22} + \cdots + a_{nn} \\ & \mathrm{d} (\textbf{tr} \mathbf{A}) = \mathrm{d} a_{11} + \mathrm{d} a_{22} + \cdots + \mathrm{d} a_{nn} \\ & \mathrm{d} \mathbf{A} = \begin{bmatrix} \mathrm{d} a_{11} & \cdots & \mathrm{d} a_{1n} \\ \vdots & \ddots & \vdots \\ \mathrm{d} a_{n1} & \cdots & \mathrm{d} a_{nn} \end{bmatrix} \\ & \textbf{tr}(\mathrm{d}\mathbf{A}) = \mathrm{d} a_{11} + \mathrm{d} a_{22} + \cdots + \mathrm{d} a_{nn} = \mathrm{d}(\textbf{tr}\mathbf{A}) \tag{36} \end{align*}

证毕。

公式(四)


$$ \nabla_\mathbf{A} \textbf{tr} (\mathbf{A} \mathbf{B} \mathbf{A}^T \mathbf{C}) = \mathbf{C}^T \mathbf{A} \mathbf{B}^T + \mathbf{C} \mathbf{A} \mathbf{B} \tag{37} $$

证明:

设 $u = u(\mathbf{A}) = \mathbf{AB} $,$v = v(\mathbf{A}^T) = \mathbf{A}^T \mathbf{C}$,则:


\begin{align*} \nabla_\mathbf{A} \textbf{tr} (\mathbf{A} \mathbf{B} \mathbf{A}^T \mathbf{C}) &= \nabla_{\mathbf{A}} \textbf{tr} (uv) \\ &= \dfrac{\partial \textbf{tr} (uv)}{\partial u} \dfrac{\partial u}{\partial \mathbf{A}} + \dfrac{\partial \textbf{tr} (uv)}{\partial v} \dfrac{\partial v}{\partial \mathbf{A}} \\ &= v^T \mathbf{B}^T + \left(\dfrac{\partial \textbf{tr} (uv)}{\partial v} \dfrac{\partial v}{\partial \mathbf{A}^T}\right)^T \\ &= \mathbf{C}^T \mathbf{A} \mathbf{B}^T + \left(u^T \mathbf{C}^T \right)^T \\ &= \mathbf{C}^T \mathbf{A} \mathbf{B}^T + \left( \mathbf{B}^T \mathbf{A}^T \mathbf{C}^T \right)^T \\ &= \mathbf{C}^T \mathbf{A} \mathbf{B}^T + \mathbf{C} \mathbf{A} \mathbf{B} \tag{38} \end{align*}

证毕。

求解

(如果你已经很熟悉如何求解最小二乘法或者想直接看结论,可以跳过这段,从这里继续阅读。)

下面求 $(8)$ 式的最小值。首先将 $(4)$ 式带入 $(8)$ 式,得到:


$$ \mathcal{L} = \frac{1}{m} \sum_{i=1}^m (\mathbf{y}_i - f(\mathbf{x}_i;\boldsymbol{\theta}))^2 = \frac{1}{m} \sum_{i=1}^m (y_i - \boldsymbol{\theta}^T \mathbf{x}_i)^2 \tag{39} $$

为了简化计算,令 $\mathbf{X} = \left[ \mathbf{x}_1^T \ \mathbf{x}_2^T \ \cdots \ \mathbf{x}_m^T \right]^T$, $ \mathbf{y} = \left[ y_1 \ y_2 \ \cdots \ y_m \right]^T $,则可将上式写成矩阵形式:


$$ \mathcal{L} = \frac{1}{m} (\mathbf{y} - \mathbf{X} \boldsymbol{\theta})^T (\mathbf{y} - \mathbf{X} \boldsymbol{\theta}) \tag{40} $$

将上面的式子展开,合并同类项:


\begin{align*} \mathcal{L} &= \frac{1}{m} (\mathbf{y} - \mathbf{X} \boldsymbol{\theta})^T (\mathbf{y} - \mathbf{X} \boldsymbol{\theta}) \\ &= \frac{1}{m} [\mathbf{y}^T - (\mathbf{X} \boldsymbol{\theta})^T] (\mathbf{y} - \mathbf{X} \boldsymbol{\theta}) \\ &= \frac{1}{m} (\mathbf{y}^T - \boldsymbol{\theta}^T \mathbf{X}^T) (\mathbf{y} - \mathbf{X} \boldsymbol{\theta}) \\ &= \frac{1}{m} (\mathbf{y}^T \mathbf{y} - \mathbf{y}^T \mathbf{X} \boldsymbol{\theta} - \boldsymbol{\theta}^T \mathbf{X}^T \mathbf{y} + \boldsymbol{\theta}^T \mathbf{X}^T \mathbf{X} \boldsymbol{\theta}) \\ &= \frac{1}{m} (\mathbf{y}^T \mathbf{y} - 2 \mathbf{y}^T \mathbf{X} \boldsymbol{\theta} + \boldsymbol{\theta}^T \mathbf{X}^T \mathbf{X} \boldsymbol{\theta}) \tag{41} \end{align*}

最后一步 $\mathbf{y}^T \mathbf{X} \boldsymbol{\theta}$ 与 $\boldsymbol{\theta}^T \mathbf{X}^T \mathbf{y}$ 可以合并的原因同 $(5)$,也就是说:互为转置的两个矩阵表达式,如果它们的结果为标量,那么它们相等。

好了,现在已经把 $\mathcal{L}$ 化简得不能再化简了,那么如何求 $\mathcal{L}$ 的最小值呢?现在请回顾一下高等数学的知识,在高等数学中,要求一个函数的最小值点,有两个步骤:

  1. 对函数求一阶导数并令其等于 $0$,求得的点就是可能的极值点。
  2. 对函数求二阶导数,如果恒大于 $0$,则函数图形为向上凹的,1 中求得的点即是最小值点。

比如 $ y = x^2 $,一阶导数为 $y’ = 2x$,令它等于 $0$ 解得 $x = 0$,那么 $0$ 就是一个可能的极值点,而二阶导数 $y’’ = 2$ 恒大于 $0$,那么函数 $ y = x^2 $ 的图像就是向上凹的,$x = 0$ 即为最小值点。

现在来求 $\mathcal{L}$ 的一阶导数:


\begin{align*} \frac{\partial{\mathcal{L}}}{\partial{\boldsymbol{\theta}}} &= \frac{\partial}{\partial{\boldsymbol{\theta}}} \left[\frac{1}{m} (\mathbf{y}^T \mathbf{y} - 2 \mathbf{y}^T \mathbf{X} \boldsymbol{\theta} + \boldsymbol{\theta}^T \mathbf{X}^T \mathbf{X} \boldsymbol{\theta}) \right] \\ &= \frac{\partial}{\partial{\boldsymbol{\theta}}} \left(- \frac{2}{m} \mathbf{y}^T \mathbf{X} \boldsymbol{\theta} + \frac{1}{m} \boldsymbol{\theta}^T \mathbf{X}^T \mathbf{X} \boldsymbol{\theta} \right) \\ &= - \frac{2}{m} \frac{\partial{ \mathbf{y}^T \mathbf{X} \boldsymbol{\theta} }}{\partial{ \boldsymbol{\theta} }} + \frac{1}{m} \frac{\partial{ \boldsymbol{\theta}^T \mathbf{X}^T \mathbf{X} \boldsymbol{\theta} }}{ \partial{ \boldsymbol{\theta} }} \\ &= - \frac{2}{m} \left( \mathbf{y}^T \mathbf{X} \right)^T + \frac{1}{m} 2 \mathbf{X}^T \mathbf{X} \boldsymbol{\theta} \\ &= - \frac{2}{m} \mathbf{X}^T \mathbf{y} + \frac{2}{m} \mathbf{X}^T \mathbf{X} \boldsymbol{\theta} \tag{42} \end{align*}

上式也可以通过矩阵的迹求得:


\begin{align*} \frac{\partial{\mathcal{L}}}{\partial{\boldsymbol{\theta}}} &= \frac{\partial}{\partial{\boldsymbol{\theta}}} \left[\frac{1}{m} (\mathbf{y}^T \mathbf{y} - 2 \mathbf{y}^T \mathbf{X} \boldsymbol{\theta} + \boldsymbol{\theta}^T \mathbf{X}^T \mathbf{X} \boldsymbol{\theta}) \right] \\ &= \dfrac{1}{m} \dfrac{\partial}{\partial \boldsymbol{\theta}} \textbf{tr} \left( \mathbf{y}^T \mathbf{y} - 2 \mathbf{y}^T \mathbf{X} \boldsymbol{\theta} + \boldsymbol{\theta}^T \mathbf{X}^T \mathbf{X} \boldsymbol{\theta} \right) \\ &= - \dfrac{2}{m} \dfrac{\partial}{\partial \boldsymbol{\theta}} \textbf{tr} (\mathbf{y}^T\mathbf{X}\boldsymbol{\theta}) + \dfrac{1}{m} \left( \dfrac{\partial}{\partial \boldsymbol{\theta}^T} \textbf{tr} (\boldsymbol{\theta}^T \mathbf{X}^T \mathbf{X} \boldsymbol{\theta} \mathbf{I}) \right)^T \\ &= - \dfrac{2}{m} \mathbf{X}^T \mathbf{y} + \dfrac{1}{m} \left(\boldsymbol{\theta}^T \mathbf{X}^T \mathbf{X} + \boldsymbol{\theta}^T \mathbf{X}^T \mathbf{X}\right)^T \\ &= - \frac{2}{m} \mathbf{X}^T \mathbf{y} + \frac{2}{m} \mathbf{X}^T \mathbf{X} \boldsymbol{\theta} \tag{43} \end{align*}

为了求 $\mathcal{L}$ 的极值点,令上式等于 $0$,并求解 $\boldsymbol{\theta}$:


\begin{align*} - \frac{2}{m} \mathbf{X}^T \mathbf{y} + \frac{2}{m} \mathbf{X}^T \mathbf{X} \boldsymbol{\theta} &= \mathbf{0} \\ \mathbf{X}^T \mathbf{X} \boldsymbol{\theta} &= \mathbf{X}^T \mathbf{y} \\ \boldsymbol{\theta} &= \left( \mathbf{X}^T \mathbf{X} \right)^{-1} \mathbf{X}^T \mathbf{y} \tag{44} \end{align*}

上式也叫作正规方程。注意上式最后一步成立的条件是 $\mathbf{X}^T \mathbf{X}$ 必须可逆,也就是说假设 $\mathbf{X}$ 是一个 $m \times n$ 的矩阵,则必须保证 $R(\mathbf{X}) = n$( $R$ 表示矩阵的秩)且 $n \le m$。即 $(44)$ 成立的条件为:


$$ R(\mathbf{X}) = n \quad \text{且} \quad n \le m \tag{45} $$

注意 $n \le m$ 就是说样本特征的数量必须小于等于样本的数量。因而在特征数量大于样本数量的场合,不能使用最小二乘法!

现在得到了 $\mathcal{L}$ 可能的极值点,为了验证该点是否最小值点,对 $\mathcal{L}$ 求二阶导数:


$$ \frac{\partial^2 \mathcal{L}}{\partial \boldsymbol{\theta} \partial \boldsymbol{\theta}^T} = \frac{2}{m} \mathbf{X}^T \mathbf{X} \tag{46} $$

上式就是 $\mathcal{L}$ 在任意一点 $\mathbf{X}$ 的海森矩阵。我们的目的是要判断函数是否取得最小值,因此下面验证 $\mathbf{X}^T \mathbf{X}$ 是否为正定矩阵($\dfrac{2}{m}$ 为正的常数,暂时忽略)。这里假设 $\mathbf{z}$ 为任意非零向量,则:


\begin{align*} \mathbf{z}^T \left( \mathbf{X}^T \mathbf{X} \right) \mathbf{z} &= \left( \mathbf{z}^T \mathbf{X}^T \right) \left( \mathbf{X} \mathbf{z} \right) \\ &= \left( \mathbf{X} \mathbf{z} \right)^T \left( \mathbf{X} \mathbf{z} \right) \\ &= \left \| \mathbf{X} \mathbf{z} \right \|^2 \geqslant 0 \tag{47} \end{align*}

到这里我们证明了 $\mathbf{X}^T \mathbf{X}$ 是一个半正定矩阵,离正定矩阵还有一步之遥。接下来我们看看 $\mathbf{z}^T ( \mathbf{X}^T \mathbf{X} ) \mathbf{z}$ 为 $0$ 的情况:

令 $ \mathbf{z}^T ( \mathbf{X}^T \mathbf{X} ) \mathbf{z} = 0 $,根据 $(47)$ 式,有:$ \mathbf{X} \mathbf{z} = \mathbf{0} $,这就转化为了求 $n$ 元线性方程组的解的问题。我们知道 $\mathbf{X} \mathbf{z} = \mathbf{0}$ 有非零解的充要条件是 $R(\mathbf{X}) < n$,也就是说只有 $R(\mathbf{X}) < n$ 时才能保证对任意非零向量 $\mathbf{z}$ 有 $ \mathbf{X} \mathbf{z} = \mathbf{0} $。那么反过来就可以说,当 $R(\mathbf{X}) = n$ 时,对任意非零向量 $\mathbf{z}$ 有 $ \mathbf{X} \mathbf{z} \ne \mathbf{0} $。结合 $(47)$ 式得:当 $R(\mathbf{X}) = n$ 时,有 $ \mathbf{z}^T ( \mathbf{X}^T \mathbf{X} ) \mathbf{z} > 0 $,即当 $R(\mathbf{X}) = n$ 时,$\mathbf{X}^T \mathbf{X}$ 为正定矩阵。

注意上面讨论的讨论中 $\mathbf{X}$ 可以是任意取值的,说明海森矩阵在最小二乘函数取值范围内的任意一点上都为正定矩阵,因此最小二乘函数是一个凸函数。所谓凸函数,可以简单的理解为在其图形上方的区域中任取两点,这两点的连线仍在该区域内。

好了,说结论:当 $R(\mathbf{X}) = n$ 时,$\mathbf{X}^T \mathbf{X}$ 为正定矩阵,函数取得最小值。这样我们就证明了 $(44)$ 式确实为 $\mathcal{L}$ 的最小值。

总结

终于,我们得到了最小二乘法的结论:

假设样本集 $\mathbf{X}$ 中样本的个数为 $m$,样本的特征个数为 $n$,即 $\mathbf{X}$ 是一个 $m \times n$ 的矩阵,并且 $R(\mathbf{X}) = n$,则回归参数的最优解 $\widehat{\boldsymbol{\theta}}$ 为:


$$ \widehat{\boldsymbol{\theta}} = \left( \mathbf{X}^T \mathbf{X} \right)^{-1} \mathbf{X}^T \mathbf{y} \tag{48} $$

最后,我们来看看 $R(\mathbf{X}) = n$ 这个条件,它说了两件事:

  1. 样本数量必须大于特征数量, 这在实际应用中很容易达成。
  2. 样本集中不能存在特征完全相同的两个或两个以上样本,这也很容易达成,只要把特征相同的样本剔除就行了。

到这里最小二乘法就讲完了,下面我们来看最大似然法。

最大似然法

最大似然法的核心思想是:假设真实值是由服从某个特定概率分布的随机变量产生的,算法的目标则是调整模型参数使模型产生这些真实值的概率达到最大。

比如去买彩票,为了预测中奖率我们训练了一个模型,这个模型预测的中奖率为 $99\%$,我们又问了 100 位朋友是否中奖,他们说都没中奖,这就说明模型预测结果是非常不准确的,需要进一步调整模型参数。

最大似然法并不像最小二乘法那样假定样本集中的数据有完美的线性规律,而是假设它们有一定的误差。也就是说,最小二乘法假定样本数据是这样计算的:


$$ y = f(\mathbf{x}; \boldsymbol{\theta}) = \boldsymbol{\theta}^T \mathbf{x} \tag{49} $$

而最大似然法假定样本数据是这样计算的:


$$ y = f(\mathbf{x}; \boldsymbol{\theta}, \epsilon) = \boldsymbol{\theta}^T \mathbf{x} + \epsilon \tag{50} $$

其中 $y$ 表示样本值,$\epsilon$ 表示(服从某个概率分布的)误差。也就是说最大似然法考虑到样本在采集的过程中会出现一定的误差,因此用带有误差的算式来模拟这种情况。

也许你已经迫不及待地想要求解最大似然法了,不过请稍等一下,我们还需要复习一些概率论的基础知识。

基础知识

概率和条件概率

首先从最简单的开始,介绍什么是概率和条件概率。

概率就是一件事情发生的可能性,比如抛一枚硬币,得到正面的概率为 $50\%$ 。概率通常记为 $P(x)$ 或 $p(x)$。

条件概率就是在一件事情已经发生的情况下,另一件事情发生的概率。比如 “打雷喽下雨收衣服啊”,此时要计算的就是在打雷的条件下,下雨的概率。条件概率通常记为 $P(x \mid y)$ 或 $p(x \mid y)$。比如 $p \left( \text{下雨} \mid \text{打雷} \right)$ 就表示在打雷发生的情况下,下雨的概率。

在两个事件毫不相关的时候,有 $p(x \mid y) = p(x) p(y)$ 。比如抛两次硬币,第一次出现正面的概率是 $50\%$,第二次出现正面的概率还是 $50\%$,不会因为第一次出现或不出现正面而影响到第二次出现或不出现正面。

随机变量

下面介绍随机变量,随机变量也是一种实值变量,但它的取值是有一定概率的。比如我们抛两次硬币,用变量 $X$ 表示出现正面的次数($H$ 表示正面,$T$ 表示反面),则 $X$ 的取值为:


\begin{equation} X = \left\{ \begin{array}{ll} 2, & HH \\ 1, & HT, TH \\ 0, &TT \end{array} \right. \tag{51} \end{equation}

可以看到 $X$ 的每个取值是有一定概率的,比如 $X = 2$ 的概率为 $25\%$,这样的变量就叫做随机变量。

随机变量分为离散型随机变量和连续性随机变量,上面例子里的随机变量就是离散型随机变量。而 “每年的降雨量” 就属于连续型随机变量。

概率密度

对于连续型随机变量,通常要考察的是变量小于某个实数值的概率,为此引入分布函数的概念。分布函数用来描述随机变量 $X$ 小于某个实数值 $x$ 的概率,表示为:


$$ F(x) = P \{X \le x \} , \ - \infty \lt x \lt \infty \tag{52} $$

其中 $F(x)$ 即为分布函数,$P \{X \le x \}$ 表示 $\{X \le x \}$ 这一事件发生的概率。

通常分布函数 $F(x)$ 可以写成一个非负函数 $f(x)$ 的积分形式:


$$ F(X) = \int_{-\infty}^{x} f(t) \mathrm{d} t \tag{53} $$

$f(x)$ 就叫做 $X$ 的概率密度函数,简称概率密度。

正态分布

如果连续性随机变量 $X$ 的概率密度为


$$ f(x) = \frac{1}{\sqrt{2 \pi} \sigma} \mathrm{e}^{-\dfrac{(x - \mu)^2}{2 \sigma^2}}, \ - \infty \lt x \lt \infty \tag{54} $$

其中 $\mu$, $\sigma (\sigma \gt 0)$ 为常数,则称 $X$ 服从参数为 $\mu$, $\sigma$ 的正态分布或高斯分布,记为 $X \thicksim \mathcal{N} \left( \mu, \sigma^2 \right)$ . 这里的 $\mu$ 和 $\sigma$ 有着特殊的含义,它们分别表示均值和方法(这里不予证明)。

在实际生活中很多现象都符合正态分布,比如我们通常所说的误差就是如此。

下图为 $X \thicksim \mathcal{N} \left( 0, 1 \right)$ 的概率密度图像:

高斯分布

可以看到图像呈钟形,也就是说越往中间的区间, $x$ 取得其中的值的概率越大,而越往两边的区间,$x$ 取得其中的值的概率越小。

从公式和图像中很容易得到 $X \thicksim \mathcal{N} \left(\mu, \sigma^2 \right)$ 的图像以 $\mu$ 为对称轴,当 $x = \mu$ 时取得最大值: $\dfrac{1}{\sqrt{2 \pi} \sigma}$ .

如果 $X$ 服从高斯分布,那么 $X + C$ (其中 $C$ 是任意常数)也服从高斯分布,并且图像向左(当 $C \lt 0$ 时)或向右(当 $C \gt 0$ 时)平移 $C$ 个单位,如下图:

常数对分布的影响

以上就是我们需要复习的基础,接下来就可以求解最大似然法了。

求解

上面我们提到最大似然法假定样本数据有一个误差,并且给出了计算单个样本输出值的函数:


$$ y = f(\mathbf{x}; \boldsymbol{\theta}, \epsilon) = \boldsymbol{\theta}^T \mathbf{x} + \epsilon \tag{55} $$

假设 $\epsilon$ 服从正态分布(这是最大似然法求解线性回归问题的第一个假设),即:


$$ \epsilon \thicksim \mathcal{N} \left( \mu, \sigma^2 \right) \tag{56} $$

那么这里的 $\mu$ 和 $\sigma$ 如何取值呢?由于误差通常都是有正有负的,比如测量实际为10厘米的零件,测出来的数据通常有的大于10厘米一点,有的小于10厘米一点,这样一平均,就可以认为误差的均值为 $0$ 。因此我们这里令 $\mu = 0$(这是最大似然法求解线性回归问题的第二个假设)。而 $\sigma$ 我们是没办法事先得知的,因此仍把它作为模型参数。把 $\mu = 0$ 代入 $(56)$ 式,得:


$$ \epsilon \thicksim \mathcal{N} \left( 0, \sigma^2 \right) \tag{57} $$

由于 $\boldsymbol{\theta}^T \mathbf{x}$ 是一个与 $\epsilon$ 无关的常量(误差是采集数据时产生的,和数据的本质没有任何关系),根据上面的基础知识,可知输出值也服从正态分布:


$$ \left( \boldsymbol{\theta}^T \mathbf{x} + \epsilon \right) \thicksim \mathcal{N} \left( \boldsymbol{\theta}^T \mathbf{x}, \sigma^2 \right) \tag{58} $$

输出值服从高斯分布,也就是输出值的概率等于高斯分布的概率密度:


$$ p \left( y \mid \mathbf{x}, \boldsymbol{\theta}, \sigma^2 \right) = \mathcal{N} \left( \boldsymbol{\theta}^T \mathbf{x}, \sigma^2 \right) \tag{59} $$

其中 $y$ 表示输出值,$p \left( y \mid \mathbf{x}, \boldsymbol{\theta}, \sigma^2 \right)$ 表示在 $\mathbf{x}$, $\boldsymbol{\theta}$ 以及 $\sigma^2$ 同时取定时,$y$ 值的概率。

最大似然法的思想,就是固定 $(59)$ 式中的 $y$ 和 $\mathbf{x}$(因为它们是已知数据,我们无法改变它们),调整模型参数 $\boldsymbol{\theta}$ 和 $\sigma^2$,以使模型产生 $y$ 的概率 $p \left( y \mid \mathbf{x}, \boldsymbol{\theta}, \sigma^2 \right)$ 达到最大。比如下图中蓝色曲线是我们当前的模型,可以看到模型输出 $5$ 的概率最大(A点),而输出实际数据 $5.25$(B点)的概率不是最大,因此模型需要进一步调整,比如向右调整到红色虚线的位置。

参数调整

$(59)$ 式描述的是单个样本的概率,为了描述整个样本集的概率,将所有样本的概率整合到一起:


$$ \mathcal{L} = p \left( y_1, y_2, \cdots, y_m \mid \mathbf{x}_1, \mathbf{x}_2, \cdots, \mathbf{x}_n, \boldsymbol{\theta}, \sigma^2 \right) = p \left(\mathbf{y} \mid \mathbf{X}, \boldsymbol{\theta}, \sigma^2 \right) \tag{60} $$

假设每个样本的误差都是相互独立的,那么 $(60)$ 式就可以拆成单个样本概率的积的形式:


$$ \mathcal{L} = p \left(\mathbf{y} \mid \mathbf{X}, \boldsymbol{\theta}, \sigma^2 \right) = \prod\limits_{i=1}^m p \left( y_i \mid \mathbf{x}_i, \boldsymbol{\theta}, \sigma^2 \right) = \prod\limits_{i=1}^m \mathcal{N} \left( \boldsymbol{\theta}^T \mathbf{x}_i, \sigma^2 \right) \tag{61} $$

上式就叫做似然函数。求解使似然函数达到最大值时回归参数的方法就叫做最大似然法。也就是说我们需要求解的是使似然函数 $\mathcal{L}$ 达到最大值时的参数 $\widehat{\boldsymbol{\theta}}$ 和 $\widehat{\sigma}$,即:


$$ \widehat{\boldsymbol{\theta}}, \widehat{\sigma} = \underset{\boldsymbol{\theta}, \sigma}{\arg\max}\ p \left(\mathbf{y} \mid \mathbf{X}, \boldsymbol{\theta}, \sigma^2 \right) \tag{62} $$

$(61)$ 式是一个由 $m$ 个子式连乘的式子,研究起来很困难,而加法相对来说就简单多了,又由于取对数不影响一个函数的最值,因此我们对 $(61)$ 式求对数,得到一个连加的式子,然后整理,得:


\begin{align*} \log{\mathcal{L}} &= \log{ \prod\limits_{i=1}^m \mathcal{N} \left( \boldsymbol{\theta}^T \mathbf{x}_i, \sigma^2 \right) } \\ &= \log{ \prod\limits_{i=1}^m \frac{1}{\sqrt{2 \pi} \sigma} \mathrm{e}^{-\dfrac{(y_i - \boldsymbol{\theta}^T \mathbf{x}_i)^2}{2 \sigma^2}} } \\ &= \sum\limits_{i=1}^m \log \frac{1}{\sqrt{2 \pi} \sigma} \mathrm{e}^{-\dfrac{(y_i - \boldsymbol{\theta}^T \mathbf{x}_i)^2}{2 \sigma^2}} \\ &= \sum\limits_{i=1}^m \left[ - \dfrac{1}{2} \log{ 2 \pi } - \log \sigma - \dfrac{1}{2 \sigma^2} \left( y_i - \boldsymbol{\theta}^T \mathbf{x}_i \right)^2 \right] \\ &= - \dfrac{m}{2} \log{ 2 \pi } - m \log \sigma - \dfrac{1}{2 \sigma^2} \sum\limits_{i=1}^m \left( y_i - \boldsymbol{\theta}^T \mathbf{x}_i \right)^2 \tag{63} \end{align*}

$(63)$ 式就叫做对数似然函数。上式中最后一个子式是不是非常熟悉呢?如果没想起来,请参看 $(39)$ 式。

根据 $(39)$ $(40)$ $(41)$ 式,可以进一步整理:


\begin{align*} \log{\mathcal{L}} &= - \dfrac{m}{2} \log{ 2 \pi } - m \log \sigma - \dfrac{1}{2 \sigma^2} \left( \mathbf{y} - \mathbf{X} \boldsymbol{\theta} \right)^T \left( \mathbf{y} - \mathbf{X} \boldsymbol{\theta} \right) \\ &= - \dfrac{m}{2} \log{ 2 \pi } - m \log \sigma - \dfrac{1}{2 \sigma^2} \left(\mathbf{y}^T \mathbf{y} - 2 \mathbf{y}^T \mathbf{X} \boldsymbol{\theta} + \boldsymbol{\theta}^T \mathbf{X}^T \mathbf{X} \boldsymbol{\theta} \right) \tag{64} \end{align*}

为了求最大值,我们需要对该式的参数 $\boldsymbol{\theta}$ 和 $\sigma$ 分别求导,并令导数等于 $0$。首先来求 $\dfrac{\partial{\log \mathcal{L}}}{\partial \boldsymbol{\theta}}$,求解这个导数与最小二乘法差不多,这里仅列出推导过程,如果有不明白的地方,请参考最小二乘法的推导过程。


\begin{align*} \dfrac{\partial{\log \mathcal{L}}}{\partial \boldsymbol{\theta}} &= \dfrac{1}{\sigma^2} \frac{\partial{ \mathbf{y}^T \mathbf{X} \boldsymbol{\theta} }}{\partial{ \boldsymbol{\theta} }} - \frac{1}{2 \sigma^2} \frac{\partial{ \boldsymbol{\theta}^T \mathbf{X}^T \mathbf{X} \boldsymbol{\theta} }}{ \partial{ \boldsymbol{\theta} }} \\ &= \frac{1}{\sigma^2} \mathbf{X}^T \mathbf{y} - \frac{1}{\sigma^2} \mathbf{X}^T \mathbf{X} \boldsymbol{\theta} \tag{65} \end{align*}

可以注意到,上式和最小二乘法对 $\boldsymbol{\theta}$ 的导数非常像,为了比较,把最小二乘法对 $\boldsymbol{\theta}$ 的导数 $(42)$ 式再次列出:


$$ \frac{\partial{\mathcal{L}}}{\partial{\boldsymbol{\theta}}} = - \frac{2}{m} \mathbf{X}^T \mathbf{y} + \frac{2}{m} \mathbf{X}^T \mathbf{X} \boldsymbol{\theta} \tag{66} $$

可以看到,它们的形式一模一样,但是它们的符号相反

形式一模一样说明它们都在同一点取得极值,即令 $(65)$ 式等于 $0$,可以解得与 $(42)$ 式一样的 $\widehat{\boldsymbol{\theta}}$:


$$ \widehat{\boldsymbol{\theta}} = \left( \mathbf{X}^T \mathbf{X} \right)^{-1} \mathbf{X}^T \mathbf{y} \tag{67} $$

而符号相反说明,二乘函数在 $\widehat{\boldsymbol{\theta}}$ 取得最小值,而似然函数在 $\widehat{\boldsymbol{\theta}}$ 处可能取得最大值。由于二乘函数对 $\boldsymbol{\theta}$ 的二阶偏导为正定矩阵,因此似然函数对 $\boldsymbol{\theta}$ 的二阶导数为负定矩阵,这里就不再证明了。

为什么说似然函数在 $\widehat{\boldsymbol{\theta}}$ 处可能取得最大值呢?因为我们还需要考虑另外一个模型参数 $\sigma$ 的导数。

$(64)$ 式中,对 $\sigma$ 求导并令导数等于 $0$,解出 $\widehat{\sigma}$:


\begin{align*} \dfrac{\partial{\log \mathcal{L}}}{\partial \sigma} &= - \dfrac{m}{\sigma} + \dfrac{1}{\sigma^3} \left(\mathbf{y}^T \mathbf{y} - 2 \mathbf{y}^T \mathbf{X} \boldsymbol{\theta} + \boldsymbol{\theta}^T \mathbf{X}^T \mathbf{X} \boldsymbol{\theta} \right) = 0 \\ \widehat{\sigma^2} &= \dfrac{1}{m} \left(\mathbf{y}^T \mathbf{y} - 2 \mathbf{y}^T \mathbf{X} \boldsymbol{\theta} + \boldsymbol{\theta}^T \mathbf{X}^T \mathbf{X} \boldsymbol{\theta} \right) \tag{68} \end{align*}

$(67)$ 式代入 $(68)$ 式,得:


\begin{align*} \widehat{\sigma^2} &= \dfrac{1}{m} \left[ \mathbf{y}^T \mathbf{y} - 2 \mathbf{y}^T \mathbf{X} \left( \mathbf{X}^T \mathbf{X} \right)^{-1} \mathbf{X}^T \mathbf{y} + \mathbf{y}^T \mathbf{X} \left( \mathbf{X}^T \mathbf{X} \right)^{-1} \mathbf{X}^T \mathbf{X} \left( \mathbf{X}^T \mathbf{X} \right)^{-1} \mathbf{X}^T \mathbf{y} \right] \\ &= \dfrac{1}{m} \left[ \mathbf{y}^T \mathbf{y} - 2 \mathbf{y}^T \mathbf{X} \left( \mathbf{X}^T \mathbf{X} \right)^{-1} \mathbf{X}^T \mathbf{y} + \mathbf{y}^T \mathbf{X} \left( \mathbf{X}^T \mathbf{X} \right)^{-1} \mathbf{X}^T \mathbf{y} \right] \\ &= \dfrac{1}{m} \left[ \mathbf{y}^T \mathbf{y} - \mathbf{y}^T \mathbf{X} \left( \mathbf{X}^T \mathbf{X} \right)^{-1} \mathbf{X}^T \mathbf{y} \right] \\ &= \dfrac{1}{m} \left( \mathbf{y}^T \mathbf{y} - \mathbf{y}^T \mathbf{X} \widehat{\boldsymbol{\theta}} \right) \tag{69} \end{align*}

对 $\sigma$ 的二阶导数为:


\begin{align*} \dfrac{\partial^2 \log {\mathcal{L}}}{\partial \sigma^2} &= \dfrac{m}{\sigma^2} - \dfrac{3}{\sigma^4} \left(\mathbf{y}^T \mathbf{y} - 2 \mathbf{y}^T \mathbf{X} \boldsymbol{\theta} + \boldsymbol{\theta}^T \mathbf{X}^T \mathbf{X} \boldsymbol{\theta} \right) \tag{70} \end{align*}

$(68)$ 式代入,得:


$$ \dfrac{\partial^2 \log {\mathcal{L}}}{\partial \sigma^2} = \dfrac{m}{\widehat{\sigma^2}} - \dfrac{3} {\left( \widehat{\sigma^2} \right)^2} m \widehat{\sigma^2} = - \frac{2m}{\widehat{\sigma^2}} \tag{71} $$

可以看到对 $\sigma$ 的二阶导数恒为负数。结合上面对 $\boldsymbol{\theta}$ 的二阶导数为负定矩阵的结论,即可知似然函数在 $\widehat{\boldsymbol{\theta}}$ 和 $\widehat{\sigma^2}$ 处取得最大值。

到这里我们求得了最大似然法的最优解,值得注意的是最大似然法是从实际样本的分布规律来求解线性回归参数的,再来看一下我们是如何模拟实际数据的:


$$ y = f(\mathbf{x}; \boldsymbol{\theta}, \epsilon) = \boldsymbol{\theta}^T \mathbf{x} + \epsilon \tag{72} $$

实际数据包含了两个部分,一部分是内在的、不变的线性规律:$\boldsymbol{\theta}^T \mathbf{x}$,另一部分是变化的,不确定的误差:$\epsilon$。注意在实际求解时,$x$ 和 $y$ 分别表示实际样本的特征值与目标值,它们是常量。

在预测时,我们更关心内在的、不变的线性规律,而不在乎实际数据的误差是多少,因此仅用第一部分进行预测:


$$ y = f(\mathbf{x}; \widehat{\boldsymbol{\theta}}) = \widehat{\boldsymbol{\theta}}^T \mathbf{x} \tag{73} $$

这里的 $\mathbf{x}$ 是自变量,表示输入的样本特征,$y$ 是因变量,表示输出的预测值,$\widehat{\boldsymbol{\theta}}$ 则是通过最大似然法求得的 $\boldsymbol{\theta}$ 最优解,为常量。

关于最大似然法我们就讲到这里了,最后做一个总结。

总结

最大似然法将样本的误差考虑在内,以达到输出实际数据集的概率最大为目标对回归参数进行求解。

样本数据的模型:


$$ y = f(\mathbf{x}; \boldsymbol{\theta}, \epsilon) = \boldsymbol{\theta}^T \mathbf{x} + \epsilon \tag{74} $$

最优解为:


\begin{equation} \begin{array}{ll} \widehat{\boldsymbol{\theta}} &= \left( \mathbf{X}^T \mathbf{X} \right)^{-1} \mathbf{X}^T \mathbf{y} \\ \widehat{\sigma^2} &= \dfrac{1}{m} \left( \mathbf{y}^T \mathbf{y} - \mathbf{y}^T \mathbf{X} \widehat{\boldsymbol{\theta}} \right) \tag{75} \end{array} \end{equation}

用于预测的模型:


$$ y = f(\mathbf{x}; \widehat{\boldsymbol{\theta}}) = \widehat{\boldsymbol{\theta}}^T \mathbf{x} \tag{76} $$

示例

(本例完整的代码可以到我的码云下载。)

现在就可以很容易地解决开始的那个问题了,比如有下面的数据集:

所有样本点

这些数据集是由下面的函数产生的:

1
2
3
4
5
6
7
8
9
10
11
m = 500  # 样本个数
theta = np.array([[2], [3]]) # 实际的模型参数

def data_set(m, theta):
np.random.seed(42)
x = np.linspace(-10, 10, m).reshape(-1, 1)
X = np.c_[np.ones([m, 1]), x]
y = np.dot(X, theta) + np.random.normal(0, 10, m).reshape(-1, 1)
return X, y

X, y = data_set(m, theta)

把 $\mathbf{X}$ 和 $\mathbf{y}$ 代入 $(48)$$(75)$ 式(分别是最小二乘法和最大似然法求得的 $\widehat{\boldsymbol{\theta}}$,结果都是一样的):

1
2
best_theta = np.linalg.inv(X.T.dot(X)).dot(X.T).dot(y)
print(best_theta)

最后我们把所得的直线绘制在散点图上:

1
2
3
4
plot_scatter(X, y)
y_pred = np.dot(X, best_theta)
plt.plot(X[:,1], y_pred, color='r', linewidth=2)
plt.show()

最优解直线

模型的评估

线性模型通常使用均方差(Mean Square Error, MSE)来评估,均方差就是实际值与真实值的差的平方的平均数,与 $(8)$ 式一样,只不过把 $(8)$ 式应用于验证集罢了:


$$ \mathbf{C} = \frac{1}{m} \sum_{i=1}^m (y_i - \boldsymbol{\theta}^T \mathbf{x}_i)^2 \tag{77} $$

上式中的 $\mathbf{C}$ 表示损失函数(代价函数),$m$ 为验证集的样本数量,$\mathbf{x}_i$ 表示第 $i$ 个验证集样本的特征,$y_i$ 则表示第 $i$ 个验证集样本的真实值。

均方差的大小说明了模型的拟合度,均方差越大,说明模型的拟合度(模型与真实数据的贴合程度)越低,均方差越小,说明模型的拟合度越高。

将示例中的预测结果带入 $(77)$ 式就得到均方差:

1
2
mse = np.mean(np.square(y_pred - y))
print(mse)

输出结果为:

1
96.04753271883293

多项式回归

在前面的推导过程中,始终将 $\mathbf{X}$ 看做是常数,因此使用线性回归也可以求解如下的多项式回归问题:


$$ f(x; \theta_0, \theta_1, \dots, \theta_n) = \theta_0 x^0 + \theta_1 x^1 + \theta_2 x^2 + \dots + \theta_n x^n = \boldsymbol{\theta}^T \mathbf{x} \tag{78} $$

只要保证参数 $\theta_0, \theta_1, \dots, \theta_n$ 是线性(幂为1)的,就可以通过与线性回归相同的推导过程得到上式的最优解为:


$$ \widehat{\boldsymbol{\theta}} = \left( \mathbf{X}^T \mathbf{X} \right)^{-1} \mathbf{X}^T \mathbf{y} \tag{79} $$

其中,


$$ \mathbf{X} = \begin{bmatrix} 1 & x_1 & x_1^2 & \cdots & x_1^n \\ 1 & x_2 & x_2^2 & \cdots & x_2^n \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_n & x_n^2 & \cdots & x_n^n \end{bmatrix} \tag{80} $$

示例

假设有如下的数据集:

1
2
3
4
5
6
7
8
9
10
[[ 1.49816048 17.02188475]
[ 3.80285723 64.4104382 ]
[ 2.92797577 41.73590281]
[ 2.39463394 31.08275619]
[ 0.62407456 9.02764883]
[ 0.62397808 5.7940379 ]
[ 0.23233445 3.50389893]
[ 3.46470458 54.02792058]
[ 2.40446005 30.44287786]
[ 2.83229031 40.86613746]]

多项式数据集

假设数据集表示 $x$ 的二次函数,则用 $y = \theta_0 + \theta_1 x + \theta_2 x^2 $ 的多项式来拟合,代码如下:

1
2
3
X = np.c_[np.ones([m, 1]), x, x**2]
best_theta_poly = np.linalg.inv(X.T.dot(X)).dot(X.T).dot(y)
print(best_theta_poly)

输出为:

1
2
3
[[2.9937327 ]
[4.56883305]
[2.99927292]]

画出图形:

多项式回归求解

当然也可以假设数据为 $x$ 的 $7$ 次函数,则用 $ y = \theta_0 + \theta_1 x + \theta_2 x^2 + \theta_3 x^3 + \theta_4 x^4 + \theta_5 x^5 + \theta_6 x^6 + \theta_7 x^7 $ 的多项式来拟合,代码如下:

1
2
3
X = np.c_[np.ones([m, 1]), x, x**2, x**3, x**4, x**5, x**6, x**7]
best_theta_poly = np.linalg.inv(X.T.dot(X)).dot(X.T).dot(y)
print(best_theta_poly)

得到最优解为:

1
2
3
4
5
6
7
8
[[  116.21141935]
[ -896.17156581]
[ 2305.72649171]
[-2649.00289328]
[ 1581.09144523]
[ -509.86317057]
[ 84.37417362]
[ -5.6214629 ]]

图形为:

多项式回归-过拟合

那么上面两个多项式哪个更符合实际规律呢?或者说哪个模型在测试集中表现得更好呢?我们将在过拟合与欠拟合中揭晓答案。

总结

  • 线性回归问题可以使用最小二乘法和最大似然法来解决
  • 最小二乘法的思想是将真实值与预测值的差距最小化
  • 最大似然法的思想是将模型输出真实值的概率最大化
  • 使用 MSE 评估线性回归模型
  • 使用线性回归模型可以解决多项式回归问题

以上就是线性回归的内容,感谢阅读!如果有任何疑问或建议,请在评论区留言!祝你度过愉快的一天!

版权声明:本文为原创文章,转载请注明出处。http://cynhard.com/2018/06/23/ML-Linear-Regression/

推荐文章