MENU

Numerical Analysis: 最小二乘原理

• September 20, 2019 • Read: 3404 • Knowledge,Numerical Analysis

数据的拟合、回归分析最为常用的一个原理就是最小二乘原理,这一原理不仅在数据的处理中得到广泛的应用,而且在深度学习中,线性回归中OLS损失函数是最为常用的模型,这里的OLS就是最小二乘原理。最小二乘原理应用的目的是为了控制最优化问题,即以什么样的标准才能确定问题的最优解。对于这一原理(方法)的认识,更为有价值的是为什么要研究这个原理,它能够解决什么样的问题,从这一角度去理解学习,可以避免枯燥的数学公式推导以及证明。

现实中遇到了什么问题?

现实中,对于一个理工科学生来说,在实验室进行实验研究是很重要的科学研究手段,而实验研究过程数据的记录以及对于数据的处理分析是十分重要的环节。假设,我做了一个实验,得到了一系列的数据(二维的数据,x和y值),可能几百个(也可能上千个),假设然后我发现这些数据的分布接近线性,也就是说x与y之间的变化规律比较符合线性关系,那么我自然会依据数学函数的相关知识,将数据的模型建立为$y = a x + c$,其中$a$和$b$这个一次函数的系数,它们都是未知的;既然我已经观测到了那么多的数据点,将每一个数据点都带入事先设定的这个线性模型中,一定会得到很多个方程(有多少点就有多少个方程),对于我而言,只需要确定两个未知的系数即可,即$a$和$c$;但问题在于,现在方程的数量远远大于未知数的数量,可以取任意两个方程组成方程组进行求解得到一组结果,而这种结果可以有很多个,那么究竟选择哪一个才最合适?我想这个结果,一定是要最优的,最能反映出每一个方程的要求的,这就需要利用最小二乘原理来进行解答。

总之,现实中我遇到的问题,现在转为了一个方程组求解的问题,并且具体情况是这样的:得到的方程组中方程的个数远远大于未知数的个数

从向量空间来看待最小二乘原理

前述已经阐明了面临的问题,以及这个问题具体的情形是什么。不从线性代数中线性方程组解的特性来理解,因为这样往往过于抽象和公理化,很难透彻地理解最小二乘原理的意义。这里从向量空间的角度去理解,因为从根本上来说,线性代数研究的本质就是研究$N$维向量空间中的变换问题

前面已经说了,这个问题依然是一个线性方程组求解的问题,那么自然可以用以下的式子来表示:

$$ A x = b \tag{1} $$

比较特别的是,这里的矩阵$A$不是一个方阵,它的行数实际上就是我们实测的数据点的个数,由于我们之前的模型为线性模型,只有两个待确定的系数,因此这个矩阵的列数只有2列。$x$就是我们要求的未知系数,此时它是一个长度为2列向量,具有两个元素即$a$和$c$;向量$b$就是实测数据点的$y$值。这里不用纠结用什么变量表示,因为问题的本质特性不会因为变量名字的改变而发生变化。为了更好地从向量空间的角度理解,我可以将式(1),依据矩阵乘法的运算法则进行展开,得到:

$$ \begin{bmatrix} a_{11}\\ a_{21}\\ a_{31}\\ \vdots\\ a_{n1} \end{bmatrix} x_{1} + \begin{bmatrix} a_{12}\\ a_{22}\\ a_{32}\\ \vdots\\ a_{n2} \end{bmatrix} x_{2} = \begin{bmatrix} b_{1}\\ b_{2}\\ b_{3}\\ \vdots\\ b_{n} \end{bmatrix} \tag{2} $$

从这个形式来看,实际上就是n维向量之间的线性关系的研究。上述形式中,第一个列向量和第二个列向量组成的n维向量空间进行线性运算后得到$b$向量。有一点需要说明的是,对于空间来说,人只能认识到三维空间,再高维度的空间无法被可视化,列向量中元素的个数表明了这个向量的空间维度,为了可视化表示这里假定图中的向量是$R^2$中的向量,如下图 。

最小二乘原理

这两个向量各自都是一条线,可以确定一个平面,那么这个平面中任何其他的向量可以用它们的线性组合来表示,也就是$Ax$,也必定是与他们处于同一平面上的(或者说是同以向量子空间中),而真实的列向量$b$肯定不在这个平面上,那么这个问题所要寻求的最佳结果,就是在这个平面上找到一个$b'$向量,它由$x'$来确定,即我们寻找的一个近似的解。使得其与向量$b$最接近!这一要求,自然是向量$b$在这个平面上的投影就是,这里向平面做一个垂线,由此得到投影$b'$向量;由于$Ax$与垂线$b - Ax'$垂直,利用向量的正交性质,内积必定为0,所以可以写为

