2004/06/22

先週 : Gauss の消去法
今週 : LU 分解法

LU 法 : 教科書には、Gauss 法と異るとかいてあるが、本質的には、同じ物

LU 分解法
	L : 左下三角行列
	U : 右上三角行列

	行列 A が与えられたときに、
	     A = L U
	となる L と U を求め、これを用いて、連立方程式を解く

前回の話
	Ax = b
を解くと
	\frac{1}{3}n^3 + \frac{1}{2}n^2
だけの計算量が必要
	この内、前の項目 ( n^3 ) が重い
		これは、「前進消去」の部分
		後の部分は、「後退代入」の計算量

# ここまでを復習しておいて...

応用として、複数の連立を解きたい場合がある
	Ax_1 = b_1
	Ax_2 = b_2
	..
	Ax_m = b_m

# ポイントは、係数行列 A が同じで、定数項の b_i だけが変化するという点

もしこれを、これまでと同じ方法で解くと、計算量は、

	m (\frac{1}{3}n^3 + \frac{1}{2}n^2)

しかし、
	前進消去	: 定数項の内容は無関係
	後退代入	: 定数項を利用している
		=> むだな計算 (毎回同じ計算) をしていることになる
			効率よく計算するにはどうすればよいか ?

効率化の発想
	# まず、単純なものを考えると..
		ax = b
	を解くには..
		x = a^{-1} b
	一般の行列では..
		Ax = b
	より、
		x = A^{-1} b	Where A^{-1} は、A の逆行列
	とすればよい
		x_i を計算する場合も予め A^{-1} を計算しておけば、A^{-1} b_i で Okey
			一度計算した A^{-1} を何度も利用できる点がポイント

# 発想としては、正しい、しかしながら.. 数値計算的には、やってはいけないこと。
## 本当に A^{-1} の要素が欲しい場合は別。しかし、連立方程式を解くために使ってはだめ

## Why ?
##	逆行列の計算に \frac{2}{3}n^3 だけかかる
##	精度が良くない ( せいぜい、同じか、悪くなるだけ.. )
##		# LU ならば、計算量が少い ( \frac{1}{3}n^3 ) し精度も良い。

## 定数項 (b_1, .., b_m) があらかじめ判っている場合は、逆行列を求める
## 方法と同じやりかたでやっても無駄がない。

LU 分解法が利用される場合

	始めに、
		Ax_1 = b_1
	を解いた後に、ここで、初めて b_2 = x_1 が解り、これを利用して、
		Ax_2 = b_2
	を解くような場合。

==

LU 分解の論理的な構造

	A = L U


               a_11      a_1n

	A = (	               )

               a_n1      a_nn

                  1      
               l_21 1      0
	L = (	      ..        )
                           1
               l_n1          1


               u_11 u_12    u_1n
                    u_22
	U = (	        ..       )
                 0        ..
                            u_nn


	# LU の分解を一意にするために、 L の対角要素は 1 にする。
	# もし、U の対角を 1 にすれば、L の対角要素は 1 にならない。
	# LU ではなく LDU ( D は対角要素だけ ) とし、L も U も対角が 1 になる

と分解できたとする。

左辺の要素を計算すると..

	a_ij = l_i_1 u_1_j + l_i_2 u_2_j + .. + l_i_{i-1} u_{i-1}_i + u_i_j
		(i<=j)

	a_ij = l_i_1 u_1_j + l_i_2 u_2_j + .. + l_i_{i-1} u_{i-1}_i l_i_j u_j_j
		(i<=j)

となる。

# 実は、LU は、「Gauss の消去法の手順」を記録している行列

	今
		A x = b
	とすると、
		LU x = b		( where A = LU )
	なので、
		U x = L^{-1} b
	ここで、
		U x = y
	と置けば
		y = L^{-1} b
	であることに注意

	最後の
		y = L^{-1} b
	は、
		L y = b
	ですが、L は、左下行列なので、前進消去で解くことができる
		オーダーは、\frac{1}{2}{n^2}
	一方、( 上記の方法で、y が求まったと考えれば.. )
		U x = y
	は、Gauss の後退代入そのもの
		# したがって、b の値に無関係であることに注意 !!

	では、L は何を表しているのか ???

	Gauss の前進消去をもう一度考えてみる


               a_11      a_1n        x_1        b_1

	    (	               ) (        ) = (     )

               a_n1      a_nn        x_n        b_n

	で、a_21 を消去するために、両辺に、基本行列 B_1 = B(1,2,-a21/a11) を掛けている
	同様にして、消去を繰り返すと..

		A x = b
		B_1 A x = B_1 b
		B_2 B_1 A x = B_2 B_1 b
			..
		B_m .. B_2 B_1 A x = B_m .. B_2 B_1 b
		(B_m .. B_2 B_1) A x = (B_m .. B_2 B_1) b
		U x =  (B_m .. B_2 B_1) b

	ここで、
		U = (L^{-1}L)U = L^{-1}(LU) = L^{-1}A 
	なので、結局、
		L^{-1} =  (B_m .. B_2 B_1)
	更に、L の要素をみると、消去を行うためには、
		-\frac{a_21}{a_11}
	を計算するが、それが l_21 の値となっていることがわかる。
	同様にして、a_ij を消去するための因子
		-\frac{a_ij}{a_ii}
	が l_ij となる。

	# L は消去法で、要素を消去するための過程をそこで利用される因子
        # を記憶することによって、手順を記録している行列となっている。

