お気楽 Haskell プログラミング入門

●LU 分解

連立一次方程式 \(Ax = b\) を解くとき、あらかじめ係数行列 \(A\) を下三角行列 \(L\) (lower) と上三角行列 \(U\) (Upper) に分解しておくと、簡単に解くことができます。\(Ax = b\) は \(LUx = b\) と表すことができるので、\(Ux = y\) とおくと \(Ly = b\) を解くことで \(y\) を求めることができます。次に、求めた \(y\) を使って \(Ux = y\) を解けば、\(Ax = b\) の解 \(x\) を求めることができます。簡単な例を示しましょう。

  連立方程式

  x +  y +  z = 10
 2x + 4y + 6z = 38
 2x      + 4z = 14

 [ 1.  1.  1.       [ 1.  0.  0.      [1.  1.  1. 
   2.  4.  6.    =    2.  1.  0.   *   0.  2.  4. 
   2.  0.  4.]        2. -1.  1. ]     0.  0.  6. ]

               (1) LU 分解
  
 [ 1.  0.  0.  10.     y1 = 10
   2.  1.  0.  38.     y2 = 38 - 2 * 10 = 18 
   2. -1.  1.  14. ]   y3 = 14 - 2 * 10 + 18 = 12

               (2) 前進代入

 [ 1.  1.  1.  10.     x1 = 10 - x2 - x3 = 3
   0.  2.  4.  18.     x2 = (18 - 4 * x3) / 2 = 5
   0.  0.  6.  12. ]   x3 = 12 / 6 = 2

               (3) 後退代入

LU 分解はガウスの消去法と同じアルゴリズムで求めることができます。\(U\) は前進消去で得た行列で、\(L\) の要素は前進消去するとき係数を 0 にするための倍率になります。実際に \(L \times U\) を計算すると元の行列に戻ります。\(Ly = b\) を解く場合、\(L\) は下三角行列なので前進代入で \(y\) を求めることができます。次に \(Ux = y\) を解く場合、\(U\) は上三角行列なので後退代入で \(x\) を求めることができます。

ガウスの消去法の計算量は \(O(N^3)\) ですが、前進代入と後進代入の計算量は \(O(N^2)\) です。同じ係数行列を何度も使う場合、あらかじめ係数行列を LU 分解しておけば、ガウスの消去法よりも効率的に連立方程式を解くことができるわけです。

hmatrix には LU 分解を行う関数 lu があります。もうひとつ luPacked という関数があり、それを用いて連立方程式を解く関数 luSolve が用意されています。実用的にはそれらを使えばいいのですが、Haskell とアルゴリズムの学習ということで、あえてプログラムを作ってみましょう。

●プログラムの作成

LU 分解のプログラムは次のようになります。

リスト : LU 分解

-- 上三角行列を作る
makeUp :: [Vector R] -> Int -> [Vector R]
makeUp [v] _ = [v]
makeUp (v:vs) n =
  v : makeUp (map (\w -> let temp = (w ! n) / (v ! n)
                         in w - (v * scalar temp))
                  vs)
             (n + 1)

-- 下三角行列を作る
makeLow :: [Vector R] -> Int -> [[R]]
makeLow [_] _ = [[1]]
makeLow (v:vs) n =
  let
    a = map (\w -> let temp = (w ! n) / (v ! n)
                   in (w - (v * scalar temp), temp))
        vs
    b = makeLow (map fst a) (n + 1)
  in
    (1 : map snd a) : map (0:) b

-- LU 分解
makeLU :: Matrix R -> (Matrix R, Matrix R)
makeLU xs =
  let
    a = toRows xs
    b = makeUp a 0
    c = makeLow a 0
  in
    (fromColumns $ map vector c, fromRows b)

関数 makeLU は引数の係数行列 xs を LU 分解します。実際の処理は関数 makeUp と makeLow で行います。makeUp は上三角行列を作ります。makeLow は下三角行列を作ります。どちらの関数もアルゴリズムはガウスの消去法と同じですが、makeLow は係数を 0 にするための倍率 temp をリストに格納して返すところが異なります。

