线性回归的解法、参数求解方法比较线性回归是常用的机器学习模型，本文以线性回归的解法为例，对比了矩阵解法，最小二乘法、梯度

线性回归是常用的机器学习模型，本文以线性回归的解法为例，对比了矩阵解法，最小二乘法、梯度下降法以及牛顿法等优化算法的差异，有助于更好的理解各种优化算法。

矩阵解法

numpy求解线性方程组

对于线性方程组：

$A\bold x=\bold b$ （1）

求解矢量 $\bold x$ 。

其中A为系数矩阵

$\left[ \begin{matrix} 2&1&-2 \ 3&0&1\ 1&1&-1 \end{matrix} \right]$

b为：

$\left[ \begin{matrix} -3\ 5\ -2 \end{matrix} \right]$

通过numpy求解的具体过程如下，首先构造矩阵

import numpy as np
A = np.array([[2, 1, -2], [3, 0, 1], [1, 1, -1]])
b = np.array([[-3], [5], [-2]])

利用numpy的线性代数模块进行求解

x = np.linalg.solve(A, b)

numpy用于多元线性回归

线性回归主要目的在于找出输入与输出之间的映射函数，往往数据量远大于自变量个数。每个输入数据可以用特征向量进行表示 (x_1 , x_2 , …, x_m) ，m代表自变量的维度大小，当m为1时，即为单元线性回归。

对于多元线性回归的求解，通常是利用最小二乘法进行求解。这里先利用矩阵方法求解，后续详解最小二乘法。

输入的全部数据可以利用矩阵X表征，

$X=\left[ \begin{matrix} 1 & x_{1, 1}& x_{1, 2} &\cdots & x_{1, m}\ 1& x_{2, 1}& x_{2, 2} &\cdots & x_{2, m} \ \vdots &\vdots &\ddots &\vdots &\vdots\ 1& x_{n, 1} &x_{n, 2} &\cdots &x_{n, m} \end{matrix} \right]$

输出数据可以利用向量表征如下

$\bold y =\left[ \begin{matrix} y_1\ y_2\ \vdots\ y_n \end{matrix} \right]$ （3）

其中n代表样本个数，即输入数据的数据量。

线性回归的出发点在于找到映射函数 $\boldsymbol\beta$ ，使得

$X\boldsymbol \beta = \bold y$ （4）

其中的维度为 (n, m+1) ， $\boldsymbol \beta$ 维度为 (m+1, 1) 。

对比（4）和（1）式，二者形式相似，似乎可以利用相同的算法求解。但是仔细对比会发现，式（1）中的系数为方阵，方程个数与 $\bold x$ 的维度相同，可以直接利用进行求解。但是（4）式中，系数不为方阵，其中数据样本个数n，往往多于变量维度m+1。不能直接系利用np.linalg.solve进行求解，需要做适当变换。

利用矩阵乘法，可以得出

$(X^{T}X)\boldsymbol \beta = X^{T}\bold y$ （5）

$\boldsymbol \beta = (X^{T} X)^{-1}X^{T}\bold y$ （6）

利用numpy的线性代数模块，可以求出 $\beta$ ，如下：

X_t = np.transpose(X)
X_tX = np.dot(X_t, X)
X_ty = np.dot(X_t, y)
beta = np.linalg.solve(X_tX, X_ty)

函数np.linalg.solve

np.linalg.solve(a, b)用于线性矩阵方程（线性标量方程组）的解析解求解。

其中a必须为方阵和满秩矩阵（矩阵的所有行（列）都线性独立），如果不满足上述的任何一个条件，则必须使用最小二乘法求解。

上述对于多元线性回归，通过矩阵变换将表达式（4）变换为形式（5），从而可以利用np.linalg.solve进行求解。

线性回归的最小二乘解法

最小二乘法，又称最小平方法，通过误差的平方和最小寻找数据的最佳函数匹配。

对应于线性回归中，寻找映射函数 $\boldsymbol \beta$ ，

$\bold{\hat {y}} = X\boldsymbol {\beta}$ （7）

使得误差的平方和最小：

