难舍的机器学习(1)—最小二乘法和极大似然来解线性回归线性回归最小二乘法几何角度来看最小二乘法概率角度看最小二乘法

这是我参与8月更文挑战的第29天，活动详情查看：8月更文挑战

文章有点，而且大部分都是公式推导，不过推导过程中已经尽量把每一步解释清楚，通过阅读，大家对如何使用最小二乘法和极大似然法来求解线性模型的

线性回归

最小二乘法
几何角度来看最小二乘法
概率角度看最小二乘法

什么是回归问题

今天我们来讨论一下线性模型，之前已经了解到线性模型来做回归问题，回归问题多用来预测一个具体的数值，如预测房价、股价(这个应用用随机过程来解)、未来的天气情况等等。例如我们根据一个地区的若干年的 PM2.5 数值变化来估计某一天该地区的 PM2.5 值大小，预测值与当天实际数值大小越接近，回归分析算法的可信度越高。

线性模型用于回归问题，就是根据给定样本数据训练出一个线性模型(拟合现有样本点)，然后用这个线性模型来对新样本预估一个值。通常模型

h_{\theta} = \theta_0 + \theta_1 x_1 + \theta_2 x_2 + \dots + \theta_n x_n

使用一元线性模型去做回归问题，那么该任务就称之为一元线性回归。一元线性回归的求解可采用最小二乘法，也可以采用迭代的方法去求解。

准备数据

监督学习中，样本 D 通常都是这样表示的，每一个样是一个 n 维特征的 X 向量和表示标签的 y 组成的。在回归问题中 y 是具体数值，而在分类问题中 y 是表示样本所属的类别。

D = \{ (x^{(1)},y^{(1)}),(x^{(2)},y^{(2)}),\dots, (x^{(n)},y^{(n)}) \} x^{(i)} \in \mathbb{R}^n

有了样本之后，就可以开始定义模型了，模型其实就是用数学的语言描述一些事物以及事物之间关系，这是我个人对模型的理解，就是用数学语言描述事物以及事物之间的关系。算法也就是如何基于模型求解出优化问题的解。

定义模型

那么这里模型就是表示样本中 X 和 y 之间的关系，下面来定义线性函数来表示 x 和 y 之间关系，函数参数用 $\theta$ 来表示，可以写成下面的样子

h_{\theta} = \theta_0 + \theta_1 x_1 + \theta_2 x_2 + \dots + \theta_n x_n

\theta = \{ \theta_0, \theta_1, \dots, \theta_n \}^T

这里 $\theta_0$ 表示偏置。

设定目标函数

通过计算样本点到线性模型(直线或超平面)距离来评估我们模型对数据拟合程度，也就是我们要优化的目标。因为通过衡量真实值估计值之间距离最小作为目标。我们目标就是让目标函数(也叫损失函数)取值最小，这里用 $\hat{y} - y$ 的作为点到线距离，大家可能会有疑问点到直线距离不是从点引一条直线的垂线，其实这是等价的问题，并不需要去求距离，较少预测值到真实值之间平方差就是等于样本点到线性方差直线之间的距离。

如果数据量在 100 万之内，我们无需用梯度下降，可以使用最小二乘法来直接求解问题。

最小二乘法(LMS)

在线性问题中，我们找到一条直线来拟合样本，然后计算点到直线距离来来评估找到模型好坏。

D = \{ (x^{(1)},y^{(1)}),(x^{(2)},y^{(2)}),\dots, (x^{(n)},y^{(n)}) \} x^{(i)} \in \mathbb{R}^n

h_{\theta} = \theta_0 x_0 + \theta_1 x_1 + \dots + \theta_n x_n

接下来用线性方程组来表示模型，之所以矩阵，可以将一些线性方程组数学问题转换到空间上问题

\begin{aligned} \theta_0 x_0 + \theta_1 x_1^{(1)} + \theta_2 x_2^{(1)} + \cdots + \theta_m x_m^{(1)} = y^{(1)} \\ \theta_0 x_0 + \theta_1 x_1^{(2)} + \theta_2 x_2^{(2)} + \cdots + \theta_m x_m^{(2)} = y^{(2)} \\ \vdots \\ \theta_0 x_0 + \theta_1 x_1^{(N)} + \theta_2 x_2^{(N)} + \cdots + \theta_m x_m^{(N)} = y^{(N)} \\ \end{aligned}

现在尝试用矩阵形式来表示这个线性方程组，矩阵 A 表示我们样本，每一行表示一个样本，每一列表示样本的特征，这里有 N 个样本，小写 n 表示每一个样本的特征数量。也就是每一个样本的增广矩阵，

A = \begin{bmatrix} 1 & x_1^{(1)} & x_2^{(1)} & \dots & x_m^{(1)} \\ 1 & x_1^{(2)} & x_2^{(2)} & \dots & x_m^{(2)} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_1^{(N)} & x_1^{(N)} & \dots & x_m^{(N)} \\ \end{bmatrix}

