固有値/固有ベクトル の変更点

Top/固有値/固有ベクトル

[ トップ | 新規 | 一覧 | 検索 | 最終更新 | ヘルプ ] [ 編集 | 差分 | 履歴 | 添付 | リロード ]  

• 追加された行はこの色です。
• 削除された行はこの色です。
• 固有値/固有ベクトル へ行く。
• 固有値/固有ベクトル の差分を削除

---

----
#contents
----
*ヤコビ(Jacobi)法 [#he7a28c3]
対称行列A(n×n)のn個の固有値・固有ベクトルを求める方法．

#code(C){{
/*!
 * 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;
}
}}

*ハウスホルダー変換 [#e8b8a85a]
実対称行列を三重対角行列に変換する．さらに三重対角行列から固有値/固有ベクトルを求めるまでをハウスホルダー法と呼ぶ．
ハウスホルダー変換はQR法におけるヘッセンベルグ行列への変換にも用いられる．

ハウスホルダー変換では固有値を求めたい行列Aから対称な直交行列Pを定義し，
1行1列目が三重対角行列になった行列PAPを計算，これを繰り返す．
Pが対称な直交行列(P=P^-1=P^T)なのでAとPAPの固有値は等しい．

#code(C){{
/*!
 * ハウスホルダー変換で実対称行列を三重対角行列に変換
 * @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 = k+1; j < n; ++j){
				a[i*n+j] = a[i*n+j]-2.0*u[i]*v[j]-v[i]*u[j];
			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法 [#h5d379a2]
正方行列A(n×n)の固有値と固有ベクトルを求める．Aは複素行列でも可．

QR法ではまず，非対称行列をハウスホルダー変換によりヘッセンベルグ行列に変換する
(対称行列ではないので三重対角行列ではないことに注意)．
ハウスホルダー変換によるヘッセンベルグ行列への変換は以下．

#code(C){{
/*!
 * ハウスホルダー変換で行列をヘッセンベルグ行列に変換
 * @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に渡す．

#code(C){{
/*!
 * 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
の固有値を求める．

#code(C){{
	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で上記結果と一致する．

       

Site admin: fujis

PukiWiki 1.5.4 © 2001-2022 PukiWiki Development Team. Powered by PHP 8.1.29. HTML convert time: 0.002 sec.