非线性最小二乘法 是非线性 形式的最小二乘法 ,用包含n 个未知参数的非线性模型拟合m 个观测值(
m
≥
n
{\displaystyle m\geq n}
),可用于某些形式的非线性回归 。该方法的基础是使用线性模型近似并通过连续迭代来优化参数。它与线性最小二乘法既有相同之处、也有一些显著差异。
考虑一组
(
x
1
,
y
1
)
,
(
x
2
,
y
2
)
,
…
,
(
x
m
,
y
m
)
{\displaystyle (x_{1},y_{1}),(x_{2},y_{2}),\dots ,(x_{m},y_{m})}
共
m
{\displaystyle m}
个数据点以及曲线(模型函数)
y
^
=
f
(
x
,
β
)
{\displaystyle {\hat {y}}=f(x,{\boldsymbol {\beta }})}
。该曲线同时取决于x 与
β
=
(
β
1
,
β
2
,
…
,
β
n
)
{\displaystyle {\boldsymbol {\beta }}=(\beta _{1},\beta _{2},\dots ,\beta _{n})}
共n 个参数(满足
m
≥
n
{\displaystyle m\geq n}
)。目标是找到在最小二乘意义上与数据点拟合最好的曲线所对应的参数
β
{\displaystyle {\boldsymbol {\beta }}}
,即最小化平方和
S
=
∑
i
=
1
m
r
i
2
,
{\displaystyle S=\sum _{i=1}^{m}r_{i}^{2},}
其中残差 ri 的定义为
r
i
=
y
i
−
f
(
x
i
,
β
)
,
(
i
=
1
,
2
,
…
,
m
)
.
{\displaystyle r_{i}=y_{i}-f(x_{i},{\boldsymbol {\beta }}),\qquad (i=1,2,\dots ,m).}
S 取最小值 时的梯度 为零。由于模型包含n 个参数,因此可得到n 个梯度方程:
∂
S
∂
β
j
=
2
∑
i
r
i
∂
r
i
∂
β
j
=
0
(
j
=
1
,
…
,
n
)
.
{\displaystyle {\frac {\partial S}{\partial \beta _{j}}}=2\sum _{i}r_{i}{\frac {\partial r_{i}}{\partial \beta _{j}}}=0\quad (j=1,\ldots ,n).}
在非线性系统中,偏导数
∂
r
i
∂
β
j
{\textstyle {\frac {\partial r_{i}}{\partial \beta _{j}}}}
同时是自变量x 和参数
β
{\displaystyle {\boldsymbol {\beta }}}
的函数,因此这些梯度方程通常没有封闭解。因而必须为参数选择初始值用以迭代求解。迭代表达式为
β
j
≈
β
j
k
+
1
=
β
j
k
+
Δ
β
j
.
{\displaystyle \beta _{j}\approx \beta _{j}^{k+1}=\beta _{j}^{k}+\Delta \beta _{j}.}
其中,k 是迭代次数,
Δ
β
{\displaystyle \Delta {\boldsymbol {\beta }}}
则是偏移向量。每次迭代时,使用关于
β
k
{\displaystyle {\boldsymbol {\beta }}^{k}}
的一阶泰勒级数 展开以线性化模型:
f
(
x
i
,
β
)
≈
f
(
x
i
,
β
k
)
+
∑
j
∂
f
(
x
i
,
β
k
)
∂
β
j
(
β
j
−
β
j
k
)
=
f
(
x
i
,
β
k
)
+
∑
j
J
i
j
Δ
β
j
.
{\displaystyle f(x_{i},{\boldsymbol {\beta }})\approx f(x_{i},{\boldsymbol {\beta }}^{k})+\sum _{j}{\frac {\partial f(x_{i},{\boldsymbol {\beta }}^{k})}{\partial \beta _{j}}}\left(\beta _{j}-\beta _{j}^{k}\right)=f(x_{i},{\boldsymbol {\beta }}^{k})+\sum _{j}J_{ij}\,\Delta \beta _{j}.}
雅可比矩阵 J 是常数、自变量与参数的函数,因此每次迭代时的J 并不固定。对线性化模型而言,
∂
r
i
∂
β
j
=
−
J
i
j
,
{\displaystyle {\frac {\partial r_{i}}{\partial \beta _{j}}}=-J_{ij},}
残差的表达式则为
Δ
y
i
=
y
i
−
f
(
x
i
,
β
k
)
,
{\displaystyle \Delta y_{i}=y_{i}-f(x_{i},{\boldsymbol {\beta }}^{k}),}
r
i
=
y
i
−
f
(
x
i
,
β
)
=
(
y
i
−
f
(
x
i
,
β
k
)
)
+
(
f
(
x
i
,
β
k
)
−
f
(
x
i
,
β
)
)
≈
Δ
y
i
−
∑
s
=
1
n
J
i
s
Δ
β
s
.
{\displaystyle r_{i}=y_{i}-f(x_{i},{\boldsymbol {\beta }})=\left(y_{i}-f(x_{i},{\boldsymbol {\beta }}^{k})\right)+\left(f(x_{i},{\boldsymbol {\beta }}^{k})-f(x_{i},{\boldsymbol {\beta }})\right)\approx \Delta y_{i}-\sum _{s=1}^{n}J_{is}\Delta \beta _{s}.}
将上述表达式代入梯度方程,可以得到
−
2
∑
i
=
1
m
J
i
j
(
Δ
y
i
−
∑
s
=
1
n
J
i
s
Δ
β
s
)
=
0
,
{\displaystyle -2\sum _{i=1}^{m}J_{ij}\left(\Delta y_{i}-\sum _{s=1}^{n}J_{is}\ \Delta \beta _{s}\right)=0,}
以上方程可化简为n 个联立的线性方程,称为正规方程 (normal equations):
∑
i
=
1
m
∑
s
=
1
n
J
i
j
J
i
s
Δ
β
s
=
∑
i
=
1
m
J
i
j
Δ
y
i
(
j
=
1
,
…
,
n
)
.
{\displaystyle \sum _{i=1}^{m}\sum _{s=1}^{n}J_{ij}J_{is}\ \Delta \beta _{s}=\sum _{i=1}^{m}J_{ij}\ \Delta y_{i}\qquad (j=1,\dots ,n).}
正规方程可用矩阵表示法写成
(
J
T
J
)
Δ
β
=
J
T
Δ
y
.
{\displaystyle \left(\mathbf {J} ^{\mathsf {T}}\mathbf {J} \right)\Delta {\boldsymbol {\beta }}=\mathbf {J} ^{\mathsf {T}}\ \Delta \mathbf {y} .}
上述方程是使用高斯-牛顿算法 求解非线性最小二乘问题的的基础。
需要注意的是雅可比矩阵定义中导数的符号约定。某些文献中的J 可能与此处的定义相差一个负号。
不同数据点(观测结果)的可靠性并不一定相同,此时可使用加权平方和
S
=
∑
i
=
1
m
W
i
i
r
i
2
.
{\displaystyle S=\sum _{i=1}^{m}W_{ii}r_{i}^{2}.}
权重矩阵W 是一个对角矩阵 ,理想情况下每个权重系数应等于观测误差方差 的倒数。[ 1] 此时,正规方程可扩展为
(
J
T
W
J
)
Δ
β
=
J
T
W
Δ
y
.
{\displaystyle \left(\mathbf {J} ^{\mathsf {T}}\mathbf {WJ} \right)\Delta {\boldsymbol {\beta }}=\mathbf {J} ^{\mathsf {T}}\mathbf {W} \ \Delta \mathbf {y} .}
^ 此处假定所有观测点是相互独立的。如果观测点之间相关 时,加权平方和可表示为
S
=
∑
k
∑
j
r
k
W
k
j
r
j
.
{\displaystyle S=\sum _{k}\sum _{j}r_{k}W_{kj}r_{j}.}
此时权重矩阵的理想值应为观测误差协方差矩阵 的逆。