玩不转scikit-learn（一）：LinearRegression

为了能更好的掌握机器学习，自己边学习边自己写一写算法。

为了验证自己算法实现的准确性，我会和scikit-learn的对应算法进行对比。恩，水平有限，尽力分析下为什么scikit-learn的对应算法会比自己实现的朴素算法好。

希望可以对于算法本身和scikit-learn的实现有一定的深度理解。。希望吧

代码用Python3在jupyter notebook(ipython notebook)上写的。代码+图的完整实现我放在GitHub上了：地址如下

https://github.com/lixintong1992/Machine_Learning

求求各位在GitHub上star我。。。

---

线性回归

用最简单的最小二乘（Least Squares）作为开始吧

线性回归形式

回归就是预测一个数，线性回归就是用输出值的线性组合预测这个数，公式如下：

如果只有一个输入$w$，那么再加上一个偏置 $b$ ， $x_i$ 的回归结果$f(x_i)$如下

$$f(x_i)=wx_i+b$$

这就是线性回归的最简单形式，如何衡量预测结果的准确性？我们需要一个loss function来定义预测结果的准确性。

---

最小二乘

最直观的衡量方法，就是预测值与实际值的差的平方。这个就是最小二乘，形式如下

$$(f(x_i)-y_i)^2$$
机器学习的过程就是通过调整参数减小loss function来确定参数。

---

求解

给你$m$个训练样本，loss function如下所示：

$$E_{(w,b)}=\sum_{i=1}^{m}(y_i-wx_i-b)^2$$

如何通过调整$w,b$的值来求$E_{(w,b)}$的最小值，分别对$w,b$求偏导，当导数为0的时候就是啦（因为是凸函数的说，导数为0就是最小值了）

$$\frac{\partial{E_{(w,b)}}}{\partial{w}}=2\left({w\sum_{i=1}^{m}{x^2_i}-\sum_{i=1}^{m}(y_i-b)x_i}\right)=0$$

$$\frac{\partial{E_{(w,b)}}}{\partial{w}}=2\left(mb-\sum_{i=1}^m(y_i-wx_i)\right)=0$$

结果是这样子的

$$w=\frac{\sum_{i=1}^{m}y_i(x_i-\overline{x})}{\sum_{i=1}^{m}x_i^2-\frac1m(\sum_{i=1}^mx_i)^2}$$

$$b=\frac1m\sum_{i=1}^m(y_i-wx_i)$$

$$\overline{x}=\frac1m\sum_{i=1}^mx_i$$

很难得的解，因为他是闭式(closed-form)解，就是可以直接用公式求出来的解。以上是只有一个输入变量情况下的线性回归，下面考虑下多面量输入，嘛，也向量化（vectorized）下。

---

多元线性回归

如果是输入$d$的变量，预测一个变量，公式如下

$$f(x)=w_1x_1+w_2x_2+w_3x_3+...+w_dx_d + b$$

$$x=(x_1,x_2,...,x_d)$$

都向量化了外面飘一个$b$肯定不合适，我们把$b$放进去：$\hat{w}=(w;b)$

同时，我们用$X$表示有$m$个样本的训练集

$$X= \begin{pmatrix} x_11 & x_12 &\cdots& x_1d&1 \\ x_21 & x_22 &\cdots& x_2d&1 \\ \vdots & \vdots &\cdots& \ddots &\vdots \\ x_m1 & x_m2 &\cdots& x_md&1 \\ \end{pmatrix}$$
这样子损失函数就是
$$E_\hat{w}=(y-X\hat{w})^T(y-X\hat{w})$$

老样子求导为0

$$\frac{\partial{E_\hat{w}}}{\partial{\hat{w}}}=2X^T(X\hat{w}-y)$$

$$\frac{\partial{E_\hat{w}}}{\partial{\hat{w}}}=0$$

解下面的式子就成了

$$X^TX\hat{w}=X^Ty$$

结果如下

$$f(\widehat{x_i})=\widehat{x_i}^T(X^TX)^{-1}X^Ty$$

其实也就是求

$$(X^TX)^{-1}X^T$$

---

实现

朴素人工版

X = np.insert(X,X.shape[1],values=1,axis=1)  #For constant term
y = y[:,np.newaxis]

这个是在X中加入1的列，以及把y的维数搞好

if np.linalg.det(X.T.dot(X)) == 0:
    print("singular error!")
else:
    w = np.linalg.inv(X.T.dot(X)).dot(X.T).dot(y)

因为要求$(X^TX)^{-1}X^T$，要求逆矩阵，这里通过np.linalg.det求行列式，行列式为0的矩阵没有逆矩阵

np.linalg.inv就可以求逆矩阵了

结果如下（完整代码https://github.com/lixintong1992/Machine_Learning）

scikit-learn高端版

scikit-learn算法实现

model = linear_model.LinearRegression()
model.fit(X, y)

调包真是安逸啊，就两句话。结果如下

一样的，成功~

scikit-learn实现分析

其实scikit-learn是通过scipy.linalg.lstsq来实现的

因为逆矩阵有可能不存在，一般都用SVD奇异值分解，对不起大家这个我不太懂······不细说了
通过上面的代码发现scipy.linalg.lstsq会返回四个值，第一个就是解啦。后面的依次是Sums of residues;Effective rank of matrix;Singular values of matrix
有需要可以直接查看，如下图：

需要注意的是Sums of residues要被移除了，反正都很少用到就是了

这里有个可有可无的小知识:
numpy.linalg.lstsq 和 scipy.linalg.lstsq是不一样的(Numpy 1.8.2, Scipy 0.14.1 )
numpy.linalg.lstsq() 用的是LAPACK 惯例的xGELSD ，scipy.linalg.lstsq()用的是xGELSS，翻译成人话就是Numpy的速度更快但需要更多内存

---

LinearRegression分析

scikit-learn上的相关叙述

普通最小二乘的系数估计依赖于模型的独立性（就是特征间的独立性）。如果特征之间是相关的，就是说$X$的列与列之间近似线性相关，那么矩阵就会近似奇异(singular)。那么，最小二乘估计会对随机噪音敏感，就是模型的方差会很大(variance)。这个特征之间是相关的情况叫多元共线性(multicollinearity)

通俗的理解：当线性回归模型中存在多个相关变量时，它们的系数确定性变差并呈现高方差。比如说，在一个变量上的一个很大的正系数可能被在其相关变量上的类似大小的负系数抵消

上边的方差，应该就是机器学习模型中的泛化误差里边的

$$泛化误差=bias^2+var+\epsilon$$

$bias$就是偏差，$var$是方差，$\epsilon$是噪声

补充下，多重共线性使参数估计值的方差增大，可以用杜宾-瓦特森（DW）检验，计量经济，统计分析中常用的一种检验序列一阶 自相关 最常用的方法。

算法复杂度主要集中在奇异值分解(SVD),$X$为大小为$(n,p)$的矩阵，复杂度是$O(np^2)$,$n \ge p$

---

结尾

https://github.com/lixintong1992/Machine_Learning

GitHub上求star~写blog不容易呀~

公式大部分来自《机器学习》周志华