3.1 长度
首先,以一个简单的线性回归问题(一条直线拟合数据的问题),来直观的体会“长度”。
a)对数据(z,d)的最小二乘直线拟合;b)
e
i
=
d
i
o
b
s
−
d
i
p
r
e
e_i = d_{i}^{obs} - d_{i}^{pre}
ei=diobs−dipre
那么,最佳的拟合直线所具有模型参数(截距和斜率)将使总误差E最小,即
E
=
∑
i
=
1
n
e
i
2
=
e
T
e
=
m
i
n
E = \sum\limits_{i=1}^{n}e_{i}^{2}=\bf {e^Te} = min
E=i=1∑nei2=eTe=min
式中,总误差
E
E
E恰恰使向量
e
\bf e
e的欧几里得长度。
从长度方法的观点来看,最小二乘法是通过寻找预测误差的最小长度所对应的模型参数(截距和斜率)来估计反问题的解。通常来讲在,求解反问题的过程中使用长度方法是最简单的方法,也符合人的直观感受。
3.2 范数==》长度的推广
所谓范数是指对某些长度或大小的度量,其公式如下:
∥
e
∥
p
=
(
∑
i
=
1
n
e
p
)
1
p
{\|e\|}_{p} = (\sum\limits_{i=1}^{n}e^p)^{\frac{1}{p}}
∥e∥p=(i=1∑nep)p1
以下三种范数最为常用:
L
1
L_{1}
L1 norm:
∥
e
∥
1
=
[
∑
i
∣
e
i
∣
1
]
\quad\|\mathbf{e}\|_{1}=\left[\sum_{i}\left|e_{i}\right|^{1}\right]
∥e∥1=[∑i∣ei∣1]
L
2
L_{2}
L2 norm:
∥
e
∥
2
=
[
∑
i
∣
e
i
∣
2
]
1
/
2
\quad\|\mathbf{e}\|_{2}=\left[\sum_{i}\left|e_{i}\right|^{2}\right]^{1 / 2}
∥e∥2=[∑i∣ei∣2]1/2
L
∞
L_{\infty}
L∞ norm:
∥
e
∥
∞
=
max
i
∣
e
i
∣
\quad\|\mathbf{e}\|_{\infty}=\max _{i}\left|e_{i}\right|
∥e∥∞=maxi∣ei∣
幂次越高的范数给
e
\mathbf{e}
e中最大元素更大的权重,
L
∞
L_{\infty}
L∞表明仅对
e
\mathbf{e}
e的最大元素施加权重
从上图上看,
∣
e
∣
|e|
∣e∣的显著性最明显,而
∣
e
10
∣
|e^{10}|
∣e10∣最差
那么,最小二乘法究竟该选择哪种长度(即采取哪种范数),这个问题的答案涉及对远离平均趋势的离群数据进行加权的方式,如下图所示:
图中,
L
1
L_1
L1范数给离群点的权重最小。
3.3 线性反问题的最小二乘解
最小二乘能够以一种非常直接的方式推广到一般线性反问题。
E
=
e
T
e
=
(
d
−
G
m
)
T
(
d
−
G
m
)
=
∑
i
=
1
N
[
d
i
−
∑
j
=
1
M
G
i
j
m
j
]
[
d
i
−
∑
k
=
1
M
G
i
k
m
k
]
E=\mathbf{e}^{\mathrm{T}} \mathbf{e}=(\mathbf{d}-\mathbf{G} \mathbf{m})^{\mathrm{T}}(\mathbf{d}-\mathbf{G} \mathbf{m})=\sum_{i=1}^{N}\left[d_{i}-\sum_{j=1}^{M} G_{i j} m_{j}\right]\left[d_{i}-\sum_{k=1}^{M} G_{i k} m_{k}\right]
E=eTe=(d−Gm)T(d−Gm)=i=1∑N[di−j=1∑MGijmj][di−k=1∑MGikmk]
为避免混乱,对上式进行重新整理:
E
=
∑
j
=
1
M
∑
k
=
1
M
m
j
m
k
∑
i
=
1
N
G
i
j
G
i
k
−
2
∑
j
=
1
M
m
j
∑
i
=
1
N
G
i
j
d
i
+
∑
i
=
1
N
d
i
d
i
E=\sum_{j=1}^{M} \sum_{k=1}^{M} m_{j} m_{k} \sum_{i=1}^{N} G_{i j} G_{i k}-2 \sum_{j=1}^{M} m_{j} \sum_{i=1}^{N} G_{i j} d_{i}+\sum_{i=1}^{N} d_{i} d_{i}
E=j=1∑Mk=1∑Mmjmki=1∑NGijGik−2j=1∑Mmji=1∑NGijdi+i=1∑Ndidi
令
∂
E
∂
m
p
=
0
\frac{\partial E}{\partial m_{p}}=0
∂mp∂E=0
第一项为:
∂
∂
m
q
[
∑
j
=
1
M
∑
k
=
1
M
m
j
m
k
∑
i
=
1
N
G
i
j
G
i
k
]
=
∑
j
=
1
M
∑
k
=
1
M
[
δ
j
q
m
k
+
m
j
δ
k
q
]
∑
i
=
1
N
G
i
j
G
i
k
=
2
∑
k
=
1
M
m
k
∑
i
=
1
N
G
i
q
G
i
k
\begin{aligned} \frac{\partial}{\partial m_{q}}\left[\sum_{j=1}^{M} \sum_{k=1}^{M} m_{j} m_{k} \sum_{i=1}^{N} G_{i j} G_{i k}\right] &=\sum_{j=1}^{M} \sum_{k=1}^{M}\left[\delta_{j q} m_{k}+m_{j} \delta_{k q}\right] \sum_{i=1}^{N} G_{i j} G_{i k} \\ &=2 \sum_{k=1}^{M} m_{k} \sum_{i=1}^{N} G_{i q} G_{i k} \end{aligned}
∂mq∂[j=1∑Mk=1∑Mmjmki=1∑NGijGik]=j=1∑Mk=1∑M[δjqmk+mjδkq]i=1∑NGijGik=2k=1∑Mmki=1∑NGiqGik
注意:模型参数是相互独立的变量, ∂ m i / ∂ m j \partial m_{i} / \partial m_{j} ∂mi/∂mj形式的导数只有在 i = j i=j i=j时候等于 1 1 1,而 i ≠ j i \neq j i=j是为0,因此 ∂ m i / ∂ m j \partial m_{i} / \partial m_{j} ∂mi/∂mj克罗内克函数(Kronecker delta) δ i j \delta_{i j} δij,包含它的公式可以明显的简化
第二项为:
−
2
∂
∂
m
q
[
∑
j
=
1
M
m
j
∑
i
=
1
N
G
i
j
d
i
]
=
−
2
∑
j
=
1
M
δ
j
q
∑
i
=
1
N
G
i
j
d
i
=
−
2
∑
i
=
1
N
G
i
q
d
i
-2 \frac{\partial}{\partial m_{q}}\left[\sum_{j=1}^{M} m_{j} \sum_{i=1}^{N} G_{i j} d_{i}\right]=-2 \sum_{j=1}^{M} \delta_{j q} \sum_{i=1}^{N} G_{i j} d_{i}=-2 \sum_{i=1}^{N} G_{i q} d_{i}
−2∂mq∂[j=1∑Mmji=1∑NGijdi]=−2j=1∑Mδjqi=1∑NGijdi=−2i=1∑NGiqdi
第三项为:
∂
∂
m
q
[
∑
i
=
1
N
d
i
d
i
]
=
0
\frac{\partial}{\partial m_{q}}\left[\sum_{i=1}^{N} d_{i} d_{i}\right]=0
∂mq∂[i=1∑Ndidi]=0
结合三项,可得:
∂
E
∂
m
q
=
0
=
2
∑
k
=
1
M
m
k
∑
i
=
1
N
G
i
q
G
i
k
−
2
∑
i
=
1
N
G
i
q
d
i
\frac{\partial E}{\partial m_{q}}=0=2 \sum_{k=1}^{M} m_{k} \sum_{i=1}^{N} G_{i q} G_{i k}-2 \sum_{i=1}^{N} G_{i q} d_{i}
∂mq∂E=0=2k=1∑Mmki=1∑NGiqGik−2i=1∑NGiqdi
重新写成矩阵形式:
G
T
G
m
−
G
T
d
=
0
\mathbf{G}^{\mathrm{T}} \mathbf{G m}-\mathbf{G}^{\mathrm{T}} \mathbf{d}=0
GTGm−GTd=0
假设
[
G
T
G
]
−
1
\left[\mathbf{G}^{\mathrm{T}} \mathbf{G}\right]^{-1}
[GTG]−1存在(有不存在的时候),那么有如下解:
m
e
s
t
=
[
G
T
G
]
−
1
G
T
d
\mathbf{m}^{\mathrm{est}}=\left[\mathbf{G}^{\mathrm{T}} \mathbf{G}\right]^{-1} \mathbf{G}^{\mathrm{T}} \mathbf{d}
mest=[GTG]−1GTd
对于维度较大的情况, G T G \mathbf{G}^{\mathrm{T}} \mathbf{G} GTG的计算成本可能很高,而且 G T G \mathbf{G}^{\mathrm{T}} \mathbf{G} GTG极少像 G \mathbf G G那样稀疏。在这种情况下,优先考虑迭代的矩阵求解方案,如双共轭梯度算法(biconjugate gradient method)
3.4 一些最小二乘法的例子
1) 直线拟合问题:
模型是
d
i
=
m
1
+
m
2
z
i
d_i =m_1+m_2z_i
di=m1+m2zi,方程
G
m
=
d
\mathbf{Gm=d}
Gm=d形式如下:
[
1
z
1
1
z
2
⋮
⋮
1
z
N
]
[
m
1
m
2
]
=
[
d
1
d
2
⋮
d
N
]
\left[\begin{array}{cc} 1 & z_{1} \\ 1 & z_{2} \\ \vdots & \vdots \\ 1 & z_{N} \end{array}\right]\left[\begin{array}{c} m_{1} \\ m_{2} \end{array}\right]=\left[\begin{array}{c} d_{1} \\ d_{2} \\ \vdots \\ d_{N} \end{array}\right]
⎣⎢⎢⎢⎡11⋮1z1z2⋮zN⎦⎥⎥⎥⎤[m1m2]=⎣⎢⎢⎢⎡d1d2⋮dN⎦⎥⎥⎥⎤
==》
G
T
G
=
[
1
1
⋯
1
z
1
z
2
⋯
z
N
]
[
1
z
1
1
z
2
⋮
⋮
1
z
N
]
=
[
N
∑
i
=
1
N
z
i
∑
i
=
1
N
z
i
∑
i
=
1
N
z
i
2
]
\mathbf{G}^{\mathrm{T}} \mathbf{G}=\left[\begin{array}{cccc} 1 & 1 & \cdots & 1 \\ z_{1} & z_{2} & \cdots & z_{N} \end{array}\right]\left[\begin{array}{cc} 1 & z_{1} \\ 1 & z_{2} \\ \vdots & \vdots \\ 1 & z_{N} \end{array}\right]=\left[\begin{array}{cc} N & \sum_{i=1}^{N} z_{i} \\ \sum_{i=1}^{N} z_{i} & \sum_{i=1}^{N} z_{i}^{2} \end{array}\right]
GTG=[1z11z2⋯⋯1zN]⎣⎢⎢⎢⎡11⋮1z1z2⋮zN⎦⎥⎥⎥⎤=[N∑i=1Nzi∑i=1Nzi∑i=1Nzi2]
G
T
d
=
[
1
1
⋯
1
z
1
z
2
⋯
z
N
]
[
d
1
d
2
⋮
d
N
]
=
[
∑
i
=
1
N
d
i
∑
i
=
1
N
d
i
z
i
]
\mathbf{G}^{\mathrm{T}} \mathbf{d}=\left[\begin{array}{cccc} 1 & 1 & \cdots & 1 \\ z_{1} & z_{2} & \cdots & z_{N} \end{array}\right]\left[\begin{array}{c} d_{1} \\ d_{2} \\ \vdots \\ d_{N} \end{array}\right]=\left[\begin{array}{c} \sum_{i=1}^{N} d_{i} \\ \sum_{i=1}^{N} d_{i} z_{i} \end{array}\right]
GTd=[1z11z2⋯⋯1zN]⎣⎢⎢⎢⎡d1d2⋮dN⎦⎥⎥⎥⎤=[∑i=1Ndi∑i=1Ndizi]
G
T
d
=
[
1
1
⋯
1
z
1
z
2
⋯
z
N
]
[
d
1
d
2
⋮
d
N
]
=
[
∑
i
=
1
N
d
i
∑
i
=
1
N
d
i
z
i
]
\mathbf{G}^{\mathrm{T}} \mathbf{d}=\left[\begin{array}{cccc} 1 & 1 & \cdots & 1 \\ z_{1} & z_{2} & \cdots & z_{N} \end{array}\right]\left[\begin{array}{c} d_{1} \\ d_{2} \\ \vdots \\ d_{N} \end{array}\right]=\left[\begin{array}{c} \sum_{i=1}^{N} d_{i} \\ \sum_{i=1}^{N} d_{i} z_{i} \end{array}\right]
GTd=[1z11z2⋯⋯1zN]⎣⎢⎢⎢⎡d1d2⋮dN⎦⎥⎥⎥⎤=[∑i=1Ndi∑i=1Ndizi]
m e s t = [ G T G ] − 1 G T d = [ N ∑ i = 1 N z i ∑ i = 1 N z i ∑ i = 1 N z i 2 ] − 1 [ ∑ i = 1 N d i ∑ i = 1 N d i z i ] \mathbf{m}^{\mathrm{est}}=\left[\mathbf{G}^{\mathrm{T}} \mathbf{G}\right]^{-1} \mathbf{G}^{\mathrm{T}} \mathbf{d}=\left[\begin{array}{cc} N & \sum_{i=1}^{N} z_{i} \\ \sum_{i=1}^{N} z_{i} & \sum_{i=1}^{N} z_{i}^{2} \end{array}\right]^{-1}\left[\begin{array}{c} \sum_{i=1}^{N} d_{i} \\ \sum_{i=1}^{N} d_{i} z_{i} \end{array}\right] mest=[GTG]−1GTd=[N∑i=1Nzi∑i=1Nzi∑i=1Nzi2]−1[∑i=1Ndi∑i=1Ndizi]
2 抛物线拟合
模型是
d
i
=
m
1
+
m
2
z
i
+
m
3
z
i
2
d_{i}=m_{1}+m_{2} z_{i}+m_{3} z_{i}^{2}
di=m1+m2zi+m3zi2,方程
G
m
=
d
\mathbf{Gm=d}
Gm=d形式如下:
[
1
z
1
z
1
2
1
z
2
z
2
2
⋮
⋮
⋮
1
z
N
z
N
2
]
[
m
1
m
2
m
3
]
=
[
d
1
d
2
⋮
d
N
]
\left[\begin{array}{ccc} 1 & z_{1} & z_{1}^{2} \\ 1 & z_{2} & z_{2}^{2} \\ \vdots & \vdots & \vdots \\ 1 & z_{N} & z_{N}^{2} \end{array}\right]\left[\begin{array}{c} m_{1} \\ m_{2} \\ m_{3} \end{array}\right]=\left[\begin{array}{c} d_{1} \\ d_{2} \\ \vdots \\ d_{N} \end{array}\right]
⎣⎢⎢⎢⎡11⋮1z1z2⋮zNz12z22⋮zN2⎦⎥⎥⎥⎤⎣⎡m1m2m3⎦⎤=⎣⎢⎢⎢⎡d1d2⋮dN⎦⎥⎥⎥⎤
==》
G
T
G
=
[
1
1
⋯
1
z
1
z
2
⋯
z
N
z
1
2
z
N
2
⋯
z
N
2
]
[
1
z
1
z
1
2
1
z
2
z
2
2
⋮
⋮
⋮
1
z
N
z
N
2
]
\mathbf{G}^{\mathrm{T}} \mathbf{G}=\left[\begin{array}{cccc} 1 & 1 & \cdots & 1 \\ z_{1} & z_{2} & \cdots & z_{N} \\ z_{1}^{2} & z_{N}^{2} & \cdots & z_{N}^{2} \end{array}\right]\left[\begin{array}{ccc} 1 & z_{1} & z_{1}^{2} \\ 1 & z_{2} & z_{2}^{2} \\ \vdots & \vdots & \vdots \\ 1 & z_{N} & z_{N}^{2} \end{array}\right]
GTG=⎣⎡1z1z121z2zN2⋯⋯⋯1zNzN2⎦⎤⎣⎢⎢⎢⎡11⋮1z1z2⋮zNz12z22⋮zN2⎦⎥⎥⎥⎤
G
T
d
=
[
1
1
⋯
1
z
1
z
2
⋯
z
N
z
1
2
z
N
2
⋯
z
N
2
]
[
d
1
d
2
⋮
d
N
]
=
[
∑
i
=
1
N
d
i
∑
i
=
1
N
z
i
d
i
∑
i
=
1
N
z
i
2
d
i
]
\mathbf{G}^{\mathrm{T}} \mathbf{d}=\left[\begin{array}{cccc} 1 & 1 & \cdots & 1 \\ z_{1} & z_{2} & \cdots & z_{N} \\ z_{1}^{2} & z_{N}^{2} & \cdots & z_{N}^{2} \end{array}\right]\left[\begin{array}{c} d_{1} \\ d_{2} \\ \vdots \\ d_{N} \end{array}\right]=\left[\begin{array}{c} \sum_{i=1}^{N} d_{i} \\ \sum_{i=1}^{N} z_{i} d_{i} \\ \sum_{i=1}^{N} z_{i}^{2} d_{i} \end{array}\right]
GTd=⎣⎡1z1z121z2zN2⋯⋯⋯1zNzN2⎦⎤⎣⎢⎢⎢⎡d1d2⋮dN⎦⎥⎥⎥⎤=⎣⎢⎡∑i=1Ndi∑i=1Nzidi∑i=1Nzi2di⎦⎥⎤
m
e
s
t
=
[
G
T
G
]
−
1
G
T
d
=
[
N
∑
i
=
1
N
z
i
∑
i
=
1
N
z
i
2
∑
i
=
1
N
z
i
∑
i
=
1
N
z
i
2
∑
i
=
1
N
z
i
3
∑
i
=
1
N
z
i
2
∑
i
=
1
N
z
i
3
∑
i
=
1
N
z
i
4
]
−
1
[
∑
i
=
1
N
d
i
∑
i
=
1
N
z
i
d
i
∑
i
=
1
N
z
i
2
d
i
]
\mathbf{m}^{\mathrm{est}}=\left[\mathbf{G}^{\mathrm{T}} \mathbf{G}\right]^{-1} \mathbf{G}^{\mathrm{T}} \mathbf{d}=\left[\begin{array}{ccc} N & \sum_{i=1}^{N} z_{i} & \sum_{i=1}^{N} z_{i}^{2} \\ \sum_{i=1}^{N} z_{i} & \sum_{i=1}^{N} z_{i}^{2} & \sum_{i=1}^{N} z_{i}^{3} \\ \sum_{i=1}^{N} z_{i}^{2} & \sum_{i=1}^{N} z_{i}^{3} & \sum_{i=1}^{N} z_{i}^{4} \end{array}\right]^{-1}\left[\begin{array}{c} \sum_{i=1}^{N} d_{i} \\ \sum_{i=1}^{N} z_{i} d_{i} \\ \sum_{i=1}^{N} z_{i}^{2} d_{i} \end{array}\right]
mest=[GTG]−1GTd=⎣⎢⎡N∑i=1Nzi∑i=1Nzi2∑i=1Nzi∑i=1Nzi2∑i=1Nzi3∑i=1Nzi2∑i=1Nzi3∑i=1Nzi4⎦⎥⎤−1⎣⎢⎡∑i=1Ndi∑i=1Nzidi∑i=1Nzi2di⎦⎥⎤
上图为开普勒第三定律的验证,它阐述了行星轨道半径的立方正比于轨道周期的平方。a)对于太阳系,数据(红圈)用二次方公式
d
i
=
m
1
+
m
2
z
i
+
m
3
z
i
2
d_i=m_1+m_2z_i+m_3z_{i}^{2}
di=m1+m2zi+m3zi2实现了最小二乘拟合,其中
d
i
d_i
di为半径的立方,
z
i
z_i
zi为周期。b)拟合误差,将数据和误差独立显示的目的是将误差以恰当的尺度画出。数据源于维基百科(Wikipedia)。
3、平面拟合
模型:
d
i
=
m
1
+
m
2
x
i
+
m
3
y
i
d_{i}=m_{1}+m_{2} x_{i}+m_{3} y_{i}
di=m1+m2xi+m3yi,方程
G
m
=
d
\mathbf{Gm=d}
Gm=d形式如下:
[
1
x
1
y
1
1
x
2
y
2
⋮
⋮
⋮
1
x
N
y
N
]
[
m
1
m
2
m
3
]
=
[
d
1
d
2
⋮
d
N
]
\left[\begin{array}{ccc} 1 & x_{1} & y_{1} \\ 1 & x_{2} & y_{2} \\ \vdots & \vdots & \vdots \\ 1 & x_{N} & y_{N} \end{array}\right]\left[\begin{array}{c} m_{1} \\ m_{2} \\ m_{3} \end{array}\right]=\left[\begin{array}{c} d_{1} \\ d_{2} \\ \vdots \\ d_{N} \end{array}\right]
⎣⎢⎢⎢⎡11⋮1x1x2⋮xNy1y2⋮yN⎦⎥⎥⎥⎤⎣⎡m1m2m3⎦⎤=⎣⎢⎢⎢⎡d1d2⋮dN⎦⎥⎥⎥⎤
==》
G
T
G
=
[
1
1
⋯
1
x
1
x
2
⋯
x
N
y
1
y
2
⋯
y
N
]
[
1
x
1
y
1
1
x
2
y
2
⋮
⋮
⋮
1
x
N
y
N
]
\mathbf{G}^{\mathrm{T}} \mathbf{G}=\left[\begin{array}{cccc} 1 & 1 & \cdots & 1 \\ x_{1} & x_{2} & \cdots & x_{N} \\ y_{1} & y_{2} & \cdots & y_{N} \end{array}\right]\left[\begin{array}{ccc} 1 & x_{1} & y_{1} \\ 1 & x_{2} & y_{2} \\ \vdots & \vdots & \vdots \\ 1 & x_{N} & y_{N} \end{array}\right]
GTG=⎣⎡1x1y11x2y2⋯⋯⋯1xNyN⎦⎤⎣⎢⎢⎢⎡11⋮1x1x2⋮xNy1y2⋮yN⎦⎥⎥⎥⎤
=
[
N
∑
i
=
1
N
x
i
∑
i
=
1
N
y
i
∑
i
=
1
N
x
i
∑
i
=
1
N
x
i
2
∑
i
=
1
N
x
i
y
i
∑
i
=
1
N
y
i
∑
i
=
1
N
x
i
y
i
∑
i
=
1
N
y
i
2
]
=\left[\begin{array}{ccc} N & \sum_{i=1}^{N} x_{i} & \sum_{i=1}^{N} y_{i} \\ \sum_{i=1}^{N} x_{i} & \sum_{i=1}^{N} x_{i}^{2} & \sum_{i=1}^{N} x_{i} y_{i} \\ \sum_{i=1}^{N} y_{i} & \sum_{i=1}^{N} x_{i} y_{i} & \sum_{i=1}^{N} y_{i}^{2} \end{array}\right]
=⎣⎢⎡N∑i=1Nxi∑i=1Nyi∑i=1Nxi∑i=1Nxi2∑i=1Nxiyi∑i=1Nyi∑i=1Nxiyi∑i=1Nyi2⎦⎥⎤
G T d = [ 1 1 ⋯ 1 x 1 x 2 ⋯ x N y 1 y 2 ⋯ y N ] [ d 1 d 2 ⋮ d N ] = [ ∑ i = 1 N d i ∑ i = 1 N x i d i ∑ i = 1 N y i d i ] \mathbf{G}^{\mathrm{T}} \mathbf{d}=\left[\begin{array}{cccc} 1 & 1 & \cdots & 1 \\ x_{1} & x_{2} & \cdots & x_{N} \\ y_{1} & y_{2} & \cdots & y_{N} \end{array}\right]\left[\begin{array}{c} d_{1} \\ d_{2} \\ \vdots \\ d_{N} \end{array}\right]=\left[\begin{array}{c} \sum_{i=1}^{N} d_{i} \\ \sum_{i=1}^{N} x_{i} d_{i} \\ \sum_{i=1}^{N} y_{i} d_{i} \end{array}\right] GTd=⎣⎡1x1y11x2y2⋯⋯⋯1xNyN⎦⎤⎣⎢⎢⎢⎡d1d2⋮dN⎦⎥⎥⎥⎤=⎣⎢⎡∑i=1Ndi∑i=1Nxidi∑i=1Nyidi⎦⎥⎤
m e s t = [ G T G ] − 1 G T d = [ N ∑ i = 1 N x i ∑ i = 1 N y i ∑ i = 1 N x i ∑ i = 1 N x i 2 ∑ i = 1 N x i y i ∑ i = 1 N y i ∑ i = 1 N x i y i ∑ i = 1 N x i y i 2 ] − 1 [ ∑ i = 1 N d i ∑ i = 1 N x i d i ∑ i = 1 N y i d i ] \mathbf{m}^{\mathrm{est}}=\left[\mathbf{G}^{\mathrm{T}} \mathbf{G}\right]^{-1} \mathbf{G}^{\mathrm{T}} \mathbf{d}=\left[\begin{array}{ccc} N & \sum_{i=1}^{N} x_{i} & \sum_{i=1}^{N} y_{i} \\ \sum_{i=1}^{N} x_{i} & \sum_{i=1}^{N} x_{i}^{2} & \sum_{i=1}^{N} x_{i} y_{i} \\ \sum_{i=1}^{N} y_{i} & \sum_{i=1}^{N} x_{i} y_{i} & \sum_{i=1}^{N} x_{i} y_{i}^{2} \end{array}\right]^{-1}\left[\begin{array}{c} \sum_{i=1}^{N} d_{i} \\ \sum_{i=1}^{N} x_{i} d_{i} \\ \sum_{i=1}^{N} y_{i} d_{i} \end{array}\right] mest=[GTG]−1GTd=⎣⎢⎡N∑i=1Nxi∑i=1Nyi∑i=1Nxi∑i=1Nxi2∑i=1Nxiyi∑i=1Nyi∑i=1Nxiyi∑i=1Nxiyi2⎦⎥⎤−1⎣⎢⎡∑i=1Ndi∑i=1Nxidi∑i=1Nyidi⎦⎥⎤
上图是使用平面拟合来检测地质断层上发生地震的例子,(圆圈)西北太平洋千岛群岛俯冲带发生的地震。
x
x
x坐标轴指向北,
y
y
y坐标轴指向东。地震散布在一个由最小二乘确定的倾斜平面附近。数据来自美国地质调查局(United States Geological Survey, USGS)。
最小二乘法的是为了求解没有精确解的反问题。该方法是为了寻找最佳近似解,从数学公式来看: E = e T e \mathbf {E = e^Te} E=eTe,寻找的是 ∂ E ∂ m = 0 \mathbf{\frac{\partial E}{\partial m}=0} ∂m∂E=0的解( m e s t m^{est} mest)即最小化 L 2 L_2 L2误差长度(此时,或许 e ≠ 0 e\neq0 e=0),而不是在求 E = 0 \mathbf{E = 0} E=0的解(即误差 e = 0 e=0 e=0),所以说最小二乘法的解是近似解(最佳的),而不是精确解
3.5 应用最小二乘的一些缺陷(存在性、唯一性问题)
公式 m e s t = [ G T G ] − 1 G T d \begin{array}{c} \mathbf{m}^{\mathrm{est}}=\left[\mathbf{G}^{\mathrm{T}} \mathbf{G}\right]^{-1} \mathbf{G}^{\mathrm{T}} \mathbf{d} \end{array} mest=[GTG]−1GTd中,隐含地假设仅存在唯一的最佳近似解。如果反问题存在多个解,并且这些解具有相同的最小预测误差,那么最小二乘法将会失效。
反例
如上图,考虑仅具有一个数据点的直线拟合问题,显然该问题的解是非唯一的,从图上来看,可以画出可以穿过该点的无穷条直线,且各直线都具有零预测误差(e=0)。
用公式
m
e
s
t
=
[
G
T
G
]
−
1
G
T
d
\mathbf{m}^{\mathrm{est}}=\left[\mathbf{G}^{\mathrm{T}} \mathbf{G}\right]^{-1} \mathbf{G}^{\mathrm{T}} \mathbf{d}
mest=[GTG]−1GTd去求解的话,则有
[
G
T
G
]
−
1
=
[
N
∑
i
=
1
N
z
i
∑
i
=
1
N
z
i
∑
i
=
1
N
z
i
2
]
−
1
→
[
1
z
1
z
1
z
1
2
]
−
1
\left[\mathbf{G}^{\mathrm{T}} \mathbf{G}\right]^{-1}=\left[\begin{array}{cc} N & \sum_{i=1}^{N} z_{i} \\ \sum_{i=1}^{N} z_{i} & \sum_{i=1}^{N} z_{i}^{2} \end{array}\right]^{-1} \rightarrow\left[\begin{array}{cc} 1 & z_{1} \\ z_{1} & z_{1}^{2} \end{array}\right]^{-1}
[GTG]−1=[N∑i=1Nzi∑i=1Nzi∑i=1Nzi2]−1→[1z1z1z12]−1
由于
A
−
1
=
A
∗
∣
A
∣
A^{-1}=\frac{A^{*}}{|A|}
A−1=∣A∣A∗
然而
∣
G
T
G
∣
=
1
z
1
2
−
z
1
2
=
∞
\mathbf{|G^TG|} = \frac{1}{z_1^2-z_1^2}=\infty
∣GTG∣=z12−z121=∞
显然,
[
G
T
G
]
−
1
\left[\mathbf{G}^{\mathrm{T}} \mathbf{G}\right]^{-1}
[GTG]−1是奇异的,也就意味着最小二乘法的公式失效了!!!!!
从反例来看,显然最小二乘法不能解决所有反问题,从这个角度上看,等式 G m = d Gm=d Gm=d没有提供了足够信息来唯一的确定模型参数 m \mathbf m m,基于这一点(也可以说,基于最小二乘法的解是否唯一),可以对反问题进行分类。
1) 恰定问题
G m = d \mathbf{Gm=d} Gm=d恰好存在足够的信息来 唯一的 确定模型参数,且具有零预测误差(e=0)。
2 )超定问题
当
G
m
=
d
\mathbf{Gm=d}
Gm=d拥有过多的信息来确定一个精确解时,则为超定问题(overdetermined)。这时,最小二乘法可以找到一个最佳近似解。
超定问题拥有的数据通常多于未知量(
N
>
M
N>M
N>M)。但并非是绝对的,有些问题即使当数据少于未知量(
N
<
M
N<M
N<M),在一定程度上依然是超定的。
3)欠定问题
当
G
m
=
d
\mathbf{Gm=d}
Gm=d不能提供足够的信息来唯一地确定模型参数是,则为欠定问题(underdetermined)。就如上如的反例,欠定发生在几个解通常具有零预测误差的情况。这时,最小二乘法是失效的。欠定问题通常未知量多于数据(
M
>
N
M>N
M>N)。但也存在一些问题,当
M
<
N
M<N
M<N时,它们在一定程度上也是欠定的,**在很多情况下,数据唯一的确定了某些模型参数,但不能确定其他模型参数。**如下图所示:
上图为一个声学实验,由于对第二块砖没有测定慢度,这个模型参数完全不受数据的约束。而第一块砖因多次测量慢度,是超定的。
h
[
1
0
1
0
⋮
⋮
1
0
]
[
s
1
s
2
]
=
[
d
1
d
2
⋮
d
N
]
h\left[\begin{array}{cc} 1 & 0 \\ 1 & 0 \\ \vdots & \vdots \\ 1 & 0 \end{array}\right]\left[\begin{array}{c} s_{1} \\ s_{2} \end{array}\right]=\left[\begin{array}{c} d_{1} \\ d_{2} \\ \vdots \\ d_{N} \end{array}\right]
h⎣⎢⎢⎢⎡11⋮100⋮0⎦⎥⎥⎥⎤[s1s2]=⎣⎢⎢⎢⎡d1d2⋮dN⎦⎥⎥⎥⎤
式中:
s
i
s_i
si为第
i
i
i块砖的慢度;
h
h
h为砖的宽度;
d
n
d_n
dn为第n次测量的旅行时。在用最小二乘法求解时,
[
G
T
G
]
−
1
\mathbf{{[G^TG]}^{-1}}
[GTG]−1是奇异的,此时即使
M
<
N
M<N
M<N,依然是欠定问题。这时常见的仅有部分模型参数是欠定的情况,在实际的情况中,还有很多欠定的形式更加不易察觉。
另外,并非所有的欠定问题都具有零预测误差。
为此,对欠定问题进一步分类:
- 将具有非零预测误差的欠定问题称为混定问题(mixed-determined problems);
- 将具有零预测误差的欠定问题称为纯欠定问题(purely undetermined problems);
从以上分类来看,单纯的最小二乘法无法应对欠定问题,需要其他的求解方法来应对欠定问题
3.6 纯欠定问题(purely undetermined problems)
纯欠定问题, G m = d \mathbf{Gm=d} Gm=d方程的数据的数量少于未知模型参数的数量,有多个解使预测误差为0( E = 0 \mathbf{E}=0 E=0)。对于这样的问题,我们必须使用某些手段,从无穷多个解中挑选出一个合适的解( m e s t m^{est} mest)。为此,我们必须添加一些没有包含在 G m = d \mathbf{Gm=d} Gm=d中的额外信息,即先验信息。
3.6.1 先验信息举例
先验信息可以采取任何形式,但在任何形式下,它仅描述了期望解有什么样的特征,它并非基于实际数据(即与实际数据无关)
例子1:
在拟合通过一个数据点的直线例子中,加入"期望直线穿过远点"的先验信息。此时,这个先验信息提供了足够的信息来找到唯一的解(因为,两点确定一条直线)。
例子2:
期望模型参数拥有给定的符号(如均为正/负)或位于一个给定的范围。
如假设模型参数代表了地球内部不同点的密度。即使没有测量,人们也可以肯定密度是正值。而且,由于地球内部岩石的密度必然处于一个密度范围之内,如
1000
∼
10000
k
g
/
m
3
1000 \sim 10000 kg/m^3
1000∼10000kg/m3。使用这样的先验信息可以极大的缩小解的范围,甚至使之唯一。
对于欠定问题,不得不将先验信息添加到反演问题,从而才能寻找到一个合适的解。但这种方法也存在着一些不确定性。
这些信息来自何处?它们的确定程度有多大?
这些问题没有固定答案,一些情况下,人们也许能找到合适的先验信息,有些时候则不能。如果人们想要论证估计值的唯一性,那么先验信息的有效性至关重要。
3.6.2 最小长度解
最为普遍的先验信息就是,希望反问题的解是简单的,这里所谓的“简单”是通过解的长度在某种度量下来量化的。其中采用欧几里得空间长度最为常见,即解的长度可定义为
L
=
m
T
m
=
∑
m
i
2
\mathbf{L}=\mathbf{m^Tm}=\sum{m_{i}^2}
L=mTm=∑mi2
我们所期望的解,其解的长度要最小。
虽然,这个度量也许不是对简单程度的实际度量,但它有时是有用的。在一些特定的问题,
L
\mathbf{L}
L可以是一个物理系统的能量(离散化后则代表功率–王家映)
使用先验信息来约束纯欠定反问题可表述如下:
寻找
m
e
s
t
\mathbf{m}^{est}
mest,以在
e
=
G
m
−
d
=
0
\mathbf{e=Gm-d=0}
e=Gm−d=0的约束下,使
L
=
m
T
m
\mathbf{L=m^Tm}
L=mTm最小化。
该问题可通过拉格朗日乘数法来求解:
Φ
(
m
)
=
L
+
∑
i
=
1
N
λ
i
e
i
=
∑
i
=
1
M
m
i
2
+
∑
i
=
1
N
λ
i
[
d
i
−
∑
j
=
1
M
G
i
j
m
j
]
\Phi(\mathbf{m})=L+\sum_{i=1}^{N} \lambda_{i} e_{i}=\sum_{i=1}^{M} m_{i}^{2}+\sum_{i=1}^{N} \lambda_{i}\left[d_{i}-\sum_{j=1}^{M} G_{i j} m_{j}\right]
Φ(m)=L+i=1∑Nλiei=i=1∑Mmi2+i=1∑Nλi[di−j=1∑MGijmj]
式中,
λ
i
\lambda_i
λi为拉格朗日乘数,求导可得
∂
Φ
∂
m
q
=
∑
i
=
1
M
2
∂
m
i
∂
m
q
m
i
−
∑
i
=
1
N
λ
i
∑
j
=
1
M
G
i
j
∂
m
j
∂
m
q
=
2
m
q
−
∑
i
=
1
N
λ
i
G
i
q
\frac{\partial \Phi}{\partial m_{q}}=\sum_{i=1}^{M} 2 \frac{\partial m_{i}}{\partial m_{q}} m_{i}-\sum_{i=1}^{N} \lambda_{i} \sum_{j=1}^{M} G_{i j} \frac{\partial m_{j}}{\partial m_{q}}=2 m_{q}-\sum_{i=1}^{N} \lambda_{i} G_{i q}
∂mq∂Φ=i=1∑M2∂mq∂mimi−i=1∑Nλij=1∑MGij∂mq∂mj=2mq−i=1∑NλiGiq
令
∂
Φ
∂
m
q
=
0
\frac{\partial \Phi}{\partial m_{q}}=0
∂mq∂Φ=0,并将其重新写成矩阵形式,得到方程
2
m
=
G
T
λ
2\mathbf{m}=\mathbf{G^T\lambda}
2m=GTλ
将其带入
d
=
G
m
\mathbf{d}=\mathbf{Gm}
d=Gm中,得到
d
=
G
m
=
G
[
G
T
λ
/
2
]
\mathbf{d}=\mathbf{Gm}=\mathbf{G[G^T\lambda/2]}
d=Gm=G[GTλ/2]
式中
G
G
T
\mathbf{GG^T}
GGT是一个
N
×
N
N \times N
N×N的方阵,若其存在逆矩阵,那么有
λ
=
2
[
G
G
T
]
−
1
d
\mathbf{\lambda}=2\mathbf{[GG^T]^{-1}d}
λ=2[GGT]−1d
将其回代到
2
m
=
G
T
λ
2\mathbf{m}=\mathbf{G^T\lambda}
2m=GTλ中,可得到
m
e
s
t
=
G
T
[
G
G
T
]
−
1
d
\mathbf{m}^{\mathrm{est}}=\mathbf{G}^{\mathrm{T}}\left[\mathbf{G} \mathbf{G}^{\mathrm{T}}\right]^{-1} \mathbf{d}
mest=GT[GGT]−1d
注意,上式的一个条件是 G m = d \mathbf{Gm=d} Gm=d是纯欠定的。
3.7 混定问题(mixed-determined problems)
实际中的大多数反问题,既不是完全超定的,也不是完全欠定的。如下例所示。
在X射线层析成像中,模型离散为多个小方格,有些方格有多条射线穿过(上图A),显然就这一个网格来说X射线不透明度(模型参数),是超定的。也有一些网格,完全没有X射线穿过(上图B),这些方格是完全欠定的。另外,还存在一些方格,因为射线从两个(或多个)方格中穿过(上图C),所以这些方格不能单独求取,这些方格也是欠定的,因为仅能确定它们的平均不透明度。
方法1
如果可能(理想情况下),我们希望将未知的模型参数分为两组:超定的和欠定的。这样将
G
m
=
d
=
=
>
G
′
m
′
=
d
′
\mathbf{Gm=d} ==>\mathbf{G^\prime m^\prime=d^\prime}
Gm=d==>G′m′=d′,其中
m
′
\mathbf{m^\prime}
m′分为超定的部分
m
′
o
\mathbf{m^{\prime o}}
m′o和欠定部分
m
′
u
\mathbf{m^{\prime u}}
m′u:
[
G
′
0
0
0
G
′
u
]
[
m
′
0
m
′
u
]
=
[
d
′
0
d
′
u
]
\left[\begin{array}{cc} \mathbf{G}^{\prime 0} & 0 \\ 0 & \mathbf{G}^{\prime \mathrm{u}} \end{array}\right]\left[\begin{array}{c} \mathbf{m}^{\prime 0} \\ \mathbf{m}^{\prime \mathrm{u}} \end{array}\right]=\left[\begin{array}{l} \mathbf{d}^{\prime 0} \\ \mathbf{d}^{\prime \mathrm{u}} \end{array}\right]
[G′000G′u][m′0m′u]=[d′0d′u]
这样上半部分求最小二乘解,下半部分求最小长度解。这样我们将找到一个解,而且给反问题添加了尽量少的先验信息。这个分割过程可以通过数据核的奇异值分解来实现(矩阵
G
G
G的奇异值分解就是把
m
′
o
\mathbf{m^{\prime o}}
m′o和
m
′
u
\mathbf{m^{\prime u}}
m′u分开的一种途径)
方法2
在不分割
m
\mathbf{m}
m的情况下,阻尼最小二乘法是最常用的。这时我们不仅要求预测误差(E)最小化,还要求解的长度(L)最小化,所以目标函数为二者的某种组合:
Φ
(
m
)
=
E
+
ε
2
L
=
e
T
e
+
ε
2
m
T
m
\Phi(\mathbf{m})=E+\varepsilon^{2} L=\mathbf{e}^{\mathrm{T}} \mathbf{e}+\varepsilon^{2} \mathbf{m}^{\mathrm{T}} \mathbf{m}
Φ(m)=E+ε2L=eTe+ε2mTm
式中,
ε
2
\varepsilon^2
ε2决定了预测误差和解的长度的相对重要性。若
ε
\varepsilon
ε足够大,则极小化的过程主要针对上式中第二项进行,这样可以明显的削弱解的欠定部分。但第一项权重相对减少,解估计的预测误差加大(不是最小的E)。结果就是,获得的解既不是最小二乘解,也不是最小长度解。若
ε
\varepsilon
ε设定为0,则最小化过程仅针对上式中的第一项,没有加入先验信息来挑选欠定部分的解。
与最小二乘的推导过程相似,可以给出阻尼最小二乘解的形式:
[
G
T
G
+
ε
2
I
]
m
e
s
t
=
G
T
d
or
m
e
s
t
=
[
G
T
G
+
ε
2
I
]
−
1
G
T
d
\left[\mathbf{G}^{\mathrm{T}} \mathbf{G}+\varepsilon^{2} \mathbf{I}\right] \mathbf{m}^{\mathrm{est}}=\mathbf{G}^{\mathrm{T}} \mathbf{d} \quad \text { or } \quad \mathbf{m}^{\mathrm{est}}=\left[\mathbf{G}^{\mathrm{T}} \mathbf{G}+\varepsilon^{2} \mathbf{I}\right]^{-1} \mathbf{G}^{\mathrm{T}} \mathbf{d}
[GTG+ε2I]mest=GTd or mest=[GTG+ε2I]−1GTd
阻尼最小二乘法可以挑选出
ε
\varepsilon
ε的一些折中值,以期在近似最小化解中欠定部分长度的基础上,近似最小化预测误差。然而,并没有直接的方法来求取折中值,它必须通过试错法来确定(实验/经验)。另外,值得注意的是,误差中不仅包含预测误差,还包含了拟合先验信息误差。
有时候,用
L
=
m
T
m
\mathbf{L=m^Tm}
L=mTm来衡量解的简单程度并不合理,例如,假设求解海洋中海水密度扰动这一反问题。我们不希望找到一个最接近0的解,而是想要找到一个接近于某一特定数值(如海水的平均密度)情况下的解。那么,对
L
\mathbf{L}
L自然可以推广为:
L
=
(
m
−
<
m
>
)
T
(
m
−
<
m
>
)
\mathbf{L}=\mathbf{(m-<m>)}^{T}\mathbf{(m-<m>)}
L=(m−<m>)T(m−<m>)
式中,
<
m
>
<m>
<m>为模型参数的先验值(在例子中为已知海水的平均值)。
3.8 利用先验信息加权
用来衡量解的简单程度并不只有长度这一度量,平直度(flat)和光滑度(smooth)也是经常用到的两种方式。
一个空间连续函数的平直度可以通过一阶导数(即陡度(stepness)(平直度的反义词))的范数来量化。对于离散化的模型参数,可以用模型空间(物理)上相邻的模型参数的差值作为一阶导数的近似:
那么向量
m
\mathbf{m}
m的陡度
I
\mathbf{I}
I是
I
=
1
Δ
x
[
−
1
1
−
1
1
⋱
⋱
−
1
]
[
m
1
m
2
⋮
m
M
]
=
D
m
\mathbf{I}=\frac{1}{\Delta x}\left[\begin{array}{cccc} -1 & 1 & & \\ & -1 & 1 & \\ & & \ddots & \ddots \\ & & & -1 \end{array}\right]\left[\begin{array}{c} m_{1} \\ m_{2} \\ \vdots \\ m_{M} \end{array}\right]=\mathbf{D}{\mathbf{m}}
I=Δx1⎣⎢⎢⎡−11−11⋱⋱−1⎦⎥⎥⎤⎣⎢⎢⎢⎡m1m2⋮mM⎦⎥⎥⎥⎤=Dm
解的光滑度(smoothness)可以通过使用二阶导数量化粗糙度(roughness)(光滑度的反义词)来实现。
于是,解的总陡度/总粗糙度是如下长度:
L
=
I
T
I
=
[
D
m
]
T
[
D
m
]
=
m
T
D
T
D
m
=
m
T
W
m
m
L=\mathbf{I}^{\mathrm{T}} \mathbf{I}=[\mathbf{D m}]^{\mathrm{T}}[\mathbf{D m}]=\mathbf{m}^{\mathrm{T}} \mathbf{D}^{\mathrm{T}} \mathbf{D} \mathbf{m}=\mathbf{m}^{\mathrm{T}} \mathbf{W}_{\mathrm{m}} \mathbf{m}
L=ITI=[Dm]T[Dm]=mTDTDm=mTWmm
矩阵
W
m
=
D
T
D
\mathbf{W_m=D^TD}
Wm=DTD可以当作参与计算模型参数向量
m
\mathbf{m}
m长度的权重因子。
因此,解的简单程度的度量
L
L
L可以推广为:
L
=
[
m
−
⟨
m
⟩
]
T
W
m
[
m
−
⟨
m
⟩
]
L=[\mathbf{m}-\langle\mathbf{m}\rangle]^{\mathrm{T}} \mathbf{W}_{\mathrm{m}}[\mathbf{m}-\langle\mathbf{m}\rangle]
L=[m−⟨m⟩]TWm[m−⟨m⟩]
另外,预测误差的加权度量也是有用的。有一些观测经常比另一些做的更准确,在这种情况下,我们希望较精确观测的预测误差
e
i
e_i
ei在总误差
E
E
E的量化过程中比不精确观测具有更大的权重。为了实现这样的权重,我们定义一个广义的预测误差:
E
=
e
T
W
e
e
E=\mathbf{e}^{\mathrm{T}} \mathbf{W}_{\mathrm{e}} \mathbf{e}
E=eTWee
式中,矩阵
W
e
\mathbf{W_e}
We定义了各单独测量误差对总预测误差的相对贡献。通常,该矩阵为一对角阵。
例如,共有5个测量数据(N=5)且第三次观测精度是其他几次观测精度的2倍,那么可以使用:
W
e
=
[
1
1
2
1
1
]
\mathbf{W}_{\mathrm{e}}=\left[\begin{array}{cccc} 1 & & & \\ & 1 & & \\ & & 2 & \\ & & 1 & \\ & & & 1 \end{array}\right]
We=⎣⎢⎢⎢⎢⎡11211⎦⎥⎥⎥⎥⎤
这样,可以通过加权的方式,进一步改进反问题的解:
-
1)加权最小二乘==》完全超定问题
通过最小化广义观测误差 E = e T W e e \mathbf{E = e^TW_ee} E=eTWee来求解:
m e s t = [ G T W e G ] − 1 G T W e d \mathbf{m}^{\mathrm{est}}=\left[\mathbf{G}^{\mathrm{T}} \mathbf{W}_{\mathrm{e}} \mathbf{G}\right]^{-1} \mathbf{G}^{\mathrm{T}} \mathbf{W}_{\mathrm{e}} \mathbf{d} mest=[GTWeG]−1GTWed -
2)加权最小长度==》纯欠定问题
通过最小化广义长度 L = [ m − ⟨ m ⟩ ] T W m [ m − ⟨ m ⟩ ] T L=[\mathbf{m}-\langle\mathbf{m}\rangle]^{\mathrm{T}} \mathbf{W}_{\mathrm{m}}[\mathbf{m}-\langle\mathbf{m}\rangle]^{\mathrm{T}} L=[m−⟨m⟩]TWm[m−⟨m⟩]T来求解:
m e s t = ⟨ m ⟩ + W m − 1 G T [ G W m − 1 G T ] − 1 [ d − G ⟨ m ⟩ ] \mathbf{m}^{\mathrm{est}}=\langle\mathbf{m}\rangle+\mathbf{W}_{\mathrm{m}}^{-1} \mathbf{G}^{\mathrm{T}}\left[\mathbf{G} \mathbf{W}_{\mathrm{m}}^{-1} \mathbf{G}^{\mathrm{T}}\right]^{-1}[\mathbf{d}-\mathbf{G}\langle\mathbf{m}\rangle] mest=⟨m⟩+Wm−1GT[GWm−1GT]−1[d−G⟨m⟩] -
3)加权阻尼最小二乘法==》混定问题(稍微欠定)
通过最小化预测误差和解的长度的组合 Φ ( m ) = E + ε 2 L \Phi(\mathbf{m})=E+\varepsilon^{2} L Φ(m)=E+ε2L来求解:
[ G T W e G + ε 2 W m ] m e s t = G T W e d + ε 2 W m ⟨ m ⟩ \left[\mathbf{G}^{\mathrm{T}} \mathbf{W}_{\mathrm{e}} \mathbf{G}+\varepsilon^{2} \mathbf{W}_{\mathrm{m}}\right] \mathbf{m}^{\mathrm{est}}=\mathbf{G}^{\mathrm{T}} \mathbf{W}_{\mathrm{e}} \mathbf{d}+\varepsilon^{2} \mathbf{W}_{\mathrm{m}}\langle\mathbf{m}\rangle [GTWeG+ε2Wm]mest=GTWed+ε2Wm⟨m⟩
或
m e s t = [ G T W e G + ε 2 W m ] − 1 [ G T W e d + ε 2 W m ⟨ m ⟩ ] \mathbf{m}^{\mathrm{est}}=\left[\mathbf{G}^{\mathrm{T}} \mathbf{W}_{\mathrm{e}} \mathbf{G}+\varepsilon^{2} \mathbf{W}_{\mathrm{m}}\right]^{-1}\left[\mathbf{G}^{\mathrm{T}} \mathbf{W}_{\mathrm{e}} \mathbf{d}+\varepsilon^{2} \mathbf{W}_{\mathrm{m}}\langle\mathbf{m}\rangle\right] mest=[GTWeG+ε2Wm]−1[GTWed+ε2Wm⟨m⟩]