具体的には、makeLow の第 1 引数のリストの要素が一つになったら [[1]] を返します。そうでなければ、リストを v と vs に分解し、vs に格納されたベクトルの n 番目の要素を 0 にします。これはガウスの消去法の前進消去と同じですが、ベクトルと一緒に倍率 temp もタプルに格納して返すところに注意してください。そして、map で第 1 要素 (ベクトル) を取り出して makeLow に渡して再帰呼び出しします。

次に、makeLow の返り値 (リスト) の要素の先頭に 0 を追加して、そこに倍率 temp を格納したリスト (map snd a) を追加します。このとき、先頭に 1 を追加することをお忘れなく。あとは、makeLU で makeUp と makeLow を呼び出して、返り値 (リスト) を行列に変換するだけです。

なお、makeUp と makeLow はひとつの関数にまとめて処理することもできます。ただし、ピボット選択を行うと処理が複雑になるので、今回は 2 つにわけてプログラムしました。

それでは、実際に試してみましょう。

ghci> a = matrix 2 [1,1,2,4]
ghci> a
(2><2)
 [ 1.0, 1.0
 , 2.0, 4.0 ]
ghci> makeLU a
((2><2)
 [ 1.0, 0.0
 , 2.0, 1.0 ],(2><2)
 [ 1.0, 1.0
 , 0.0, 2.0 ])
ghci> b = matrix 3 [1,1,1,2,4,6,2,0,4]
ghci> b
(3><3)
 [ 1.0, 1.0, 1.0
 , 2.0, 4.0, 6.0
 , 2.0, 0.0, 4.0 ]
ghci> makeLU b
((3><3)
 [ 1.0,  0.0, 0.0
 , 2.0,  1.0, 0.0
 , 2.0, -1.0, 1.0 ],(3><3)
 [ 1.0, 1.0, 1.0
 , 0.0, 2.0, 4.0
 , 0.0, 0.0, 6.0 ])
ghci> c = matrix 4 [1,1,1,1,-1,1,-1,1,8,4,2,1,-8,4,-2,1]
ghci> c
(4><4)
 [  1.0, 1.0,  1.0, 1.0
 , -1.0, 1.0, -1.0, 1.0
 ,  8.0, 4.0,  2.0, 1.0
 , -8.0, 4.0, -2.0, 1.0 ]
ghci> makeLU c
((4><4)
 [  1.0,  0.0,  0.0, 0.0
 , -1.0,  1.0,  0.0, 0.0
 ,  8.0, -2.0,  1.0, 0.0
 , -8.0,  6.0, -1.0, 1.0 ],(4><4)
 [ 1.0, 1.0,  1.0,  1.0
 , 0.0, 2.0,  0.0,  2.0
 , 0.0, 0.0, -6.0, -3.0
 , 0.0, 0.0,  0.0, -6.0 ])
ghci> d = matrix 5 [1,-1,1,-1,1,12,-6,2,0,0,1,1,1,1,1,12,6,2,0,0,4,3,2,1,0]
ghci> d
(5><5)
 [  1.0, -1.0, 1.0, -1.0, 1.0
 , 12.0, -6.0, 2.0,  0.0, 0.0
 ,  1.0,  1.0, 1.0,  1.0, 1.0
 , 12.0,  6.0, 2.0,  0.0, 0.0
 ,  4.0,  3.0, 2.0,  1.0, 0.0 ]
ghci> makeLU d
((5><5)
 [  1.0,                0.0,               0.0,                 0.0, 0.0
 , 12.0,                1.0,               0.0,                 0.0, 0.0
 ,  1.0, 0.3333333333333333,               1.0,                 0.0, 0.0
 , 12.0,                3.0, 6.000000000000001,                 1.0, 0.0
 ,  4.0, 1.1666666666666667, 2.900000000000001, 0.26666666666666655, 1.0 ],(5><5)
 [ 1.0, -1.0,               1.0,                   -1.0,                    1.0
 , 0.0,  6.0,             -10.0,                   12.0,                  -12.0
 , 0.0,  0.0, 3.333333333333333,                   -2.0,                    4.0
 , 0.0,  0.0,               0.0,    -11.999999999999998, -3.552713678800501e-15
 , 0.0,  0.0,               0.0, -4.440892098500626e-16,    -1.6000000000000023 ])
