共轭梯度法解求解大规模稀疏矩阵，对比最速梯度法（C++）-Toy模板网

这篇具有很好参考价值的文章主要介绍了共轭梯度法解求解大规模稀疏矩阵，对比最速梯度法（C++）。希望对大家有所帮助。如果存在错误或未考虑完全的地方，请大家不吝赐教，您也可以点击"举报违法"按钮提交疑问。

记录计算方法大作业，练习C++，欢迎指正。

1，共轭梯度法介绍

共轭梯度法（Conjugate Gradient）是介于最速下降法与牛顿法之间的一个方法，它仅需利用一阶导数信息，但克服了最速下降法收敛慢的缺点，又避免了牛顿法需要存储和计算Hesse矩阵并求逆的缺点。共轭梯度法不仅是解决大型线性方程组最有用的方法之一，也是解大型非线性最优化最有效的算法之一。

在实际应用中，共轭梯度法不仅可以去求解方程组，还可以推广到非二次目标函数的极小值求解。在各种优化算法中，共轭梯度法是非常重要的一种。其优点是所需存储量小，具有步收敛性，稳定性高，而且不需要任何外来参数。

2，共轭梯度法原理

求解 Ax=b时，最简单粗暴的方式为 x=A-1 *b。但是这种方法的问题很明显：矩阵求逆的计算复杂度非常高。即使我们考虑用矩阵分解的方式，仍然会很慢。因此，我们尽可能考虑用迭代的方式，而不是直接求逆的方式来解这个问题。如果构造一个二次函数：f(x)=xT Ax/2-b T x,求其最小值，即其导数为0时正好是方程Ax=b的解。因此，我们可以将线性方程组求解问题转化为二次函数求极小值问题

在寻优过程中利用当前点xk 大规模稀疏矩阵,c++,算法,矩阵处的残向量rk 和前一迭代点xk-1 处的搜索方向dk-1 对搜索方向进行如下修正：

𝑑𝑘 = 𝑟𝑘 + 𝛽𝑘−1𝑑𝑘−1

其修正系数𝛽𝑘−1的取值有一个约束条件，即要确保𝑑𝑘与𝑑𝑘 , 𝑑𝑘−1 ,⋯, 𝑑0之间满足关于 A 的共轭关系。这就是共轭梯度法的基本思想。可以看出共轭梯度法的搜索方向𝑑𝑘的计算只需要梯度向量，不需要海森矩阵 H，可以推广到非二次目标函数的极小值求解。

下面给出共轭梯度法求解方程组的伪代码：

(1) 给定初始近似向量𝑥 (0)以及精度𝜖 > 0;

(2) 计算𝑟 (0) = 𝑏 − 𝐴𝑥 (0)，𝑑 (0) = 𝑟 (0)；

(3) for k=0 to n-1 do

(i) 𝛼𝑘 = 𝑟 (𝑘)𝑇𝑟 (𝑘) 𝑑(𝑘)𝑇𝐴𝑑(𝑘) ;

(ii) 𝑥 (𝑘+1) = 𝑥 (𝑘) + 𝛼𝑘𝑑 (𝑘) ;

(iii) 𝑟 (𝑘+1) = 𝑏 − 𝐴𝑥 (𝑘+1) ;

(iv) 如果|| 𝑟 (𝑘+1) || < 𝜖或者 k+1=n，则输出近似解𝑥 (𝑘+1)，停止；否则转(v)

(v) 𝛽𝑘 = || 𝑟 (𝑘+1) ||2 2 || 𝑟 (𝑘) ||2 2 ;

(vi) 𝑑 (𝑘+1) = 𝑟 (𝑘+1) + 𝛽𝑘𝑑 (𝑘) ;

end do

3，共轭梯度法和最速下降法的比较

本次计算的实例采用的是课本第 133 页中的例 3.2，分别分析了 n 的取值为 100，200，400 时的情况，并同时用共轭梯度法和最速下降法去求解线性方程，设置两种方法的精度为 0.0001。

3.1程序说明

该程序只包含一个 cpp 文件。请将要求解的方程组的系数矩阵保存在对应文件夹下的 data.txt 中。程序中有共轭梯度法和最速下降法的函数，两种方法求得的结果和对应的迭代误差保存在对应文件夹下的 ConjugateGradientMethod.txt 和 SteepestDescentMethod.txt 中。文件夹保存和读取的路径都放在程序开头的全局变量中，可以自行修改。

3.2，迭代误差

大规模稀疏矩阵,c++,算法,矩阵

图 1-1 共轭梯度法的迭代误差变化图 1-1 是共轭梯度法迭代过程中的误差变化，可以看出该方法最终收敛，且共轭梯度解出来的结果是一个全 1 的向量，也就是所有的未知数都是 1。通过最速下降法得到的解为[0.5,0,0,0,0,0,……,0,0,0,0.5],即 x1 和 x100 是 0.5，其他 x 的值为 0，带入到原方程组中发现符合解特征。而且最速下降法只迭代一次就得到最终解

