本文是关于最小二乘问题的一些总结。
问题引入
最小二乘法最常见于曲线拟合的问题。比如有一组样本点
(
x
1
,
y
1
)
(x_1, y_1)
(x1,y1),
(
x
2
,
y
2
)
(x_2, y_2)
(x2,y2),
(
x
3
,
y
3
)
(x_3, y_3)
(x3,y3), …,
(
x
m
,
y
m
)
(x_m, y_m)
(xm,ym)。每个样本点的
x
x
x都是n-1维列向量。然后根据用户定义的模型可以列出m个方程组成的方程组来求解模型中的参数。
以
x
∈
C
3
x \in { {\rm{C} }^3}
x∈C3为例,每个样本是3维的数据,共5组样本数据,列出一次多项式方程组,求解n个参数
a
,
b
,
c
,
d
a, b, c, d
a,b,c,d:
{ a x 11 + b x 12 + c x 13 + d = y 1 a x 21 + b x 22 + c x 23 + d = y 2 a x 31 + b x 32 + c x 33 + d = y 3 a x 41 + b x 42 + c x 43 + d = y 4 a x 51 + b x 52 + c x 53 + d = y 5 \left\{ \begin{array}{l} ax_{11} + bx_{12} + c{x_{13} } + d = {y_1}\\ ax_{21} + bx_{22} + c{x_{23} } + d = {y_2}\\ ax_{31} + bx_{32} + c{x_{33} } + d = {y_3}\\ ax_{41} + bx_{42} + c{x_{43} } + d = {y_4}\\ ax_{51} + bx_{52} + c{x_{53} } + d = {y_5} \end{array} \right. ⎩ ⎨ ⎧ax11+bx12+cx13+d=y1ax21+bx22+cx23+d=y2ax31+bx32+cx33+d=y3ax41+bx42+cx43+d=y4ax51+bx52+cx53+d=y5
也可以采用更加复杂的模型,例如:
{ a e x 11 + b x 12 + c ln x 13 + d = y 1 a e x 21 + b x 22 + c ln x 23 + d = y 2 a e x 31 + b x 32 + c ln x 33 + d = y 3 a e x 41 + b x 42 + c ln x 43 + d = y 4 a e x 51 + b x 52 + c ln x 53 + d = y 5 \left\{ \begin{array}{l} a{e^{ {x_{11} } } } + \frac{b}{ { {x_{12} } } } + c\ln {x_{13} } + d = {y_1}\\ a{e^{ {x_{21} } } } + \frac{b}{ { {x_{22} } } } + c\ln {x_{23} } + d = {y_2}\\ a{e^{ {x_{31} } } } + \frac{b}{ { {x_{32} } } } + c\ln {x_{33} } + d = {y_3}\\ a{e^{ {x_{41} } } } + \frac{b}{ { {x_{42} } } } + c\ln {x_{43} } + d = {y_4}\\ a{e^{ {x_{51} } } } + \frac{b}{ { {x_{52} } } } + c\ln {x_{53} } + d = {y_5} \end{array} \right. ⎩ ⎨ ⎧aex11+x12b+clnx13+d=y1aex21+x22b+clnx23+d=y2aex31+x32b+clnx33+d=y3aex41+x42b+clnx43+d=y4aex51+x52b+clnx53+d=y5
求解的时候所有的样本点都是会代入方程组中的,所以只要待求的参数最后能组成线性的方程组就行。
在线性代数中我们学过线性方程组解的结构。
[ x 11 x 12 x 13 1 x 21 x 22 x 23 1 x 31 x 32 x 33 1 x 41 x 42 x 43 1 x 51 x 52 x 53 1 ] [ a b c d ] = A [ a b c d ] = [ y 1 y 2 y 3 y 4 ] {\begin{bmatrix} {x_{11} }&{x_{12} }&{ {x_{13} } }&1\\ {x_{21} }&{x_{22} }&{ {x_{23} } }&1\\ {x_{31} }&{x_{32} }&{ {x_{33} } }&1\\ {x_{41} }&{x_{42} }&{ {x_{43} } }&1\\ {x_{51} }&{x_{52} }&{ {x_{53} } }&1 \end{bmatrix} } {\begin{bmatrix} a\\ b\\ c\\ d \end{bmatrix} } = A {\begin{bmatrix} a\\ b\\ c\\ d \end{bmatrix} } = {\begin{bmatrix} { {y_1} }\\ { {y_2} }\\ { {y_3} }\\ { {y_4} } \end{bmatrix} } ⎣ ⎡x11x21x31x41x51x12x22x32x42x52x13x23x33x43x5311111⎦ ⎤⎣ ⎡abcd⎦ ⎤=A⎣ ⎡abcd⎦ ⎤=⎣ ⎡y1y2y3y4⎦ ⎤
可以把前面的第一个方程组写成矩阵的形式,系数用矩阵A表示,A的秩
r
a
n
k
(
A
)
\rm rank(A)
rank(A)如果小于n(n是参数的个数),那么就有无穷解;
r
a
n
k
(
A
)
=
n
\rm rank(A)=n
rank(A)=n,则有唯一解;
r
a
n
k
(
A
)
>
n
\rm rank(A)>n
rank(A)>n则无解。这和我们熟知的,要解n个变量,就需要n个方程是一个道理。
最小二乘法主要用在“无解”的情况,在实际应用中,采集大量样本数据后,往往是得不到唯一解的。但我们想要得到一个曲线能最好地描述这一批数据,最小二乘法就是能得到令估计误差的平方和最小的一种估计方法。
而最佳体现在哪里?一般很少会有人提到最佳的最小二乘解是因为在“无解”的情况下,最小二乘解是唯一的,而只有在“无穷解”的情况下才有“最佳”的说法。而在“无穷解”的情况下相当于是在用极少的样本估计模型的参数,这在实际应用中一般也是没有意义的,所以“最佳”这个概念基本没人提。
问题的数学描述
目标函数可以这样表示:
arg min ω ∥ X ω − Y ∥ 2 2 \mathop {\arg \min }\limits_\omega {\left\| {X\omega - Y} \right\|_2^2} ωargmin∥Xω−Y∥22
ω
\omega
ω是一个n维的列向量,是我们需要估计的参数。
X
X
X是一个
m
×
n
m\times n
m×n的矩阵,每一行对应一个样本数据,共m行,也就是说有m个方程;
Y
Y
Y是每个样本的标签值组成的
m
×
1
m\times 1
m×1的列向量。
几何意义就是样本点与估计的曲线的偏差
d
y
dy
dy越小越好,目标函数取的是误差向量的二范数,忽略了偏差的正负。
解析推导
我之前在关于CSK与KCF算法推导(二)中写了岭回归的推导,岭回归就是在最小二乘的基础上加上了正则项,保证了
ω
\omega
ω一定有解。这里的推导和岭回归几乎一致。
类似二次多项式求最小值的问题,目标函数对
ω
\omega
ω求偏导,令其偏导为零向量即可求出
ω
\omega
ω在什么情况下能令目标函数最小化。
∂ ∥ X ω − Y ∥ 2 2 ∂ ω = 0 \frac{ {\partial \left\| {X\omega - Y} \right\|_2^2} }{ {\partial \omega } } = 0 ∂ω∂∥Xω−Y∥22=0
将目标函数展开,向量二范数的平方就是向量内积:
∥ X ω − Y ∥ 2 2 = ( X ω − Y ) T ( X ω − Y ) = ω T X T X ω − ω T X T Y − Y T X ω + Y T Y \left\| {X\omega - Y} \right\|_2^2 = {\left( {X\omega - Y} \right)^T}\left( {X\omega - Y} \right) = {\omega ^T}{X^T}X\omega - {\omega ^T}{X^T}Y - {Y^T}X\omega + {Y^T}Y ∥Xω−Y∥22=(Xω−Y)T(Xω−Y)=ωTXTXω−ωTXTY−YTXω+YTY
每一项分别对 ω \omega ω求偏导:
∂ ω T X T X ω ∂ ω = X T X ω + ( ω T X T X ) T = 2 X T X ω \frac{ {\partial {\omega ^T}{X^T}X\omega } }{ {\partial \omega } } = {X^T}X\omega + {\left( { {\omega ^T}{X^T}X} \right)^T} = 2{X^T}X\omega ∂ω∂ωTXTXω=XTXω+(ωTXTX)T=2XTXω
∂ ω T X T Y ∂ ω = X T Y \frac{ {\partial {\omega ^T}{X^T}Y} }{ {\partial \omega } } = {X^T}Y ∂ω∂ωTXTY=XTY
∂ Y T X ω ∂ ω = ( Y T X ) = X T Y \frac{ {\partial {Y^T}X\omega } }{ {\partial \omega } } = \left( { {Y^T}X} \right) = {X^T}Y ∂ω∂YTXω=(YTX)=XTY
∂ Y T Y ∂ ω = 0 \frac{ {\partial {Y^T}Y} }{ {\partial \omega } } = 0 ∂ω∂YTY=0
整合起来令偏导为零向量:
2 ( X T X ω − X T Y ) = 0 2\left( { {X^T}X\omega - {X^T}Y} \right) = 0 2(XTXω−XTY)=0
也就是:
X T X ω = X T Y {X^T}X\omega = {X^T}Y XTXω=XTY
如果没有额外的条件的话是无法保证 X T X X^TX XTX可逆的,如果 X T X X^TX XTX可逆,那么:
ω = ( X T X ) − 1 X T Y \omega = {\left( { {X^T}X} \right)^{ - 1} }{X^T}Y ω=(XTX)−1XTY
X T X X^TX XTX可逆也不难,只要 X X X是列满秩的就行。 X X X是 m × n m\times n m×n的矩阵,一般情况下样本个数m远大与待求的变量个数n时, X X X列满秩是很容易满足的。
从正交投影矩阵的角度推导
参考了华中科技大学出版社的《矩阵论》(第二版)教材,书中关于最佳的最小二乘解的内容得从M-P广义逆开始说起。求最佳的最小二乘解是M-P广义逆的一个应用,想要详细了解的小伙伴建议直接看书,我这里写的东西难免有些不严谨的地方,而且没法真正地把这个思路讲明白。书中的结论是更加一般化的,而且推导步骤也很严谨,前面的推导的结果是其中的一个特例。
M-P广义逆
书中介绍了M-P广义逆的一些性质,这里就不写了,主要提一下怎么计算M-P广义逆。任意矩阵
A
∈
C
m
×
n
A \in { {\rm{C} }^{m \times n} }
A∈Cm×n都存在M-P广义逆
A
+
A^+
A+。先对
A
A
A进行满秩分解,
A
=
B
C
,
B
∈
C
m
×
r
,
C
∈
C
r
×
n
,
r
a
n
k
(
B
)
=
r
a
n
k
(
C
)
=
r
A = BC,\;B \in { {\rm{C} }^{m \times r} },\;C \in { {\rm{C} }^{r \times n} },\;{\mathop{\rm rank}\nolimits} (B) = {\mathop{\rm rank}\nolimits} (C) = r
A=BC,B∈Cm×r,C∈Cr×n,rank(B)=rank(C)=r
M-P广义逆就可以表示为,
A
+
=
C
H
(
C
C
H
)
−
1
(
B
H
B
)
−
1
B
H
{A^ + } = {C^H}{\left( {C{C^H} } \right)^{ - 1} }{\left( { {B^H}B} \right)^{ - 1} }{B^H}
A+=CH(CCH)−1(BHB)−1BH
特别的,当
A
A
A为列满秩的矩阵时,
C
=
I
r
×
r
C=I_{r\times r}
C=Ir×r,
A
+
=
(
A
H
A
)
−
1
A
H
{A^ + } = {\left( { {A^H}A} \right)^{ - 1} }{A^H}
A+=(AHA)−1AH
正交投影矩阵
正交投影矩阵
P
P
P满足:
P
2
=
P
,
P
H
=
P
P^2=P,\;P^H=P
P2=P,PH=P。正交投影矩阵有什么用呢?假如
P
∈
C
n
×
n
P \in { {\rm{C} }^{n \times n} }
P∈Cn×n,是一个正交投影矩阵,
v
=
P
u
,
u
∈
C
n
×
1
,
v
∈
C
n
×
1
v=Pu, u \in \rm C^{n\times 1}, v \in \rm C^{n\times 1}
v=Pu,u∈Cn×1,v∈Cn×1
u
u
u经过
P
P
P变换以后可以得到一个
R
(
P
)
R(P)
R(P)空间中的一个向量
v
v
v,
∥
v
−
u
∥
2
≤
∥
x
−
u
∥
2
,
∀
x
∈
R
(
P
)
{\left\| {v - u} \right\|_2} \le \left\| {x - u} \right\|_2,\;\forall x \in R(P)
∥v−u∥2≤∥x−u∥2,∀x∈R(P)
向量
v
v
v与
u
u
u之间的欧式距离小于
R
(
P
)
R(P)
R(P)空间中的任何向量
x
x
x
考察矩阵
A
A
+
AA^+
AA+和矩阵
A
+
A
A^+A
A+A可以发现它们都是正交投影矩阵。书中有证明
R
(
A
A
+
)
=
R
(
A
)
R(AA^+)=R(A)
R(AA+)=R(A),这为我们带来了很大的方便。如果我们想把向量
x
x
x正交投影到
R
(
A
)
R(A)
R(A)空间,只需要计算
A
A
+
x
AA^+x
AA+x就行了。
回到问题
arg min ω ∥ X ω − Y ∥ 2 2 \mathop {\arg \min }\limits_\omega {\left\| {X\omega - Y} \right\|_2^2} ωargmin∥Xω−Y∥22
再回到原来的问题,我们想找的就是 X ω X\omega Xω如果能是 Y Y Y在 R ( X ) R(X) R(X)空间中的正交投影就好了,前面的结论就告诉了我们 X X + Y XX^+Y XX+Y就是 Y Y Y在 R ( X ) R(X) R(X)空间中的正交投影。那么自然,
ω = X + Y \omega = {X^ + }Y ω=X+Y
与前面的解析推导不一样的地方是,这里没有对 X X X的形式做限制。如果同样加上 X X X是列满秩的限制,那么可以得到与前面一样的结论。
ω = ( X T X ) − 1 X T Y \omega = {\left( { {X^T}X} \right)^{ - 1} }{X^T}Y ω=(XTX)−1XTY
ω = X + Y \omega = {X^ + }Y ω=X+Y是一个更加一般的结论,它被称为“最佳的最小二乘解”。这个解不光满足最小二乘,而且这个解本身的二范数也是最小的,所以是所有最小二乘解中最佳的解。而“最佳”也只会在有无穷解的情况下体现,下面举个最简单的例子。
所谓“最佳”的举例
用样本点(1,2)拟合一条直线,
y
=
a
x
+
b
y=ax+b
y=ax+b。虽然用一个点去拟合一条直线听着有点离谱,但这就是体现最佳的最小二乘解所谓的“最佳”的地方。
[
1
1
]
[
a
b
]
=
2
{\begin{bmatrix} 1&1 \end{bmatrix} } {\begin{bmatrix} a\\ b \end{bmatrix} } = 2
[11][ab]=2
X
=
[
1
1
]
,
X
+
=
X
T
(
X
X
T
)
−
1
=
[
0.5
0.5
]
X = {\begin{bmatrix} 1&1 \end{bmatrix} },\;{X^ + } = {X^T}{\left( {X{X^T} } \right)^{ - 1} } = {\begin{bmatrix} {0.5}\\ {0.5} \end{bmatrix} }
X=[11],X+=XT(XXT)−1=[0.50.5]
最佳的最小二乘解:
[
a
b
]
=
X
+
Y
=
[
1
1
]
{\begin{bmatrix} a\\ b \end{bmatrix} } = {X^ + }Y = {\begin{bmatrix} 1\\ 1 \end{bmatrix} }
[ab]=X+Y=[11]
过点(1,2)的直线有无数条,这无数条直线的参数 a , b a,b a,b都是这一问题的最小二乘解,它们都能令 X ω − Y = 0 X\omega-Y=0 Xω−Y=0,而这些解中 a = 1 , b = 1 a=1,b=1 a=1,b=1是最佳的,因为此时 a 2 + b 2 a^2+b^2 a2+b2最小。