http://www.ma-xy.com
目录
第一章 非线性最小二乘优化 1
1.1 理论基础 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Gauss-Newton 法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.3 Levenerg-Marquardt . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.4 MATLAB 应用实例 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1
http://www.ma-xy.com
第一章 非线性最小二乘优化
1.1 理论基础
下面来看一个特别的无约束非线性规划 - 非线性最小二乘优化。此优化问题在数据拟合、
参数估计和函数逼近等问题中经常遇到，有较高的应用价值。设 r
i
(x) : R
n
→ R 是 x 的函数
(i = 0, 1, 2, . . . , m)，最小二乘问题描述为
min
x
f(x) =
1
2
r(x)
T
r(x) =
1
2
m
∑
i=1
[r
i
(x)]
2
m ⩾ n
(1) 如果 r
i
(x) 是 x 的线性函数
r
i
(x) = a
T
i
x − b
i
其中：a
i
∈ R
n
, b
i
∈ R，则称为线性最小二乘规划问题。容易证明，线性最小二乘问题是一
个凸二次规划。
(2) 如果 r
i
(x) 是 x 的非线性函数，则称为非线性最小二乘规划。由于最小二乘优化是无约
束非线性规划的一个特例，所以前面介绍的方法也可以适用，但由其特殊性，因此，它会有些适
用于自身的特殊的方法。下面，我们将介绍一些求解非线性最小二乘的算法。我们先来给出 r(x)
的 Jacobi 矩阵的定义：
定义 (Jacobi 矩阵) 连续函数 r : R
n
→ R
m
在 x ∈ R
n
连续可微，如果其每一个分量 r
i
(x)
在 x 连续可微。r 在 x 的导数 r
′
(x) ∈ R
m×n
叫做 r 在 x 的 Jacobi 矩阵，它的转置叫做 r 在 x
的梯度，即
r
′
(x) = J(x) = ∇r(x)
T
Jacobi 矩阵的第 i, j 元素为
[r
′
(x)]
ij
= [J(x)]
ij
=
∂r
i
∂x
j
(x) i = 1, . . . , m, j = 1, . . . , n
设 J(x) 为 r(x) 的 Jacobi 矩阵，则目标函数 f 的梯度为
g(x) =
m
∑
i=1
r
i
(x)∇r
i
(x) = J(x)
T
r(x)
1
http://www.ma-xy.com
1.2 GAUSS-NEWTON 法 第一章 非线性最小二乘优化
f 的 Hesse 矩阵为
G(x) =
m
∑
i=1
(∇r
i
(x)∇r
i
(x)
T
+ r
i
(x)∇
2
r
i
(x))
= J(x)
T
J(x) + s(x)
其中：
s(x) =
m
∑
i=1
r
i
(x)∇
2
r
i
(x)
上面，我们给出了目标函数 f 的梯度 g(x) 和 Hesse 矩阵 G(x)。我们写出目标函数 f 的二
次模型
m
k
(x) = f(x
k
) + g(x
k
)
T
(x − x
k
) +
1
2
(x − x
k
)
T
G(x
k
)(x − x
k
)
=
1
2
r(x
k
)
T
r(x
k
) + (J(x
k
)
T
r(x
k
))
T
(x − x
k
)
+
1
2
(x − x
k
)
T
(J(x
k
)
T
J(x
k
) + s(x
k
))(x − x
k
)
从而，解决非线性最小二乘的牛顿法为
x
k+1
= x
k
− (J(x
k
)
T
J(x
k
) + s(x
k
)
−1
J(x
k
)r(x
k
))
我们知道牛顿法具有二阶收敛速度，但是，上述牛顿迭代格式的主要问题是 Hesse 矩阵 G(x)
中的二阶信赖域 s(x) 通常难以计算。而如果仅对 G(x) 近似 (拟牛顿) 又有些浪费，毕竟，我们在
计算 g(x) 时已经得到 J(x)，而 J
T
(x)J(x) 是 G(x) 的一阶信息项。鉴于此，我们或者忽略 s(x)，
或者用一阶导数信息逼近 s(x)。
1.2 Gauss-Newton 法
下面介绍的 Gauss-Newton 法相当于目标函数的二次模型 m
k
(x) 中忽略 G(x) 中的二阶信
息项 s(x)，这样 m
k
(x) 变为
¯m
k
(x) =
1
2
r(x
k
)
T
r(x
k
) + (J(x
k
)
T
r(x
k
))
T
(x − x
k
)
+
1
2
(x − x
k
)
T
(J(x
k
)
T
J(x
k
))(x − x
k
) (1.1)
由此得到的牛顿迭代公式为
x
k+1
= x
k
− (J(x
k
)
T
J(x
k
))
−1
J(x
k
)r(x
k
)
= x
k
+ s
k
其中：s
k
= −(J(x
k
)
T
J(x
k
))
−1
J(x
k
)r(x
k
)。
http://www.ma-xy.com 2 http://www.ma-xy.com
http://www.ma-xy.com
第一章 非线性最小二乘优化 1.3 LEVENERG-MARQUARDT
而模型 (1.1) 相当于 r(x) 在 x
k
附近的仿射模型
˜
M
k
(x) = r(x
k
) + J(x
k
)(x − x
k
)
从而求下面线性最小二乘问题的解
min
1
2
∥
˜
M
k
(x)∥
2
的解。
从 Gauss-Newton 法的迭代公式中可以看出，该方法仅需要残差函数 r(x) 的一阶导数信息，
并且 J(x)
T
J(x) 至少是正半定的。
如果 s(x
∗
) = 0，则 G-N 方法是二阶收敛的。如果 s(x
∗
) 相当于 J(x
∗
)
T
J(x
∗
) 是小的，则
G(x) 方法是局部 Q 线性收敛的。但如果 s(x
∗
) 太大，则 G-N 方法可能不收敛。下面，我们给出
G-N 方法的优缺点：
(1) 当 r(x
∗
) = 0 时，有局部二阶收敛速度；
(2) 当 r(x
∗
) 较小时，有快的局部收敛速度；
(3) 当 r(x
∗
) 不是很大时，有较慢的局部收敛速度；
(4) 当 r(x
∗
) 很大时，有不收敛；
(5) 如果 J(x
k
) 不满秩，方法没有定义；
(6) G-N 不一定总体收敛。
1.3 Levenerg-Marquardt
在 Gauss-Newton 方法中，我们要求 J(x
∗
) 是满秩的。遗憾的是，J(x
∗
) 不满秩的情况是经
常发生的。一旦 J (x
∗
) 奇异，则在距离解点的某处，s
k
与 g
k
便数值上直交（正交）。这样，由
线搜索就得不到进一步下降，为了克服这种困难，考虑采用信赖域策略。其理由是：通常 r(x)
是非线性函数，而 Gauss-Newton 法用线性化模型
˜
M
k
(x) 代替 r(x)，但这种线性化并不对所有
(x − x
k
) 都成立，因此，我们考虑约束线性最小二乘问题，即考虑信赖域模型：
min ∥r(x
k
) + J(x
k
)(x − x
k
)∥
2
s.t. ∥x − x
k
∥
2
⩽ h
k
由前面的信赖域算法，我们知道这个模型的解可以由解方程组
(J(x
k
)
T
J(x
k
) + µ
k
I)s = −J(x
k
)
T
r(x
k
)
来表示，从而
x
k+1
= x
k
− (J(x
k
)
T
J(x
k
) + µ
k
I)
−1
J(x
k
)
T
r(x
k
)
如果 ∥J(x
k
)
T
J(x
k
))
−1
J(x
k
)
T
r(x
k
)∥ ⩽ h
k
，则 µ
k
= 0，否则 µ
k
> 0。由于 J(x
k
)
T
J(x
k
) + µ
k
I 正
定，所以上面信赖域模型产生的方向 s 是下降的，此方向由 Levenberg(1944) 和 Marqurdt(1963)
提出，所以又称 L-M 方法。
http://www.ma-xy.com 3 http://www.ma-xy.com
http://www.ma-xy.com
1.4 MATLAB 应用实例 第一章 非线性最小二乘优化
1.4 MATLAB 应用实例
MATLAB 中用 lsqnonlin 函数来求解非线性最小二乘优化问题，其调用格式为
[x,resnorm,residual,exitag,output,lambda,jacobian]=lsqnonlin(fun,x0,lb,ub,options)
其中：resnorm 为残差平方和，也即为最小值 r(x)
T
r(x)；residual 为残差 r(x)；lambda 返
回最优解 x 处的拉格朗日乘子；jacobian 为最优解 x 处的雅克比矩阵。
我们用 lsqnonlin 求解如下非线性最小二乘问题
min
x
f
2
(x) =
m
∑
i=1
f
2
i
(x)
其中:
f(x) =
sin(x
1
+ x
2
− 2)
1
−(x
1
− 3)
2
+ 2
e
2x
1
+ e
2x
2
x
2
1
+ x
2
2
− x
1
x
2
+ x
1
+ 1
求解程序为
1 x0 = [ 0 , 0 ] ;
2 fun = @(x) [
3 sin (x (1 )+x (2 )−2) ;
4 1/(2−(x (1 )−3)^2) ;
5 exp (2 (x (1) )+exp(2−x(2 ) ) ;
6 x( 1)^2+x( 2)^2−x (1) *x (2 )+x (1 ) +1];
7 o ption s = optimoptions ( ’ ls q no nl in ’ , ’ Display ’ , ’ i t e r ’ ) ;
8 o ption s . Algorithm = ’ Levenberg−Marquardt ’
9 [ x , resnorm , r esidu al , e x i t fl a g , output , lambda , jac ob ian ] = lsq no nl i n ( fun , x0 , IJ , IJ , op ti ons )
10
http://www.ma-xy.com 4 http://www.ma-xy.com