ghci> disp 3 $ fst $ makeLU d
5x5
 1.000  0.000  0.000  0.000  0.000
12.000  1.000  0.000  0.000  0.000
 1.000  0.333  1.000  0.000  0.000
12.000  3.000  6.000  1.000  0.000
 4.000  1.167  2.900  0.267  1.000
ghci> disp 3 $ snd $ makeLU d
5x5
1.000  -1.000    1.000   -1.000    1.000
0.000   6.000  -10.000   12.000  -12.000
0.000   0.000    3.333   -2.000    4.000
0.000   0.000    0.000  -12.000   -0.000
0.000   0.000    0.000   -0.000   -1.600

次は LU 分解を使って連立方程式を解く関数 luSolver を作ります。

リスト : LU 分解による連立方程式の解法

-- 前進代入
lu1 :: [Vector R] -> Int -> [R] -> [R]
lu1 [] _ zs =  zs
lu1 (v:vs) 0 [] =
  lu1 vs 1 [(v ! (size v - 1)) / (v ! 0)]
lu1 (v:vs) n xs =
  let
    x = ((v ! (size v - 1)) - (subVector 0 n v <.> vector xs)) / (v ! n)
  in
    lu1 vs (n + 1) (xs ++ [x])

-- 後退代入
lu2 :: [Vector R] -> Int -> [R]
lu2 [v] n = [(v ! (n + 1)) / (v ! n)]
lu2 (v:vs) n =
  let m = n + 1
      x = lu2 vs m
      y = subVector m (size v - 1 - m) v <.> vector x
      z = ((v ! (size v - 1)) - y)  / (v ! n)
  in z : x

luSolver :: (Matrix R, Matrix R) -> Vector R -> [R]
luSolver (ls, us) ys =
  let
    a = ls ||| asColumn ys
    b = lu1 (toRows a) 0 []
    c = us ||| asColumn (vector b)
  in
    lu2 (toRows c) 0

関数 lu1 が前進代入を、lu2 が後進代入を行います。lu1 の第 1 引数が拡大係数行列を行単位で分解したリスト、第 2 引数が求める変数のインデックス、第 3 引数が変数値を格納するリストです。第 1 引数が空リストになったら引数 zs をそのまま返します。第 2 引数が 0 の場合、ベクトル v の先頭要素と末尾の値から変数の値を求め、リストに格納して lu1 を再帰呼び出しします。

それ以外の場合、リストを v : vs に分解し、ベクトル v の n 番目の変数の値を求めます。v の先頭から n 個の要素を取り出し、それと変数の値を格納したリスト xs の内積を計算します。これで変数に値を代入して計算することができます。その値を v の末尾要素から引き算し、v ! n の値で割り算すれば変数の値を求めることができます。あとは xs の末尾に値を追加して lu1 を再帰呼び出しします。

後進代入はガウスの消去法で作成した関数 gauss2 をほとんど同じです。luSolver は上下の三角行列をタプル (ls, us) で受け取ります。ls と引数 ys で拡大係数行列を生成し、lu1 を呼び出して前進代入を行います。次に、その結果 (変数 b) と us で拡大係数行列を生成し、lu2 を呼び出して後進代入を行います。

それでは実際に試してみましょう。

ghci> a
(2><2)
 [ 1.0, 1.0
 , 2.0, 4.0 ]
ghci> luSolver (makeLU a) $ vector [100, 272]
[64.0,36.0]
ghci> b
(3><3)
 [ 1.0, 1.0, 1.0
 , 2.0, 4.0, 6.0
 , 2.0, 0.0, 4.0 ]