$\sum_{i=1}^n(\hat {y}_i-y_i )^2=||\bold{\hat y} - \bold y||^2$ （8）

其中 $X, \bold y$ 的表达式如式（2）和式（3）所示， $\boldsymbol \beta$ 维度为(m+1, 1).

要使得（8）式的平方和最小，即要找出合适的 $\boldsymbol \beta$ ，使如下式子和最小

$S = ||X\boldsymbol{\beta} - \bold y||^2$ (9)

既要使得和最小，即要找到极值点，即对 $\boldsymbol\beta$ \boldsymbol\beta 的导数 $\frac{\partial S}{\partial \boldsymbol \beta}$ 为0。

利用矩阵求导方法，首先将（9）的向量模平方改写为向量与自身的內积，即

$S=(X\boldsymbol \beta - \bold y)^{T}(X\boldsymbol \beta - \bold y)$ （10）

求微分：

$\begin{align} dS &= (X d\boldsymbol\beta)^T(X\boldsymbol\beta-\bold y)+(X\boldsymbol\beta-\bold y)^{T}(Xd\boldsymbol\beta)\ &= 2(X\boldsymbol \beta -\bold y)^T X d\boldsymbol\beta \end{align}$

其中利用到 $X\boldsymbol\beta - \bold y$ 与 $X d\boldsymbol \beta$ 均为向量，向量的內积满足 $\bold v^T\bold u=\bold u^T \bold v$ 。

利用微分与导数的联系：

$dS = \frac{\partial S}{\partial \boldsymbol\beta}^T d\boldsymbol\beta$

可以得出

$\frac{\partial S}{\partial\boldsymbol \beta}= X^T(X\boldsymbol\beta-\bold y)$

令上式为0，可得

$X^TX\boldsymbol \beta = X^T \bold y$ (11)

比较（10）式和（5）式可得，利用最小二乘法得出了与矩阵变换相同的形式。至此，可以利用numpy的线性算法模块np.linalg.solve，即可求得映射函数 $\boldsymbol \beta$ 。

最大似然概率

最小二乘法，是从使残差（误差）平方和为最小为出发点，得出参数 $\boldsymbol \beta$ 的值。用深度学习的语言来讲，残差平方和为最小，相当于目标函数。即最小二乘法是以残差平方和最小为目标函数，得出参数 $\boldsymbol \beta$ 的值。

残差平方和最小值：

$\sum_{i=1}^n (\hat{y_i} - y_i)^2$

与深度学习中的MSE（均方误差），形式相同，只是相差一个系数（1/n）：

$\frac{1}{n}\sum_{i=1}^n(\hat{y_i}-y_i)^2$

MSE可以根据高斯分布和最大似然概率推出。

同样的，对于线性回归问题，也可以从概率角度，利用最大似然概率方法和高斯分布假设得出参数 $\boldsymbol \beta$ 。具体过程如下：

样本组合表示为 ${{X_i, y_i}}_{i=1}^n$

对于多元线性回归，满足以下几个基本假设：

（1）线性性质以及自变量之间相互独立： $E[Y|X_1=x_1, X_2=x_2, ..., X_m=x_m] = \beta_0 + \beta_1 x_1+\beta_2 x_2+..+\beta_m x_m$

（2）误差项的方差为常数: $Var[\epsilon|X_1=x_1, X_2=x_2, ..., X_m=x_m] = \sigma^2$ ，称为同方差性。

（3）误差满足正态分布： $\epsilon \sim N(0, \sigma ^2)$

（4）误差相互独立：即 $\epsilon_1, \epsilon_2, ..., \epsilon_n$ 之间相互独立。

样本 ${{X_i}}_{i=1}^n$ 可以用矩阵即式子（2）表示。

多元线性回归可以表示为如下：

$\begin{align} \bold y & = \bold {\hat y} + \boldsymbol \epsilon\ & = X\boldsymbol \beta + \boldsymbol \epsilon \end{align}$

对于每个样本，可以表示为

