最小二乘数值优化方法--梯度法、牛顿法、高斯牛顿法理论与c++实现

最小二乘常见求解方法

定义目标函数：
$$\text{F(x)}=\frac{1}{n}\sum _{i=0}^n(y_i-\widetilde{y_i}{)}^2=\frac{1}{n}\sum _{i=0}^n(h(x_i)-\widetilde{y_i}{)}^2\text{\hspace{1em}(1)}$$
目的是使得目标函数最小的情况下，求得未知量。

1梯度下降法

梯度的本意是一个向量（矢量），表示某一函数在该点处的方向导数沿着该方向取得最大值，即函数在该点处沿着梯度方向变化最快，变化率最大（为该梯度的模）。根据梯度的含义，要是$F(x)$最小，很自然地想到，沿着负梯度方法不断变化，能最快的到达极小值点，这就是梯度下降的算法思路。因此很容易总结出梯度下降迭代公式：
$$x_{n+1}=x_n-a.\nabla F(x_n)\text{\hspace{1em}(2)}$$
$a为学习率（决定了下降的速度），\nabla F(x_n)为梯度$。

线性回归案例：

有一个一元一次函数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)=ax+b$,(注意：这里a,b是待求量)所以式(1)变为：
$$F(a,b)=\frac{1}{n}\sum _{i=0}^n(a{x}_i+b-\widetilde{y_i}{)}^2\text{\hspace{1em}(3)}$$
求偏导$\nabla F(a)=\frac{2}{n}{\sum }_{i=0}^n(a{x}_i+b-\widetilde{y_i})x_i,\nabla F(b)=\frac{2}{n}{\sum }_{i=0}^n(a{x}_i+b-\widetilde{y_i})$，得到$F(a,b)的梯度为(\nabla F(a),\nabla 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\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\text{\hspace{1em}(4)}$$
对于$F(x)$而言，求最小值的必要条件是让其一阶偏导为0.将其转换为求根过程：
$$F^{\prime }(x)=F^{\prime }\left(x_0\right)+F^{\prime \prime }\left(x_0\right)\left(x-x_0\right)=0\text{\hspace{1em}(5)}$$
移向，可以得到
$$x=x_0-\frac{F^{\prime }(x_0)}{F^{\prime \prime }(x_0)}\text{\hspace{1em}(6)}$$
得到牛顿法迭代公式：
$$x_{n+1}=x_n-\frac{F^{\prime }(x_n)}{F^{\prime \prime }(x_n)}\text{\hspace{1em}(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^{\prime \prime }(a)=\frac{2}{n}\sum x_i,F^{\prime \prime }(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高斯牛顿法

令$\text{F(x)}=\frac{1}{2}||f(x)|{|}_2^2$，常系数不影响优化，为了方便计算将常数$\frac{1}{n}$换成$\frac{1}{2}$，$f(x)=h(x)-\tilde{y}$.

那么，我们对 $f(x+\Delta x)$ 进行一阶泰勒展开。
$$f(x+\Delta x)\approx f(x)+J(x{)}^T\Delta x+o(\Delta x)\text{\hspace{1em}(8)}$$
上式也可以利用导数的定义来推导：
$$f^{\prime }(x)=\frac{f(x+\Delta x)-f(x)}{\Delta x}\text{\hspace{1em}(9)}$$
可以得到：
$$f(x+\Delta x)\approx f(x)+J(x{)}^T\Delta x\text{\hspace{1em}(10)}$$
$J(x)$为一阶偏导矩阵形式，又称为雅可比矩阵。

我们需要求 $\Delta x$ 使得上面的式子 $\Vert f(x+\Delta x){\Vert }_2^2$ 有最小值，所以，我们可以得到最小二乘问题为:
$$\Delta x^{\ast }=\arg \min \frac{1}{2}\Vert f(x+\Delta x){\Vert }_2^2\approx \arg \min \frac{1}{2}{\Vert f(x)+J(x{)}^T\Delta x\Vert }_2^2\text{\hspace{1em}(11)}$$
展开:
$$\begin{array}{rl}& m(\Delta x)=\frac{1}{2}{\Vert f(x)+J(x{)}^T\Delta x\Vert }^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}(\Vert f(x){\Vert }^2+2f(x)J(x{)}^T\Delta x+\Delta x^{T}J(x)J(x{)}^T\Delta x\end{array}\text{\hspace{1em}(12)}$$
故，对$\Delta x$求导可以得到:
$$m^{\prime }(\Delta x)=J(x)f(x)+J(x)J(x{)}^T\Delta x\text{\hspace{1em}(13)}$$
则，此时可以转化为线性求解问题:
$$m^{\prime }(\Delta x)=0\to J(x)J(x{)}^T\Delta x=-J(x)f(x)\text{\hspace{1em}(14)}$$
所以有：
$$\Delta x=(J(x)J(x{)}^T{)}^{-1}(-J(x)f(x))\text{\hspace{1em}(15)}$$
所以高斯牛顿法过程如下：

(1)给定初始值$x_0$

(2)求解$J(x_0),f(x_0)$

(3)求解$\Delta x=(J(x)J(x{)}^T{)}^{-1}(-J(x)f(x))$

(4)如果$\Delta x$足够小，则终止迭代

(5)否则，令$x_{n+1}=x_n+\Delta 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)=ax+b$,(注意：这里a,b是待求量)所以式(1)变为：
$$F(a,b)=\frac{1}{2}\sum _{i=0}^n(a{x}_i+b-\widetilde{y_i}{)}^2$$
其中$f(a,b)=h(x)-\tilde{y}=a{x}_i+b-\widetilde{y_i}$.

所以偏导(雅可比矩阵)
$$J=\left[\begin{array}{c}{f}^{\prime }(a)\\ f^{\prime }(b)\end{array}\right]={\left[\begin{array}{c}{x}_0...x_i...x_n\\ \mathrm{1..1..1}\end{array}\right]}_{2\times n}$$

$$J{J}^T=\left[\begin{array}{cc}{\sum }_{i=0}^n{x}_i^2& {\sum }_{i=0}^n{x}_i\\ {\sum }_{i=0}^n{x}_i& n+1\end{array}\right]$$

$$G=-J(a,b)f(a,b)=\left[\begin{array}{c}-{\sum }_{i=0}^n{x}_i(a{x}_i+b-y_i)\\ -{\sum }_{i=0}^n(a{x}_i+b-y_i)\end{array}\right]$$

对于二阶矩阵而言
$${\left(\begin{array}{ll}a& b\\ c& d\end{array}\right)}^{-1}=\frac{1}{ad-bc}\left(\begin{array}{cc}d& -b\\ -c& a\end{array}\right)$$

$$H=(J{J}^T{)}^{-1}=\frac{1}{(n+1){\sum }_{i=0}^n{x}_i^2-({\sum }_{i=0}^n{x}_i{)}^2}\left[\begin{array}{cc}n+1& -{\sum }_{i=0}^n{x}_i\\ -{\sum }_{i=0}^n{x}_i& {\sum }_{i=0}^n{x}_i^2\end{array}\right]$$
这时候就可以用解析方法求得$\Delta a,\Delta 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)$进行优化从而求解未知量的。梯度下降和高斯牛顿法属于一阶方法，梯度下降收敛较慢，高斯牛顿法需要保证矩阵可逆，同时雅可比矩阵会占用大量内存空间。牛顿法属于二阶方法，需要保证二阶偏导（黑森矩阵）存在。
LM算法后续可能会更新。

参考链接：

1

2

3

4