ghci> luSolver (makeLU b) $ vector [10, 38, 14]
[3.0,5.0,2.0]
ghci> c
(4><4)
 [  1.0, 1.0,  1.0, 1.0
 , -1.0, 1.0, -1.0, 1.0
 ,  8.0, 4.0,  2.0, 1.0
 , -8.0, 4.0, -2.0, 1.0 ]
ghci> luSolver (makeLU c) $ vector [-5,-7,-31,-35]
[0.0,-9.0,1.0,3.0]
ghci> d
(5><5)
 [  1.0, -1.0, 1.0, -1.0, 1.0
 , 12.0, -6.0, 2.0,  0.0, 0.0
 ,  1.0,  1.0, 1.0,  1.0, 1.0
 , 12.0,  6.0, 2.0,  0.0, 0.0
 ,  4.0,  3.0, 2.0,  1.0, 0.0 ]
ghci> luSolver (makeLU d) $ vector [1,0,8,0,1]
[0.3124999999999991,0.0,-1.874999999999997,3.5,6.062499999999997]

正常に動作していますね。

●ピボット選択を使った LU 分解

次は LU 分解にピボット選択を適用してみましょう。プログラムは次のようになります。

リスト : LU 分解 (ピボット選択バージョン)

selectPivot' :: [Vector R] -> Int -> [Int] -> ([Vector R], [Int])
selectPivot' zs n idx =
  let i = findMaxRow zs n
  in if i == 0
     then (zs, idx)
     else (swapRow zs i 0, swapRow idx (n + i) n)

makeUp' :: [Vector R] -> Int -> [Int] -> ([Vector R], [Int])
makeUp' [v] _ idx = ([v], idx)
makeUp' (v:vs) n idx =
  let
    a = map (\w -> let temp = (w ! n) / (v ! n)
                   in w - (v * scalar temp))
            vs
    (b, d) = selectPivot' a (n + 1) idx
    (xs, ys) = makeUp' b (n + 1) d
  in
    (v : xs, ys)

transfer :: [a] -> [Int] -> [a]
transfer vs xs = map (\k -> vs !! k) xs

-- LU 分解
makeLU' :: Matrix R -> (Matrix R, Matrix R, [Int])
makeLU' xs =
  let
    a = toRows xs
    (b, d) = selectPivot' a 0 [0..((rows xs) - 1)]
    (ys, d1) = makeUp' b 0 d
    zs = makeLow (transfer a d1) 0
  in
    (fromColumns $ map fromList zs, fromRows ys, d1)

関数 selectPivot' の引数 idx に交換した行の情報をセットします。係数行列の大きさを n とすると、idx の初期値は 0 から n - 1 までの整数を格納したリストになります。selectPivot' は引数 zs の 0 番目と最大値の行 (i 番目) を交換します。このとき、idx は n 番目と n + i 番目の値を交換することに注意してください。関数 makeUp' はピボット選択しながら上三角行列を生成します。返り値は上三角行列を表すリストと行の交換を記録したリストです。

関数 makeLU' は selectPivot' で 0 行目と最大値の行を交換してから makeUp' を呼び出します。変数 ys が上三角行列を表すリスト、変数 d1 が行の交換位置を記録したリストになります。そして、関数 transfer で係数行列 a をリスト d1 の順序で並べ替えてから makeLow を呼び出します。makeLow でピボット選択する必要はありません。あとは、リストを行列に変換して d1 と一緒にタプルに格納して返します。

それでは実際に試してみましょう。

ghci> a
(2><2)
 [ 1.0, 1.0
 , 2.0, 4.0 ]
ghci> (x, y, z) = makeLU' a
ghci> x
(2><2)
 [ 1.0, 0.0
 , 0.5, 1.0 ]
ghci> y
(2><2)
 [ 2.0,  4.0
 , 0.0, -1.0 ]
ghci> z
[1,0]
ghci> x <> y
(2><2)
 [ 2.0, 4.0
 , 1.0, 1.0 ]

ghci> b
(3><3)
 [ 1.0, 1.0, 1.0
 , 2.0, 4.0, 6.0
 , 2.0, 0.0, 4.0 ]