$\begin{align} y_i &= \hat{y_i} + \epsilon_i \ &=\beta_0+X_{i1}\beta_1+X_{i2}\beta_2+...+X_{im}\beta_m + \epsilon_i \end{align}$ (12)

其中误差 $\epsilon_1, \epsilon_2, ..., \epsilon_n$ 满足独立同分布 $N(0, \sigma^2)$ ，因此 y_i 满足

$y_i|(X_{i1}, X_{i2},..., X_{im}) \sim N(\beta_0+\beta_1X_{i1}+\beta_2X_{i2}+...+\beta_mX_{im}, \sigma^2)$ (13)

由于 y_i 之间满足条件独立分布，上式可以写为矩阵形式为：

$\bold y|X \sim \boldsymbol N_n (X\boldsymbol \beta, \sigma^2\boldsymbol I)$

具体表达式为：

$\begin{align} \bold y|X & \sim \prod_{i=1}^n \frac{1}{ \sqrt {2\pi} \sigma}e^{-\frac{1}{2}(\frac{y_i-\hat{y_i}}{\sigma})^2} \ &= (\frac{1}{\sqrt{2\pi}\sigma})^n e^{-\frac{1}{2}\sum_{i=1}^n(\frac{y_i-\hat{y_i}}{\sigma})^2} \ &=(\frac{1}{\sqrt{2\pi}\sigma})^n e^{-\frac{1}{2}(\bold y- \bold {\hat y})^T(\sigma ^2 \boldsymbol I)^{-1}(\bold y-\bold {\hat y})} \end{align}$

对上式的概率分布求对数可得，

$\begin{align} \boldsymbol l(\boldsymbol \beta) &= -log((2\pi)^{\frac{n}{2}}\sigma^n)-\frac{1}{2\sigma^2}(\bold y-\bold {\hat y})^T(\bold y-\bold {\hat y}) \ &=--log((2\pi)^{\frac{n}{2}}\sigma^n)-\frac{1}{2\sigma^2}(\bold y - X\boldsymbol \beta)^T(\bold y -X\boldsymbol \beta) \end{align}$ （14）

上式中去除常数与 $\boldsymbol \beta$ 相关的表达式为：

$-\frac{1}{2\sigma^2}(\bold y - X\boldsymbol \beta)^T(\bold y -X\boldsymbol \beta)$

与式（10）相比，只相差一个常数 $-\frac{1}{2\sigma^2}$ ，因此对式（14）进行求导 $\frac{\partial \boldsymbol l}{\partial \boldsymbol \beta}$ 的过程与对（10）式求导过程相同，因此会得出与式（11）相同的结果。即在线性回归情形下，利用最大似然概率和高斯分布假设，求得 $\boldsymbol \beta$ 的结果与最小二乘法的结果相同。

通过上述的推导过程可以看出，在线性回归中，最大似然概率和最小二乘法等价，基于不同的角度得出了相同的结果。最大似然概率，是从概率角度出发，最小二乘法是从误差平方和最小作为出发点。

最小二乘法

最小二乘法的目标为：

$\begin{align} arg\ min_{\boldsymbol \beta \in R^{m+1}}\ S & =\sum_{i=1}^n(f(x_i, \boldsymbol \beta)-y_i)^2 \ & = \sum_{i=1}^n r_i^2 \end{align}$

其中 $r_i = f(x_i, \boldsymbol \beta)-y_i$ ， (x_i, y_i) 对应第 i 个样本， $f(x_i, \boldsymbol \beta)$ 为需要拟合的函数。

要使上式最小，需要满足 $\frac{\partial S}{\partial \beta_i}=0$ ，即

$\begin{align} \frac{\partial S}{\partial \beta_j} &= \sum_{i=1}^n 2r_i \frac{\partial r_i}{\partial \beta_j} \ &=\sum_{i=1}^n 2r_i\frac{\partial f(x_i, \boldsymbol \beta)}{\partial \beta_j}\ &=0 \end{align}$

即为：

$\sum_{i=1}^n 2r_i\frac{\partial f(x_i, \boldsymbol \beta)}{\partial \beta_j} = 0$ (15)