（2）N=200

大规模稀疏矩阵,c++,算法,矩阵

（3）N=400

大规模稀疏矩阵,c++,算法,矩阵

3.3两种方法比较

（1）迭代结果

大规模稀疏矩阵,c++,算法,矩阵

（2）迭代速度

大规模稀疏矩阵,c++,算法,矩阵

图 1-4 是两种方法解方程组迭代收敛的过程，可以看出共轭梯度在进行四次迭代时就已经收敛，最速下降法在 9 次迭代之后才收敛。从迭代误差的变化中可以发现此时梯度下降法的每次迭代的误差总是比上次小，但是最速下降法每次迭代的误差还可能比前一次大，说明共轭梯度在每次迭代都是向着准确解的方向移动。

最速下降法相邻两次迭代的方向是正交的，这就说明最速下降法每次迭代的方向并不是都朝着最优解的方向，而且这种正交性会导致最速下降法向极小点逼近是曲折前进的，这种现象会影响收敛速度。同时这种迭代方式很可能使迭代陷入局部最优解，影响优化的结果。

4，代码

主函数：

int main()
{
	//vector<vector<double>> A = { {2,0,1},{0,1,0},{1,0,2} };
	//vector<double> B = { 3,1,3 };
	// 
	//定义矩阵
	vector<vector<double>> A;
	vector<double> B;


	//加载算例3.2
	//SuanLi32(A, B,400);
	//加载数据，数据保存路径在全局变量S中
	Load_data(A, B);

	//定义精度
	double Accuracy = 0.0001;
	
	//定义结果x，和迭代的误差erro
	vector<double> x, erro;
	vector<double> x1, erro1;

	//梯度下降法，结果保存在路径S的ConjugateGradientMethod.txt
	ConjugateGradientMethod(A, B, Accuracy, x, erro);

	//显示结果
	for_each(x.begin(), x.end(), PrintV);
	cout << endl;
	for_each(erro.begin(), erro.end(), PrintV);

		//最速下降法，结果保存在路径S的ConjugateGradientMethod.txt
	SteepestDescentMethod(A, B, Accuracy, x1, erro1);

	//显示结果
	for_each(x1.begin(), x1.end(), PrintV);
	cout << endl;
	for_each(erro1.begin(), erro1.end(), PrintV);

	system("pause");
	return 0;
}

共轭梯度法

void ConjugateGradientMethod(vector<vector<double>>& A, vector<double>& B, 
	double& Accuracy, vector<double>& x, vector<double>& erro)
{
	int n = A.size();
	vector<double> r(n), d(n),Ad(n),xx(n),ad(n),Ax,r1,betad,d1;
	double rr,r1r,dAd,a,beta=0,err=1;
	//赋初值
	x.assign(n, 0);
	r.assign(B.begin(), B.end());
	d.assign(B.begin(), B.end());
	while (err> Accuracy)
	{
		rr = Vect_vec(r, r);
		Ad = Mat_vector(A, d);
		dAd = Vect_vec(d, Ad);
		a = rr / dAd;//i
		ad = Constant_vector(a, d);
		xx = Add_vector(ad, x);//ii
		Ax = Mat_vector(A, xx);
		r1 = Sub_vector(B, Ax);//iii
		r1r = Vect_vec(r1, r1);
		beta = r1r / rr;
		betad = Constant_vector(beta, d);
		d1 = Add_vector(r1, betad);
		err = sqrt(r1r);
		//迭代
		r.assign(r1.begin(), r1.end());
		d.assign(d1.begin(), d1.end());
		x.assign(xx.begin(), xx.end());
		erro.push_back(err);
	}

	//保存结果
	ofstream ofs;
	ofs.open("ConjugateGradientMethod.txt", ios::out | ios::app);
	if (!ofs.is_open())
	{
		cout << "储存文件打开失败！" << endl;
	}
	else
	{
		for (int i = 0; i < n; i++)
		{
			ofs << x[i] << "\t";
		}
		ofs << endl;
		int l = erro.size();
		for (int i = 0; i < l; i++)
		{
			ofs << erro[i] << "\t";
		}
	}
	ofs.close();

	return;
}

void PrintV(const double& a)
{
	cout << a << " ";
}

最速梯度法：

