最小二乘拟合直线真的没问题么?

一幅封面插画

偶然间我看到下图中的插画(下图)

一下子想起很多事情, 记得有次需要求解某个椭圆前后位置的旋转角, 已知了椭球边上点的坐标,于是我就想: 我用最小二乘不就可以得到椭圆的长轴了吗, 然后再由最小二乘拟合的直线斜率一下子就可以得到椭圆旋转的角度. 听起来似乎很完美, 结果最小二乘拟合后我得到上图中的实线(其实我想得到的是虚线). 另外, 在本科学习最小二乘拟合直线的时候, 一直存有疑惑为什么不是最小化点到直线的距离, 而是直接 $kx+b-y$ . 直到后来仔细分析过才发现里面其实大有乾坤! 下文将慢慢道来….

问题的提出&最小二乘的疑问

已知 N 个点数据坐标 $\left(\begin{array}{c} x_i \\ y_i \end{array}\right),~~i=1,\cdots,N$ 利用矩阵表示 $\boldsymbol{D}= \left[\begin{array}{cc} x_1 & y_1 \\ x_2 & y_2 \\ \vdots&\vdots\\ x_N & y_N \end{array}\right]$ ,

这些数据集大致呈现出一条直线的形态, 甚至这些数据本身就是呈一条直线, 只不过在测量、观察或计算过程中带入大量噪声, 使得这些数据并不严格呈现一条直线状. 那么如何找到一条直线 $y=kx+b$

使得这些数据点集和该直线最为”贴合“?

其实用泛泛的文字很难准确描述很多东西, 此时我们需要严格定义何为”贴合”. 首先“贴合”即点和直线的某种”距离”(范数)之和最小. Fig 2中给出了两种”距离”的定义: 左图为竖直距离, 右图为垂直距离. 显然最小二乘法是使得点和直线在”竖直”距离下最贴合, 然而不知大家在学习的时候有没有这个疑问, 明明是右边的”垂直”距离更直观有意义!

某一种解释: 实际问题中， $x$ 是没有误差的固定值，只有 $y$ 才是有误差的观测值，所以在只考虑 $x$ 的误差情况下, 竖直的距离更有意义一些.

但不知道有没有人想过这个问题, 如果严格按照Fig 2右图中的距离来拟合直线会怎样, 以及那又该怎么求? 毕竟没有了最小二乘那么方便的求解公式, 我将在下文中仔细解答, 首先回顾最小二乘方法, 再来给出究竟如何求解”垂直”距离下的拟合直线.

最小二乘(least square)拟合直线

最小二乘拟合直线应该是大学里必修的知识点, 不过很多同学应该都是为了考试背一下公式、做大作业网上Download一份MATLAB代码交差而已, 下文就帮助大家回忆一下究竟公式怎么求来的, 以及最小二乘拟合直线的一个彩蛋: 拟合的直线过数据点集的重心! 有兴趣可以自己先推到一下看看.

首先由上文所述, 我们定义”贴合“所采用的距离为”竖直”距离, 即如Fig 3所示, 即

$d_i=kx_i+b-y_i,\qquad~i=1,2,\cdots,N$

则直线拟合可写成如下优化问题:

优化问题: 找到斜率 $k$ 和截距 $b$

使得 $\min_{k,b}{(d_1^2+d_2^2+\cdots+d_N^2)}$ 成立.

由距离 $d_i$ 的定义

$d_1^2+d_2^2+\cdots+d_N^2 = (kx_1+b-y_1)^2+(kx_2+b-y_2)^2+\cdots+(kx_N+b-y_N)^2$

优化问题等价于

$k,b = arg\min_{k,b} \sum_{i=1}^{N}{d_i^2} = arg\min_{k,b} \sum_{i=1}^{N}{(kx_i+b-y_i)^2}$

换成矩阵描述形式, 并定义系数矩阵 $\boldsymbol{A}$

$\boldsymbol{d}= \left(\begin{array}{c} d_1 \\ d_2\\ \vdots\\ d_N \end{array}\right) = \left[\begin{array}{cc} x_1 & 1 \\ x_2 & 1 \\ \vdots&\vdots\\ x_N & 1 \end{array}\right] \left(\begin{array}{c} k \\ b \end{array}\right) -\left(\begin{array}{c} y_1 \\ y_2\\ \vdots\\ y_N \end{array}\right)=:\boldsymbol{A}\boldsymbol{k}-\boldsymbol{y}$

