lassoの話からだんだん離れて言ってるような気がしますが・・・(汗
QR 分解が求められた後どうしたらいいんでしょう?
例えば gsl では gsl_linalg_QR_Rsolve という関数があって gsl_linalg_QR_decomp で求めた QR 分解を元に方程式を解きます。
行列 A の QR 分解を Q, R として
は
となります。ここで なので
を解くことになります。 gsl_linalg_QR_Rsolve はこの部分の計算を行います。
実際はこの関数内で (cblas の) dtrsv という上三角行列による方程式を後退代入で解く関数を呼び出します。
今 A の行数を m、列数を n として
m >= n の場合 Q は m x m の正方行列、R は m x n の縦長の行列になります。
しかし R は n x n の正方行列の下に m - n x n のゼロ行列を加えた形になるので、冗長部分を除いた
Q1 = Q(1 : m , 1 : n)
R1 = R(1 : n , 1 : n)
でも Q1 * R1 = A が成り立ちます。ここで R は n x n の正方上三角行列になります。
これは octave で
[ q , r ] = qr (A , 'econ' )
とした場合に相当し、この R1 を dtrsv に渡して
が求められます。
lars、elastic net では常に m > n ですが、逆に m < n の場合 Q は m x m、R は m x n の横長の行列になります。
この場合 R は m x m の上三角の正方行列の右に m x n - m の密行列を加えたものになるので dtrsv には渡せません。
この場合は A の転置の QR 分解を用います。
とすれば方程式
は
即ち
となります。 は A の LQ 分解に相当します。
Q は n x n、 R は n x m の縦長になるので m >= n の場合と同様冗長部分を除くと
Q2' = Q' (1 : m , 1 : n)
R2' = R' (1 : m , 1 : m)
と R2' は m x m の正方下三角行列になります。
ただ今度の場合は
なので dtrsv で得られる解は Q2' * x となり、それに後から Q を掛けて x が求められます。
以上から、qrinsert、qrdelete で Q、R を逐次更新し dtrsv で等角ベクトルを求める、とすれば計算時間はかなり短縮されそうです。
ところで lapack には dtrsv に相当する dtrtrs という関数がありますが blas の dtrsv の方が速いです。というのも、dtrtrs内では dtrsv (dtrsm) を呼び出しますが、その前にR が特異で無いか、つまり R の対角成分が 0 で無いか調べているからです。