通过求解上述方程组，得到参数 $\boldsymbol \beta$ 的值。

根据模型函数 $f(x, \boldsymbol \beta)$ 的形式，可以将最小二乘法分为线性最小二乘法和非线性最小二乘法。

线性最小二乘法

在线性回归模型中，如式子（7）所示，

$\bold {\hat{y}} = X\boldsymbol \beta$

对于分量 $\hat{y_i}$ 而言：

$\hat{y_i} = \beta_0+X_{i1}\beta_1+X_{i2}\beta_2+...+X_{im}\beta_m$

即在线性回归模型中，是参数 $\beta_0, \beta_1, ..., \beta_m$ 的线性组合表达式。

参数的线性组合表达式，写成更一般的形式为：

$f(\bold x, \boldsymbol \beta) = \sum_{i=0}^m \beta_i \phi(\bold x)$

更一般的讲，只要 $\frac{\partial f}{\partial \beta_i }$ 值为常数或者非 $\boldsymbol \beta$ 的相关变量，即为参数 $\boldsymbol \beta$ 的线性函数。

对于线性最小二乘法，求解的过程类似上述线性回归模型的求解，往往有解析解（closed-form solution)。

非线性最小二乘法

非线性最小二乘法，往往对应 $f(x, \boldsymbol \beta)$ 不能表示为 $\boldsymbol \beta$ 的线性表达式，对于这种情况往往没有解析解（closed-form solution）。通过迭代的方式，求解 $\boldsymbol \beta$ 的最优解。

$\beta_{j}^{k+1} = \beta_{j}^{k} + \Delta \beta_j$

其中 $\Delta \beta_j$ 的更新通过下述方法求得。

一般的做法是对 $f(x, \boldsymbol \beta)$ 做一阶泰勒展开，

$\begin{align} f(x_i, \boldsymbol \beta)&=f^k(x_i, \boldsymbol \beta)+\sum_{j=1}^{m+1} \frac{\partial f(x_i, \boldsymbol \beta)}{\partial \beta_j}(\beta_j -\beta_j^k)\ &=f^k(x_i, \boldsymbol\beta) + \sum_{j}J_{ij}(\beta_j-\beta_j^k) \end{align}$

可以得出

$r_i = f^k(x_i, \boldsymbol \beta) + \sum_{j}J_{ij}(\beta_j -\beta_j^k)-y_i$

代入（15）式，可得

$\begin{align} \sum_{i}r_i \frac{\partial r_i}{\partial \beta_j} & =\sum_i r_iJ_{ij} \ &= \sum_i (f(x_i, \boldsymbol \beta)-y_i+\sum_jJ_{ij}(\beta_j - \beta_j^k))J_{ij} \ & = 0 \end{align}$

记 $\Delta \beta_j = \beta_j -\beta_j^k$ ， $\Delta y_i = f(x_i, \boldsymbol\beta)-y_i$ ,则上式可以写为：

$\sum_i \Delta y_i J_{ij} = -\sum_i\sum_t J_{ij} J_{it}\Delta \beta_t$

其中 j = (1, 2, ..., m+1) ，即上述对应m+1个方程，写为矩阵形式为：

$\bold J^T \Delta \bold y = -(\bold J^T \bold J)\Delta \boldsymbol \beta$

上述方程组称为 Gauss–Newton algorithm 的定义方程，通过上述方程即可求出 $\Delta \boldsymbol \beta$ ，即可以利用迭代法求出参数 $\boldsymbol \beta$ 。

其中 $\bold J$ 矩阵。

线性最小二乘法 vs 非线性最小二乘法

对比上述线性和非线性最小二乘法的求解过程，可知对于非线性最小二乘法，需要给出初始值然后利用迭代法进行求解。

无论是线性还是非线性，都需要计算Jacobian矩阵，如果Jacobian矩阵没有解析解，需要进行数值求解。

对于线性最小二乘法，属于凹函数，不存在不收敛问题，而且值为全局最小且唯一。

对于非线性最小二乘法，会存在不收敛问题，而且会存在多个极小值。

梯度下降算法

