有名な数学の未解決問題に”リーマン予想”と言うものがあるが、その証明は数学的に重要な内容を含んでおり、多くの数学者が取り組んでいるとのこと。最近その解説本(文献1)を読んでいると、まったく関係の無い様に思える『素数』と『行列の固有値』の間に何らかのつながりが有りそうだとのことだ。そもそも、その本は素数と量子力学(サイト内リンク) のつながりの説明に興味を引かれて読み始めたのだが、これは意外な発見だ。そこで『固有値ってどんな奴だったっけ?』程度に、知識を整理してみよう。文中のグラフや図、計算結果等は数式処理のMathematicaを利用した。興味がある方は計算機代数ノート『 [CAS-Lab] 行列の固有値(座標変換と図形変形) Mathematica[編]』もの方もどうぞ。
尚、用語に関しては自分の理解したイメージや感覚でつぶやいているので、その場で浮かんだ用語を使っている。あまり正確じゃないと思うので信用しないでほしい。
1. ベクトル空間
線形代数の試験で『次に示す行列Wの固有値と固有ベクトルを示せ』とか出されて苦しめられた嫌な思い出は無いだろうか?
Wは正則な正方行列でvは固有ベクトル、λは固有値でスカラーだ。Wがn次なら、vとλの組もn個存在する。いろいろなサイトで解き方の詳しい説明がなされている。なので解き方については(何せ嫌な思い出なので)ここでは省略する。ここでは主に固有値・固有ベクトルとはいったい何者で、どういう所で、どの様に働いているか?を中心に実感が得られる程度に理解を掘り下げてみよう。
今まで調べた自分の知識で言うと、固有値の物理実態は”フィルター”とか”歪み”に近いイメージだ。音楽を聴くときに”イコライザー”を操作したことは無いだろうか?例えばイヤホンが貧弱なので低音を持ちあげたいとか、特定の音域を強調したいとかに利用する。私の思いつく限りで、この”イコライザー”という一種のフィルターが、身の回りで実際に触れて、かつ、いじることができる固有値の数少ない実例だと思う。イヤホンを電気信号から音響信号への変換器と考えるとその周波数特性は固有値と言えると思うのだが、こちらの方はイヤホンの値段固有で買い替えるまで固定されている。なのでイコライザーをいじって補正している?
ということで数学の試験問題でしかないようにに思える『固有値』がいかに『固有ベクトル』と連携して、何をどのように歪ませるかに注目してほしい。
1.1 基底とベクトルの成分の関係
まずは座標変換で図形が歪む例をみてみよう。そのためにはベクトル空間で考えを整理する必要がある。手始めに3次元ベクトル空間の表し方を思い出そう。いわゆる3次元ユークリッド空間の任意の座標点に対応する位置ベクトル p は、直交する3本の基底ベクトル b1, b2, b3 を使って、以下の様に表すことができる。
x、y、zは基底ベクトルを何倍するかを示すスカラー量で位置ベクトルの成分となる。この仕組みは3次元物理空間にこだわらなければ何次元に増やしてもよい。
次に図形の変形を考察するので2次元に落としてみよう。
となり、以下の図の様ないわゆるx-y座標だ。もちろん基底ベクトルの要素数も3から2に減っている。この例の通り基底ベクトルに単位行列 が設定された場合は、基本ベクトルと呼び、eで表そう。
Fig1.1-1
この基底ベクトル b1(=e1), b2(=e2) をそれぞれ何倍かすることにより座標上の任意の点Pを指定できる。この座標を、X系と呼ぼう。例えば半径2の円をを描く場合は以下の様になる。いま在る矢印の位置ベクトルPxは b1, b2 のそれぞれ√2倍して作ったたものである。
Fig1.1-2
逆に任意の点 Px が与えられた場合、基底ベクトルb1、b2それぞれ何倍すればその位置ベクトルを再現できるのか?つまり、成分を求める方法だが、この場合は単純に与えられた位置ベクトルと基底ベクトルの内積をとれば良い。
Fig1.1-3
任意の基底ベクトルを使って位置を指定する方法を数式で表してみよう。二つの基底ベクトルを横に束ねて、2行2列の行列B=[b1 b2]を基底行列と定義し、倍率の方をベクトルSとする。これによりR=B.Sと書ける。これはベクトル空間そのものの定義といえよう。なので、Rを与えてその成分を算出するには、S=B-1.Rとすれば良いことになる。数式上はその通りなのだが、なかなか実感がわかない。直観的に納得できる様にその計算の仕組みを座標変換と併せて説明しよう。
まず傾いた基底ベクトルが張る緑の座標を表示すると、以下のFig-1.1-4の様になる。この緑の座標をU系と呼ぼう。この傾いた系のままでも点Rを緑の基底ベクトルを使って表すことができる。b1は1.5倍、b2は2倍だ。
Fig.1-4
ここで少し座標変換について触れておこう。座標変換とは今、X系の任意の点を {x1, x2}と表した場合、U系の座標 {u1, u2} を2組の変換関数 u1=f1(x1, x2) と u2=f2(x1, x2) で表すことだ。この変換を順方向の変換とした場合、逆変換とは逆関数の x1=f1(u1, u2)-1 と x2=f2(u1, u2)-1 によりU系の座標からX系の座標を求めることを指す。
上記Fig.1-4は青のX系を基準系として、その座標{x1, x2}を指定すると、対応するU系の座標 {u1, u2}が読み取れるようになっている。逆に緑のU系の座標を指定してX系の座標を読み取る場合は逆変換 x1=f1(u1, u2)-1, x2=f2(u1, u2)-1 を行っていることになる。
以降、本稿で出てくる例では、変換関数の中身が以下の様に線形変換の形式で表すことができる
と、逆変換では
のみを扱う。尚、座標変換にさらに興味がある場合は本サイト:AcademiaNotes/物理学ノート/相対性理論における線形座標変換関連項目や一般座標変換県連項目をご参照願いたい。
それでは、緑の座標系に視点を移してみよう。『視点をU系に移す』とは、緑の直立した直交座標を用意し、Eq1.1-3aで算出されたU系成分を使って、基底ベクトルや図形を描き直すということだ。つまり表示されるすべてのベクトルに基底の逆行列B-1を掛ければ良い。
結果、以下の様にR点を指す位置ベクトルがRからSに変わり少し変化すると同時に、緑の斜交基底が直交し、基本ベクトルに化けて直交座標と重なっていることが分かる。さらに元のX座標にもB-1を掛け、薄い青で表示を残した。見ての通り逆方向に傾いているが、矢印先頭点の青座標は {2, 2} を正しく指している。これがEq1.1-3aの座標変換の視覚化例だ。
Fig1.1-5
見ての通り単純な直交系なので、b1の成分が1.5、b2の成分が2であることが簡単に分かる。
このB-1で注出した成分Sを元の座標系の基底ベクトルBに掛ける(つまり、R’=B.Sの計算をする)と、以下の様にRが再現できるというからくりだ。
Fig1.1-6
基底ベクトルと内積を取って正しい成分が得られるのは。基底ベクトルそれぞれの長さが1で互いに直交している場合の特例と考えよう。
話をまとめると、成分を参照するベクトルを R、参照した成分で再現されるベクトルをR’ とすると成分の抽出は
と一般化できる。そして抽出成分 S で R を再現しようとする場合は
となる。従ってEq1.1-4aと4bを続けて書くと
1.2 座標変換と固有値の関係
と書ける。一方 Eq1.1-4 参照ベクトルの再現式:R’=B.B-1.R は固有値行列でλ1=1, λ2=1 と置くと
と書ける。
かたや、Eq1.2-1の変換行列Wを使ったベクトルの変換操作の式を書くと、
あれ、こいつらは同じ仲間かい?両者比較すると、
ベクトル成分の抽出→再現 | R’=[B.Λ.B-1].R | λ1=1, λ2=1 | R’=R |
固有値の定義 | R’=[V.Λ.V-1].R | λ1,λ2は任意 | Vn方向に R’=λnR |
基底ベクトルは目的の座標系を張るために、意図して選んで設定する。固有ベクトルVの方は、何かの変換行列が与えられ、その固有値と固有ベクトルを計算しろと言われる。しかし両者は同じ物っぽい。
ならば、基底ベクトルを固有ベクトルとみなして、座標変換の過程で図形の位置ベクトルが固有ベクトル方向(つまり基底ベクトル方向)にΛで圧縮伸張する様子を確認することができる。最初は直観的に分かりやすい回転変換を試してみよう。変換行列 W の中身を以下の様に3段階に分けて変換の様子を見てみる。
変換1 | W=B-1 |
変換2 | W=Λ.B–1 |
変換3 | W=B.Λ.B-1 |
回転変換の例 W=B-1
以下の様に現在のX系基底ベクトル e に対し左にπ/4 回した基底ベクトルを新たに作り、U系の基底行列とする。
このU系を緑で描画し、元の青の座標系Xと重ねてみよう。
Fig1.2-1
この様に視点を青のX系とすると、U系は緑の基底ベクトルとともに π/4 左に回転している。X系とU系とでの位置座標の読み替えはPu=B-1.Px、 Px=B.Pu となる。
次に以下の通り視点を緑の座標系Uに移してみる。そのためには基底ベクトルや位置ベクトルに式B-1を掛け描画すれよい。結果は見ての通り、基底ベクトルは緑の基準系の座標と重なり、X系で描かれた図形は45度右に回転した。
Fig1.2-2
回転変換の例 W=Λ.B-1
(λ1=1, λ2=0.5) を持つ固有値行列Λを追加して、位置ベクトルPxにW=Λ.B-1を乗じると以下の様に図形が縦方向(λ2方向)への1/2圧縮となる。
Fig1.2-3
正確にはここで Λ(図中ではLと表示) を加えることは、さらにもう一段の縦横への圧縮伸張変換が起きるので黄色のV系座標でも用意して、(図形の位置ベクトルだけではなく)基底ベクトル、つまり座標系込みで、変換後を表示すべきだ。が、見ての通り縦横にしか圧縮伸張されないので理解は容易だ。黄色い座標系では図形は変形しない代わりに、緑と青の座標が伸びる様子がうかがえるのだが、ここでは省略する。
回転変換の例 W=B.Λ.B-1
最後に全ベクトルにW=B.Λ.B-1を掛けると緑のU系と図形は以下の通り最初の角度に戻ってくる。
Fig1.2-4
しかし、図形は基底ベクトルBのb2方向にλ2=0.5だけ圧縮されている。これが変換行列 W=B.Λ.B-1の固有値の図形を歪める効果であり、b2はその固有ベクトルとなる。もちろんλ1に何か倍率を設定すればb1方向にも同時に圧縮伸張させる事ができる。
斜交座標変換の例
このWによる変換を利用すると、好きな方向に図形を引き延ばしたり縮めたりできる。もう少し図形をいじってみよう。横幅を1/2してみる。さらにPxをそれ自身の矢印の方向(右45度方向)に2/√2引き延ばすと、丁度 {2, 2} の位置に届くはずだ。そのためには固有ベクトル役の基底ベクトルと固有値を以下の様に設定する。
基底ベクトル | 固有値 | |
b1=[0, 1] | λ1=1/2 | 水平方向に1/2圧縮 |
b2=[1, 1] | λ2=2/√2 | 右上45度方向に√2伸張 |
Fig1.2-5
W=B.Λ.B-1変換を実行する。背景の緑座標はW=B-1変換時のU系を表示している。
Fig1.2-6
以上が座標変換で固有値、固有ベクトルがどのように図形の変形に働くかを明らかにした。
いままでのまとめ
行列Wを基底行列Bから以下のように構成した場合の基底ベクトル(基底行列Bの各列ベクトル)は固有ベクトルvとなる。
なぜなら固有値、固有ベクトルの定義の式は
であり、WにEq1.2-4を代入しは変形すると、以下の様にEq1.1-4と同じく機能別に二つのカッコで分けて考えることができる。右側のカッコ内で与えられたベクトルvの成分を抽出し左側カッコに注入する。左側カッコは、その注入成分に対応する基底ベクトルを対応する固有値[λ]倍して出力する。
固有値、固有ベクトルの定義の通りだ。従って変換行列Wの基底ベクトルを特定できたら固有値問題が解けたと言るだろう。
ただし、図形の変形時に理解できたと思うが、圧縮伸張には基底ベクトルの大きさは問われず、方向のみが必要だった。また左右に圧縮とか上下に伸張などの様に方向のプラスマイナスを問わない緩さがある。なので基底ベクトルは固有ベクトルではあるが、固有ベクトルが求まったとしても必ずしも基底ベクトルそのものと言うわけでは無い。
尚、Fig.2-6などの図形の変形例で分かる通り、固有ベクトルと一致していないベクトルが入力されたの場合でも、右カッコ内では各基底成分に分解されるので、ちゃんと線形加算された変形図形がが生成・出力される。
この先、固有値の話は基底を周波数領域に取るフーリエ変換につながり、さらにラプラス変換により複素領域まで広がる。参考文献には複素領域での議論も含まれているのでもう少し掘り返してみよう。
1.3 時間-周波数領域変換
続いて時間領域から周波数領域への変換を見てみよう。代表的なものに有名なフーリエ変換がある。ちょっと難しそうに聞こえるが、非常に大雑把に言うと、基底BにSin/Cos波を使っただけだ。
前節からのつながりの良さから、離散時間及び離散周波数の領域から始める。尚、本節の詳細は本サイト:AcademiaNotes/計算機代数ノート/信号処理の基礎数学における関連項目、文献2をベースにしている。御参照願いたい。
アダマール変換
基底ベクトルを以下の様にとってみよう。この基底行列Bはアダマール行列(ウォルシュ型)とよばれている。
nは基底サイズで決まり、この場合4となるが、つじつま合わせの様な係数なのであまり気にする必要はない。重要なことは基本ベクトルと同様、各基底が直交していることだ。例えば b2.b3= 0とかb3.b4= 0の通りだ。
変換に際しては前出の線形変換と全く同じ仕組みが使える。例えば時間→周波数の変換はS=B-1.Rであり、参照ベクトル成分の抽出や基底注入によるベクトルの再現では、下記の式(Eq1.1-4)が同様に使える。
さらに言うと固有値の下記の式(Eq1.2-1)も同様に使える事を確認するのが本章の目的だ。
どちらも行列とベクトルの掛け算だけで定義された変換なので、これ等の式自身が正しく機能することは想像できるが、疑問は『何でちょっと基底が変わっただけで、周波数領域に変換できるの?』ということだろう。Eq1.3-1の基底の中身をプロットしてみよう。
Fig1.3-1
上から順に 基底ベクトル b1, b2, b3, b4 だ。少しばかり基底の波動を感じられないだろか?基底数を8に増やしてみる。上から順にゼロクロスの数を数えてみよう。
Fig1.3-2
ゼロクロス数は0、1、…7となる。表示期間を1秒とした場合、それぞれ 0Hz、0.5Hz、…3.5Hzの波だと言えないこともないか?Rに以下の1Hzの離散Cos波(3番目の基底に最も近い!)を設定し、周波数領域に変換してみよう。
Fig1.3-3 入力R(Cos波)
ベクトルの再現を示す式 Eq1-1はカッコ内が成分抽出の機能を担っていた。
従って周波数成分の抽出はカッコ内をそのまま使えるはずだ。以下の様にカッコ内を取り出してRを入力信号とすると、出力Sをに周波数成分が得られるはずだ。結果は以下の通り。
Fig1.3-4 アダマール周波数成分
期待通り、3番目の基底成分(1Hz相当)が一番多く検出された。つまりこれがアダマール変換だ。Eq1-1の示す通り、この抽出された成分(周波数成分)を基底Bに掛けると入力された時間領域の信号が再現されるはずだ。(もちろん理由は、間に固有値を含ませていないからだが…)
計算の結果は以下の通り。
Fig1.3-5 再現波R’
時間信号Rを完全に再現したといえる。つまりこれが周波数領域から時間領域への変換であり、アダマールの逆変換となる。尚、アダマール行列は要素が±1だけで出来ているばかりか、なんと B=B-1 であり、変換行列と逆行列が共通に使えると言う非常にエコな変換行列だ。
離散コサイン変換(DCT)
もう少し波動っぽい基底に離散コサイン行列がある。以下の様にちゃんと波打っているのが分かる。フーリエ変換ではSin波Cos波を合わせて使うが、Cos波だけ使った半欠け変換といったところか?
Fig1.3-6
このDCTで画像信号を変換すると低周波側に成分がうまく集まるとのことで画像圧縮に向いており、JPEGや、MPEGなどで利用されているとのこと。
Fig1.3-2と同じ1HzのCos波を入れてみよう。周波数成分の結果は以下のFig1.3-6 の通り。
Fig1.3-7 DCT周波数成分
参照:アダマール周波数成分
アダマール変換での周波数成分に比べて、より1Hz近辺への周波数成分が集中していようだ。やはりDCT の方が手間が掛かっている分だけノイズが少なく優秀なのか?
入力信号の再現はアダマール変換と同様、正確に行われており結果は割愛する。
離散フーリエ変換(DFT)
DFTはフーリエ変換(数学的関数)のデジタル版(離散時間・離散周波数)といったところか。フーリエ変換・フーリエ級数との関連は割愛する。DFTの基底は以下の通りだ。
この基底は以下の様にSin波Cos波を合体させた、複素指数関数 [ Aeiθ ] の形式で表されている。
Fig1.3-8
基底行列BをN=8で生成してみよう。
Fig1.3-9a DFT基底行列B(複素数表示)
なじみのあるSin,Cosと違って、複素指数関数と言うと少しばかり小難しい響きがするが、Fig1.3-8で分かる通り、このくるくる回っている棒の先端の水平方向の成分がCosで垂直方向がSinだと読み替えれば良い。各基底ベクトルは各列に相当し、縦に目で追ってみると、くるくると左回転しており、担当周波数としては左第1列から第8列までを
[+0Hz,+1Hz,+2Hz,+3Hz,-4Hz,-3Hz,-2Hz,-1Hz]
と見なせる。
変換行列B-1の方は以下の通りで、基底行列とは右回りの逆回転(共役関係)をしている。
Fig1.3-9b DFT変換行列B-1(複素数表示)
DCTと同様に参照信号Rを1HzのCos波とし、上記変換行列を使って周波数成分を抽出してみよう。Fig1.3-10はRの波形をDCTの場合と同じく実数表示したグラフだ。
Fig1.3-10 入力R Cos波
しかしフーリエ変換では複素数で表示する方が理解が早い場合が多い。以下Fig1.3-11に複素数で表示してみる。入力信号はあくまで実数信号となるので複素数表示をした場合でも実数軸、つまり横軸に振幅が現れる。
Fig1.3-11 入力R Cos波(複素数表示)
横方向にうねうね波打っているのが分かるだろうか?
続いて S=B-1.R でフーリエ変換した結果が以下のFig1.3-12となる。
Fig1.3-12 Cos波の周波数成分
フーリエ変換で抽出された周波数成分は通常複素数となる。しかし入力が純粋にCos波をの場合は上記の通り実数部のみの成分を持つ。基底の+1Hzに相当する2列目成分と、-1Hzに相当する8列目成分に全て集約されており、変換でのノイズは見当たらない。
次に、この成分を基底行列Bに掛けてRを再現してみよう。基底行列Bと周波数成分ベクトルSの掛け算は、通常の行列演算と同様だか、複素数どうしの演算となるので少しトリッキーだ。
以下の図では基底行列Bの上にベクトルSを重ねてみた。
Fig1.3-12a B.S演算の内側 Cos波
赤い枠で囲った様に、基底Bの2列目の全ての要素に、周波数成分ベクトルSの2列目の要素が掛かる。Sの方は実数なので向きは変わらずSの振幅に応じて大きくなる。8列目も同様だ。2列目と8列目以外はSの要素はゼロなので、演算の対象にならない。
その後、演算結果は2列目と8列目それぞれ行の向きに[r1,r2,…r8]に足しこまれて、R‘の演算結果として出力される。
以上がR’=B.Sの裏側で起きている演算内容であり、以下がその結果だ。2列目と8列目の各行で複素数は、見ての通り実数軸を挟んで上下対称になっているので、足し合わせると結局以下の様に実数成分しか残らない。
Fig1.3-13 再現出力 R’ (複素数表示)
御覧の通り、正しく入力信号Rが再現されている。
それではSin波の場合はどうだろうか?実数表示では以下の通りだ。
Fig1.3-14 入力R Sin波
B.Sの演算結果は、以下基底行列Bの上に重ねたSの通りで、2列目は右45度回転、8列目は左45度回転を示している。その結果Bの各行で、2列目は右に45度回転しSの振幅分増やされ、8列目は逆に左に45度回転させられ、やはりSの振幅分増やされてR’の要素として出力される。この様に2列目と8列目は±45度回転が加わったことにより位相がずらされたが、実数軸を挟んでの上下に対称性は保たれており、その結果、実数成分のみの残ってR’として出力されている。
Fig1.3-15 B.S演算の内側 Sin波
R’の結果は以下の通り、正しく入力参照信号RのSin波が再現された。
Fig1.3-16 再現出力 R’ (複素数表示)
固有値を加えた場合
続いてFig1.2-6の図形変形と同じことを周波数領域で試してみよう。DFTで構成した変換行列Wに固有値を設定し、抽出された周波数成分Sのレベルを操作する。それによりどのような変換が行われるかを見てみる。抽出された周波数成分に手を加えるというこは、言い換えるとフィルターの係数を設定するという事だ。
具体的にはDFT の基底サイズは32×32としRを矩形波として、その応答を見てみる。
固有値(フィルターの係数)を以下の様に左端[0Hz~+3Hz]=1とし、右端 [-3Hz~-1Hz]=1とする(本来フーリエ変換で周波数領域を図示する場合は、0Hz中心にする方が見やすいのだが…)。これは3Hz以下の成分を素通し、それ以上を阻止する低域通過フィルターと言えよう。
Fig1.3-17 固有値Λ(低域通過フィルター)
5サンプル時間分の矩形波を入力信号としてRに設定。
Fig1.3-18 入力信号R
以下の変換W(=B.Λ.B-1)を実施する。
R’=B.Λ.B-1.R
結果は以下の様に矩形波の入力信号Rがなまって出力された。固有値Λに低域通過型のフィルターを設定したので、高域の成分が制限されてためだ。尚、DFTの場合は基底サイズで繰り返す周期波を想定しているので、右端まで行くと左の頭につながっていると見てほしい。
Fig1.3-19 変換出力信号R’
もしこれが音声信号なら、きっともやもやして聞き取りにくい音になっているだろう。
これが固有値の働きの視覚化であり、図形の変形と同様、歪具合を決定(あるいは表現)していると言える。もう少し大きな目で見ると対象の持つ『特性』と言えるのかもしれない。
以上で固有値に関するおさらいは一旦終わりにする。続きは仮称『複素平面とラプラス変換』としよう。
付録稿
本稿使用のMathematicaプログラムの解説
本稿のおまけ章節
参考文献
文献1『素数からゼータへ、そしてカオスへ』 小山信也著 日本評論社
文献2『[特講]工業数学 ようこそフーリエ・ラプラスの世界へ!』 前田裕 工学社
本サイト:AcademiaNotes/計算機代数ノート/信号処理の基礎数学をベースに参考書として再編したものです。本稿の内容に興味があり、フーリエ変換やラプラス変換などさらに詳しく知りたい方は是非ご購入をご検討ください。
コメント