原优化问题又等价于 $k,b = arg\min_{\boldsymbol{k}} {\|\boldsymbol{d}\|} = arg\min_{\boldsymbol{k}}{\boldsymbol{d}^T\boldsymbol{d}} = arg\min_{\boldsymbol{k}}{(\boldsymbol{k}^T\boldsymbol{A}^T\boldsymbol{A}\boldsymbol{k} -2\boldsymbol{y}^T\boldsymbol{A}^T\boldsymbol{k}+\boldsymbol{y}^T\boldsymbol{y} )}$

求导可求得最优化问题等价于方程组

$\Leftrightarrow \boldsymbol{A}^T\boldsymbol{A}\boldsymbol{k}= \boldsymbol{A}^T\boldsymbol{y}$

其中左端矩阵和右端向量分别为

$\boldsymbol{A}^T\boldsymbol{A}=\left[\begin{array}{cc} x_1^2+\cdots+x_N^2 & x_1+\cdots+x_N\\ x_1+\cdots+x_N & 1+\cdots+1\\ \end{array}\right] =\left[\begin{array}{cc} \sum_{i=1}^{N}{x_i^2} & \sum_{i=1}^{N}{x_i} \\ \sum_{i=1}^{N}{x_i} & N\\ \end{array}\right]$

$\boldsymbol{A}^T\boldsymbol{y}=\left(\begin{array}{c} x_1y_1+\cdots+x_Ny_N \\ y_1+\cdots+y_N \\ \end{array}\right)=\left(\begin{array}{c} \sum_{i=1}^{N}{x_iy_i} \\ \sum_{i=1}^{N}{y_i} \\ \end{array}\right)$

最小二乘法拟合直线的待求方程组如下

$\left[\begin{array}{cc} \sum_{i=1}^{N}{x_i^2} & \sum_{i=1}^{N}{x_i} \\ \sum_{i=1}^{N}{x_i} & N\\ \end{array}\right]\left(\begin{array}{c} k \\ b \end{array}\right)=\left(\begin{array}{c} \sum_{i=1}^{N}{x_iy_i} \\ \sum_{i=1}^{N}{y_i} \\ \end{array}\right)$

有意思的事是上面方程组的第二个方程

$k \sum_{i=1}^{N}{x_i}+Nb= \sum_{i=1}^{N}{y_i} 即 \left(\frac{1}{N}\sum_{i=1}^{N}{x_i}\right)k+b= \frac{1}{N}\sum_{i=1}^{N}{y_i}$

换言之, 重心 $x_0=\frac{1}{N}\sum_{i=1}^{N}{x_i} , y_0=\frac{1}{N}\sum_{i=1}^{N}{y_i}$ 必定经过该拟合直线(彩蛋).

主成分分析(Principal Component Analysis)拟合直线

在”垂直”距离(Fig 2右图)下求解拟合直线, 我在下文中称为pca方法拟合直线, 至于为什么, 大家可以先看看pca相关的知识, 或者以后有空我在专栏里介绍.

由上一节最小二乘法的分析, 我们不妨假设数据集的重心就在原点, 或者做数据的平移使得该条件成立

$x_i=x_i-\frac{1}{N}\sum_{i=1}^{N}{x_i},\qquad y_i=y_i-\frac{1}{N}\sum_{i=1}^{N}{y_i},\qquad i=1,2,\cdots,N$

因此在拟合该过原点直线的时候, 我们仅仅需要考虑斜率 $k$ (如Fig 4.)

换个角度, 如果我们知道了该直线的单位法向量 $\boldsymbol{n}=(n_1,n_2)^T$ ,则斜率也就知道了, 由此问题转化为如何确定单位向量 $\boldsymbol{n}$ . 下一步就是如何求得各个数据点到假定直线之间的最短距离, 这也是为什么需要引入直线的法向量.

如Fig 4所示, 对于点 $(x_i,y_i)^T$ , 该点的坐标向量与直线法向量 $\boldsymbol{n}$ 的内积正好等于该点到直线的距离(同向为正, 反向为负), 即

$h_i = \left(\begin{array}{c} x_i \\ y_i \end{array}\right)\cdot \boldsymbol{n} =\left(\begin{array}{c} x_i \\ y_i \end{array}\right)\cdot \left(\begin{array}{c} n_i \\ n_i \end{array}\right)= \left(\begin{array}{cc} x_1&y_1 \end{array}\right)\left(\begin{array}{c} n_i \\ n_i \end{array}\right),\qquad i=1,2,\cdots,N$

优化问题: 找到单位法向量 $\boldsymbol{n}$ ,

使得 $\min_{\|\boldsymbol{n}\|=1}{(h_1^2+h_2^2+\cdots+h_N^2)}$ 成立.

带入距离 $h_i$ 的计算公式, 得

