雅可比迭代法法
在图形图像中很多地方用到求矩阵的特征值和特征向量,比如主成分分析、OBB包围盒等。编程时一般都是用数值分析的方法来计算,这里介绍一下雅可比迭代法求解特征值和特征向量。雅可比迭代法的原理,网上资料很多,详细可见参考资料1。这里我们简单介绍求解矩阵S特征值和特征向量的步骤:
- 初始化特征向量为对角阵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示例对比。
复制代码
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159#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
最后
以上就是悲凉小蘑菇最近收集整理的关于雅可比(Jacobi)计算特征值和特征向量雅可比迭代法法代码实现参考资料的全部内容,更多相关雅可比(Jacobi)计算特征值和特征向量雅可比迭代法法代码实现参考资料内容请搜索靠谱客的其他文章。
本图文内容来源于网友提供,作为学习参考使用,或来自网络收集整理,版权属于原作者所有。
发表评论 取消回复