Jacobi 方法计算特征值特征向量

在求解协方差矩阵的特征值和特征向量时用到了Jacobi方法

原理

设A是n阶实对称矩阵，则必有正交矩阵P，使

PTAP=⎡⎣⎢⎢λ1⋮0⋯⋱⋯0⋮λn⎤⎦⎥⎥=∧;
P^TAP=
\begin{bmatrix}
\lambda_1 &\cdots &0 \\
\vdots &\ddots &\vdots\\
0 &\cdots &\lambda_n
\end{bmatrix}
=\wedge;

其中∧\wedge的对角线元素是A的n个特征值，正交阵PP的第ii列是A的对应于特征值λi\lambda_i的特征向量。

所以对于一个任意的n阶实对称矩阵A，只要能求得一个正交阵P，使PTAP=∧P^TAP=\wedge，则可得到A的全部特征值和特征向量。

具体

我们选择旋转矩阵

P=[cosθsinθ−sinθcosθ]
P=
\begin{bmatrix}
\cos\theta&-\sin\theta \\
\sin\theta&\cos\theta
\end{bmatrix}

P是正交矩阵，之前看的资料介绍了P是坐标的逆时针旋转，这样我就不明白PTAPP^TAP的含义了，还不如不把它理解为旋转。。。。。

PTAPP^TAP只改变A的第i，j行和i，j列，可以计算出改变的值，程序实现时只要改变这些值就可以了。

P=P1P2..PnP=P_1P_2..P_n就是特征向量

程序

//雅克比方法求实对称矩阵的特征值和特征向量

#include <iostream>
#include<math.h>
#include <iomanip>
#define PI 3.1415926
using namespace std;

/**
/* Input：
a:待求解矩阵
p:初始化为0
n:矩阵维数
eps:迭代停止条件
T:最大迭代次数
Output:
a:对角线为特征值
p：a中的每个特征值对应的特征向量，每一列为一个向量
*/
bool Jacobi(double a[][3], double p[][3], int n,double eps,int T)
{
int i, j, k;
double max = a[0][1];
int row = 0;
int col = 0;
double Ptemp[3][3];
int ite_num = 0;
for (i = 0; i < n; i++)
p[i][i] = 1;
while (1)
{
max = fabs(a[0][1]);
row = 0;
col = 1;
for (i = 0; i < n; i++)
for (j = 0; j < n; j++)
if (i != j && fabs(a[i][j])>max)
{
max = fabs(a[i][j]);
row = i;
col = j;
}
if (max < eps)
{
cout << "max:" << max << endl;
cout << "t:" << ite_num << endl;
return true;
}

if (ite_num>T)
{
cout << "max:" << max << endl;
cout << "t:" << ite_num << endl;
return false;
}

double theta;
if (a[row][row] == a[col][col])
theta = PI / 4;
else
theta = 0.5*atan(2 * a[row][col] / (a[row][row] - a[col][col]));
double aii = a[row][row];
double ajj = a[col][col];
double aij = a[row][col];
double sin_theta = sin(theta);
double cos_theta = cos(theta);
double sin_2theta = sin(2 * theta);
double cos_2theta = cos(2 * theta);
a[row][row] = aii*cos_theta*cos_theta + ajj*sin_theta*sin_theta + aij*sin_2theta;
a[col][col] = aii*sin_theta*sin_theta + ajj*cos_theta*cos_theta - aij*sin_2theta;
a[row][col] = 0.5*(ajj - aii)*sin_2theta + aij*cos_2theta;
a[col][row] = a[row][col];
for (k = 0; k < n; k++)
{
if (k != row && k != col)
{
double arowk = a[row][k];
double acolk = a[col][k];
a[row][k] = arowk * cos_theta + acolk * sin_theta;
a[k][row] = a[row][k];
a[col][k] = acolk * cos_theta - arowk * sin_theta;
a[k][col] = a[col][k];
}
}
if (ite_num == 0)
{
p[row][row] = cos_theta;
p[row][col] = -sin_theta;
p[col][row] = sin_theta;
p[col][col] = cos_theta;
}
else
{
for (k = 0; k < n; k++)
{
double pki = p[k][row];
double pkj = p[k][col];
p[k][row] = pki*cos_theta + pkj*sin_theta;
p[k][col] = pkj*cos_theta - pki*sin_theta;
}
}
ite_num++;
}
}

void main()
{
double a[3][3] = { 2, -1, 0, -1, 2, -1, 0, -1, 2 };
int n = 3;
double eps = 1e-10;
int T = 10000;
double p[3][3] = { 0 };
bool re = Jacobi(a, p, n, eps, T);
if (re)
{
cout << "a:" << endl;
for (int i = 0; i < n; i++)
{
for (int j = 0; j < n; j++)
cout << setw(15) << a[i][j];
cout << endl;
}
cout << "p:" << endl;
for (int i = 0; i < n; i++)
{
for (int j = 0; j < n; j++)
cout << setw(15) << p[i][j];
cout << endl;
}

}

else
cout << "false" << endl;
}