雅可比(Jacobi)计算特征值和特征向量

雅可比迭代法法

在图形图像中很多地方用到求矩阵的特征值和特征向量，比如主成分分析、OBB包围盒等。编程时一般都是用数值分析的方法来计算，这里介绍一下雅可比迭代法求解特征值和特征向量。雅可比迭代法的原理，网上资料很多，详细可见参考资料1。这里我们简单介绍求解矩阵S特征值和特征向量的步骤：

1. 初始化特征向量为对角阵V，即主对角线的元素都是1.其他元素为0。
2. 在S的非主对角线元素中，找到绝对值最大元素 Sij。
3. 用下 式计算tan2θ，求 cosθ、sinθ 及旋转矩阵Gij 。
   
4. 用下面公式求S‘；用当前特征向量矩阵V乘以矩阵Gij得到当前的特征向量V。
   
5. 若当前迭代前的矩阵A的非主对角线元素中最大值小于给定的阈值e时，停止计算；否则, 令S =S‘, 重复执行(2) ~ (5)。 停止计算时，得到特征值 li≈(S‘) ij ，i,j= 1,2,…,n.以及特征向量V。
6. 这一步可选。根据特征值的大小从大到小的顺序重新排列矩阵的特征值和特征向量。

代码实现

用C++实现，并与参考资料1示例对比。

#include <iostream>
#include <map>
#include<math.h>
#include <iomanip>
using namespace std;

/**
* @brief                        Jacobi eigenvalue algorithm
* @param matrix				    n*n array
* @param dim					dim represent n
* @param eigenvectors			n*n array
* @param eigenvalues			n*1 array
* @param precision   			precision requirements
* @param max					max number of iterations
* @return
*/
bool Jacobi(double* matrix, int dim, double* eigenvectors, double* eigenvalues, double precision, int max)
{
	for (int i = 0; i < dim; i++) {
		eigenvectors[i*dim + i] = 1.0f;
		for (int j = 0; j < dim; j++) {
			if (i != j)
				eigenvectors[i*dim + j] = 0.0f;
		}
	}

	int nCount = 0;		//current iteration
	while (1) {
		//find the largest element on the off-diagonal line of the matrix
		double dbMax = matrix[1];
		int nRow = 0;
		int nCol = 1;
		for (int i = 0; i < dim; i++) {			//row
			for (int j = 0; j < dim; j++) {		//column
				double d = fabs(matrix[i*dim + j]);
				if ((i != j) && (d > dbMax)) {
					dbMax = d;
					nRow = i;
					nCol = j;
				}
			}
		}

		if (dbMax < precision)     //precision check 
			break;
		if (nCount > max)       //iterations check
			break;
		nCount++;

		double dbApp = matrix[nRow*dim + nRow];
		double dbApq = matrix[nRow*dim + nCol];
		double dbAqq = matrix[nCol*dim + nCol];
		//compute rotate angle
		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);
		matrix[nRow*dim + nRow] = dbApp*dbCosTheta*dbCosTheta +
			dbAqq*dbSinTheta*dbSinTheta + 2 * dbApq*dbCosTheta*dbSinTheta;
		matrix[nCol*dim + nCol] = dbApp*dbSinTheta*dbSinTheta +
			dbAqq*dbCosTheta*dbCosTheta - 2 * dbApq*dbCosTheta*dbSinTheta;
		matrix[nRow*dim + nCol] = 0.5*(dbAqq - dbApp)*dbSin2Theta + dbApq*dbCos2Theta;
		matrix[nCol*dim + nRow] = matrix[nRow*dim + nCol];

		for (int i = 0; i < dim; i++) {
			if ((i != nCol) && (i != nRow)) {
				int u = i*dim + nRow;	//p  
				int w = i*dim + nCol;	//q
				dbMax = matrix[u];
				matrix[u] = matrix[w] * dbSinTheta + dbMax*dbCosTheta;
				matrix[w] = matrix[w] * dbCosTheta - dbMax*dbSinTheta;
			}
		}

		for (int j = 0; j < dim; j++) {
			if ((j != nCol) && (j != nRow)) {
				int u = nRow*dim + j;	//p
				int w = nCol*dim + j;	//q
				dbMax = matrix[u];
				matrix[u] = matrix[w] * dbSinTheta + dbMax*dbCosTheta;
				matrix[w] = matrix[w] * dbCosTheta - dbMax*dbSinTheta;
			}
		}

		//compute eigenvector
		for (int i = 0; i < dim; i++) {
			int u = i*dim + nRow;		//p   
			int w = i*dim + nCol;		//q
			dbMax = eigenvectors[u];
			eigenvectors[u] = eigenvectors[w] * dbSinTheta + dbMax*dbCosTheta;
			eigenvectors[w] = eigenvectors[w] * dbCosTheta - dbMax*dbSinTheta;
		}
	}

	//sort eigenvalues
	std::map<double, int> mapEigen;
	for (int i = 0; i < dim; i++) {
		eigenvalues[i] = matrix[i*dim + 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*dim + 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;
		}
	}
	memcpy(eigenvectors, pdbTmpVec, sizeof(double)*dim*dim);
	delete[]pdbTmpVec;
	return true;
}


int main()
{
	double a[4][4] = { {4, -30, 60, -35}, {-30, 300, -675, 420},{ 60, -675, 1620, -1050},{ -35, 420, -1050, 700}};
	int n = 4;
	double eps = 1e-10;
	int T = 10000;
	double p[4][4] = { 0 };
	double v[4] = { 0 };
	bool re = Jacobi(&a[0][0], n, &p[0][0], v, eps, T);
	if (re)	{
		cout << "Matrix:" << endl;
		for (int i = 0; i < n; i++)	{
			for (int j = 0; j < n; j++)
				cout << setw(15) << a[i][j];
			cout << endl;
		}
		cout << "eigenvectors:" << endl;
		for (int i = 0; i < n; i++)	{
			for (int j = 0; j < n; j++)
				cout << setw(15) << p[i][j];
			cout << endl;
		}
		cout << "eigenvalues:" << endl;
		for (int i = 0; i < n; i++)	{
				cout << setw(15) << v[i];
			cout << endl;
		}
	}
	else
		cout << "false" << endl;
	system("pause");
}

运行结果：

参考资料1示例结果：

参考资料

【Jacobi eigenvalue algorithm】https://en.wikipedia.org/wiki/Jacobi_eigenvalue_algorithm