最小2乗法 - その7

しばらく間があいてしまって何をやってたのかまた例によって忘れそうになる。前回までは行列GがQR分解できたとしたら、という話だった。ここからはどうやってQR分解するか、という話に入る。

7  QR分解の実行

QR分解にはいろいろなやり方があって、やりかたによって見た目異なる結果が得られることがあるが、列ベクトルの相対関係が変わらないと言う意味では同じである。

考え方としては、もとのm ×n行列Gの列ベクトルを回転させたり鏡に映したり座標軸を入れ替えたりして、上三角行列にすればいい。上三角行列Rって

こんなだから、幾何学的なイメージを考えてみる。列ベクトルの空間R^mの軸にそった単位ベクトルをe_iとする。最初の列ベクトルは回転させたり鏡に映したり座標軸を入れ替えたりして、e₀の向きにあわせる、ということになる。その次のベクトルは最初の列ベクトルを変えずにこんどはe₀とe₁がはる２次元平面内にはいるようにする。その次はe₀とe₁とe₂がはる3次元超平面内に入るようにする。これを最後の列ベクトルが残るところまでやればいい。

めんどくさそうだけど、3次元で考えれば簡単にわかる。

こんな3 ×3の行列G_eを考えてみる。なんだか恣意的な数字が並んだ行列だけど、これはなんでもいい。後の操作がわかりやすいように選んだのと、列ベクトルが適度に独立で長さが1になるようにした。この列ベクトルを絵に描いてみたのが図-10である。

立体的に見たいので、公差法の裸眼立体視で描いた。列ベクトルの順に緑青赤である。

これをQR分解しよう。まず第1例ベクトルのe₃成分が0になるように、つまりe₁-e₂平面内に入るようにe₂を軸に（つまりe₁-e₃平面内で）回転させる。とりあえず今はベクトルの成分からヒューリスティックにπ/6だけ向かって右回りに（図中の矢印のように）回転させる。そうすると図-11のようになる。

こうすると行列は

となって、第1列のe₃要素は0になる。

さらに今度は図-12のようにe₃を軸に（つまりe₁-e₂平面内で）−π/6だけ回転させると、第1列ベクトルの方向はe₁と一致して

となる。

また、さらにこんどは図-13のように2番目の列ベクトルのe₃成分を0にするためにe₁を軸に（つまりe₂-e₃平面内で）−π/3だけ回転させると、

となって、めでたく上三角行列になった。回転操作を整理するために、e_i-e_j平面内でθだけ回転させる行列をΦ_i,j(θ)とすると、

だから、

と分解できたと言うことになる。

これがいわゆるGivens回転によるQR分解のエッセンスである。

Givens 回転は、ベクトル幾何なんかで使う回転の行列を多次元に拡張したもので、n次元のGivens 回転の行列Φ_i,j(θ)は

というようなものである（都合によりWikipediaと符号が違うが、回転方向が逆なだけで同じである）。3次元の場合でもいちいち回転軸をいうついでに○○平面内で、と断っていたのは、3次元以上の場合には回転軸の指定よりも、「どの面内で」回転させるかを指定した方が簡単（n次元での軸を指定しようとするとn−2本の軸を指定する必要がある）で、Givens回転はその拡張と見やすいと思ったからである。

i,j面内のGivens回転によって、どの列ベクトルもi番目要素とj番目要素以外は変わらないことは実際に計算してみればわかる。列ベクトルeのj番目の要素e_jを0にするようなGivens回転Φ_i,j(θ)は

となるようなθを持っていればいい（複号は、お互い異なる符号を持っていればいいと言う意味）。これも実際にかけ算してみればすぐわかる。

このGivens回転を使って、まず第1列の第1要素e_1,1以外を全部0にできる。次に第2列の第3要素までを同じように0にできるが、第2列の第2要素e_2,2を0にしようとすると、第1列の第2要素e_1,2が0でなくなってしまうのでこれはできない。これをずっと繰り返せば上三角行列になる、と言う次第である。