1. 概述

最小二乘法Least Square Method，是分类回归算法的基础
主要思想
- 通过最小化误差的平方和寻找数据的最佳函数匹配
- 利用最小二乘法可以简便地求得未知的数据，并使得这些求得的数据与实际数据之间误差的平方和为最小

2. 最小二乘法原理

假设我们现在有一系列数据点$(x_i, y_i), i = 1, 2, \dots, m$，那么，给定拟合函数$h(x)$，我们可以得到拟合出来的数据点。但是，给出的拟合函数与实际待求解的函数的拟合程度比较高呢？因此，提出残差来衡量，计算公式如下：

$r_i = h(x_i) - y_i$，

其中，$r_i$为第$i$个数据点的残差，$y_i$为第$i$个数据点的真值，$h(x_i)$就是第$i$个数据点由拟合函数计算得到的数值。

评估拟合程度的主要思想就是真值和拟合函数计算的数据距离越小越相似。相似性计算方法有如下几种：

$\infty-$范数
- 残差绝对值的最大值$\max_{1\leq i\leq m}|r_i|$，即所有数据点中残差距离的最大值
$1-$范数
- 绝对残差和$\sum_{i=1}^{m}|r_i|$，即所有数据点残差距离之和
$2-$范数
- 残差平方和$\sum_{i=1}^{m}r_i^2$

由于前两种评价方式不好计算微分，因此，较为广泛的采用二范数。

由此，我们可以写出最小二乘法的定义了：

对于给定的数据$(x_i, y_i), i = 1, 2, \dots, m$，在取定的假设空间H中，求解$h(x)\in H$，使得残差$r_i = h(x_i) - y_i$的二范数最小，即

$h(x)\ s.t.\ \min\sum_{i=1}^m(h(x_i) - y_i)^2$.

从几何上讲，就是寻找与给定点$(x_i, y_i), i = 1, 2, \dots, m$距离平方和最小的曲线$y = h(x)$。$h(x)$称为拟合函数或者最小二乘解，求解拟合函数$h(x)$的方法称为曲线拟合的最小二乘法。

这里的$h(x)$一般情况下是一种多项式曲线，即$h(x, w) = w_0 + w_1x + \dots + w_nx^n$，其中$w$是参数。也就是说，最小二乘就是要找到一组$w$，使得$\sum_{i=1}^m(h(x_i) - y_i)^2$最小。

如何求解得到$w$呢？

下面以线性最小二乘举例，更高的非线性函数原理相同，即$h(x, w) = w_0 + w_1x$。优化目标为

$w = (w_0, w_1),\ s.t.\ \min\sum_{i=1}^m(w_0 + w_1x - y_i)^2$。

代数解法

令$Q(w) = \sum_{i=1}^m(w_0 + w_1x - y_i)^2$为样本的平方损失函数，这是一个典型的求解极值的问题，只需要分别对$w_0$和$w_1$求偏导数，令偏导数为0即可求出极值点，即

$\left\{\begin{array}{rc1} \frac{\partial Q}{\partial w_0} = 2\sum_{i=1}^m(w_0 + w_1x_i - y_i) = 0 \\ \frac{\partial Q}{\partial w_1} = 2\sum_{i=1}^m(w_0 + w_1x_i - y_i) = 0 \end{array}\right.$.

接下来只需要求解这个方程组。

矩阵解法

这里用多元线性回归例子来描。假设函数$h_{w}(x_1, x_2, \dots, x_n) = w_0 + w_1x_1 + \dots + w_nx_n$的矩阵表示为

$h_w(x) = Xw$，其中，$X$为$m\times n$为矩阵，即m个样本数和n个特征。

损失函数定义为：

$J(w) = \frac{1}{2}(Xw - Y)^T(Xw - Y)$，其中Y是样本的真实值。根据最小二乘法的原理，我们要对这个损失函数对w求导取0，即

$\frac{\partial J}{\partial w} = X^T(Xw - Y) = 0$，得$w = (X^TX)^{-1}X^TY$。

3. 代码实现

import numpy as np  # 引入numpy
import scipy as sp
import pylab as pl
from scipy.optimize import leastsq  # 引入最小二乘函数
 
n = 9  # 多项式次数
 
 
# 目标函数
def real_func(x):
    return np.sin(2 * np.pi * x)
 
 
# 多项式函数
def fit_func(p, x):
    f = np.poly1d(p)
    return f(x)
 
 
# 残差函数
def residuals_func(p, y, x):
    ret = fit_func(p, x) - y
    return ret
 
 
x = np.linspace(0, 1, 9)  # 随机选择9个点作为x
x_points = np.linspace(0, 1, 1000)  # 画图时需要的连续点
 
y0 = real_func(x)  # 目标函数
y1 = [np.random.normal(0, 0.1) + y for y in y0]  # 添加正太分布噪声后的函数
 
p_init = np.random.randn(n)  # 随机初始化多项式参数
 
plsq = leastsq(residuals_func, p_init, args=(y1, x))
 
print 'Fitting Parameters: ', plsq[0]  # 输出拟合参数
 
pl.plot(x_points, real_func(x_points), label='real')
pl.plot(x_points, fit_func(plsq[0], x_points), label='fitted curve')
pl.plot(x, y1, 'bo', label='with noise')
pl.legend()
pl.show()

4. 参考

最小二乘法及其python实现