ヤコビ(Jacobi)法†
対称行列A(n×n)のn個の固有値・固有ベクトルを求める方法．
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
| |
int EigenJacobiMethod(float *a, float *v, int n, float eps = 1e-8, int iter_max = 100)
{
float *bim, *bjm;
float bii, bij, bjj, bji;
bim = new float[n];
bjm = new float[n];
for(int i = 0; i < n; ++i){
for(int j = 0; j < n; ++j){
v[i*n+j] = (i == j) ? 1.0 : 0.0;
}
}
int cnt = 0;
for(;;){
int i, j;
float x = 0.0;
for(int ia = 0; ia < n; ++ia){
for(int ja = 0; ja < n; ++ja){
int idx = ia*n+ja;
if(ia != ja && fabs(a[idx]) > x){
i = ia;
j = ja;
x = fabs(a[idx]);
}
}
}
float aii = a[i*n+i];
float ajj = a[j*n+j];
float aij = a[i*n+j];
float alpha, beta;
alpha = (aii-ajj)/2.0;
beta = sqrt(alpha*alpha+aij*aij);
float st, ct;
ct = sqrt((1.0+fabs(alpha)/beta)/2.0); st = (((aii-ajj) >= 0.0) ? 1.0 : -1.0)*aij/(2.0*beta*ct);
for(int m = 0; m < n; ++m){
if(m == i || m == j) continue;
float aim = a[i*n+m];
float ajm = a[j*n+m];
bim[m] = aim*ct+ajm*st;
bjm[m] = -aim*st+ajm*ct;
}
bii = aii*ct*ct+2.0*aij*ct*st+ajj*st*st;
bij = 0.0;
bjj = aii*st*st-2.0*aij*ct*st+ajj*ct*ct;
bji = 0.0;
for(int m = 0; m < n; ++m){
a[i*n+m] = a[m*n+i] = bim[m];
a[j*n+m] = a[m*n+j] = bjm[m];
}
a[i*n+i] = bii;
a[i*n+j] = bij;
a[j*n+j] = bjj;
a[j*n+i] = bji;
for(int m = 0; m < n; ++m){
float vmi = v[m*n+i];
float vmj = v[m*n+j];
bim[m] = vmi*ct+vmj*st;
bjm[m] = -vmi*st+vmj*ct;
}
for(int m = 0; m < n; ++m){
v[m*n+i] = bim[m];
v[m*n+j] = bjm[m];
}
float e = 0.0;
for(int ja = 0; ja < n; ++ja){
for(int ia = 0; ia < n; ++ia){
if(ia != ja){
e += fabs(a[ja*n+ia]);
}
}
}
if(e < eps) break;
cnt++;
if(cnt > iter_max) break;
}
delete [] bim;
delete [] bjm;
return cnt;
}
|
ハウスホルダー変換†
実対称行列を三重対角行列に変換する．さらに三重対角行列から固有値/固有ベクトルを求めるまでをハウスホルダー法と呼ぶ．
ハウスホルダー変換はQR法におけるヘッセンベルグ行列への変換にも用いられる．
ハウスホルダー変換では固有値を求めたい行列Aから対称な直交行列Pを定義し，
1行1列目が三重対角行列になった行列PAPを計算，これを繰り返す．
Pが対称な直交行列(P=P^-1=P^T)なのでAとPAPの固有値は等しい．
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
| |
void HouseholderTransformation(float *a, int n)
{
float *u, *v, *q;
u = new float[n];
v = new float[n];
q = new float[n];
for(int j = 0; j < n; ++j){
u[j] = 0.0;
v[j] = 0.0;
q[j] = 0.0;
}
for(int k = 0; k < n-2; ++k){
float s = 0;
for(int i = k+1; i < n; ++i){
s += a[i*n+k]*a[i*n+k];
}
s = -((a[(k+1)*n+k] >= 0) ? 1 : -1)*sqrt(s);
float alpha = sqrt(2.0*s*(s-a[(k+1)*n+k]));
if(fabs(alpha) < 1e-8) continue;
u[k+1] = (a[(k+1)*n+k]-s)/alpha;
for(int i = k+2; i < n; ++i){
u[i] = a[i*n+k]/alpha;
}
q[k] = alpha/2.0;
for(int i = k+1; i < n; ++i){
q[i] = 0.0;
for(int j = k+1; j < n; ++j){
q[i] += a[i*n+j]*u[j];
}
}
alpha = 0.0;
for(int i = k+1; i < n; ++i){
alpha += u[i]*q[i];
}
v[k] = 2.0*q[k];
for(int i = k+1; i < n; ++i){
v[i] = 2.0*(q[i]-alpha*u[i]);
}
a[k*n+(k+1)] = a[(k+1)*n+k] = s;
for(int i = k+2; i < n; ++i){
a[k*n+i] = 0.0;
a[i*n+k] = 0.0;
}
for(int i = k+1; i < n; ++i){
a[i*n+i] = a[i*n+i]-2*u[i]*v[i];
for(int j = i+1; j < n; ++j){
a[i*n+j] = a[i*n+j]-u[i]*v[j]-v[i]*u[j];
a[j*n+i] = a[i*n+j];
}
}
}
delete [] u;
delete [] v;
delete [] q;
}
|
QR法†
正方行列A(n×n)の固有値と固有ベクトルを求める．Aは複素行列でも可．
QR法ではまず，非対称行列をハウスホルダー変換によりヘッセンベルグ行列に変換する
(対称行列ではないので三重対角行列ではないことに注意)．
ハウスホルダー変換によるヘッセンベルグ行列への変換は以下．
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
| |
void HouseholderTransformationForQR(float *a, int n)
{
float *b, *p, *q;
float *u;
b = new float[n*n];
p = new float[n*n];
q = new float[n*n];
u = new float[n];
for(int k = 0; k < n-2; ++k){
float s = 0;
for(int i = k+1; i < n; ++i){
s += a[i*n+k]*a[i*n+k];
}
s = -((a[(k+1)*n+k] >= 0) ? 1 : -1)*sqrt(s);
float alpha = sqrt(2.0*s*(s-a[(k+1)*n+k]));
if(fabs(alpha) < 1e-8) continue;
u[k+1] = (a[(k+1)*n+k]-s)/alpha;
for(int i = k+2; i < n; ++i){
u[i] = a[i*n+k]/alpha;
}
for(int i = k+1; i < n; ++i){
for(int j = i; j < n; ++j){
if(j == i){
p[i*n+j] = 1.0-2.0*u[i]*u[i];
}
else{
p[i*n+j] = -2.0*u[i]*u[j];
p[j*n+i] = p[i*n+j];
}
}
}
for(int i = k+1; i < n; ++i){
for(int j = k+1; j < n; ++j){
q[i*n+j] = 0.0;
for(int m = k+1; m < n; ++m){
q[i*n+j] += p[i*n+m]*a[m*n+j];
}
}
}
for(int i = 0; i <= k; ++i){
b[i*n+k] = a[i*n+k];
}
b[(k+1)*n+k] = s;
for(int i = k+2; i < n; ++i){
b[i*n+k] = 0.0f;
}
for(int j = k+1; j < n; ++j){
for(int i = 0; i <= k; ++i){
b[i*n+j] = 0.0;
for(int m = k+1; m < n; ++m){
b[i*n+j] += a[i*n+m]*p[j*n+m];
}
}
for(int i = k+1; i < n; ++i){
b[i*n+j] = 0.0;
for(int m = k+1; m < n; ++m){
b[i*n+j] += q[i*n+m]*p[j*n+m];
}
}
}
for(int j = 0; j < n*n; ++j){
a[j] = b[j];
}
}
delete [] u;
delete [] b;
delete [] p;
delete [] q;
}
|
行列Aをヘッセンベルグ行列Hへ変換後，QR法による固有値を算出する．
ハウスホルダー変換後の行列を引数hに渡す．
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
| |
void QR(float *h, int n, float eps = 1e-8, int iter_max = 100)
{
float *r, *q, *t;
r = new float[n*n];
q = new float[n*n];
t = new float[n*n];
float *u, *v;
u = new float[n];
v = new float[n];
for(int i = 0; i < n; ++i){
for(int j = 0; j < n; ++j){
int idx = i*n+j;
r[idx] = h[idx];
}
}
int cnt = 0;
for(;;){
for(int i = 0; i < n; ++i){
for(int j = 0; j < n; ++j){
q[i*n+j] = (i == j) ? 1.0 : 0.0;
}
}
for(int k = 0; k < n-1; ++k){
float alpha = sqrt(r[k*n+k]*r[k*n+k]+r[(k+1)*n+k]*r[(k+1)*n+k]);
if(fabs(alpha) < 1e-8) continue;
float c = r[k*n+k]/alpha;
float s = -r[(k+1)*n+k]/alpha;
for(int j = k+1; j < n; ++j){
u[j] = c*r[k*n+j]-s*r[(k+1)*n+j];
v[j] = s*r[k*n+j]+c*r[(k+1)*n+j];
}
r[k*n+k] = alpha;
r[(k+1)*n+k] = 0.0;
for(int j = k+1; j < n; ++j){
r[k*n+j] = u[j];
r[(k+1)*n+j] = v[j];
}
for(int j = 0; j <= k; ++j){
u[j] = c*q[k*n+j];
v[j] = s*q[k*n+j];
}
q[k*n+(k+1)] = -s;
q[(k+1)*n+(k+1)] = c;
for(int j = 0; j <= k; ++j){
q[k*n+j] = u[j];
q[(k+1)*n+j] = v[j];
}
}
for(int i = 0; i < n; ++i){
for(int j = 0; j < n; ++j){
float rq = 0.0;
for(int m = 0; m < n; ++m){
rq += r[i*n+m]*q[j*n+m];
}
t[i*n+j] = rq;
}
}
float e = 0.0;
for(int i = 0; i < n; ++i){
e += fabs(t[i*n+i]-h[i*n+i]);
}
if(e < eps) break;
if(cnt > iter_max) break;
for(int i = 0; i < n; ++i){
for(int j = 0; j < n; ++j){
int idx = i*n+j;
r[idx] = t[idx];
h[idx] = t[idx];
}
}
cnt++;
}
delete [] r;
delete [] q;
delete [] t;
delete [] u;
delete [] v;
}
|
QR関数を使用した例を示す．
行列
3 0 0
-2 -2 4
0 -1 3
の固有値を求める．
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
| | int n = 3;
float *h = new float[n*n];
float h0[3][3] = { {3, 0, 0},
{-2, -2, 4},
{0, -1, 3} };
for(int i = 0; i < n; ++i){
for(int j = 0; j < n; ++j){
h[i*n+j] = h0[i][j];
}
}
HouseholderTransformationForQR(h, n);
QR(h, n, 1e-8, 100);
cout << "QR" << endl;
for(int i = 0; i < n; ++i){
for(int j = 0; j < n; ++j){
cout << h[i*n+j] << " ";
}
cout << endl;
}
cout << "eigen values : " << h[0] << ", " << h[4] << ", " << h[8] << endl;
delete [] h;
|
結果は，
|3 -0.333333 -0.29814|
QR : |3.66271e-007 2 -5.36656|
|0 6.50788e-011 -1|
eigen values : 3, 2, -1
数学的の求めてみる．特性方程式は，
(λ-3)((λ+2)(λ-3)+4) = 0
変形すると
(λ-3)(λ-2)(λ+1) = 0
となり，λ=3,2,-1で上記結果と一致する．
![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
