HaskellでQR分解を実装する

QR分解は与えられた\(m\times n\)行列\(A\)を\(m\times m\)のユニタリ行列（実数の場合、直交行列）\(Q\)と\(m\times n\)の上三角行列\(R\)の積、すなわち\(A=QR\)と分解する手法です。 数値的に安定な計算アルゴリズムが知られており、固有値の計算（QR法）やカルマンフィルターの安定的な計算にも応用されています。またこういった分解はより抽象的な対象で考えられることも多く、QR分解は半単純リー群の岩澤分解に一般化されることが知られています。

QR分解を実現するアルゴリズムはWikipediaにも詳しく載っており

• グラム・シュミットの正規直交化法
• ハウスホルダー変換
• ギブンス回転

を利用した手法などがあります。

Haskellでも例えば hmatrix が qr というQR分解を行う関数を提供していたり、hmatrixを使ったギブンス回転やハウスホルダー変換によるQR分解の実装を解説した記事もあります（お気楽 Haskell プログラミング入門 線形代数編）。しかし本稿ではあえて vector-sized を使って自分で実装してみようと思い、数値的にも安定しているハウスホルダー変換を利用した実行列のQR分解の実装したいと思います。

即席線形代数

まずは Haskellで実装する即席線形代数 を参考に実装に必要なベクトルと行列の型と関数の定義を行います。

import GHC.TypeLits
import Text.Printf

import Data.Vector.Sized (Vector)
import qualified Data.Vector.Sized as V

type Matrix m n a = Vector m (Vector n a)

-- | ベクトルのスカラー倍
(*^) :: Num a => a -> Vector n a -> Vector n a
(*^) a = V.map (*a)

-- | ベクトルをスカラー値で割る
(^/) :: Fractional a => Vector n a -> a -> Vector n a
(^/) v a = recip a *^ v

-- | 内積
dot :: Num a => Vector n a -> Vector n a -> a
dot = (V.sum .) . V.zipWith (*)

-- | 外積
outer :: Num a => Vector m a -> Vector n a -> Matrix m n a
outer xs ys = V.map (\x -> V.map (*x) ys) xs

-- | ユークリッドノルム
norm2V :: Floating a => Vector n a -> a
norm2V = sqrt . V.sum . V.map (^2)

-- | リストから行列を作成する
fromList :: (KnownNat m, KnownNat n) => [[a]] -> Maybe (Matrix m n a)
fromList = (=<<) V.fromList . mapM V.fromList

-- | 行列を整形して表示する
displayM :: PrintfArg a
         => Int  -- 数値の表示幅
         -> Int  -- 有効数字
         -> Matrix n m a
         -> IO ()
displayM w p = putStrLn . drop 1 . V.foldl (\x v -> x ++ '\n' : V.foldl (++) "" (V.map (printf "%*.*f" w p) v)) ""

-- | 単位行列
identity :: (KnownNat n, Num a) => Matrix n n a
identity = V.generate (\x -> V.generate (\y -> if x == y then 1 else 0))

-- | 行列のスカラー倍
(*!!) :: Num a => a -> Matrix m n a -> Matrix m n a
(*!!) a = V.map (V.map (*a))

-- | 行列の転置
transpose :: KnownNat n => Matrix m n a -> Matrix n m a
transpose = sequenceA

-- | 行列積
(!*!) :: (KnownNat r, Num a) => Matrix m n a -> Matrix n r a -> Matrix m r a
a !*! b = fmap (flip fmap (transpose b) . dot) a

ハウスホルダー変換

ハウスホルダー変換は与えられたベクトル\(x\)を単位法線ベクトル\(v\)で表された原典を通る超平面で鏡映変換する変換です。変換後のベクトルは \(x - 2 v \langle v, x \rangle\) と表すことができ、これは行列 \(I - 2vv^{\rm T}\)を\(x\)に左から掛けて変換していると考えることもできます。このハウスホルダー変換を使えば、与えられた行列の列ベクトルを左から順番に第n成分までの部分空間に射影していくことでQR分解を得ることができます。

アルゴリズムの詳しい解説は他の記事に譲るとして（例えばWikipedia）、さっそく実装を見ていきたいと思います。以下 householder として実装するのは添字\(i\)とベクトル\(x\)が与えられた時に、\(x\)の第\(i\)成分以降を第\(i\)成分までの部分空間に射影する（すなわち残りの成分を0にする）ハウスホルダー変換を表す行列を計算する関数です。

import Data.Maybe (fromJust)

import Data.Finite (Finite)
import qualified Data.Vector as V'

-- | ハウスホルダー変換
householder :: (KnownNat n, Ord a, Floating a) => Finite n -> Vector n a -> Matrix n n a
householder i' x =
  let i = fromIntegral i'
      y = V'.drop i $ V.fromSized x
      u = y V'.// [(0, V'.head y - y `V.withSized` norm2V)]
      padding = (V'.++) (V'.replicate i 0)
      u_norm = u `V.withSized` norm2V
      v = fromJust . V.toSized . padding $ V'.map (/u_norm) u
   in if abs u_norm < 1e-12 then identity else identity - 2 *!! outer v v

ベクトルと行列の型にはサイズに関する情報を持たせていますが householder では最初からその情報を捨てて素の Data.Vector で変換を行っています。理由としてはハウスホルダー変換を計算するベクトルの長さ（すなわち型）は第一引数である i' の値に依存しており、今のHaskellの依存型だと今回の様な状況では簡潔に実装できる方法がないため型からサイズの情報を削ることにしました。

