3次元アフィン変換行列のべき乗

ついつい慣例に従って、アルゴリズムの定まっているほうに堕落してしまいます。

ほんとならpixivにもアップしたほうが、というかpixivのほうを重点的にアップしたほうが宣伝効果がずっと高いとは思うのですが

ブログにアップするのとちょっとアルゴリズムが違うんですよね。


ブログだったらテキストの合間に図を挟むことが可能なのに対し
pixivのほうは、上にキャプションがあって、メインは画像という形を取るので
テキストを画像に埋め込む形も取らなければならず

ブログをはいアップしました、といってもそこからそのまんま流用してpixivに張り付けるという単純な作業ではないわけですよ。


なんかね、ブログのほうでさえ時々アップするアルゴリズムを忘れそうになるんです。

「文字を作りながらその都度計算して、その様子をキャプチャーしてアップするだけの簡単なお仕事」ってアルゴリズムに落として考えないと、最近はろくにブログも更新できなくって・・・


精神的余裕がないと、そのアルゴリズムの存在をすっかり忘れてどんよりしてしまうので困るんですよね。


pixivのアップ方法のアルゴリズムも、さっさと確立させて慣れていかないといけませんね。
といってももうpixivもだいぶ老年期に入ったようなコンテンツのようなので
何か別のSNSを探したほうがもしかしたら有意義なのかもしれませんが
僕はいつもいつも後手に回ってしまうのです。



＝＝＝＝＝＝＝＝＝
さて、今回はただの自己満足以外の何物でもないんですが

3次元のアフィン変換(でいいんでしたっけ？)行列のべき乗についての話です。


ジョルダン標準形以前に、まずこれを整えておかないとと思いまして。



こんな変換行列Aがあるとします。
x軸方向にa、y軸方向にb、z軸方向にcだけ拡大しながら、
x軸方向にx0、y軸方向にy0、z軸方向にz0だけ平行移動させる変換行列です。
この行列Aのべき乗を考えてみましょう。


平行移動を含んでいるので、3次元の変換といっても、ダミー次元のせいで4次の行列になっています。


べき乗を計算したいので、この4次の行列の固有値λをまず求めてみましょう。
(a-λ)(b-λ)(c-λ)(1-λ)=0なので、λ=a,b,c,1となりますね。
4次の行列式なのでサラスの方法が使えませんが、すでに行ごとに掃き出し終えている形なのですんなり計算できるのがわかるかと思います。




では、λにa,b,c,1を順に代入していった際の固有ベクトルを求めてみましょう。


まずはλ=a
この連立方程式(ただし永年方程式)を解けばいいので

v2=v3=0とわかります。
また、x0v1=(a-1)v4であることもわかるので
求めたい固有ベクトルvは以下のようになります。
 


次にλ=bについて解きますと
 
λ=aのときと同様に
v1=v3=0と
yv2=(b-1)v4が条件なので、固有ベクトルは
 


λ=cでは

v1=v2=0
zv3=(c-1)v4



λ=1では

v1=v2=v3=0



これらの固有ベクトル(縦ベクトル)を、λがa,b,c,1の順に横に並べると
対角化のための行列の片割れPが以下のように求められます。
 
ちなみに、固有ベクトルを4つとも規格化できるときはPはユニタリ行列になるはずです。
量子力学でいうところのケットベクトルに相当しますね。


この「対角化のための行列P」のもう片方は、Pがユニタリならエルミート共役(転置して複素共役)をとればすぐに求まるのですが

規格化してないときは、Pの逆行列を求めることになります。

その際に、Pの行列式も求めることになりますが、これは以下のようになります。



Pの逆行列は、次のように余因子展開をしてから転置して求めます。

 


ここでは、invP×A×P=Jという対角化された行列になる確認は割愛して

変換行列Aのべき乗を直接求めてみます。

A^n=P×J^n×invPとなるので、以下のようになります。



実は、この計算にはa-1とb-1とc-1という部分が分母に含まれています。
ですので、a=1またはb=1またはc=1で計算不能になってしまいます。
この「または」は排他的XORではない論理和ORなので、a=b=1やc=a=1、a=b=c=1のときも
もちろん計算不能になります。


しかし、この式をよく見ると
たとえば をa-1で割る計算だったら以下のように多項式の割り算を行って


と、aの0乗からaの2乗までの級数として表現しなおすことが可能で
ｎを一般化しても


とすることが可能です。

特に、これを使うときはa=1の場合に限られるため

とすることが可能で

の例外処理は、最悪、a=b=c=1だった場合でも

このように処理することが可能です。
もちろんa,b,cどれか1つか2つが1でないときの例外処理も大丈夫です。

まあ、等倍の平行移動というイメージを忘れていなければ、当然といえば当然の帰結です。

また、平行移動がなくそれぞれの軸方向にスカラー倍だけしたいときの計算も問題ありません。


解析的に計算すれば、これらの例外はこのままでも簡単に取り除くことができます。


しかし、もし行列のまま計算結果を得たい場合は少々異なってきます。

たとえばa=bだった場合やb=c=1だった場合などは、固有値がカブるので、
それなりの措置を取らなくてはいけません。

ランク落ちがあるかもしれないので、必ずしもジョルダン標準形になるとは限らないようなのですが。

ただ、解析的に計算した場合の結果のみを用いる際は
このように 横ベクトルの後ろに変換行列を掛け算する方式でも
このように 縦ベクトルの前に変換行列を掛け算する方式でも、議論は同じということになります。(平行移動の項の位置に注意)


このような、数値計算ではない解析計算でありがちな困ることは、汎用性がなく、「計算が閉じている」ことですね。


以下に、アフィン変換の様子の例をgifであげておきます。

 
x軸倍率→y軸倍率→z軸倍率→x軸平行移動→y軸平行移動→z軸平行移動


x軸a倍しながら平行移動→y軸b倍しながら平行移動→z軸c倍しながら平行移動


3軸同時拡大→3軸同時平行移動→3軸同時拡大しながら平行移動



あ、そうだ。忘れてましたがExcelソースファイルを置いておきますね。
ダウンロード

応援したいと思ってくれたら押してね＾＾

にほんブログ村