ghci> (x, y, z) = makeLU' b
ghci> x
(3><3)
 [ 1.0,  0.0, 0.0
 , 1.0,  1.0, 0.0
 , 0.5, 0.25, 1.0 ]
ghci> y
(3><3)
 [ 2.0,  4.0,  6.0
 , 0.0, -4.0, -2.0
 , 0.0,  0.0, -1.5 ]
ghci> z
[1,2,0]
ghci> x <> y
(3><3)
 [ 2.0, 4.0, 6.0
 , 2.0, 0.0, 4.0
 , 1.0, 1.0, 1.0 ]

ghci> c
(4><4)
 [  1.0, 1.0,  1.0, 1.0
 , -1.0, 1.0, -1.0, 1.0
 ,  8.0, 4.0,  2.0, 1.0
 , -8.0, 4.0, -2.0, 1.0 ]
ghci> (x, y, z) = makeLU' c
ghci> x
(4><4)
 [    1.0,     0.0,  0.0, 0.0
 ,   -1.0,     1.0,  0.0, 0.0
 ,  0.125, 6.25e-2,  1.0, 0.0
 , -0.125,  0.1875, -1.0, 1.0 ]
ghci> y
(4><4)
 [ 8.0, 4.0,  2.0,  1.0
 , 0.0, 8.0,  0.0,  2.0
 , 0.0, 0.0, 0.75, 0.75
 , 0.0, 0.0,  0.0,  1.5 ]
ghci> z
[2,3,0,1]
ghci> x <> y
(4><4)
 [  8.0, 4.0,  2.0, 1.0
 , -8.0, 4.0, -2.0, 1.0
 ,  1.0, 1.0,  1.0, 1.0
 , -1.0, 1.0, -1.0, 1.0 ]

ghci> d
(5><5)
 [  1.0, -1.0, 1.0, -1.0, 1.0
 , 12.0, -6.0, 2.0,  0.0, 0.0
 ,  1.0,  1.0, 1.0,  1.0, 1.0
 , 12.0,  6.0, 2.0,  0.0, 0.0
 ,  4.0,  3.0, 2.0,  1.0, 0.0 ]
ghci> (x, y, z) = makeLU' d
ghci> x
(5><5)
 [                  1.0,                    0.0,   0.0,                  0.0, 0.0
 ,                  1.0,                    1.0,   0.0,                  0.0, 0.0
 ,   0.3333333333333333,     0.4166666666666667,   1.0,                  0.0, 0.0
 , 8.333333333333333e-2, -4.1666666666666664e-2, 0.625,                  1.0, 0.0
 , 8.333333333333333e-2,                  0.125, 0.625, -0.23076923076923078, 1.0 ]
ghci> y
(5><5)
 [ 12.0, -6.0,                     2.0,    0.0,                0.0
 ,  0.0, 12.0,                     0.0,    0.0,                0.0
 ,  0.0,  0.0,      1.3333333333333335,    1.0,                0.0
 ,  0.0,  0.0, -1.1102230246251565e-16, -1.625,                1.0
 ,  0.0,  0.0, -1.3664283380001927e-16,    0.0, 1.2307692307692308 ]
ghci> z
[1,3,4,0,2]
ghci> x <> y
(5><5)
 [ 12.0, -6.0, 2.0,  0.0, 0.0
 , 12.0,  6.0, 2.0,  0.0, 0.0
 ,  4.0,  3.0, 2.0,  1.0, 0.0
 ,  1.0, -1.0, 1.0, -1.0, 1.0
 ,  1.0,  1.0, 1.0,  1.0, 1.0 ]

ghci> e = matrix 3 [0,2,4,1,1,1,4,2,6]
ghci> e
(3><3)
 [ 0.0, 2.0, 4.0
 , 1.0, 1.0, 1.0
 , 4.0, 2.0, 6.0 ]
ghci> (x, y, z) = makeLU' e
ghci> x
(3><3)
 [  1.0,  0.0, 0.0
 ,  0.0,  1.0, 0.0
 , 0.25, 0.25, 1.0 ]
