c++ 实现计算特征值和特征向量（matlab eig实现）

原文链接

#pragma once

#include <iostream>
#include <vector>
#include <map>
#include <iomanip>

using namespace std;

namespace Algorithm
{
	using namespace std;

	bool My_Jacobi(const vector<vector<double>>& matrix\
		, vector<vector<double>>& eigenvectors\
		, vector<double>& eigenvalues, double precision = 1e-10, int max = 10000)
	{
		int dim = matrix.size();
		//vector<vector<double>> mat(matrix.size(), vector<double>(matrix[0].size(), 0.0));
		vector<vector<double>> mat = matrix;
		eigenvectors = vector<vector<double>>(matrix.size(), vector<double>(matrix[0].size(), 0.0));
		// 初始化特征矩阵
		for (int i = 0; i < dim; i++)
		{
			eigenvectors[i][i] = 1.0;
		}
		eigenvalues = vector<double>(matrix.size());

		int nCount = 0;  //current iteration
		while (1)
		{
			// 搜索矩阵绝对值最大元素及下标
			double dbMax = mat[0][1];
			int nRow = 0;
			int nCol = 1;
			for (int i = 0; i < dim; i++) //row
			{
				for (int j = 0; j < dim; j++)
				{
					double d = fabs(mat[i][j]);
					if ((i != j) && (d > dbMax))
					{
						dbMax = d;
						nRow = i;
						nCol = j;
					}
				}
			}
			cout << dbMax << "," << nRow << "," << nCol << endl;
			// 阈值条件判断
			if (dbMax < precision)     //precision check
				break;
			if (nCount > max)       //iterations check
				break;
			nCount++;
			double dbApp = mat[nRow][nRow];
			double dbApq = mat[nRow][nCol];
			double dbAqq = mat[nCol][nCol];
			// 计算旋转矩阵
			double dbAngle = 0.5*atan2(-2 * dbApq, dbAqq - dbApp);
			double dbSinTheta = sin(dbAngle);
			double dbCosTheta = cos(dbAngle);
			double dbSin2Theta = sin(2 * dbAngle);
			double dbCos2Theta = cos(2 * dbAngle);
			mat[nRow][nRow] = dbApp*dbCosTheta*dbCosTheta + dbAqq*dbSinTheta*dbSinTheta + 2 * dbApq*dbCosTheta*dbSinTheta;
			mat[nCol][nCol] = dbApp*dbSinTheta*dbSinTheta + dbAqq*dbCosTheta*dbCosTheta - 2 * dbApq*dbCosTheta*dbSinTheta;
			mat[nRow][nCol] = 0.5*(dbAqq - dbApp)*dbSin2Theta + dbApq*dbCos2Theta;
			mat[nCol][nRow] = mat[nRow][nCol];
			for (int i = 0; i < dim; i++)
			{
				if ((i != nCol) && (i != nRow))
				{
					dbMax = mat[i][nRow];
					mat[i][nRow] = mat[i][nCol] * dbSinTheta + dbMax*dbCosTheta;
					mat[i][nCol] = mat[i][nCol] * dbCosTheta - dbMax*dbSinTheta;
				}
			}
			for (int j = 0; j < dim; j++) {
				if ((j != nCol) && (j != nRow)) {
					dbMax = mat[nRow][j];
					mat[nRow][j] = mat[nCol][j] * dbSinTheta + dbMax*dbCosTheta;
					mat[nCol][j] = mat[nCol][j] * dbCosTheta - dbMax*dbSinTheta;
				}
			}
			//compute eigenvector
			for (int i = 0; i < dim; i++)
			{
				dbMax = eigenvectors[i][nRow];
				eigenvectors[i][nRow] = eigenvectors[i][nCol] * dbSinTheta + dbMax*dbCosTheta;
				eigenvectors[i][nCol] = eigenvectors[i][nCol] * dbCosTheta - dbMax*dbSinTheta;
			}
		}
		// 特征值排序
		std::map<double, int> mapEigen;
		for (int i = 0; i < dim; i++)
		{
			eigenvalues[i] = mat[i][i];
			mapEigen.insert(make_pair(eigenvalues[i], i));
		}
		double *pdbTmpVec = new double[dim*dim];
		std::map<double, int>::reverse_iterator iter = mapEigen.rbegin();
		for (int j = 0; iter != mapEigen.rend(), j < dim; ++iter, ++j) {
			for (int i = 0; i < dim; i++) {
				pdbTmpVec[i*dim + j] = eigenvectors[i][iter->second];
			}
			eigenvalues[j] = iter->first;
		}
		for (int i = 0; i < dim; i++)
		{
			double dSumVec = 0;
			for (int j = 0; j < dim; j++)
				dSumVec += pdbTmpVec[j * dim + i];
			if (dSumVec < 0)
			{
				for (int j = 0; j < dim; j++)
					pdbTmpVec[j * dim + i] *= -1;
			}
		}
		for (int i = 0; i < dim; i++)
		{
			for (int j = 0; j < dim; j++)
			{
				eigenvectors[i][j] = pdbTmpVec[i * dim + j];
			}
		}
		delete[]pdbTmpVec;
		return true;
	}
};

int main(int argc, char *argv[])
{
	vector<vector<double>> matrix;
	matrix.push_back({ 161.079031576618, 208.250301732992, -64.1043482618227 });
	matrix.push_back({ 208.250301732992, 273.792923146327, -84.8950989429858 });
	matrix.push_back({ -64.1043482618228, -84.8950989429857, 26.9151793309588 });

	//vector<vector<double>> eigenvectors(3, vector<double>(3, 0));
	//vector<double> eigenvalues(4, 0);
	vector<vector<double>> eigenvectors;
	vector<double> eigenvalues;
	double eps = 1e-10;
	int T = 10000;

	bool ret = Algorithm::My_Jacobi(matrix, eigenvectors, eigenvalues, eps, T);

	qDebug() << "matrix";
	for (int i = 0; i < matrix.size(); ++i)
	{
		TemplateFun::printFun(matrix[i].data(), matrix[i].size());
	}

	qDebug() << "eigenvectors";
	for (int i = 0; i < eigenvectors.size(); ++i)
	{
		TemplateFun::printFun(eigenvectors[i].data(), eigenvectors[i].size());
	}

	qDebug() << "eigenvalues";
	TemplateFun::printFun(eigenvalues.data(), eigenvalues.size());

	return 0;
}