QR factorization

戻る topへ

ハウスホルダー(Householder)変換を用いたQR 分解(factorization) 

 

概要:QR分解とは,行列A (n×mm, rank(A) = m)を直交行列Q(n×n)と上三角行列R(n×m)の積に分解する手法である.連立方程式を解くときなどに活躍する.ここでは,グラムシュミットの直交化を用いたQR分解とHouseholder変換を用いたQR分解のアルゴリズムをまとめる

図1. QR分解

 

グラムシュミットの直交化を用いたQR分解

n×m行列A (m),rank(A) = mは,n×n直交行列Qn×m上三角行列Rの積,A=QRと分解できる(図1).

 

Ai列を列ベクトルaiで表す, A=( aa2 … am).rank(A) = mなので,グラムシュミットの直交化により,m本の直交ベクトルbiが取得できる.


上式を次のように整理すると,

 


 

次の関係は明らかである.


 


 

ここでBを直交化したいので,BXの間に,対角行列D=diag(||b1||, ||b2||, …, ||bm||), D-1=diag(1/||b1||, 1/||b2||, …, 1/||bm||),を挟む

 

n×m行列Q = (b1/||b1||, b2/||b2||, … , bm/||bm||)の列ベクトルは,n次元空間のm本の正規直行規定をなし,Rは上三角行列である(図2左).Qを直交行列にしたいので,bi/||bi||に直交する(n – m)本の単位ベクトルを適当に選び,これを並べた各列に並べた行列をQ1とする.すると,

となり目的のQR分解を得られる.(Q Q1)はn×n直交行列である(図2右).

図2, QR分解

 

グラムシュミットの直交化法を用いるとQR分解が計算できる.しかし,グラムシュミットの直交化法は,逐次的に,既存のベクトルに直交する成分を引き延ばして利用する(正規化する)ため,計算機を用いて計算する場合には,数値的に不安定になる.これから説明するHouseholder変換を用いた方が安定といわれている.

 

 

Householder 変換を用いたQR分解

 

Householder変換とは,長さ等しい2本のn次元ベクトルxyが与えられた時(||x||=||y||),次のn×n行列Hで定義される線形変換のことである.

Inn×n単位行列である.この変換は, xyに,yをに移す変換である.

 

よってHx=yであり, Hy=xも同様に示せる.上の変形において xTxyTy,yTxxTyを用いた.また,Householder変換行列は,直交行列である.

 

 

Householder変換を用いたQR分解. 上記のHouseholder変換を用いて,QR分解を行う方法を考える.n×m行列A(n>m)が与えられ,k-1個のn×n直交行列H1, H2,…,Hk-1により既にk-1列目まで上三角行列に分解されているとする;

ここで「上式に左乗すると,k-1列目までの上三角性を失わず,k列目のk+1~n成分までをゼロにする.」という性質を満たすHouseholder変換Hkを求める.

 

天下り的だが,u =( 0 u),  0∊Rk-1,  uk ∊Rn-k+1 というベクトルuが作るHouseholder行列は, 

となり,これを式(1)の左辺に左乗する.


上式より,k-1列目までの上三角性は失われない事が分かる.さらに,k列目のk+1~n行の要素をゼロにしたい.つまり以下を満たす,ukを取得したい.


上式のは,n-k+1次元のHouseholder変換であるので,ukを以下のように定義すればよい.


上の±はどちらの符号を選んでも良いが,xの第1成分と逆の符号を選んだ方が,桁落ちが少ない(u=x-yなので,同符号だと引き算による桁落ちの可能性がある).

 

このukにより定義されるhouseholder変換

を式(1)の両辺に左乗すると

k列目まで対角化できる.これを,m回繰り返すと,

 

と上三角化できる.Householder変換の積 は直交行列なので,H-1=HTを両辺に左乗すると

とQR分解が得られる(Q=HT).

 

 

Householder変換によるQR分解アルゴリズム

 

+A-1)n×m行列A, n>m, rank(A)=mが入力される.

+A-2)準備

    n×m行列Rを用意し RA

    n×n行列Qを用意し QI(用途によっては必要ない)

+ A-3)以下をm回繰り返す

    + A-3-1) k番目のHouseholder変換を次のように求める.

    A-3-2) RR  (HmHm-1…H1Aを求めている)

      + A-3-3) Q←Q  (H=Hm Hm-1…H1を求めている.用途によっては必要ない)

+ A-4) A = QTRが得られる

 

高速化のための注意.

1) 上記アルゴリズムにおいて,Qは計算する必要のない事が多い(例えばAx=bの線形システムを解くときなどはQ= Hm Hm-1…H1の代わりにHm Hm-1…H1bを求めればよい.).

2) u uTH (u∊RnH∊Rn×n)という積の計算が必要になるが,これを(u uT) Hと行列どうしの積にするとO(n3)かかるので, u (uTH)と計算する必要がある.後者だとO(n2)になる.

 

 

 

参考文献

[1] D.A.ハーヴィル (著), 伊理 正夫 (監訳), 統計のための行列代数 上.

 

井尻敬 @ 生物情報基盤, 理研 

勉強した日 2012/4/12