なぜ A^{T}A は逆行列を持つのか(1) Linear Algebra

A^{T}A という行列のかけ算は最小二乗法で頻出する．この行列はカラムが独立であれば: 正方行列，Symmetric, 逆行列を持つ というすばらしい特徴がある．これはなぜかということを説明してみたい．この説明は，Gilbert Strang のIntroduction of Linear Algebra, 4th edition の説明を基礎にしており，それは null space と二次形式の考えを使っていて，とてもすばらしい．しかし，単に Strang の方法を説明するだけでは芸がないので，geometric な私の説明も加えてみたい．ただし，証明できなかったので，Strang の説明に比較すると蛇足でしかないが，私には geometry 的ということもあって直感的である気がする．お楽しみ頂けたら幸いです．私の説明が駄目でも Strang のすばらしい説明があるので大丈夫でしょう．


A^{T}Aの性質


A^{T}A という行列のかけ算は頻出する．これは Projection の時に出てくるものであり，また，最小二乗法でも出てくるパターンである．(実際にはこれら2つは同じものであるが)

通常の問題では以下のような状況が仮定できる場合が多い．


- Aが正方行列ではない場合の最小二乗法．通常，m by n で m > n の場合，つまり，式が解よりも多い場合である．たとえば問題として，ある信号を観測する場合や統計をとるような場合がある．これは観測を考えると，同じ観測を何度も繰り返すことで実際のパラメータ数よりも観測数が多い場合に相当する． m < n の場合でも観測を増やすことは通常可能なのでこの形にすることができる．

- A のカラムが独立であることが期待できる．そうでない場合にはそういうカラムを捨てることができる．ただやみくもに捨てると重要なモデルのパラメータも捨ててしまい，間違った答えを導く可能性もあるので，これは単に技術的にできるという問題であってちょっと弱いかもしれない．


上記の条件のうち，二番目のカラムの独立性が満たされる場合には，A^{T}A は次のすばらしい性質を持つ．

- 正方行列
- Symmetric
- 逆行列を持つ

これは使える．

A のカラムが独立であることは必要条件である．そうでない場合にはA^{T}A が逆行列を持たないことは，

の例を見れば一目瞭然であろう．

正方行列にはるのは，[n by m] [m by n] が [n by n] になることからわかる．

Symmetric は，

からわかる．ここで，A^{{T}^{T}} = A, 転置の転置は元の行列であることを使った．問題は逆行列があるかどうかである．

Filter design (2)

Input actual data

では，前回設計したフィルタに実際に信号を入れてみよう．コンスタントな入力の場合を表 1 に示す．コンスタントな入力はつまらない例だが，まずは簡単なもので小手試しといこう．

Table 1 Constant input


表1 の赤の部分をどうやって計算したのか一応示しておこう．y_n で n=1の時なので，

これは最初の3つの信号をフィルタに通した場合である．入力が全部 1 なので，出力も全て同じになっている．さて，ここで transfer function を計算してみると，

まさか Jojo貴様! と驚く．(古すぎたか? それに英語版を書く時には Jojo の奇妙な冒険を知っている人がいるのかどうか．) Transfer function は入力と出力の比に等しくなった!

ところでちょっとしたことだが表でy_nが n=0,8 でも値があるのは cos の値が計算できるからで，実際に信号を測定する場合には，n=-1の値がないだろうから (n=0から測定を始めるとすると，その前のデータがないので)，y_0 は存在しないとしても良い．

コンスタントの入力では出力もコンスタントなのでちょっと退屈という人には次
の表2と3を示そう．

Table 2 cos π/3 input


Table 3 cos 2π/3 input




表2の赤の部分も先の場合と同様に計算すると，

この transfer function は

transfer function は入力がどれだけ transfer されるかであったが，1 というのはそのまま出力されるという意味である．表2を見ると，確かに入力と出力が一致している! この周波数はそのまま出力されるというわけだ．すばらしい．

表3の場合は，

この transfer function は

もう説明する必要はあまりないだろう．この周波数は設計どおりこのフィルタを通過しない．


まとめ

ここで使っているそれぞれの数学はそんなに複雑ではなく， Euler の式を受け入れることができれば，あとは高校生の数学である． Euler の式も e^x,sin(x), cos(x) 関数を無限級数展開して，比較することで関係があることは感じることができるのではないだろうか．

eigenvalue である transfer function とサンプリングした結果を linear combination するフィルタの結果が理想的な状況とはいえこれだけ一致するのは心地いい．

Filter design (1)

Filter design

Hamming の Digital Filter の話だが，3.9 節があまりにすばらしいので何か書いておきたい．

しばらく前に話題になっていたのば，サッカーのテレビ放送でゲームでブブゼラ(Vuvuzela)の音をは抑えたいが，歓声やゲームの音は通したいというようなことだった．この場合，どんな周波数がブブゼラの音にあるのかを分析し，それを選択的に通さないようにする．デジタルフィルタの設計というのはある周波数を通すとか通さないとかいうことになる．具体的な例としては以下の記事が面白かった．http://blogs.mathworks.com/loren/2010/06/30/vuvuzela-denoising-with-parametric-equalizers/

3.9 節では non recursive filter の単純な例が述べられる．これは

という3つのサンプルを使うフィルタだ．この場合，フィルタの設計はa,bの値を望みのフィルタに適するように求めることだ．まずは，transfer function を求めてみる．前回述べたように，transfer function はフィルタ関数のeigenvalue である．つまり，フィルタ F の入力に対して何が出てくるかを示している．ここで，信号を x (= 時間の関数 x(t)) として, フィルタ操作関数をF とすると，eigenvalue λを考えて，F x =λ x である．線形代数の方法になじんでいる人は，F を matrix A と考えて，λ をスカラと考えると，A x = λx になることに注意すると良いと思う．)

したがって，transfer function は


もし，望みのフィルタが π/3 の周波数をそのまま通し (=1)，2π/3 の周波数をカットしたい(=0)という場合，transfer functionをそのようにする．

これを解くと

が得られる．フィルタは

となった．

これだけでは何が起こったのかまだわからない．そこで次回はこれに実際の信号を入れてみて，フィルタがどう働くかを見てみよう．