$$ (Ax)^{T}(b - Ax') = 0 \tag{3} $$

转置运算后,得到

$$ x^{T}A^{T} \cdot (b - Ax') = 0 \tag{4} $$

由此得到,满足

$$ A^{T}Ax' = A^{T}b \tag{5} $$

即可,求解该方程组,得到最小二乘原理的解答。式(3)也被称为法线方程。

对于数据拟合的流程

从上述对于最小二乘原理的解释以及法线方程的意义,可以得出做数据拟合的基本步骤,此时不借助于数据分析软件,也能完成这一项工作。

首先,具有一系列的数据点$(t_1, y_1), (t_2, y_2),(t_3, y_3),\cdots,(t_n, y_n)$,共有$n$个点。

  1. 选择一个合适的数据模型。比如$y =a t + b$ 等。
  2. 强制数据点满足模型,即直接带入到模型中,则获得n个方程,形成$Ax=b$的形式。
  3. 利用法线方程进行求解。$A^{T}Ax = A^{T}b$

例子

利用随机数构造出拟合数据

设定我们理论上精确的模型关系式为$y = 1.25 x^2 + 0.7x + 2.1$。利用该式子构造出一些数据点,然后添加noise进行扰动,得到“实测”的数据点供拟合用

import numpy as np
import matplotlib.pyplot as plt

# x的范围设定在-2,2之间
x = np.linspace(-2, 2, 50)  # 50个x的值
# 定义出1.25 x^2 + 0.7x + 2.1这个原始的函数
def fun(x):
    return 1.25 * x**2 + 0.7 * x + 2.1

# 产生对应的y
y = fun(x)

# 产生干扰波动,利用随机数
wave = np.random.randn(y.shape[0]) * 0.7
# 构成拟合所需的ydata
ydata = y + wave

# 数据点的可视化
plt.style.use('ggplot')  # ggplot风格
fig,ax = plt.subplots()
ax.scatter(x, ydata, marker='o', label='data points')
ax.plot(x, y, label='real point')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_title('The points')
ax.legend(loc='best')
plt.show()

数据点

# 对数据进行拟合
# 1.依据数据图像,确定一个比较合适的模型
# 这里选用y = a x^2 + b x + c这个模型

# 2. 将所有数据点x,ydata带入,强制构造方程组,
# 主要任务是确定系数矩阵A
A1 = x**2
A2 = x
A3 = np.ones(x.shape[0])
A = np.c_[A1, A2, A3]
# b向量就是ydata
b = ydata
# 3. 获得法线方程,并求解
result = np.linalg.solve(np.dot(A.T, A), np.dot(A.T, b))
print(result)

结果为:

$$ \begin{bmatrix} 1.23244467 & 0.7380719 & 2.10100332 \end{bmatrix} $$

对拟合结果于原数据进行对比。

# 参数解包
a, b, c = result
# fitting result
y_fit = a * x**2 + b * x + c
plt.style.use('ggplot')  # ggplot风格
fig,ax = plt.subplots()
ax.scatter(x, ydata, marker='o', label='data points')
ax.plot(x, y, label='real point', color='tab:blue')
ax.plot(x, y_fit, linestyle='-.', color='tab:green')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_title('The points')
ax.legend(loc='best')
plt.show()

datas

在scipy库中,optimize模块的curve_fit函数可以进行数据的拟合,只不过,需要先构造一个含有未知数的函数,具体例子如下。

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

def func(x, a, b, c):
    return a * np.exp(-b * x) + c

xdata = np.linspace(0, 4, 50) 
y = func(xdata, 2.5, 1.3, 0.5) 
np.random.seed(1729) 
y_noise = 0.2 * np.random.normal(size=xdata.size)
ydata = y + y_noise
plt.style.use('ggplot')
plt.plot(xdata, ydata, color='tab:blue', label='data')
popt, pcov = curve_fit(func, xdata, ydata)
plt.plot(xdata, func(xdata, *popt), 'r-',
          label='fit: a=%5.3f, b=%5.3f, c=%5.3f' % tuple(popt))
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.show()

datas

如果你觉得这些对你还有些收获!当然可以打赏我。

请点击右下方Tip进行打赏!

Last Modified: October 6, 2019
Archives Tip
QR Code for this page
Tipping QR Code
Leave a Comment

已有 1 条评论
  1. 连特曼 连特曼

    本文要点:
    1.最小二乘法的背景是超定方程组的最优解问题;
    2.在低维向量空间中阐释了最小二乘法的几何意义,直观反映了最小二乘解的唯一性;
    3.利用随机数(离散/连续)构造拟合数据的实现。