ghci> y
(3><3)
 [ 4.0, 2.0,  6.0
 , 0.0, 2.0,  4.0
 , 0.0, 0.0, -1.5 ]
ghci> z
[2,0,1]
ghci> x <> y
(3><3)
 [ 4.0, 2.0, 6.0
 , 0.0, 2.0, 4.0
 , 1.0, 1.0, 1.0 ]

正常に動作していますね。連立方程式を解く関数 luSolver' は次のようになります。

リスト : LU 分解による連立方程式の解法 (ピボット選択)

luSolver' :: (Matrix R, Matrix R, [Int]) -> Vector R -> [R]
luSolver' (ls, us, idx) ys =
  let
    a = ls ||| (col $ transfer (toList ys) idx)
    b = lu1 (toRows a) 0 []
    c = us ||| col b
  in
    lu2 (toRows c) 0

ベクトル ys をリストに変換し、transfer で idx の順序に並べ替えます。あとは luSolver と同じです。簡単な実行例を示します。

ghci> a
(2><2)
 [ 1.0, 1.0
 , 2.0, 4.0 ]
ghci> luSolver' (makeLU' a) $ vector [100, 272]
[64.0,36.0]

ghci> b
(3><3)
 [ 1.0, 1.0, 1.0
 , 2.0, 4.0, 6.0
 , 2.0, 0.0, 4.0 ]
ghci> luSolver' (makeLU' b) $ vector [10, 38, 14]
[3.0,5.0,2.0]

ghci> c
(4><4)
 [  1.0, 1.0,  1.0, 1.0
 , -1.0, 1.0, -1.0, 1.0
 ,  8.0, 4.0,  2.0, 1.0
 , -8.0, 4.0, -2.0, 1.0 ]
ghci> luSolver' (makeLU' c) $ vector [-5,-7,-31,-35]
[0.0,-9.0,1.0,3.0]

ghci> d
(5><5)
 [  1.0, -1.0, 1.0, -1.0, 1.0
 , 12.0, -6.0, 2.0,  0.0, 0.0
 ,  1.0,  1.0, 1.0,  1.0, 1.0
 , 12.0,  6.0, 2.0,  0.0, 0.0
 ,  4.0,  3.0, 2.0,  1.0, 0.0 ]
ghci> luSolver' (makeLU' d) $ vector [1,0,8,0,1]
[0.31249999999999994,0.0,-1.8749999999999998,3.5,6.0625]

ghci> e
(3><3)
 [ 0.0, 2.0, 4.0
 , 1.0, 1.0, 1.0
 , 4.0, 2.0, 6.0 ]
ghci> luSolver' (makeLU' e) $ vector [14,10,38]
[5.0,3.0,2.0]

●hmatrix の LU 分解

ここで、Haskell (hmatrix) に用意されている LU 分解の関数について簡単に説明しておきましょう。lu A は行列 A を LU 分解します。

lu :: Field t => Matrix t -> (Matrix t, Matrix t, Matrix t, t)

返り値とタプル (l, u, p, s) を照合した場合、変数 l に下三角行列、u に上三角行列、p にピボット選択を表す転置行列、s に p のシグネチャ ^[*1] がセットされます。次の例を見てください。

ghci> a = matrix 2 [1,1,2,4]
ghci> a
(2><2)
 [ 1.0, 1.0
 , 2.0, 4.0 ]
ghci> (l,u,p,s) = lu a
ghci> l
(2><2)
 [ 1.0, 0.0
 , 0.5, 1.0 ]
ghci> u
(2><2)
 [ 2.0,  4.0
 , 0.0, -1.0 ]
ghci> p
(2><2)
 [ 0.0, 1.0
 , 1.0, 0.0 ]
ghci> s
-1.0
ghci> l <> u
(2><2)
 [ 2.0, 4.0
 , 1.0, 1.0 ]
ghci> p <> l <> u
(2><2)
 [ 1.0, 1.0
 , 2.0, 4.0 ]

