共轭梯度法（FR法）

       优化算法的学习是有一定难度的（也许有人不这么想，但即使你是天赋异禀的天之骄子，敬畏知识也是一种美德）。笔者本着敬畏知识的态度，从基本理论出发，逐步引入优化算法的殿堂。

     关于共轭梯度法理论的介绍，可以通过参考文献【1】了解

    这里我通过一个例子来介绍共轭梯度法

                  (1)

                   (2)     

   公式（1）（2）给出了2D点与2D直线的表达方法

   现有多条直线，求一个，使其到每条直线的距离的平方和最小（这是一道图像算法入门的经典例题）

  对每条直线进行归一化有

令 

  即求 

令

又可知

算法过程：

1.给定初始点

2.计算出目标函数在这一点的梯度

3求得步长 

4.沿方向

搜索，得到点

5.

,

6.

7.

8.

 即求得最优解

通过代码实现这个过程

// 如果你认为有必要打赏我，我的支付宝账号是[email protected]
//我是cclplus，一个热衷于研究算法的男人
#include <iostream>
#include <cmath>
#include <memory.h>
using namespace std;
double a[5] = { 1,2,3,4,5 };
double b[5] = { 2,1,5,3,4 };
double c[5] = { 6,7,8,9,10 };
double _A[3][3];
int main()
{
	memset(_A, 0, sizeof(_A));
	double temp;
	//进行归一化
	for (int i = 0; i < 5; i++) {
		temp = sqrt(a[i] * a[i] + b[i] * b[i]);
		a[i] /= temp;
		b[i] /= temp;
		c[i] /= temp;
	}
	//构建矩阵A(_A)
	for (int i = 0; i < 5; i++) {
		_A[0][0] += a[i] * a[i];
		_A[0][1] += a[i] * b[i];
		_A[0][2] += a[i] * c[i];
		_A[1][0] += b[i] * a[i];
		_A[1][1] += b[i] * b[i];
		_A[1][0] += b[i] * c[i];
		_A[2][0] += c[i] * a[i];
		_A[2][0] += c[i] * b[i];
		_A[2][0] += c[i] * c[i];
	}
	for (int i = 0; i < 3; i++) {
		for (int j = 0; j < 3; j++) {
			_A[i][j] *= 2;
		}
	}
	//设定初始点
	double x_pre[2] = { 1,1 };
	double x_new[2];
	//计算梯度
	double d_f1, d_f2;
	d_f1 = _A[0][0] * x_pre[0] + _A[0][1] * x_pre[1] + 0.5*_A[0][2];
	d_f2= _A[1][0] * x_pre[0] + _A[1][1] * x_pre[1] + 0.5*_A[1][2];
	//求步长
	double lambda= d_f1* d_f1+ d_f2* d_f2;
	temp = d_f1 * d_f1*_A[0][0] + 2 * d_f1*d_f2*_A[1][0] + d_f2 * d_f2*_A[1][1];
	lambda /= temp;
	x_new[0] = x_pre[0] + lambda * d_f1;
	x_new[1] = x_pre[1] + lambda * d_f2;
	//计算新的梯度
	double d_f1n, d_f2n;
	d_f1n = _A[0][0] * x_new[0] + _A[0][1] * x_new[1] + 0.5*_A[0][2];
	d_f2n = _A[1][0] * x_new[0] + _A[1][1] * x_new[1] + 0.5*_A[1][2];
	lambda = (d_f1n*d_f1n + d_f2n * d_f2n) / (d_f1*d_f1 + d_f2 * d_f2);
	x_pre[0] = x_new[0];
	x_pre[1] = x_new[1];
	double g_1, g_2;
	g_1 = -d_f1n + lambda * d_f1;
	g_2= -d_f2n + lambda * d_f2;
	 lambda = g_1 * g_1 + g_2 * g_2;
	temp = g_1 * g_1*_A[0][0] + 2 * g_1*g_2*_A[1][0] + g_2 * g_2*_A[1][1];
	lambda /= temp;
	double ans[2];
	ans[0] = x_pre[0] + lambda * g_1;
	ans[1] = x_pre[1] + lambda * g_2;
	cout << ans[0] << " " << ans[1] << endl;
	return 0;
}

           

参考文献

【1】 最优化理论与算法 陈宝林 清华大学出版社