最小2乗法 - その9

どうでもいいけどサーボモータやギヤボックスや配線があらわになった亜里栖に、かえってエロスを感じるというのはどういうことだろう。

すっかり忘れていた。新しい会社でいっぱいコードを書かなければいけなくなって、必要になった計算をgslの中から探していた。FFTや最小2乗法なんかが欲しいんだけど、そういえばなんか書いてなかったっけ？自分で書いた昔のコードをあさっていたら、どっちも自分で書いたのが出てきた。デバグはすんでObjective-Cのインターフェイスになっている。ちっとも知らなかった。

いやいや、それって自分で自分のために書いたもので、忘れてるだけ。もう他人が書いたコードと区別がつかないけど、読み返してみると結構いろんなところで親切じゃん、と思った。そりゃそうだわ。だって自分が使いやすいように書いたんだもの。

そういえば最小2乗法の話も、もうちょっとだけ残っていることを改めて知った。これは一番最近で、でも4ヶ月もあいだがあいてしまうと、何やってたのかきれいさっぱり。自分の鋭い忘却力に驚いたけど、せっかくなので思い出すことにした。

ということで気を取り直して、途中でぷっつり切れていた記事の続きを書くことにした。最小2乗法についてのこれまでの記事を読み返して、なるほどこういう方向に話を持っていきたかったのね、では後は僕が引き受けましょう．．．．．

前回はHouseholder変換という、ちょうどベクトルを鏡にうつしたような変換を考えた。Householder変換は任意の二つのベクトルの内積の値を変えない、という性質があるのでQR変換に使える。前回はHouseholder変換を使って行列の第1列の第1要素以外を0にする方法までをおさらいしてあった。今回はそのあと、残りの左下半分の要素を0にする方法を考えることから続けることにする。

3.5  2列目以降のHouseholder行列

2列目以降も前回と同じようにすればいいんだけど、g′₂はg′₂₁とg′₂₂のふたつの要素以外を0にするu₂を作らないといけない。つまり

となるようにしないといけない。なんとかの一つ覚えで前回の式-66に従って同じようにu₂を作ってしまうと、g′₁がまたあさっての方向を向いてしまって、せっかくのg₁の第1要素以外がまた0でなくなってしまう。

どうすればいいかというとu₂を、g′₁と直交する超平面内に選べばいい。3次元で言えばg′₁がミラーの面内に寝ているようにすれば、g′₁はミラーによる反射はない（法線との内積が0）、ということである。

具体的には

として、

とすると

u₂の第1要素は0になって、g′₁は影響を受けない（u₂との内積は0になる）、ということになる。

これを3次元の場合に絵に描いたのが図-15である。

この図では第1列ベクトルg₁はx軸、つまり第1行要素だけになっている。赤いベクトルである第2列ベクトルg₂の第3要素を0にしたベクトルが青のg₂で、x−y平面内にある。そしてu₂はy−z平面内（図ではz成分の値は負）、つまり第1要素は0となっている。図中の赤い点線は第2列ベクトルg₂の原点の位置をずらしてベクトルの先端が原点に来るように描いたものである。つまりこれが入射光線のベクトルで、法線がu₂であるミラーによって反射されて青のg₂になっている、というとイメージがわきやすくないだろうか？光学屋以外には関係ないか。

このとき、式-72には符号の選び方の任意性が残る。これは前回のu₁とu₁₋との関係と同じで、u₂のノルムが大きなるほうを選ぶのが安全である。

残りの列ベクトルに対しても同じこと、つまり第j列ベクトルは第1からj−1要素までを同じにして第j要素を式-72と同様に決めてやることで最終的に上三角行列に変換できる、ということになる。

3.6  2列目以降の別のみかた

結果に影響は無いんだけど、ちょっとだけ簡単になる考え方がある。

式-75のu₂を使うと、式-62で行ベクトルを作るとすべての行ベクトルの第1要素は変化しない、ということがわかる。それ以降もu₃では第2要素も変わらないし、あとずっと同じである。

つまり、2列目以降は行列G′全体ではなく、小行列G′₁

を考えて、それに対して前回書いた式-68

を考えればいい、ということになる。こうすると後ろの列ベクトルの計算ではどんどん行列が小さくなる、ということになる。これのおかげでどのくらい計算量が減るか、というとそれほどたいしたことではないんだけど1列目2列目、と区別する必要がなく、実装にとってはは一つのパターンだけを考えればいいので簡単にできる（もちろん実装のしかたに依存するので、かならず簡単になるというわけではない）。