最小二乘法(least sqaure method)

最小二乘法（least sqauremethod）专栏文章汇总文章结构如下：1：最小二乘法的原理与要解决的问题2 ：最小二乘法的矩阵法解法3：最小二乘法的几何解释4：最小二乘法的局限性和适用场景5：案例python实现6：参考文献1：最小二乘法的原理与要解决的问题最小二乘法是由勒让德在19世纪发现的，形式如下式：标函数 = \sum（观测值-理论值）^2\\观测值就是我们的多组样本，理论值就是我们的假设拟合函数。

目标函数也就是在机器学习中常说的损失函数，我们的目标是得到使目标函数最小化时候的拟合函数的模型。

举一个最简单的线性回归的简单例子，比如我们有 m 个只有一个特征的样本： (x_i, y_i)(i=1, 2, 3...,m)样本采用一般的 h_{\theta}(x) 为 n 次的多项式拟合，h_{\theta}(x)=\theta_0+\theta_1x+\theta_2x^2+...\theta_nx^n,\theta(\theta_0,\theta_1,\theta_2,...,\theta_n) 为参数最小二乘法就是要找到一组\theta(\theta_0,\theta_1,\theta_2,...,\theta_n) 使得\sum_{i=1}^n(h_{\theta}(x_i)-y_i)^2 (残差平方和) 最小，即，求 min\sum_{i=1}^n(h_{\theta}(x_i)-y_i)^22 ：最小二乘法的矩阵法解法最小二乘法的代数法解法就是对 \theta_i 求偏导数，令偏导数为0，再解方程组，得到 \theta_i 。

矩阵法比代数法要简洁，下面主要讲解下矩阵法解法，这里用多元线性回归例子来描：假设函数h_{\theta}(x_1,x_2,...x_n)=\theta_0+\theta_1x_1+...+\t heta_nx_n 的矩阵表达方式为：h_{\theta}(\mathbf{x})=\mathbf{X}\theta\\其中，假设函数 h_{\theta}(\mathbf{x})=\mathbf{X}\theta 为 m\times1 的向量, \theta 为 n\times1 的向量，里面有 n 个代数法的模型参数。

X 为 m\times n 维的矩阵。

m 代表样本的个数， n 代表样本的特征数。

损失函数定义为 J(\theta)=\frac{1}{2}(\mathbf{X}\theta-\mathbf{Y})^T(\mathbf{X}\theta-\mathbf{Y}) ，其中\mathbf{Y} 是样本的输出向量，维度为 m\times 1 。

\frac{1}{2} 在这主要是为了求导后系数为1，方便计算。

根据最小二乘法的原理，我们要对这个损失函数对 \theta 向量求导取0。