如上所述，最小二乘法是以残差平方和最小为目标，通过使得函数对参数的一阶导数为0来满足目标，进而求得参数的一种方法。对于线性最小二乘法，可以通过矩阵解法求得参数，对于非线性最小二乘法，可以利用迭代法求解。

梯度下降算法是机器学习和深度学习中，用的最广的求解参数的方法，是一种迭代算法。最小二乘法，用机器学习的语言来讲，目标函数（或者说损失函数）为残差平方和最小。相比最小二乘法，梯度下降算法没有以残差平方和最小为目标的假设，可以适用于任何目标函数。相同之处在于二者都需要计算目标函数的一阶导数（Jocabi矩阵）。

梯度下降算法的求解过程如下：

第 k+1 k+1 步的解为：

$\beta_i^{k+1} =\beta_i^k -\alpha \frac{\partial L(\boldsymbol \beta)}{\partial \beta_i}|_{\boldsymbol \beta=\boldsymbol \beta_k}$

其中 $\boldsymbol \beta$ ，为损失函数， $\alpha$ 为超参数学习率，控制更新步长大小。

对于线性回归问题，既可以利用最小二乘法求解，也可以利用梯度下降法求解。利用最小二乘法，需要求解矩阵的逆矩阵，计算相当耗时（尤其是参数较多时），比较而言，利用梯度下降算法计算量并不是很大，反而是一种更常用更好的选择。

牛顿法

如上所述，梯度下降算法是利用一阶导数求解参数值的方法。参数更新方式可以写为：

$\Delta \beta_i =-\alpha \frac{\partial L}{\partial \beta_i}$ i = 1, 2, ..., m+1

写为向量形式为：

$\Delta \boldsymbol \beta = - \alpha L^{'}(\boldsymbol \beta^{k})$ （16）

其中k代表第k步。

牛顿法同时利用了一阶导数和二阶导数信息，参数更新方式如下：

$\Delta \boldsymbol \beta = -\alpha H(\boldsymbol \beta^k)^{-1} L^{'}(\boldsymbol \beta^k)$ (17)

其中 $H(\boldsymbol \beta^k) = L^{''}(\boldsymbol \beta^k)$ ，称为Hessian矩阵

比较式（17）和式（16），牛顿法是二阶收敛，梯度算法是一阶收敛，牛顿法相比梯度算法收敛速度更快，但是计算也会更加复杂。梯度算法只需要函数一阶可导，牛顿法需要函数二阶可导。尽管牛顿法的收敛速度更快，但是计算复杂，且需要函数二阶可导的限制，使得牛顿法不如梯度下降法运用的更加广泛。

牛顿算法的推导过程

对函数 $L(\boldsymbol \beta)$ 做二阶泰勒展开，

$L(\boldsymbol \beta)\approx L(\boldsymbol \beta^k) +(\boldsymbol \beta -\boldsymbol \beta^k)^TL^{'}(\boldsymbol \beta^k)+\frac{1}{2}(\boldsymbol \beta - \boldsymbol \beta^k)^T H(\boldsymbol \beta^k)(\boldsymbol \beta - \boldsymbol \beta^k)$

对上式做微分可得：

$\begin{align} dL(\boldsymbol \beta) &= (d{\boldsymbol \beta})^TL^{'}(\boldsymbol \beta^k)+\frac{1}{2}(d\boldsymbol \beta)^T H(\boldsymbol \beta^k)(\boldsymbol \beta-\boldsymbol \beta^k)+\frac{1}{2}(\boldsymbol \beta - \boldsymbol \beta^k)^T H(\boldsymbol \beta^k)(d\boldsymbol \beta) \ &=(L^{'}(\boldsymbol \beta^k) )^T(d\boldsymbol \beta) +(H(\boldsymbol \beta^k) (\boldsymbol \beta - \boldsymbol \beta^k))^T(d\boldsymbol \beta) \ & = (L^{'}(\boldsymbol \beta^k) +H(\boldsymbol \beta^k)(\boldsymbol \beta - \boldsymbol \beta^k))^T(d\boldsymbol \beta) \end{align}$