因为在方程中的参数我们将截距用 $\theta_0$ 来表示，所在样本 A 矩阵中第一列是的值 1 因为这一列将和表示 $\theta_0$ (截距)相乘，这里 $\theta$ 是 $(n+1) \times 1$ 形状的参数矩阵，参数矩阵形状是由 A 矩阵形状和表示标签 Y $(N \times 1)$ 的向量的形状所确定的。

A_{N \times (n+1)} \theta_{(n+1) \times 1} = y_{(N \times 1)}

数据矩阵 A 是 $N\times(n+1)$ 和 $(n+1)\times1$ 得到矩阵就是 $N\times$ 1 ,这就是标签的形状。

\sum_{i=1}^N (h_{\theta}(x^{(i)}) - y^{(i)})^2

有了矩阵表示形式，我们就可以用矩阵形式来表示我们损失函数

S = \min ||A \theta - Y||^2

接下来对矩阵进行求模并进行化简。详细推导过程如下，这里我们有必要介绍一下矩阵的运算，特别是矩阵的转置操作，矩阵中加减法比较好理解，但是乘法就比较特别了在矩阵没有交换率 $AB \neq BA$

\begin{aligned} = (A\theta - Y)^T(A\theta - Y) \\ = (\theta^TA^T - Y^T)(A\theta - Y) \\ = \theta^TA^TA\theta - \theta^TA^TY - Y^TA\theta +Y^TY\\ = \theta^TA^TA\theta - 2\theta^TA^TY +Y^TY \end{aligned}

$A^T$ 形状是 $(n+1) \times N$ 而 $\theta^T$ 的形状为 $1 \times (n+1)$ 所以 $\theta^TA^T$ 相乘得到 $1\times N$ 的矩阵在乘以 $Y$ 这个 $N\times 1$ 的向量得到是一个标量，大家可以尝试推导一下 $Y^TA\theta$ 得到也是一个标量，因为他们都是标量所以可以相加得到

\theta^TA^TY - Y^TA\theta = -2\theta^TA^TY

之前我们都是通过梯度下降方式来不断更新参数。具体做法就是我们对参数对于损失函数进行求导后来更新参数，下面就是通过导数乘以一个学习率，学习率表示模型学习速度，也就是每一个更新参数幅度，然后用这个值来更新参数得到新的参数。

\theta_j = \theta_j + \alpha \frac{\partial L}{\partial \theta_j}

通过通过不断迭代来优化参数 $\theta$ 来找到最优解。这是之间我们做回归问题的一般步骤。

这是之前通过回归我们已经了解到了线性模型，今天我们来系统学习线性模型。

\frac{\partial S}{\partial \theta} = 0

\frac{\partial S}{\partial \theta} = \frac{\partial( \theta^TA^TA\theta - 2 \theta^TA^TY + Y^TY)}{\partial \theta}

有关矩阵求导大家已经了解到有以下公式

\frac{\partial X^T}{\partial X} = 1

然后我们可以开始对矩阵求 $\theta$ 的偏导，化简得到下面公式