$\boldsymbol{n} = arg\min_{\|n\|=1}{\left(\left| \left(\begin{array}{c} x_1 \\ y_1 \end{array}\right)\cdot \boldsymbol{n} \right|^2+\cdots+\left| \left(\begin{array}{c} x_N \\ y_N \end{array}\right)\cdot \boldsymbol{n} \right|^2 \right)}$

由矩阵相关的计算

$\boldsymbol{n} = arg\min_{\|\boldsymbol{n}\|=1}{\sum_{i=1}^{N}{\left| \left(\begin{array}{c} x_i \\ y_i \end{array}\right)\cdot \boldsymbol{n} \right|^2}} = arg\min_{\|\boldsymbol{n}\|=1}{\sum_{i=1}^{N}{\boldsymbol{n}^T \left(\begin{array}{c} x_i \\ y_i \end{array}\right) \left(\begin{array}{cc} x_i & y_i \end{array}\right) \boldsymbol{n}}}$

求和号里的结合律, 并定义数据集矩阵

$\boldsymbol{n}= arg\min_{\|\boldsymbol{n}\|=1}{\boldsymbol{n}^T \left( \sum_{i=1}^{N}{\left(\begin{array}{c} x_i \\ y_i \end{array}\right) \left(\begin{array}{cc} x_i & y_i \end{array}\right)}\right) \boldsymbol{n}} =:arg\min_{\|\boldsymbol{n}\|=1}{\boldsymbol{n}^T\boldsymbol{P}\boldsymbol{n}}$

其中定义的矩阵为数据集 $\boldsymbol{D}$ 的转置乘以本身:

$\boldsymbol{P} = \boldsymbol{D}^T\boldsymbol{D}= \left[\begin{array}{cc} x_1 & y_1 \\ x_2 & y_2 \\ \vdots&\vdots\\ x_N & y_N \end{array}\right]^T \left[\begin{array}{cc} x_1 & y_1 \\ x_2 & y_2 \\ \vdots&\vdots\\ x_N & y_N \end{array}\right] =\left[\begin{array}{cc} \sum_{i=1}^{N}{x_i^2} & \sum_{i=1}^{N}{x_iy_i}\\ \sum_{i=1}^{N}{x_iy_i} & \sum_{i=1}^{N}{y_i^2}\\ \end{array}\right]$

注意到矩阵 $\boldsymbol{P}$ 是对称矩阵, 因此对矩阵 $\boldsymbol{P}$ 做svd分解(也可视作特征分解):

$\boldsymbol{P} = \boldsymbol{U}\boldsymbol{\Sigma}\boldsymbol{U} ^T$

其中 $\boldsymbol{U}$ 是正交矩阵, $\boldsymbol{\Sigma}=diag(\sigma_1, \sigma_2)$ 对角矩阵 $( \sigma_1\geq\sigma_2 )$ , 则有:

$\min_{\|\boldsymbol{n}\|=1}{\boldsymbol{n}^T\boldsymbol{P}\boldsymbol{n}}~~\Leftrightarrow~~ \min_{\|\boldsymbol{n}\|=1}{\boldsymbol{n}^T\boldsymbol{U}\boldsymbol{\Sigma}\boldsymbol{U}^T\boldsymbol{n}} ~~\Leftrightarrow~~ \min_{\|\boldsymbol{n}\|=1}{\left(\boldsymbol{U}^T\boldsymbol{n}\right)^T\boldsymbol{\Sigma} \left(\boldsymbol{U}^T\boldsymbol{n}\right)}$

数据集为固定值, 则其数据集矩阵的svd分解后的 $\boldsymbol{U}$ 也为固定值, 再因为 $\boldsymbol{U}$ 是正交矩阵, 因此假设新的单位长度向量

$\boldsymbol{s} = \left(\begin{array}{c} s_1 \\ s_2 \end{array}\right):= \boldsymbol{U}^T\boldsymbol{n}$

接下来, 原优化问题等价为:

$\min_{\|\boldsymbol{n}\|=1}{\boldsymbol{n}^T\boldsymbol{P}\boldsymbol{n}}~~\Leftrightarrow~~\min_{\|\boldsymbol{s}\|=1} \left(\begin{array}{cc} s_1 & s_2 \end{array}\right) \left[\begin{array}{cc} \sigma_1 & 0 \\ 0 & \sigma_2\\ \end{array}\right] \left(\begin{array}{c} s_1 \\ s_2 \end{array}\right)$

不难得知

$\left(\begin{array}{cc} s_1 & s_2 \end{array}\right) \left[\begin{array}{cc} \sigma_1 & 0 \\ 0 & \sigma_2\\ \end{array}\right] \left(\begin{array}{c} s_1 \\ s_2 \end{array}\right) = s_1^2\sigma_1+s_2^2\sigma_2=\sigma_2+(\sigma_1-\sigma_2)s_1^2$

