![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
ヤコビ(Jacobi)法†対称行列A(n×n)のn個の固有値・固有ベクトルを求める方法． /*!
* Jacobi法による固有値の算出
* @param[inout] a 実対称行列．計算後，対角要素に固有値が入る
* @param[out] v 固有ベクトル(aと同じサイズ)
* @param[in] n 行列のサイズ(n×n)
* @param[in] eps 収束誤差
* @param[in] iter_max 最大反復回数
* @return 反復回数
*/
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); // sinθ
st = (((aii-ajj) >= 0.0) ? 1.0 : -1.0)*aij/(2.0*beta*ct); // cosθ
// A = PAPの計算
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;
// V = PVの計算
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の固有値は等しい． /*!
* ハウスホルダー変換で実対称行列を三重対角行列に変換
* @param[inout] a 元の行列(n×n)．変換された三重対角行列を格納
* @param[in] n 行列のサイズ
*/
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){
// sの計算
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);
// |x-y|の計算
float alpha = sqrt(2.0*s*(s-a[(k+1)*n+k]));
if(fabs(alpha) < 1e-8) continue;
// uの計算
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;
}
// Auの計算
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];
}
}
// v=2(Au-uu^T(Au))の計算
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 = PAP = A-uv^T-vu^Tの計算
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法ではまず，非対称行列をハウスホルダー変換によりヘッセンベルグ行列に変換する (対称行列ではないので三重対角行列ではないことに注意)． ハウスホルダー変換によるヘッセンベルグ行列への変換は以下． /*!
* ハウスホルダー変換で行列をヘッセンベルグ行列に変換
* @param[inout] a 元の行列(n×n)．変換されたヘッセンベルグ行列を格納
* @param[in] n 行列のサイズ
*/
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){
// sの計算
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);
// |x-y|の計算
float alpha = sqrt(2.0*s*(s-a[(k+1)*n+k]));
if(fabs(alpha) < 1e-8) continue;
// uの計算
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;
}
// Pの計算
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];
}
}
}
// PAの計算
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];
}
}
}
// A = PAP^Tの計算
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に渡す． /*!
* QR法による固有値の算出
* @param[inout] h 元の行列(ヘッセンベルグ行列)．計算後，対角要素に固有値が入る
* @param[in] n 行列のサイズ(n×n)
* @param[in] eps 収束誤差
* @param[in] iter_max 最大反復回数
* @return 反復回数
*/
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];
// R = H : ヘッセンベルグ行列
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(;;){
// Q=I (単位行列)
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){
// sinθ,cosθの計算
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;
// Rの計算
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];
}
// Qの計算
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];
}
}
// RQの計算
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 の固有値を求める． 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で上記結果と一致する． |