## 領域を節約するために..
##	U は左上行列なので、そのまま収める、右下には 0 が入る ( 0 の所は記憶しない )
##	L は右下行列で、対角要素は 1 ( 1 の所は記憶しない ) 右上も 0 
##	よって、n x n 配列で、二つの行列 L U の両方を記憶することができる。

	Gauss の消去法の前進消去の計算と L U 分解の計算は、ほぼ同じ
		Gauss の消去法では、U を求めるだけ、L の要素 -\frac{a_ij}{a_ii} は
		捨てている
		LU 分解では、L の要素もきちんと保存しておけばよい
			# U の行列で 0 になる ( する ) 所に l_ij =  -\frac{a_ij}{a_ii} を入れる。

# LU 分解を理解する上で、「実は Gauss の変形」というのは、解りやすい。
## このようなことがかいてある教科書はみたことがない。


==

LU 分解のプログラム例 ( Fortran による )

!
! LU 分解
!
!	A = LU
!
!		A = a(n,n)
!	を入力し
!		L = a(n,n)	( i < j )
!		U = a(n,n)	( i >= j )
!	とする
!
	do k = 1, n-1
		do i = k+1, n
			do j = k+1, n
 				a(i,j)=a(i,j)-a(i,j)*a(k,j)/a(k.k)
			end do	
		end do	
	end do	

! Gauss では a(i,j) ( i<j ) は、0 にするが、LU では、l_ij を入れるために使う
! 本来は l_ij = - a(i,j)/a(i,i) だが、a(i,i) が残っているので、l_ij は、a(i,j) のままで良い。
! その変りに、後で、a(i,j) と a(i,i) から l_ij を計算する

!
! U y = b の y を計算する ( 配列 x(n) に入れる ) 
!
	do k = 2, n
		s = 0
		do j = 1, k-1
			s = s + a(k,j) * x(j)
		end do
		x(k) = x(k) - s
	end do

	! この形は、少し、時間が余分にかかるが、精度は良い。

!
! L x = y の x を計算する ( y は 配列 x(n) に入っており、x も x(n) に入れる )
!

	x(n) = x(n) / a(nn)

	do k = n-1, 1, -1
		s = 0
		do j = k+1, n
			s = s + a(k,j)*x(j)
		end do
		x(k) = (x(k)-s) / a(k,k)
	end do

# このアルゴリズムは、福井先生オリジナル !!
## できるかぎり、元の情報 ( a_ij ) を残しておく形になっている。
### 0 を入れてしまうということは、重要な情報を「壊している」ことになる。
### 元の形の値の方が、後で使いやすいかも.. ( 計算をすると情報が落る.. )

==

LDU 分解

	A = LDU

	LDU の U は、LU の U の列要素を D の対角要素を割った数になっているだけ

	# LDU はどのような場合に利用するか ?
	#	n x n の場合は、LU で Okey
	#	n x m の場合は、LDU の場合が望ましい場合がある
	#		# 解が一意に定まらない場合は、最小二乗値を小くする値を選ぶ
	#		# 物理的には、エネルギーを最小にするような解になっている

==

直接法には、色々なバリエーションがあるが、本質的には、Gauss の消去法になっている。
	# 二百年たっても、誰も Gauss を越えられない ( やっぱり、Gauss は天才 )

	直接法に関しては、Gauss の消去法だけをしっかり学べば十分
		# 計算機の手法は、移りかわりが激しいが.. 
		# 皆なが死ぬまで、この方法だけは変らない !!


# シラバスの中で、変化が絶対にないと思われる内容は、この Gausss の消去法だけ
# ( 他のものは、数年後に変っているかも.. )
## 研究分野としては、面白くないが、使う方からいえば安心して利用できる方法と言える。

==

[課題]

LU 分解 ( 福井先生の方法 ) で、連立方程式を解く

データは、text p.38 を利用する

当然、結果も同じになるはず。