void SteepestDescentMethod(vector<vector<double>>& A, vector<double>& B,
	double& Accuracy, vector<double>& x, vector<double>& erro)
{
	int n = A.size();
	vector<double> r(n), d(n), Ar(n), xx(n), ar(n), Ax, r1, betad, d1;
	double  a=0, err = 1,a1,a2;
	//赋初值
	x.assign(n, 0);
	r.assign(B.begin(), B.end());
	while (err > Accuracy)
	{
		Ar = Mat_vector(A, r);
		a=Vect_vec(r,r)/ Vect_vec(Ar,r);
		//a2=inner_product(Ar.begin(), Ar.end(), r.begin(), 0) ;
		//a = double(inner_product(r.begin(), r.end(), r.begin(), 0) )/double(inner_product(Ar.begin(), Ar.end(), r.begin(), 0));
		ar = Constant_vector(a, r);
		xx = Add_vector(x, ar);
		Ax = Mat_vector(A, xx);
		r1 = Sub_vector(B, Ax);//iii
		err = sqrt(inner_product(r1.begin(),r1.end(),r1.begin(),0));
		//迭代
		r.assign(r1.begin(), r1.end());
		x.assign(xx.begin(), xx.end());
		erro.push_back(err);
		//i++; cout << err << endl;
	}

	//保存结果
	ofstream ofs;
	ofs.open("SteepestDescentMethod.txt", ios::out | ios::app);
	if (!ofs.is_open())
	{
		cout << "储存文件打开失败！" << endl;
	}
	else
	{
		for (int i = 0; i < n; i++)
		{
			ofs << x[i] << "\t";
		}
		ofs << endl;
		int l = erro.size();
		for (int i = 0; i < l; i++)
		{
			ofs << erro[i] << "\t";
		}
	}
	ofs.close();

	return;
}

其他矩阵向量或者函数

vector<double> Mat_vector(vector<vector<double>>& A, vector<double>& d)
{
	vector<double> Ad;
	int m = A.size();
	for (int i = 0; i < m; i++)
	{
		double temp = 0;
		for (int j = 0; j < m; j++)
		{
			temp += A[i][j] * d[j];
		}
		Ad.push_back(temp);
	}
	return Ad;
}

double Vect_vec(vector<double>& a1, vector<double>& a2)
{
	int m = a1.size();
	double a=0;
	for (int i = 0; i < m; i++)
	{
		a += a1[i] * a2[i];
	}
	return a;
}

vector<double> Sub_vector(vector<double>& a1, vector<double>& a2)
{
	int m = a1.size();
	vector<double> a;
	double temp;
	for (int i = 0; i < m; i++)
	{
		temp = a1[i] - a2[i];
		a.push_back(temp);
	}
	return a;
}

vector<double> Add_vector(vector<double>& a1, vector<double>& a2)
{
	int m = a1.size();
	vector<double> a;
	double temp;
	for (int i = 0; i < m; i++)
	{
		temp = a1[i] + a2[i];
		a.push_back(temp);
	}
	return a;
}

vector<double> Constant_vector(double &c,vector<double>& a1)
{
	int m = a1.size();
	vector<double> a;
	double temp;
	for (int i = 0; i < m; i++)
	{
		temp = a1[i]*c;
		a.push_back(temp);
	}
	return a;
}


void SuanLi32(vector<vector<double>>& A, vector<double>& B, int n)
{
	//定义B矩阵
	B.push_back(-1);
	for (int i = 1; i < n - 1; i++)
	{
		B.push_back(0);
	}
	B.push_back(-1);
	//定义A矩阵
	vector<double> a1(n - 2, 0);
	a1.insert(a1.begin(),1);
	a1.insert(a1.begin(), -2);
	A.push_back(a1);
	for (int i = 1; i < n - 1; i++)
	{
		vector<double> a;
		for (int j = 0; j < n; j++)
		{
			if (j == i)
			{
				a.push_back(-2);
			}
			else if (j == i - 1 || j == i + 1)
			{
				a.push_back(1);
			}
			else
			{
				a.push_back(0);
			}
		}
		A.push_back(a);
	}
	vector<double> a2(n - 2, 0);
	a2.push_back(1);
	a2.push_back(-2);
	A.push_back(a2);
}

加载数据

void Load_data(vector<vector<double>>& A, vector<double>& B)
{
	ifstream ifs;
	int n = 0;
	string s;
	double a,b;
	//行数n
	ifs.open(path, ios::in);
	if (ifs.fail())
	{
		cout << "文件打开失败！"<<endl;
	}
	else
	{
		while (getline(ifs, s))
		{
			n++;
		}
	}
	ifs.close();
	ifs.open(path, ios::in);

	while (!ifs.eof())
	{
		for (int i = 0; i < n; i++)
		{
			vector<double> aa;
			for (int j = 0; j < n; j++)
			{
				ifs >> a;
				aa.push_back(a);
			}
			ifs >> b;
			B.push_back(b);
			A.push_back(aa);
		}
	}
	ifs.close();
	return;
}

头文件和全局变量文章来源地址https://www.toymoban.com/news/detail-715259.html

#include<iostream>
#include<vector>
#include<cmath>
#include<algorithm>
#include<fstream>
#include<string>
#include<numeric>

using namespace std;

string path = "data.txt";

到了这里，关于共轭梯度法解求解大规模稀疏矩阵，对比最速梯度法（C++）的文章就介绍完了。如果您还想了解更多内容，请在右上角搜索TOY模板网以前的文章或继续浏览下面的相关文章，希望大家以后多多支持TOY模板网！