実装中に単位法線ベクトル\(v\)を計算するために法線ベクトル\(u\)をそのノルムで割る処理がありますが、\(u\)のノルムが非常に小さい場合この処理は不安定になります。しかし\(u\)のノルムが非常に小さいということは\(x\)と変換後のベクトルがほぼ等しいという状況を表しており、このような場合には結果となる変換行列をただの単位行列にするようにしています。

QR分解

QR分解は与えられた行列の列ベクトルを左から順番にハウスホルダー変換して上三角行列を作ることにより得ることができます。

{-# LANGUAGE ScopedTypeVariables #-}

import Data.Proxy

import Data.Finite (finite)

qr :: forall m n a. (KnownNat m, KnownNat n, Ord a, Floating a) => Matrix m n a -> (Matrix m m a, Matrix m n a)
qr a =
  transpose <$>
    foldl (\(q, r) i ->
      let p = householder (finite i) (V.index r (finite i))
       in (q !*! p, r !*! p)
    ) (identity, transpose a) [0..k-1]
  where k = fromInteger $ min (natVal (Proxy @n)) (natVal (Proxy @m))

実装上の都合で行列は行ベクトルのベクトルとなっているので、列ベクトルを扱うために最初に転置を行い\(A^{\rm T}\)、得られた\(R^{\rm T}\)を最後にもう一度転置することにより計算しています。\(Q\)については本来転置したものが計算結果になるのであえて転置をしていません。

数値実験

それでは実装した qr によって実際に行列のQR分解ができるか実験してみましょう。

まずはWikipediaに載っている例を元に実験してみます。

{-# LANGUAGE DataKinds #-}

do
  let x = fromJust $ fromList
            [ [12, -51,   4]
            , [ 6, 167, -68]
            , [-4,  24, -41]]
            :: Matrix 3 3 Double
      (q, r) = qr x
  putStrLn "Q = "
  displayM 8 3 q
  putStrLn "R = "
  displayM 8 3 r
  putStrLn "QR = "
  displayM 8 3 $ q !*! r
  putStrLn "Q^TQ ="
  displayM 8 3 $ transpose q !*! q

Q = 
   0.857  -0.394  -0.331
   0.429   0.903   0.034
  -0.286   0.171  -0.943
R = 
  14.000  21.000 -14.000
   0.000 175.000 -70.000
   0.000   0.000  35.000
QR = 
  12.000 -51.000   4.000
   6.000 167.000 -68.000
  -4.000  24.000 -41.000
Q^TQ =
   1.000   0.000   0.000
   0.000   1.000   0.000
   0.000   0.000   1.000

WolframAlphaでも同様の計算を行った結果と比べてみても値が一致していることが分かります。

次に非正則行列の場合を見てみましょう。先程の例の行ベクトルと列ベクトルを一つずつ0に変えたような行列を使って実験してみます。

do
  let x = fromJust $ fromList
            [ [ 0,   0,   0]
            , [ 6, 167,   0]
            , [-4,  24,   0]]
            :: Matrix 3 3 Double
      (q, r) = qr x
  putStrLn "Q = "
  displayM 8 3 q
  putStrLn "R = "
  displayM 8 3 r
  putStrLn "QR = "
  displayM 8 3 $ q !*! r
  putStrLn "Q^TQ ="
  displayM 8 3 $ transpose q !*! q

Q = 
  -0.000  -0.000   1.000
   0.832   0.555   0.000
  -0.555   0.832   0.000
R = 
   7.211 125.640   0.000
  -0.000 112.604   0.000
  -0.000   0.000   0.000
QR = 
  -0.000  -0.000   0.000
   6.000 167.000   0.000
  -4.000  24.000   0.000
Q^TQ =
   1.000   0.000  -0.000
   0.000   1.000   0.000
  -0.000   0.000   1.000

問題なさそうですね。

次に非正方行列の場合を見てみましょう。

do
  let x = fromJust $ fromList
            [ [12, -51]
            , [ 6, 167]
            , [-4,  24]]
            :: Matrix 3 2 Double
      (q, r) = qr x
  putStrLn "Q = "
  displayM 8 3 q
  putStrLn "R = "
  displayM 8 3 r
  putStrLn "QR = "
  displayM 8 3 $ q !*! r
  putStrLn "Q^TQ ="
  displayM 8 3 $ transpose q !*! q

Q = 
   0.857  -0.394   0.331
   0.429   0.903  -0.034
  -0.286   0.171   0.943
R = 
  14.000  21.000
   0.000 175.000
  -0.000   0.000
QR = 
  12.000 -51.000
   6.000 167.000
  -4.000  24.000
Q^TQ =
   1.000   0.000  -0.000
   0.000   1.000  -0.000
  -0.000  -0.000   1.000

do
  let x = fromJust $ fromList
            [ [12, -51,   4]
            , [ 6, 167, -68]]
            :: Matrix 2 3 Double
      (q, r) = qr x
  putStrLn "Q = "
  displayM 8 3 q
  putStrLn "R = "
  displayM 8 3 r
  putStrLn "QR = "
  displayM 8 3 $ q !*! r
  putStrLn "Q^TQ ="
  displayM 8 3 $ transpose q !*! q

Q = 
   0.894  -0.447
   0.447   0.894
R = 
  13.416  29.069 -26.833
  -0.000 172.177 -62.610
QR = 
  12.000 -51.000   4.000
   6.000 167.000 -68.000
Q^TQ =
   1.000   0.000
   0.000   1.000

行より列が多い場合でも列より行が多い場合でも問題なく計算できています。