雅可比迭代法法

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

初始化特征向量为对角阵V，即主对角线的元素都是1.其余元素为0。

在S的非主对角线元素中，找到绝对值最大元素 Sij。

用下 式计算tan2θ，求 cosθ、sinθ 及旋转矩阵Gij 。

用下面公式求S‘；用当前特征向量矩阵V乘以矩阵Gij获得当前的特征向量V。

若当前迭代前的矩阵A的非主对角线元素中最大值小于给定的阈值e时，中止计算；不然, 令S =S‘, 重复执行(2) ~ (5)。 中止计算时，获得特征值 li≈(S‘) ij ，i,j= 1,2,…,n.以及特征向量V。

这一步可选。根据特征值的大小从大到小的顺序从新排列矩阵的特征值和特征向量。

代码实现

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

#include

#include

#include

#include

using namespace std;

/**

* @brief Jacobi eigenvalue algorithm

* @param matrix n*n array

* @param dimdim represent n

* @param eigenvectorsn*n array

* @param eigenvaluesn*1 array

* @param precision precision requirements

* @param maxmax 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 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::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示例结果：

参考资料