ghci> b = matrix 3 [1,1,1,2,4,6,2,0,4]
ghci> b
(3><3)
 [ 1.0, 1.0, 1.0
 , 2.0, 4.0, 6.0
 , 2.0, 0.0, 4.0 ]
ghci> (l,u,p,s) = lu b
ghci> l
(3><3)
 [ 1.0,  0.0, 0.0
 , 1.0,  1.0, 0.0
 , 0.5, 0.25, 1.0 ]
ghci> u
(3><3)
 [ 2.0,  4.0,  6.0
 , 0.0, -4.0, -2.0
 , 0.0,  0.0, -1.5 ]
ghci> p
(3><3)
 [ 0.0, 0.0, 1.0
 , 1.0, 0.0, 0.0
 , 0.0, 1.0, 0.0 ]
ghci> s
1.0
ghci> l <> u
(3><3)
 [ 2.0, 4.0, 6.0
 , 2.0, 0.0, 4.0
 , 1.0, 1.0, 1.0 ]
ghci> p <> l <> u
(3><3)
 [ 1.0, 1.0, 1.0
 , 2.0, 4.0, 6.0
 , 2.0, 0.0, 4.0 ]

行列 a を LU 分解して変数 l, u, p, s にセットします。置換行列 p を左から掛け算すると、行を交換する動作になります。したがって、p <> l <> u を計算すると元の行列 a に戻ります。同様に、行列 b を LU 分解して p <> l <> u を計算すると b に戻ります。

関数 luSolve は LU 分解を使って連立方程式を解きます。この場合、関数 lu ではなく関数 luPacked を使って LU 分解した値を渡します。

luSolve :: Field t => LU t -> Matrix t -> Matrix t
luPacked :: Field t => Matrix t -> LU t

簡単な実行例を示します。

ghci> a
(2><2)
 [ 1.0, 1.0
 , 2.0, 4.0 ]
ghci> luPacked a
LU (2><2)
 [ 2.0,  4.0
 , 0.5, -1.0 ] [1,1]
ghci> luSolve (luPacked a) $ col [100, 272]
(2><1)
 [ 64.0
 , 36.0 ]

ghci> b
(3><3)
 [ 1.0, 1.0, 1.0
 , 2.0, 4.0, 6.0
 , 2.0, 0.0, 4.0 ]
ghci> luPacked b
LU (3><3)
 [ 2.0,  4.0,  6.0
 , 1.0, -4.0, -2.0
 , 0.5, 0.25, -1.5 ] [1,2,2]
ghci> luSolve (luPacked b) $ col [10, 38, 14]
(3><1)
 [ 3.0
 , 5.0
 , 2.0 ]

詳しい説明は hmatrix のマニュアル Numeric.LinearAlgebra をお読みください。

-- note --------
[*1] 置換行列 p のシグネチャは p の行列式のことです。置換行列は直交行列なので、行列式の値は 1 または -1 になります。行列 A が PLU と表される場合、A の行列式 |A| は |P||L||U| になります。三角行列の行列式は対角成分を乗算した値になるので、A を LU 分解すれば行列式を簡単に求めることができます。

●参考文献・URL

奥村晴彦,『Ｃ言語による最新アルゴリズム事典』, 技術評論社, 1991
ガウスの消去法 - Wikipedia

●プログラムリスト

--
-- lu.hs : LU 分解
--
--         Copyright (C) 2021 Makoto Hiroi
--
import Prelude hiding ((<>))
import Numeric.LinearAlgebra

-- 上三角行列
makeUp :: [Vector R] -> Int -> [Vector R]
makeUp [v] _ = [v]
makeUp (v:vs) n =
  v : makeUp (map (\w -> let temp = (w ! n) / (v ! n)
                         in w - (v * scalar temp))
                  vs)
             (n + 1)

-- 下三角行列
makeLow :: [Vector R] -> Int -> [[R]]
makeLow [_] _ = [[1]]
makeLow (v:vs) n =
  let
    a = map (\w -> let temp = (w ! n) / (v ! n)
                   in (w - (v * scalar temp), temp))
        vs
    b = makeLow (map fst a) (n + 1)
  in
    (1 : map snd a) : map (0:) b

