2003/10/14

2 週間、ヤコビ法を学んだ。今週は、べき乗法

ヤコビ法は、基本的に、対称行列を扱う
	# 保存則を満す、物理法則を普通に離散化すれば、対称行列になる。

非対称行列を扱う方法
	- 左下の要素を 0 にして QR 法 ( QR 分解を繰り返す ) を適用する
	- べき乗法

# べき乗法は、非対称の場合によく用いられるが、対称の場合でも適用可能。
# 以下の説明では、対称の場合で説明するが、非対称の場合も本質的に同じ

==

べき乗法

	A : n x n の対称行列
	A x = \lambda x を満す解 ( 固有値 ) を \lambda_i とし、それに
	対応した、固有ベクトルを x_i する。

	この時、
		A = \Sum_{i=1}^n \lambda x_i x_i^t
	が成立する。

	固有ベクトルの直交性より、	
		AA = \Sum_{i=1}^n \lambda_i^2 x_i x_i^t
	がいえる。したがって、
		A^k = \Sum_{i=1}^n \lambda_i^k x_i x_i^t

	今、
		|\lambda_1| > |\lambda_2| > .. > |\lambda_n|
	 となるようにしたとし、更に、\lambda_1 で括れば、
		A^k = \lambda_1^k\Sum_{i=1}^n \frac{\lambda_i}{\lambda_1}^k x_i x_i^t
	ここで、
		\lim_{k\rightarrow \infty} \frac{\lambda_i}{\lambda_1}^k = 0
	なので、\lambda_1 と x_1 以外の要素は、無視できるようになる。
	つまり、
		A^k \eq \lambda_1^k x_1 x_1^t
	となる。
	そこで、
		| A^{k+1) | / |A| \eq | \lambda_1 |
	より、最大絶対値の固有値の絶対値

	[収束性]
		\lambda_1 に対して、\lambda_2 の絶対値がどの位の比かどうか
		# 根が近接している場合は、収束が遅い

べき乗法の性質
	固有ベクトルの直交性が本質
	絶対値の最大の固有値だけがわかる
		# 応用上では、結構それだけで十分のことがある
		#	建物の振動では、最低周波を意味するので、
		#	地震対策ではこれだけで十分

		二つ目以後はどうするか ?
		\lambda_1 は求まったので、これを利用して固有ベクトル x_1
		が計算できるので、

			A' = A - \lambda_1 x_1 x_1^t

	 	とすれば、この A' は、ほぼ A と同じで、\lambda_1
		がなく 0 に近い固有値 ( \lambda_1 の誤差の大きさ ) 
		をもつ。

		これに、べき乗法を再度適用すれば、\lambda_2 が得られる。

		これを繰り返せば、理論的には、全ての固有値がえられるが、
		実際には、誤差が累積するために、後の方の固有値は、信頼
		できない。

		振動の応用などでは、最大値の大きい方の数個が欲しいので、
		( ヤコビ法で全部求めてから不必要なものを捨てるより..)
		こちらの方が計算量も少い。

		# Q:絶対値の一番小さい固有値を求めるには ?
		# A: A^{-1} に対してべき乗法を適用すれば Okey (逆反復法)
		# prof)
		#	A^{-1} の固有値は \lambda_i^{-1} なので..

	行列の積を計算するだけなので簡単

		並列化もしやすい

			Ax = b を直接法で解く .. 並列化しにくい
			Ax, AA を計算する .. 並列化しやすい

			# 今、「並列化」は流行..

		行列の積も、密行列だと大変 ( n^3 ) だが、粗行列なら、
		それほど苦にならない

			応用 ( 偏微分方程式の差分技術 )
				有限差分法 ..粗
				有限要素法 ..粗
				境界要素法 .. 密 ( ただし、次元が低い )

			# 差分では、近接領域の影響のみをかんがえるので、
			# 離れた処からの影響は小さい ( 無視して 0 とみなす.. )
			#	=> 粗行列へ..

		# ランチャフ法 ( べき乗法の拡張で、複数の固有値を同時に求める )

==

数値計算の世界では、時代が変って計算機の種類が代ると、昔の手法が
見直されるということがよくある。

	# トンカチの形が変れば、良い使い方も変る。

==

	非対称の場合には、
		A = \Sum_{i=1}^n \lambda x_i x_i^t
	の変りに
		A = \Sum_{i=1}^n \lambda x_i x_i'
			where	x_i' は A^t の固有値
	を利用する点が異なるので、本質的には全く同じ方法を利用する。

==

べき乗法のプログラム text p.69

	|v^(0)| = 1 となる v^(0) を初期値として取る
		# (1,0,..,0) でもよいが本当は、
                #  (\frac{1}{\sqrt{n}},\frac{1}{\sqrt{n}},..,\frac{1}{\sqrt{n}}がよい。	
		# 結局 v^{(k)} が固有ベクトルに収束するが、その成分が、
		# 初期値の v^{(0)} から出てくる時間に差がある。
		#	# 0 の要素はなかなか 0 以外にならない..(?)

	基本は、次の繰り返しを行い、v^{(k)} の収束すればそれが解
		# \lim | v^{(k+1)} - v^{(k)} | = 0

		u^{(k+1)} = A v^{(k)}
		\mu^{(k+1)} = v^{(k)}^tu^{(k+1)}
		v^{(k)} = \frac{u^{(k+1)}}{|u^{(k+1)}|}

		# \mu の収束をみるという方法もある。
		# 重複する固有値があれば、固有ベクトルが収束しない可能
                # 性があり、固有ベクトルで判定するより、こちらが望ましい。

		## 問題そのものに対称性がある場合には、固有値が重複する
		## 	=> 問題を変形し、対称性を無くして計算する。
		## cf. 正方形の真中を押す
		##		-> 四隅で同じ
		##	問題の領域を四分の一にして、対称性を無くしてから差分化

		### テスト問題 ( 教室.. ) では、これにひっかかる
		### 現実世界では、対称性があるような問題はほとんどない
		###	-> 固有値の重複で悩むことはない..


時間がある人は、次の \lambda_2 ( 二つ目 ) も求める
	# ヒント : 「A' = A - \lambda_1 x_1 x_1^t」の計算を倍精度で行えばよい。

==

数値計算では、答がわかると、どうやって解くと良いか ( 最適解法 ) が解る。
	=> これは矛盾 !! ( 答えが解らない状態で、最適な方法が知りたい !!)

	# 実際は、完全な答でないくても、最適法が推測できるので、通常は、
	# ラフな方法で数桁を求め、それを用いて、「最適な方法を探す」よ
	# うな方法がとられる。

==

[課題]

	べき乗法のプログラム text p.69 をやる。
	例題は、p.70 の行列を解く