结果如下式：\frac{\partial }{\partial\theta}J(\theta)=\mathbf{X}^T(\mathbf{X}\theta-\mathbf{Y})=0\\对上述求导等式整理后可得：\theta=(\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\mathbf{Y}\\3：最小二乘法的几何解释先说结论:最小二乘法的几何意义是高维空间中的一个向量在低维子空间中的投影。

考虑这样一个简单的问题，求解二元一次方程组：\left\{\begin{matrix} x_1+x_2=3\leftarrow a\\ -x_1+x_2=1\leftarrow b \end{matrix}\right.\\方程组的解也就是直线$a$与$b$的交点，并且很容易算出 x_1=1，x_2=2 .它的矩形形式：\begin{bmatrix}1\\ -1\end{bmatrix}\timesx_1+\begin{bmatrix}1\\ 1\end{bmatrix}\timesx_2=b\Leftrightarrow a_1\times x_1+a_2\times x_2=b\\表示 x_1 倍的向量 a_1 加上 x_2 倍的向量 a_2 等于向量b 。

或者说， b 是向量 a_1 与 a_2 的线性组合。

可以看到，1倍的 a_1 加上2倍的 a_2 既是 b ，而1和2正是我们的解。

而最小二乘所面临的问题远不止两个点，拿三个点来说吧。

（0,2）,（1,2）,（2,3）假设我们要找到一条直线 y=kx+b 穿过这三个点（虽然不可能），为表述方便，用 x_1 代替 k ， x_2 代替 b ：\left\{\begin{matrix}1\times k +b=2\\ 0\times k +b=2\\ 2\times k +b=3\end{matrix}\right.\Leftrightarrow\left\{\begin{matrix}1\times x_1 +x_2=2\\ 0\times x_1 +x_2=2\\ 2\times x_1+x_2=3\end{matrix}\right.\Leftrightarrow\begin{bmatrix}1 &1 \\ 0 &1 \\ 2 &1\end{bmatrix}\begin{bmatrix} x_1\\x_2\end{bmatrix}=\begin{bmatrix}2\\ 2\\3\end{bmatrix}\Leftrightarrow A\times x=b\\进一步的：\begin{bmatrix}1\\ 0\\ 2\end{bmatrix}\timesx_1+\begin{bmatrix}1\\ 1\\ 1 \end{bmatrix}\timesx_2=\begin{bmatrix}2\\2\\3\end{bmatrix}\Leftrightarrow a_1\times x_1 + a_2\times x_2=b\\向量 b 是向量 a_1 与a_2 的线性表示。

用图形表示：作图之后，我们惊讶的发现，无论我们怎样更改 a_1 和 a_2 的系数都不可能得到b，因为 a_1 与 a_2 的线性组合成的向量只能落在它们组成的子空间S里面，也就是说，向量 b 不在平面 S 上，虽然我们找不到这样的向量，但在 S 上找一个比较接近的可以吧。

很自然的想法就是将向量 b 投影在平面S 上，投影在 S 上的向量 P 就是 b 的近似向量，并且方程A\hat{x}=P $是有解的。

这个误差最小的时候就是 e 正交与平面 S ，也正交与 S 中的向量 a_1，a_2 （矩阵 A 的列向量），即点乘为0，a_1^Te=0 ， a_2^Te=0 矩阵表示：A^Te=0\\A^T(b-A\hat{x})=0\\A^TA\hat{x}=A^Tb\\所以，我们可以得出，它的几何意义就是高维空间中的一个向量在低维子空间上的投影。

4：最小二乘法的局限性和适用场景从上面可以看出，最小二乘法简单高效，似乎比梯度下降迭代法方便很多。

但是这里我们就来说说最小二乘法的局限性。

首先，最小二乘法需要计算 X^TX 的逆矩阵，有可能它的逆矩阵不存在，这样就没有办法直接用最小二乘法了，此时梯度下降法仍然可以使用。

当然，我们可以通过对样本数据进行整理，去掉冗余特征。

让 X^TX 的行列式不为0，然后继续使用最小二乘法。

第二，当样本特征 n 非常的大的时候，计算 X^TX 的逆矩阵是一个非常耗时的工作（ n\times n 的矩阵求逆），甚至不可行。

此时以梯度下降为代表的迭代法仍然可以使用。

那这个n 到底多大就不适合最小二乘法呢？如果你没有很多的分布式大数据计算资源，建议超过10000个特征就用迭代法吧。

或者通过主成分分析降低特征的维度后再用最小二乘法。

再次，如果拟合函数不是线性的，此时不能使用最小二乘法，需要通过一些技巧转换成线性。

这个时候还是可以用梯度下降。

5：案例python实现举例：我们用目标函数 y=sin2{\pi}x , 加上一个正太分布的噪音干扰，用多项式去拟合【《统计学习方法》例1.1 11页】importnumpyasnpimportscipyasspfromscipy.optimizeimport leastsqimportmatplotlib.pyplotasplt%matplotlibinline# 目标函数defreal_func(x):returnnp.sin(2*np.pi*x)# 多项式# ps: numpy.poly1d([1,2,3]) 生成$1x^2+2x^1+3x^0$*deffit_func(p,x):f=np.poly1d(p)return f(x)# 残差defresiduals_func(p,x,y):ret=fit_func(p,x)-yreturnret# 十个点x=np.linspace(0,1,10)x_points=np.linspace(0,1,1000)# 加上正态分布噪音的目标函数的值y_=real_func(x)y=[np.random.normal(0,0.1)+y1fory1iny_] deffitting(M=0):"""n 为多项式的次数"""# 随机初始化多项式参数p_init=np.random.rand(M+1)# 最小二乘法p_lsq=leastsq(residuals_func,p_init,args=(x,y))print(' Fitting Parameters:',p_lsq[0])# 可视化plt.plot(x_points,real_func(x_points),label='real')plt .plot(x_points,fit_func(p_lsq[0],x_points),label='fitt edcurve')plt.plot(x,y,'bo',label='noise')plt.legend()ret urnp_lsq# M=0p_lsq_0=fitting(M=0)# M=1p_lsq_1=fitting(M=1)# M=3p_lsq_3=fitting(M=3)# M=9p_lsq_9=fitting(M=9)当M=9时，多项式曲线通过了每个数据点，但是造成了过拟合正则化结果显示过拟合，引入正则化项(regularizer)，降低过拟合Q(x)=\sum_{i=1}^n(h(x_i)-y_i)^2+\lambda||w||^2\\回归问题中，损失函数是平方损失，正则化可以是参数向量的L2范数,也可以是L1范数。