最小二乘拟合

本博客是说明最小二乘拟合，在李航博士的《统计学习方法》中，在第11页的例1.1作为示例进行了演示。但是，文中所给出的答案明显是错误的。因此，在这里，我们试图给出正确的解法，以及利用Python语言来进行实验验证。

一、问题定义

给定一个数据训练集：
$T = \{ (x_{1},y_{2}), (x_{2},y_{2}),\cdots, (x_{N},y_{N}) \}$
其中，$x_{i}\in R$是输入$\textbf{x}$的观测值，$y_{i}\in R$是输入$\textbf{y}$的观测值，$i=1,2,\cdots,N$。多项式函数拟合的任务是假设给定数据由$M$次多项式函数生成，选择出最有可能产生这些数据的$M$次多项式函数。
设$M$次多项式为

$$f_{M}(x,\textbf{w}) = w_{0}+w_{1}x^{1}+w_{2}x^{2}+\cdots+w_{M}x^{M}=\sum_{j=0}^{M} w_{j}x^{j}$$
其中，$x$是单变量输入，$\textbf{w} = \{ w_{0},w_{1},\cdots,w_{M} \}$是$M+1$个参数。

目标函数

$$\textbf{w} = \arg\min L(\textbf{x} ,\textbf{w} ) = \frac{1}{2}*\sum_{i=1}^{N}(\sum_{j=0}^{M} w_{j}x^{j}_{i}-y_{i})^{2}$$

二、求解过程

将$(x_{i},y_{i})$看做已知值，为了求得$L(\textbf{x} ,\textbf{w} )$的最小值，对每一个$w_{j}$求偏导，即：

$$\frac{\partial L}{\partial w_{j}} = \frac{1}{2}*\sum_{i=1}^{N} [2*(\sum_{j=0}^{M} w_{j}x^{j}_{i}-y_{i})*x_{i}^{j}]\\ =\sum_{i=1}^{N}(x_{i}^{j}*\sum_{j=0}^{M}w_{j}x_{i}^{j})-\sum_{i=1}^{N}y_{i}x_{i}^{j}$$

这里需要注意的是，左式$x_{i}^{j}*\sum_{j=0}^{M}w_{j}x_{i}^{j}$并不能直接计算为$\sum_{j=0}^{M}w_{j}x_{i}^{2j}$，因为前面的$j$和后面的$j$并不是同一个概念，不能混淆。为了区分这一点，我们用矩阵来代表左式：

$$\begin{pmatrix} x_{1}^{j} & x_{2}^{j} & \cdots & x_{N}^{j} \end{pmatrix} \begin{pmatrix} x_{1}^{0} & x_{1}^{1} & \cdots & x_{1}^{M}\\ x_{2}^{0} & x_{2}^{1} & \cdots & x_{2}^{M}\\ \vdots & \vdots & \vdots & \vdots \\ x_{N}^{0} & x_{N}^{1} & \cdots & x_{N}^{M} \end{pmatrix}\begin{pmatrix} w_{0} \\ w_{1} \\ \vdots \\w_{M} \end{pmatrix}$$

则，该矩阵可以转换为

$$\begin{pmatrix} x_{1}^{j}+x_{2}^{j}+\cdots+x_{N}^{j}\\ x_{1}^{j+1}+x_{2}^{j+1}+\cdots+x_{N}^{j+1}\\ \vdots\\ x_{1}^{j+M}+x_{2}^{j+M}+\cdots+x_{N}^{j+M} \end{pmatrix}^{T}\begin{pmatrix} w_{0} \\ w_{1} \\ \vdots \\w_{M} \end{pmatrix}$$

按照求极值的步骤，先求偏导，令其为0，解出$\textbf{w}$。即，对每个$w_{j}$求偏导，令其为0，总共得到$M$个$M$元线性方程，求解该方程，得出$\textbf{w}$。

$$\begin{pmatrix} \sum_{i=1}^{N}x_{i}^{0} & \sum_{i=1}^{N}x_{i}^{1} &\cdots & \sum_{i=1}^{N}x_{i}^{M} \\ \sum_{i=1}^{N}x_{i}^{1} & \sum_{i=1}^{N}x_{i}^{2} &\cdots & \sum_{i=1}^{N}x_{i}^{M+1} \\ \vdots & \vdots & \vdots & \vdots \\ \sum_{i=1}^{N}x_{i}^{M} & \sum_{i=1}^{N}x_{i}^{M+1} &\cdots & \sum_{i=1}^{N}x_{i}^{2M} \end{pmatrix}\begin{pmatrix} w_{0} \\ w_{1} \\ \vdots \\w_{M} \end{pmatrix} = \begin{pmatrix} \sum_{i=1}^{N}y_{i}x_{i}^{0} \\ \sum_{i=1}^{N}y_{i}x_{i}^{1} \\ \vdots \\ \sum_{i=1}^{N}y_{i}x_{i}^{M} \end{pmatrix}$$

三、实验

下面就是利用python来做实验了。

实验代码

import random
import matplotlib.pyplot as plt
import numpy as np

#get 10 randomized points
x = range(10)
y = np.zeros(10)

for w in range(10) :   
   y[w] = random.randint(5,10)

z = np.zeros(5) 

#first get Y(10)
for i in range(5):
    for j in range(10):
        z[i] = z[i] + y[j]*x[j]**(i)


matrix = np.zeros((5,5))
for i in range(5):
    for j in range(5):
        for k in range(10):
            matrix[i][j] = matrix[i][j] + x[k]**(i+j)


out = np.linalg.solve(matrix,z)

output = np.zeros(10)

for i in range(10):
    for j in range(5):
        output[i] = output[i] + out[j]*x[i]**(j)


#display the points and curve
plt.plot(range(10),y,'bo',range(10),output,'r-')
plt.legend()
plt.show()

在上面的代码中，我们随机生成10个点，然后生成 $M=4$ 的多项式。

实验结果

Talking is cheap, show me the code. YOU KNOW!!!

吴小同
2016-5-16
2016-5-17