Maxima

2011年04月25日

Maximaを使ってみた その37

さて、実際の数値で直線近似をしてみたい。

例えば下のようなX−Y対になったデータがあったとして、

X Y
1 5
8 9.5
6 8
3 5.2
12 10

これらの点に対して直線の方程式

ax+by + c = 0

を適用すると、全ての点を通過する直線はありえないので誤差εが出る。

fig05

これらの式を行列の形にまとめると、

fig06

となり、

fig07

と置くと、

||AX+c1||2

が最小になるXを求めると、rank(A)=nの場合はAを特異値分解して

fig02

としたときに

=−cdiag(σ1-1,・・・,σn-1

になるのが前回の結論だった。

そこでAを定義して、

fig01

これを特異値分解していく。まずはAAを計算して

fig05

ここから固有値と固有ベクトルを求める。

fig02

固有値は大きい順に並べる必要があるので、λ=546.0、λ=14.26 とし、それに対する固有ベクトルはv1=[1.0,1.104]T 、v2=[1.0 , -0.91]Tになる。

ここで得られた固有ベクトルの内積をとってみるとほぼ0なのでこの2つのベクトルは直交している。

fig11

右特異行列Vは実直交行列で、VV=Iにならなくちゃならないので、固有ベクトルのノルムがそれぞれ1になるように正規化する。 これらのベクトルを列とする行列Vを組み立てる。

fig03

次に固有値λから特異値σ(ただしσ=√λ)を作り、それを組み立てて行列Σを作る。

fig04

次に、得られたσと右特異行列から左特異行列Uを作る。まずAVを計算して、

fig06

得られた行列の各列をベクトルと見なして対応する特異値で割ると左特異ベクトルu が得られる。それらを列ベクトルとして組み立てて左特異行列Uを得る。

fig07

本来ならUは正規直交化しなくちゃならないんだけど、Xの式では3列目以降のベクトルは必要ないのでそのまま使う。

これでUもΣもVも得られたのでこれらが本当に行列Aの特異値分解になっているか掛け合わせてみると、

fig09

となって確かにこれは行列Aだ。

そこでこれらを前回求めた誤差ベクトルεのノルムを最小化するベクトルXの式に当てはめると、

fig08

となる。これらを直線の式に代入してyの式に直すと、

fig10

となり、与えられた点を最もよく表わしている直線の法的式はy=0.48X+4.753となった。実際にプロットしてみたのが下の画像だ。確かに各ポイントの傾向を表わすような直線になっている。

fig12

それではまた次回。


take_z_ultima at 11:30|この記事のURLComments(0)TrackBack(0)

2011年04月18日

Maximaを使ってみた その36

数学は難しいなぁ。いろいろ思い込みで処理しちゃって、あとからよく考えてみるとダメじゃんって事が結構ある。このブログもいろんな事を間違いながら徐々に正解に漸近してくしかないかなぁ。

さて今回は特異値分解の応用として最小自乗法をやってみたい。

最小自乗法は前にもやったけど、測定などで得られたバラツキのあるデータを適当な関数で近似する時に関数とデータの残差の自乗の合計を最小にするように関数のパラメーターを決める方法だ。

さてそこで下のようにX−Y平面に対して点が散乱しているところへ1本の直線を引いてこれらの点が表わす傾向を直線で近似してみたい。

fig01
 

2次元空間での直線の方程式は

ax+by + c = 0

で表わされる。そこで散乱する点の1つP の座標値( x ,y )を直線の式に代入した時、パラメータa,b,cが点P を通過する線であれば、

ax+by+c = 0

となって式が成り立つ。

しかし実際には全ての点でこの式を満たす事は無理で、誤差ε が出る。

ax+by+c = ε

そこで次善の策としてこの誤差ε の合計が最小になるようにa,b,cを定めようと考えるわけだけど、εにはプラスもマイナスもあるのでそのまま足すのはうまくない(εが100でε が−100でも合計誤差が0って事になっちゃうからね)。

そこで符号が消えるように誤差εの自乗の合計を使って誤差合計の評価をしちゃおうってのが最小自乗法なわけだ。

1〜mまでの点を直線の式にあてはめると、

fig05

これを行列にまとめて

fig06

 ここでここで x ,y を行ベクトルとして持つ行列Aを定義し、パラメータa,b,cをそれぞれベクトル(a,b)T、c1(ただし1は[1,・・・,1]T)とすると、

fig07

と置き換えて本来は

+c1=0

となればいいんだけど実際は誤差ベクトルε(ε1,・・・,εn)を使って、

+c1ε

となる。ここで誤差が最小になる事をεベクトルのノルムが最小になる事とすると、それは||ε||が最小になるのと同じで

||ε||=||AX+c1||2

ここでrank(A)=nとしてAを特異値分解すると、

fig02

だから

fig03

dの部分は定数ベクトルなのでこれを最小にするのはΣ1y+dの部分を0に持って行くしかない。

fig04

としてパラメータa,bがcの式で求まる。直線の方程式は

ax+by+c=0

だからaとbの項にcが入るなら式の両辺をcで割ればcの項は消える。

それではまた次回。



take_z_ultima at 11:30|この記事のURLComments(0)TrackBack(0)

2011年04月11日

Maximaを使ってみた その35

今回は特異値分解を実際にやってみたい。

MaximaではLAPACKをロードするとdgesvd( )関数が使えてそのまま分解することが出来るのは以前に書いたけど、今回は自分でやってみたい。

まずは表示するケタを調整しておく。

fig02

次に分解する行列を用意。

fig01

次にAAを計算して変数ATAに代入する。

fig03

次にAAをシュア分解する。

シュア分解は

fig02
 

という形でQはユニタリ行列である。

今回分解するAAは1つの行列を転置して掛け合わせた行列なので対称行列になっているので、相似変換した時も対称性が保たれることが期待できて、右上三角行列を作るように相似変換すると対角要素の上の部分も0になって対角要素だけが残る。だから2X2の行列AAのシュア分解はこういう形になって、

fig05

 Qがユニタリ行列だから

fig06

であり、式の両辺にQをかけると

fig07

となる。ここで行列Qを列ベクトルに分けて[q1,q2]とすると、

fig08

となって、q1、q2は行列AAの固有ベクトルであることがわかる。そこでATAの固有値と固有ベクトルを求めると、

fig04

となる。そしてQQ = I であるためには

・q=1

・q=1

・q=0

にならないとならないので、それぞれのノルムが1になるようにベクトルをノルムで割った。

fig09

固有値は3.447と36.55で、その固有ベクトルをq、qにしたけど、特異値分解の標準形は特異値を大きい順に並べることになっているので、ここから先はq、qの順に並べて使う事にする。ここでこれらのベクトルを行列に組み立ててVベクトルとする。

fig10

上の式は転置するのを忘れてましたが、結果が同じになるので気付きませんでした。

これがユニタリ行列になっているかちょっと掛け合わせて見ると、確かに単位行列になるし、ATAを相似変換してみると、確かに対角に固有値が現れて、他の要素はほぼ0になる(数値誤差はとりあえず見て見ぬふりw)。

fig11

これが前回出てきた

VT(ATA)V=diag(σ12,・・・,σr2,0)

なわけで、Vは右特異行列になる。

σ12=36.55 σ22=3.447

だからそれぞれルートで開けば

fig12

ここでAVを縦ベクトルPi の行列として表わすと、

AV=[P1 ,・・・,Pr ,・・・,Pn

となり、左特異行列U の列ベクトルは

i = Pii

で求められ、Aが3X2の行列なのでUは3X3になって、Piは2までしか無いので残りの1つはRの正規直交基底に拡大するためにu1とu2が張る面に対して垂直な単位ベクトルu3をベクトルの外積から作って、

fig13

これらを組み立てて左特異行列Uにした。

fig14

転置行列と掛け合わせて見ると(誤差を無視すれば)単位行列になるのがわかる。

fig15

特異値分解の標準形Σも3X2で組み立てて、

fig16

特異値分解出来てるかどうか、これらをかけ合わせてみると、

fig17

となり、行列Aと一致した。ちょっと途中が怪しいけど行列Aを特異値分解できた。

それではまた次回。



take_z_ultima at 11:30|この記事のURLComments(0)TrackBack(0)

2011年04月04日

Maximaを使ってみた その34

今回は特異値分解の話だ。

まずは定義から。

与えられたm×n実行列Aは適当なm次実直交行列Uおよびn次実直交行列Vをとれば、

fig01

の形に分解できる。σ1≧・・・≧σr>σr+1=・・・=0をAの特異値といい、この分解をAの特異値分解と言う。また、Σを特異値分解の標準形という。ここでrはAのランクで太字の0は空いているところが全部0であることを表している。

行列を特異値分解するにはシュア分解を使う。シュア分解は適当なユニタリ行列Qによる相似変換がn×n行列Aの特性多項式の解(λ1〜λn)を対角要素とする上三角行列Tになるというものだ。以前にやったQR法のところで嫌と言うほど見た奴だ。

fig02

特異値分解の式においてUもVも実直交行列であるから

UTU=I

VTV=I

であり、AAT、ATAの式に特異値分解の式を代入して整理すると

AAT=(UΣVT)(UΣVT)T=(UΣVT)(VΣTUT)

  =UΣΣTUT

  =U diag(σ12,・・・,σr2,0) UT :m×m

ATA=(UΣVT)T(UΣVT)=(VΣTUT)(UΣVT)

  =VΣTΣVT

  =V diag(σ12,・・・,σr2,0) VT :n×n

となり、さらに変形すると

UT(AAT)U=diag(σ12,・・・,σr2,0) :m×m

VT(ATA)V=diag(σ12,・・・,σr2,0) :n×n

AAT、ATAは実対称行列であり、実対称行列のシュア分解は対角行列になるので、これらはAAT、ATAのシュア分解の形になっている。このAAT、ATAを使うのが特異値分解のミソで、扱うのが大変そうなm×nの行列が自乗の固有値を持つ素直な実対称行列で扱えちゃうし、対象行列はシュア分解の操作をしても対称なので対角要素の下の部分を0にすると対角要素の上も0になっていつのまにか対角行列が出来ちゃうようだ。

ここでAVを縦ベクトルPi の行列として表わすと、

AV=[P1 ,・・・,Pr ,・・・,Pn

となり、

VT(ATA)V=diag(σ12,・・・,σr2,0)

よりPiTPj は、

i = j = k なら PiTPj = σk2

i ≠ j なら PiTPj = 0

である。よってδ関数(添字iとjが等しいとき1でそれ以外は0になる関数)で表わせば、

(Pii )T(Pjj ) = δij  ( i , j = 1,・・・,r)

Pj = 0  ( r < j ≦ n)

ここで

i = Pii

と定義してこれをRの正規直交基底に拡大し、これらを列ベクトルに持つm×m行列U

U = [ u ,・・・,u ] :m×m

を定義すると、

PiTPi = σi2

i = Pii

だから

iTPi =( PiTi )Pi = PiTPi ii2i = σi

よって、

UTAV = [ u ,・・・,u ][P1 ,・・・,Pr ,0,・・・,0

    = diag(σ1,・・・,σr ,0,・・・,0)

これはまさにAの特異値分解の形である。よって行列Aの特異値分解はAAT、ATAのシュア分解から導出することが出来る。

実際の導出はまた次回。



take_z_ultima at 11:30|この記事のURLComments(0)TrackBack(0)

2011年03月28日

Maximaを使ってみた その33

前回は回転行列から回転角度を求めてみた。今回は回転軸ベクトルを求めて見たい。

そこで再びロドリゲスの公式に戻って

fig02

上の式の両辺に右から r をかけると、[ r ]x r = r × r = 0 よりsinθの項と(1−cosθ)の項は0になって

fig05

となる。これはRが r を回転軸にしているのでRで回転変換しても変化しない事をあらわしている。

また、固有値と固有ベクトルの関係

fig06

より rはRの固有値λ=1の時の固有ベクトルと捕らえることも出来る。

Rは回転行列なのでRR=Iであり、上のR r = rの式の両辺に左からRをかけると、

fig07

をかける前の式から辺々引くと、

fig01

よって、

fig02

ここから

fig03

となる。ここで

fig04

となり、

fig08

とすると、

fig09

が回転行列Rの回転軸ベクトルになる。

試しに前回同様に回転軸が[0.2,0.3,−0.4]方向で回転角度30度の回転行列Rを作って、そこから回転軸ベクトルr_hatを算出し、軸ベクトルrを単位ベクトル化して比較してみた。

fig10

回転行列Rから軸ベクトルを算出できる事が確認できた。

それではまた次回。



take_z_ultima at 11:30|この記事のURLComments(0)TrackBack(0)
Archives