-- LU 分解
makeLU :: Matrix R -> (Matrix R, Matrix R)
makeLU xs =
  let
    a = toRows xs
    b = makeUp a 0
    c = makeLow a 0
  in
    (fromColumns $ map vector c, fromRows b)

-- 前進代入
lu1 :: [Vector R] -> Int -> [R] -> [R]
lu1 [] _ zs =  zs
lu1 (v:vs) 0 [] =
  lu1 vs 1 [(v ! (size v - 1)) / (v ! 0)]
lu1 (v:vs) n xs =
  let
    x = ((v ! (size v - 1)) - (subVector 0 n v <.> vector xs)) / (v ! n)
  in
    lu1 vs (n + 1) (xs ++ [x])

-- 後退代入
lu2 :: [Vector R] -> Int -> [R]
lu2 [v] n = [(v ! (n + 1)) / (v ! n)]
lu2 (v:vs) n =
  let m = n + 1
      x = lu2 vs m
      y = subVector m (size v - 1 - m) v <.> vector x
      z = ((v ! (size v - 1)) - y)  / (v ! n)
  in z : x

-- 連立方程式の解法
luSolver :: (Matrix R, Matrix R) -> Vector R -> [R]
luSolver (ls, us) ys =
  let
    a = ls ||| asColumn ys
    b = lu1 (toRows a) 0 []
    c = us ||| asColumn (vector b)
  in
    lu2 (toRows c) 0

--
-- ピボット選択
--
findMaxRow :: [Vector R] -> Int -> Int
findMaxRow [] _ = 0
findMaxRow (z:zs) n =
  iter zs 1 (abs (z ! n)) 0
  where
    iter [] _ _ b = b
    iter (z:zs) i a b =
      let c = abs (z ! n)
      in if a < c
         then iter zs (i + 1) c i
         else iter zs (i + 1) a b

swapRow :: [a] -> Int -> Int -> [a]
swapRow zs m n =
  iter zs 0
  where
    x = zs !! m
    y = zs !! n
    iter [] _ = []
    iter (z:zs) i
      | i == m    = y : iter zs (i + 1)
      | i == n    = x : iter zs (i + 1)
      | otherwise = z : iter zs (i + 1)

selectPivot' :: [Vector R] -> Int -> [Int] -> ([Vector R], [Int])
selectPivot' zs n idx =
  let i = findMaxRow zs n
  in if i == 0
     then (zs, idx)
     else (swapRow zs i 0, swapRow idx (n + i) n)

transfer :: [a] -> [Int] -> [a]
transfer vs xs = map (\k -> vs !! k) xs

-- 上三角行列 (ピボット選択)
makeUp' :: [Vector R] -> Int -> [Int] -> ([Vector R], [Int])
makeUp' [v] _ idx = ([v], idx)
makeUp' (v:vs) n idx =
  let
    a = map (\w -> let temp = (w ! n) / (v ! n)
                   in w - (v * scalar temp))
            vs
    (b, d) = selectPivot' a (n + 1) idx
    (xs, ys) = makeUp' b (n + 1) d
  in
    (v : xs, ys)

-- LU 分解 (ピボット選択)
makeLU' :: Matrix R -> (Matrix R, Matrix R, [Int])
makeLU' xs =
  let
    a = toRows xs
    (b, d) = selectPivot' a 0 [0..((rows xs) - 1)]
    (ys, d1) = makeUp' b 0 d
    zs = makeLow (transfer a d1) 0
  in
    (fromColumns $ map fromList zs, fromRows ys, d1)

-- 連立方程式の解法 (ピボット選択)
luSolver' :: (Matrix R, Matrix R, [Int]) -> Vector R -> [R]
luSolver' (ls, us, idx) ys =
  let
    a = ls ||| (col $ transfer (toList ys) idx)
    b = lu1 (toRows a) 0 []
    c = us ||| col b
  in
    lu2 (toRows c) 0