\frac{\partial S}{\partial \theta} = \frac{\partial( \theta^TA^TA\theta}{\theta} - 2 \frac{\partial (\theta^TA^TY )}{\partial \theta} + \frac{\partial Y^TY}{\partial \theta}

\frac{\partial \theta^T}{\partial \theta} = 1

\frac{\partial S}{\partial \theta} = \frac{\partial (\theta^TA^TA\theta) }{\partial \theta} - 2A^TY

这一部分 $\frac{\partial (\theta^TA^TA\theta) }{\partial \theta}$ 是求导的难点，接下来我们重点就是攻克这部分求导。

有关向量求导我们先介绍一个公式

\frac{d(u^Tv)}{d(x)} = \frac{du^T}{dx} v + \frac{dv^T}{dx} u

这里 u v 是关 x 的函数，我们在线性代数已经学习过有关矩阵求导，上面就是对矩阵进行求导中一个公式，在接下来在推导过程中，我们会用这个公式，如果大家还不了解可以回去看看线性代数。

\frac{d(x^Tx)}{dx} = \frac{dx^T}{dx} x + \frac{dx^T}{dx} x = 2x

进一步来看我们引入方阵 B ，注意因为这里 B 是方阵，就可以 BX 看成一个整体来解决下面问题。

\frac{d(X^TBX)}{dx} = \frac{dX^T}{dX} BX + \frac{dX^TB^T}{dX} X = (B + B^T)X

我们知道 $A^TA$ 是一个方阵，这样我们问题就好解决了。我们通过上面推导等式

\frac{\partial \theta^TA^TA\theta}{\partial \theta} = (A^TA + AA^T)\theta = 2 A^TA\theta

\frac{\partial S}{\partial \theta} = 2 A^TA\theta - 2A^TY = 0

\theta = (A^TA)^{-1}A^TY

对于数据量不是很大情况下，可以选择最小二乘法，因为对于矩阵的逆是带来大量运算，求矩阵的逆运算，对于梯度下降问题，收敛比较慢。

最小二乘几何概念

\begin{aligned} x_1 + x_2 = 3\\ -x_1 + x_2 = 1 \end{aligned}

\begin{cases} x_1 =1\\ x_2 = 2 \end{cases}

\begin{bmatrix} 1 & 1\\ -1 & 1 \end{bmatrix} \begin{bmatrix} x_1\\ x_2\\ \end{bmatrix} = \begin{bmatrix} 3\\ 1 \end{bmatrix}

a_1x_1 + a_2x_2 = y

\left( \begin{matrix} 1\\ -1 \end{matrix} \right) x_1 + \left( \begin{matrix} 1\\ 1 \end{matrix} \right) x_2 =\left( \begin{matrix} 3\\ 1 \end{matrix} \right)

$D = \{(x^1,y^1),(x^2,y^2),(x^3,y^3) \}$ $D = \{(0,2),(1,2),(2,3) \}$

$\theta_1 x + \theta_0 = y$

\left( \begin{matrix} 0\\ 1\\ 2 \end{matrix} \right) \theta_1 + \left( \begin{matrix} 1\\ 1\\ 1 \end{matrix} \right) \theta_0 =\left( \begin{matrix} 2\\ 2\\ 3 \end{matrix} \right)

$a_1\theta_1 + a_2 \theta_0 = \hat{y}$ $L_{\min} = (y - \hat{y})$

找到 $\hat{y}$ 是向量构成平面上投影就是 y

$\hat{y} = A\theta^{*}$

$e = y - \hat{y} = y - A\theta^{*}$

\begin{cases} e\cdot a_1 = 0\\ e\cdot a_2 = 0\\ \end{cases} \rightarrow \begin{cases} a_1^Te = 0\\ a_2^Te = 0\\ \end{cases} \rightarrow A^Te = 0

\begin{aligned} A^T(y - A\theta^{*}) = 0\\ A^Ty - A^TA\theta^{*} = 0\\ A^Ty = A^TA\theta^{*}\\ \theta^{*} = (A^TA)^{-1}A^Ty \end{aligned}

概率方式来理解最小二乘

$y^i = \theta^Tx^i + \epsilon^i$

概率角度(极大似然)

从概率角度来看一下线性回归问题，以及最小二乘

D = \{ (x^{(1)},y^{(1)}),(x^{(2)},y^{(2)}),\dots, (x^{(n)},y^{(n)}) \} x^{(i)} \in \mathbb{R}^n

最小二乘法

L(\theta) = \sum_{i=1}^N ||\theta^Tx^i - y^i||^2

P(\epsilon^i) = \frac{1}{\sqrt{2\pi}\sigma}\exp\left( - \frac{(\epsilon^i)^2}{2\sigma^2} \right)

数据带有一定噪声，也就是数据有一定随机性 $y = \theta^Tx + \epsilon$

极大似然估计(MLE)

极大似然估计，通俗理解来说，就是利用已知的样本结果信息，反推最具有可能（最大概率）导致这些样本结果出现的模型参数值

$L(\theta) = \log P(Y|X;\theta) = \log \prod_i^{N} P(y^i|x^i;\theta) = \sum_{i=1}^N \log P(y^i|x^i;\theta)$

现在是模型已定，我们用样本数据去估计参数 $\theta$ , 也就是什么参数和样本数据组合后得到 $Y$ 。其实我们首先知道 y 和 x 存在已定关系，这种关系是由 $\theta$ 所确定的。

也可以这样理解，假设我们可以找到 $f(x^i) = \theta_i x^i + \theta_0$ 那么 $y^i$ 与这个函数之间差距可以看成服从正态分布。

P(y^i|x^i,\theta) = \frac{1}{\sqrt{2\pi}\sigma} \exp \left( - \frac{(y - \theta^Tx^i)^2}{2\sigma^2} \right)

L(\theta) = \sum_{i=1}^N \log \frac{1}{\sqrt{2\pi}\sigma} \exp \left( - \frac{(y - \theta^Tx^i)^2}{2\sigma^2} \right)

L(\theta) = \sum_{i=1}^N \left( \log \frac{1}{\sqrt{2\pi} \sigma} - \frac{1}{2\sigma^2}(y^i - \theta^Tx^i)^2 \right)

其实我们就是找到一个 $\theta$ 让上面的目标函数值最大，因为 $\sigma$ 是不变的值，那么想要目标函数越大，我们就需要让 $(y^i - \theta^Tx^i)$ 越小

\theta = \arg \max_{\theta} - \frac{1}{2\sigma^2}(y^i - \theta^Tx^i)^2

\theta = \arg \min_{\theta} (y^i - \theta^Tx^i)^2

最小二乘估计

$\theta = \arg \min_{\theta} (y^i - \theta^Tx^i)^2$

$L(\theta) = \sum_{i=1}^N ||\theta^Tx^i - y^i||^2$