注意到其中优化问题的变量 $s_1^2+s_2^2=1$ 且 $\sigma_1\geq\sigma_2$ 只与原始数据相关的固定值, 因此当且仅当 $s_1=0,~~s_2=1$ 时取得最小值, 即

$\boldsymbol{n} = \boldsymbol{U}\left(\begin{array}{c} 0 \\ 1 \end{array}\right) 时 \min_{\|n\|=1}{\left(\| \left(\begin{array}{c} x_1 \\ y_1 \end{array}\right)\cdot \boldsymbol{n} \|^2+\cdots+\| \left(\begin{array}{c} x_N \\ y_N \end{array}\right)\cdot \boldsymbol{n} \|^2 \right)}$ 成立取最小值

即直线的法向方向应该为数据集矩阵 $\boldsymbol{D}^T\boldsymbol{D}= \boldsymbol{U}\boldsymbol{\Sigma}\boldsymbol{U} ^T$ 的最小特征值对应的特征方向!

最后, 不难得知直线的斜率 $k=-n_1/n_2$ , 此外再将数据点的重心并不在原点, 即该直线通过数集的重心, 得到 $b= \frac{1}{N}\sum_{i=1}^{N}{y_i}-\left(\frac{1}{N}\sum_{i=1}^{N}{x_i}\right)k$ .

这其实就是在数据降维中常用的主成分分析方法!

数值测试对比

为验证上文的方法和结论等, 我测试了三组数据:

数据基本在一条直线附近, 带有轻微扰动
数据扰动较大, 略呈斜扁椭圆状
数据扰动很大, 呈斜椭圆状

以下图

不难发现, 在绝大多数的直线拟合问题下(数据1), 两种”距离”定义下拟合的结果(ls和pca)几乎重合, 均很好拟合了原数据.然而, 在斜椭圆数据(可以联系协方差相关知识)下, “垂直”距离求解(pca)比传统的”竖直”距离(最小二乘)求解要好很多, 如上图以及本文最初的封面插画所示.

因此对于标题的以为, 我想回答是: 嗯, 没有问题, 有问题的是数据……….

测试代码

测试平台为Anaconda, 环境为python2

"""
Created on Thu May  3 21:48:26 2018

@author: HSSheng
"""
import numpy as np
import matplotlib.pyplot as plt

# creat test data===========================================
# data0: random in [-1,1]x[-1,1]
data0 = np.random.rand(1000,2)
data0 = 2*data0-1
# data1: data0 cutted in ellipse
data1 = np.empty((1000,2))
count = 0
for i in range(np.size(data0,0)):
    if data0[i,0]**2+4.*data0[i,1]**2<1:
        data1[count,:] = data0[i,:]
        count = count+1
data1 = data1[:count,:]
# data: data1 rotation and displacement
theta = np.pi/5
rotMatrix = np.array([[np.cos(theta),np.sin(theta)],
                      [-np.sin(theta),np.cos(theta)]])
data = np.dot(data1,rotMatrix)
data[:,0] += 2
data[:,1] += 3
# Least square Method======================================
N = np.size(data,0) # total number of data
coeMatrix = np.vstack((data[:,0],np.ones(N))).transpose()
coeRhs = data[:,1]
A = np.dot(coeMatrix.transpose(),coeMatrix)
f = np.dot(coeMatrix.transpose(),coeRhs)
kb = np.linalg.solve(A, f)
k = kb[0]
b = kb[1]
# data for plot
x_ls = np.linspace(1,3)
y_ls = x_ls*k+b
# PCA method===============================================
# move center to zero point
dataHomo = data.copy()
dataHomo[:,0] -= np.sum(data[:,0])/N
dataHomo[:,1] -= np.sum(data[:,1])/N
# data matrix
dataMatrix = np.dot(dataHomo.transpose(),dataHomo)
u, s, vh = np.linalg.svd(dataMatrix, full_matrices=True)
n = u[:,-1]
k2 = -n[0]/n[1]
b2 = np.sum(data[:,1])/N-k2*np.sum(data[:,0])/N
# data for plot
x_pca = np.linspace(1,3)
y_pca = x_pca*k2+b2


plt.plot(data[:,0],data[:,1],'.')
plt.legend('data', shadow=True,)
plt.plot(x_ls,y_ls,linewidth=3)
plt.legend('least square', shadow=True,)
plt.plot(x_pca,y_pca,linewidth=3)
plt.legend(('data', 'least square','pca'), shadow=True,)

plt.axis('equal')
plt.show()

来源：知乎 www.zhihu.com

作者：Dr.Sheng

【知乎日报】千万用户的选择，做朋友圈里的新鲜事分享大牛。
点击下载