最小二乘数值优化方法--梯度法、牛顿法、高斯牛顿法理论与c++实现
最小二乘常见求解方法
定义目标函数：
F(x)
=
1
n
∑
i
=
0
n
(
y
i
?
y
i
~
)
2
=
1
n
∑
i
=
0
n
(
h
(
x
i
)
?
y
i
~
)
2
(1)
\text{F(x)}=\frac{1}{n}\sum_{i = 0}^{n}(y_i-\tilde{y_i})^2=\frac{1}{n}\sum_{i = 0}^{n}(h(x_i)-\tilde{y_i})^2\tag{1}
F(x)=n1?i=0∑n?(yi??yi?~?)2=n1?i=0∑n?(h(xi?)?yi?~?)2(1)
目的是使得目标函数最小的情况下，求得未知量。
1梯度下降法
梯度的本意是一个向量（矢量），表示某一函数在该点处的方向导数沿着该方向取得最大值，即函数在该点处沿着梯度方向变化最快，变化率最大（为该梯度的模）。根据梯度的含义，要是
F
(
x
)
F(x)
F(x)最小，很自然地想到，沿着负梯度方法不断变化，能最快的到达极小值点，这就是梯度下降的算法思路。因此很容易总结出梯度下降迭代公式：
x
n
+
1
=
x
n
?
a
.
?
F
(
x
n
)
(2)
x_{n+1}=x_n-a.\nabla F({x_n})\tag{2}
xn+1?=xn??a.?F(xn?)(2)
a
为学习率（决定了下降的速度），
?
F
(
x
n
)
为梯度
a为学习率（决定了下降的速度），\nabla F({x_n})为梯度
a为学习率（决定了下降的速度），?F(xn?)为梯度。
线性回归案例：
有一个一元一次函数y=ax+b有如下观测值：
x=[1.0,2.1,2.9,3.03,5.01,8.093]
y=[3.02,4.97,7.1,6.88,10.88,17.06]
求解a,b。
对于式(1)有
h
(
a
,
b
)
=
a
x
+
b
h(a,b)=ax+b
h(a,b)=ax+b,(注意：这里a,b是待求量)所以式(1)变为：
F
(
a
,
b
)
=
1
n
∑
i
=
0
n
(
a
x
i
+
b
?
y
i
~
)
2
(3)
F(a,b)=\frac{1}{n}\sum_{i = 0}^{n}(ax_i+b-\tilde{y_i})^2\tag{3}
F(a,b)=n1?i=0∑n?(axi?+b?yi?~?)2(3)
求偏导
?
F
(
a
)
=
2
n
∑
i
=
0
n
(
a
x
i
+
b
?
y
i
~
)
x
i
,
?
F
(
b
)
=
2
n
∑
i
=
0
n
(
a
x
i
+
b
?
y
i
~
)
\nabla F({a})=\frac{2}{n}\sum_{i = 0}^{n}(ax_i+b-\tilde{y_i})x_i,\nabla F({b})=\frac{2}{n}\sum_{i = 0}^{n}(ax_i+b-\tilde{y_i})
?F(a)=n2?∑i=0n?(axi?+b?yi?~?)xi?,?F(b)=n2?∑i=0n?(axi?+b?yi?~?)，得到
F
(
a
,
b
)
的梯度为
(
?
F
(
a
)
,
?
F
(
b
)
)
F(a,b)的梯度为(\nabla F({a}),\nabla F({b}))
F(a,b)的梯度为(?F(a),?F(b)).
求和可以写成向量点乘，因此偏导也可以写成矩阵的形式，具体的可以看[1](梯度下降算法 线性回归拟合（附Python/Matlab/Julia源代码） - 知乎 (zhihu.com))。
codes：
double a = 0.0,b=0.0;
double rate = 0.01;
const int size = 6;
double x[size]{ 1.0,2.1,2.9,3.03,5.01,8.093 };
double y[size]{ 3.02,4.97,7.1,6.88,10.88,17.06 };
//double x[size]{ 1,2,3,5.01,8.093 };
//double y[size]{ 3,5,7,10.88,17.06 };
while (true)
{
//求偏导
double da1 = 0.0, db1 = 0.0;
for (size_t j = 0; j < size; j++)
{
da1 += (a * x[j] + b - y[j]) * x[j];
db1 += (a * x[j] + b - y[j]);
}
//偏导
double d_a = da1 / double(size);
double d_b = db1 / double(size);
if (std::abs(d_a)<1e-3&&std::abs(d_b)<1e-3)
{
break;
}
//更新未知量
a = a - rate * d_a;
b = b - rate * d_b;
std::cout << "a:" << a << "b:" << b << std::endl;
}
得到：a:1.98363 b:1.00005，效果如下：
2牛顿法
将
F
(
x
)
F(x)
F(x)按照二阶泰勒多项式展开：
F
(
x
)
=
F
(
x
0
)
+
F
′
(
x
0
)
(
x
?
x
0
)
+
F
′
′
(
x
0
)
2
(
x
?
x
0
)
2
(4)
F(x)=F\left(x_0\right)+F^{\prime}\left(x_0\right)\left(x-x_0\right)+\frac{F^{\prime \prime}\left(x_0\right)}{2}\left(x-x_0\right)^2\tag{4}
F(x)=F(x0?)+F′(x0?)(x?x0?)+2F′′(x0?)?(x?x0?)2(4)
对于
F
(
x
)
F(x)
F(x)而言，求最小值的必要条件是让其一阶偏导为0.将其转换为求根过程：
F
′
(
x
)
=
F
′
(
x
0
)
+
F
′
′
(
x
0
)
(
x
?
x
0
)
=
0
(5)
F^{\prime}(x)=F^{\prime}\left(x_0\right)+F^{\prime\prime}\left(x_0\right)\left(x-x_0\right)=0\tag{5}
F′(x)=F′(x0?)+F′′(x0?)(x?x0?)=0(5)
移向，可以得到
x
=
x
0
?
F
′
(
x
0
)
F
′
′
(
x
0
)
(6)
x=x_0-\frac{F^{\prime}(x_0)}{F^{\prime\prime}(x_0)}\tag{6}
x=x0??F′′(x0?)F′(x0?)?(6)
得到牛顿法迭代公式：
x
n
+
1
=
x
n
?
F
′
(
x
n
)
F
′
′
(
x
n
)
(7)
x_{n+1}=x_n-\frac{F^{\prime}(x_n)}{F^{\prime\prime}(x_n)}\tag{7}
xn+1?=xn??F′′(xn?)F′(xn?)?(7)
同样使用上面的案例：
一个一元一次函数y=ax+b有如下实际值：
x=[1.0,2.1,2.9,3.03,5.01,8.093]
y=[3.02,4.97,7.1,6.88,10.88,17.06]
求解a,b。
显然对于牛顿法，需要求解二阶偏导，根据上面一阶偏导的可以得到二阶偏导：
F ′ ′ ( a ) = 2 n ∑ x i , F ′ ′ ( b ) = 2 F^{\prime\prime}(a)=\frac{2}{n}\sum{x_i},F^{\prime\prime}(b)=2 F′′(a)=n2?∑xi?,F′′(b)=2.
codes
double a = 0.0, b = 0.0;
double rate = 0.01;
const int size = 6;
double x[size]{ 1.0,2.1,2.9,3.03,5.01,8.093 };
double y[size]{ 3.02,4.97,7.1,6.88,10.88,17.06 };
//double x[size]{ 1,2,3,5.01,8.093 };
//double y[size]{ 3,5,7,10.88,17.06 };
//for (size_t i = 0; i < iter_num; i++)
while (true)
{
//求偏导
double da1 = 0.0, db1 = 0.0;
double dd_a = 0.0, dd_b = 1.0;
for (size_t j = 0; j < size; j++)
{
da1 += (a * x[j] + b - y[j]) * x[j];
db1 += (a * x[j] + b - y[j]);
dd_a += x[j] * x[j];
}
//偏导
double d_a = da1 / double(size);
double d_b = db1 / double(size);
//二阶偏导
dd_a=dd_a/ double(size);
if (std::abs(d_a) < 1e-3 && std::abs(d_b) < 1e-3)
{
break;
}
//
a = a - d_a / dd_a;
b = b - d_b / dd_b;
std::cout << "a:" << a << "b:" << b << std::endl;
}
结果：
a:1.98291 b:1.00392，相比于梯度下降，牛顿法收敛更快，但是计算量更大，要求存在二阶偏导。
3高斯牛顿法
令 F(x) = 1 2 ∣ ∣ f ( x ) ∣ ∣ 2 2 \text{F(x)}=\frac{1}{2}||f(x)||^2_2 F(x)=21?∣∣f(x)∣∣22?，常系数不影响优化，为了方便计算将常数 1 n \frac{1}{n} n1?换成 1 2 \frac{1}{2} 21?， f ( x ) = h ( x ) ? y ~ f(x)=h(x)-\tilde{y} f(x)=h(x)?y~?.
那么，我们对
f
(
x
+
Δ
x
)
f(x+\Delta x)
f(x+Δx) 进行一阶泰勒展开。
f
(
x
+
Δ
x
)
≈
f
(
x
)
+
J
(
x
)
T
Δ
x
+
o
(
Δ
x
)
(8)
f(x+\Delta x) \approx f(x)+J(x)^T \Delta x+o(\Delta x)\tag{8}
f(x+Δx)≈f(x)+J(x)TΔx+o(Δx)(8)
上式也可以利用导数的定义来推导：
f
′
(
x
)
=
f
(
x
+
Δ
x
)
?
f
(
x
)
Δ
x
(9)
f^{\prime}(x)=\frac{f(x+\Delta x)-f(x)}{\Delta x}\tag{9}
f′(x)=Δxf(x+Δx)?f(x)?(9)
可以得到：
f
(
x
+
Δ
x
)
≈
f
(
x
)
+
J
(
x
)
T
Δ
x
(10)
f(x+\Delta x) \approx f(x)+J(x)^T \Delta x\tag{10}
f(x+Δx)≈f(x)+J(x)TΔx(10)
J
(
x
)
J(x)
J(x)为一阶偏导矩阵形式，又称为雅可比矩阵。
我们需要求
Δ
x
\Delta x
Δx 使得上面的式子
∥
f
(
x
+
Δ
x
)
∥
2
2
\|f(x+\Delta x)\|_2^2
∥f(x+Δx)∥22? 有最小值，所以，我们可以得到最小二乘问题为:
Δ
x
?
=
arg
?
min
?
1
2
∥
f
(
x
+
Δ
x
)
∥
2
2
≈
arg
?
min
?
1
2
∥
f
(
x
)
+
J
(
x
)
T
Δ
x
∥
2
2
(11)
\Delta x^*=\arg \min \frac{1}{2}\|f(x+\Delta x)\|_2^2 \approx \arg \min \frac{1}{2}\left\|f(x)+J(x)^T \Delta x\right\|_2^2\tag{11}
Δx?=argmin21?∥f(x+Δx)∥22?≈argmin21?
?f(x)+J(x)TΔx
?22?(11)
展开:
m
(
Δ
x
)
=
1
2
∥
f
(
x
)
+
J
(
x
)
T
Δ
x
∥
2
=
1
2
(
f
(
x
)
+
J
(
x
)
T
Δ
x
)
T
(
f
(
x
)
+
J
(
x
)
T
Δ
x
)
=
1
2
(
∥
f
(
x
)
∥
2
+
2
f
(
x
)
J
(
x
)
T
Δ
x
+
Δ
x
T
J
(
x
)
J
(
x
)
T
Δ
x
(12)
\begin{aligned}& m(\Delta x)=\frac{1}{2}\left\|f(x)+J(x)^T \Delta x\right\|^2=\frac{1}{2}\left(f(x)+J(x)^T \Delta x\right)^T\left(f(x)+J(x)^T \Delta x\right) \\& =\frac{1}{2}\left(\|f(x)\|^2+2 f(x) J(x)^T \Delta x+\Delta x^T J(x) J(x)^T \Delta x\right.\end{aligned}\tag{12}
?m(Δx)=21?
?f(x)+J(x)TΔx
?2=21?(f(x)+J(x)TΔx)T(f(x)+J(x)TΔx)=21?(∥f(x)∥2+2f(x)J(x)TΔx+ΔxTJ(x)J(x)TΔx?(12)
故，对
Δ
x
\Delta x
Δx求导可以得到:
m
′
(
Δ
x
)
=
J
(
x
)
f
(
x
)
+
J
(
x
)
J
(
x
)
T
Δ
x
(13)
m^{\prime}(\Delta x)=J(x) f(x)+J(x) J(x)^T \Delta x\tag{13}
m′(Δx)=J(x)f(x)+J(x)J(x)TΔx(13)
则，此时可以转化为线性求解问题:
m
′
(
Δ
x
)
=
0
→
J
(
x
)
J
(
x
)
T
Δ
x
=
?
J
(
x
)
f
(
x
)
(14)
m^{\prime}(\Delta x)=0 \rightarrow J(x) J(x)^T \Delta x=-J(x) f(x)\tag{14}
m′(Δx)=0→J(x)J(x)TΔx=?J(x)f(x)(14)
所以有：
Δ
x
=
(
J
(
x
)
J
(
x
)
T
)
?
1
(
?
J
(
x
)
f
(
x
)
)
(15)
\Delta x=(J(x) J(x)^T)^{-1}(-J(x) f(x))\tag{15}
Δx=(J(x)J(x)T)?1(?J(x)f(x))(15)
所以高斯牛顿法过程如下：
(1)给定初始值 x 0 x_0 x0?
(2)求解 J ( x 0 ) , f ( x 0 ) J(x_0),f(x_0) J(x0?),f(x0?)
(3)求解 Δ x = ( J ( x ) J ( x ) T ) ? 1 ( ? J ( x ) f ( x ) ) \Delta x=(J(x) J(x)^T)^{-1}(-J(x) f(x)) Δx=(J(x)J(x)T)?1(?J(x)f(x))
(4)如果 Δ x \Delta x Δx足够小，则终止迭代
(5)否则，令 x n + 1 = x n + Δ x x_{n+1}=x_n+\Delta x xn+1?=xn?+Δx,返回(2)
同样使用上面的案例：
一个一元一次函数y=ax+b有如下实际值：
x=[1.0,2.1,2.9,3.03,5.01,8.093]
y=[3.02,4.97,7.1,6.88,10.88,17.06]
求解a,b。
对于式(1)有
h
(
a
,
b
)
=
a
x
+
b
h(a,b)=ax+b
h(a,b)=ax+b,(注意：这里a,b是待求量)所以式(1)变为：
F
(
a
,
b
)
=
1
2
∑
i
=
0
n
(
a
x
i
+
b
?
y
i
~
)
2
F(a,b)=\frac{1}{2}\sum_{i = 0}^{n}(ax_i+b-\tilde{y_i})^2
F(a,b)=21?i=0∑n?(axi?+b?yi?~?)2
其中
f
(
a
,
b
)
=
h
(
x
)
?
y
~
=
a
x
i
+
b
?
y
i
~
f(a,b)=h(x)-\tilde{y}=ax_i+b-\tilde{y_i}
f(a,b)=h(x)?y~?=axi?+b?yi?~?.
所以偏导(雅可比矩阵)
J
=
[
f
′
(
a
)
f
′
(
b
)
]
=
[
x
0
.
.
.
x
i
.
.
.
x
n
1..1..1
]
2
×
n
J=\left[\begin{matrix}f^\prime(a)\\f^\prime(b)\end{matrix}\right]=\left[\begin{matrix}x_0...x_i...x_n\\1..1..1\end{matrix}\right]_{2\times n}
J=[f′(a)f′(b)?]=[x0?...xi?...xn?1..1..1?]2×n?
J J T = [ ∑ i = 0 n x i 2 ∑ i = 0 n x i ∑ i = 0 n x i n + 1 ] JJ^T=\left[\begin{matrix}\sum_{i = 0}^{n}x_i^2 &\sum_{i = 0}^{n}x_i\\\sum_{i = 0}^{n}x_i&n+1\end{matrix}\right] JJT=[∑i=0n?xi2?∑i=0n?xi??∑i=0n?xi?n+1?]
G = ? J ( a , b ) f ( a , b ) = [ ? ∑ i = 0 n x i ( a x i + b ? y i ) ? ∑ i = 0 n ( a x i + b ? y i ) ] G=-J(a,b)f(a,b)=\left[\begin{matrix}-\sum_{i = 0}^{n}x_i(ax_i+b-y_i)\\-\sum_{i = 0}^{n}(ax_i+b-y_i)\end{matrix}\right] G=?J(a,b)f(a,b)=[?∑i=0n?xi?(axi?+b?yi?)?∑i=0n?(axi?+b?yi?)?]
对于二阶矩阵而言
(
a
b
c
d
)
?
1
=
1
a
d
?
b
c
(
d
?
b
?
c
a
)
\left(\begin{array}{ll}a & b \\c & d\end{array}\right)^{-1}=\frac{1}{a d-b c}\left(\begin{array}{cc}d & -b \\-c & a\end{array}\right)
(ac?bd?)?1=ad?bc1?(d?c??ba?)
H
=
(
J
J
T
)
?
1
=
1
(
n
+
1
)
∑
i
=
0
n
x
i
2
?
(
∑
i
=
0
n
x
i
)
2
[
n
+
1
?
∑
i
=
0
n
x
i
?
∑
i
=
0
n
x
i
∑
i
=
0
n
x
i
2
]
H=(JJ^T)^{-1}=\frac{1}{(n+1)\sum_{i = 0}^{n}x_i^2-(\sum_{i = 0}^{n}x_i)^2}\left[\begin{matrix} n+1&-\sum_{i = 0}^{n}x_i\\-\sum_{i = 0}^{n}x_i&\sum_{i = 0}^{n}x_i^2\end{matrix}\right]
H=(JJT)?1=(n+1)∑i=0n?xi2??(∑i=0n?xi?)21?[n+1?∑i=0n?xi???∑i=0n?xi?∑i=0n?xi2??]
这时候就可以用解析方法求得
Δ
a
,
Δ
b
\Delta a,\Delta b
Δa,Δb，然后进行迭代了。
codes:
double a = -0.0, b = 0.0;
const int size = 6;
double x[size]{ 1.0,2.1,2.9,3.03,5.01,8.093 };
double y[size]{ 3.02,4.97,7.1,6.88,10.88,17.06 };
double sum_xi = 0.0;
double sum_xi_sqr = 0.0;
for (int i = 0; i < size; ++i)
{
sum_xi += x[i];
sum_xi_sqr += x[i] * x[i];
}
//计算逆矩阵H=(JJ^T)^(-1)
double h00 = 1.0 / (double(size) * sum_xi_sqr - sum_xi * sum_xi) * double(size);
double h01 = 1.0 / (double(size) * sum_xi_sqr - sum_xi * sum_xi) * (-1.0) * sum_xi;
double h10 = h01;
double h11 = 1.0 / (double(size) * sum_xi_sqr - sum_xi * sum_xi) * sum_xi_sqr;
while (true)
{
//计算G
double g00 = 0.0, g10 = 0.0;
for (size_t j = 0; j < size; j++)
{
g00 += x[j] * (a * x[j] + b - y[j]);
g10 += a * x[j] + b - y[j];
}
g00 *= -1.0;
g10 *= -1.0;
double delta_a = h00 * g00 + h01 * g10;
double delta_b = h10 * g00 + h11 * g10;
a += delta_a;
b += delta_b;
if (std::abs(delta_a)<1e-3&& std::abs(delta_b) < 1e-3)
{
break;
}
std::cout << "a:" << a << "b:" << b << std::endl;
}
结果：a:1.9829b:1.00373。仅用一次就收敛了。迭代速度更快。
梯度下降和牛顿法是对目标函数
F
(
x
)
F(x)
F(x)优化，高斯牛顿法是对误差函数
f
(
x
)
f(x)
f(x)进行优化从而求解未知量的。梯度下降和高斯牛顿法属于一阶方法，梯度下降收敛较慢，高斯牛顿法需要保证矩阵可逆，同时雅可比矩阵会占用大量内存空间。牛顿法属于二阶方法，需要保证二阶偏导（黑森矩阵）存在。
LM算法后续可能会更新。
参考链接：
本文来自互联网用户投稿，该文观点仅代表作者本人，不代表本站立场。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。 如若内容造成侵权/违法违规/事实不符，请联系我的编程经验分享网邮箱：chenni525@qq.com进行投诉反馈，一经查实，立即删除！