今回は、前回までの関数を利用して逆行列を行うプログラムを作成した。

  ソースコード

  追加した関数は次の 3 つ。

inverseMatrix :: [[Int]] -> [[Ratio Int]]
inverseMatrix mx = let factors = allFactor mx
                       detA = determinant mx
                   in if detA == 0 
                      then []
                      else cutArray (length mx) $ map (\x -> (1 % (fromIntegral detA)) * (x%1)) 
                       $ checkSign factors $ map determinant $ map (cofactor mx) factors
 where
  checkSign :: [(Int, Int)] -> [Int] -> [Int]
  checkSign fs xs | null xs = []
                  | True = let low = fst (head fs)
                               col = snd (head fs)
                               x = head xs
                           in if mod (low + col) 2 > 0 then (-1*x):checkSign (tail fs) (tail xs)
                                                       else x:checkSign (tail fs) (tail xs)

  -- 配列から行列への変換
  cutArray :: Int -> [a] -> [[a]]
  cutArray len xs | null xs = []
                  | True     = [(take len xs)] ++ (cutArray len $ reverse $ take ((length xs) - len) $ reverse xs)

-- 余因子行列を求める
cofactor :: [[Int]] -> (Int, Int) -> [[Int]]
cofactor mx (low, col) = map (removeColumn col) $ removeColumn low mx
 where
  removeColumn :: Int -> [a] -> [a]
  removeColumn num xs | num > length xs = xs
                      | True = (take (num - 1) xs) ++ (reverse $ take ((length xs) - num) $ reverse xs)

-- 行列の各要素のインデックスを求める
allFactor :: [a] -> [(Int, Int)]
allFactor mx = do x <- allFactor' 1 mx
                  y <- allFactor' 1 mx
                  [(x, y)]
 where
  allFactor' :: Int -> [a] -> [Int]
  allFactor' cnt mx | length mx == 0  = []
                    | True = cnt:(allFactor' (cnt+1) (tail mx))

  allFactor が行列の各要素のインデックスを求めている。行列式の各要素の符号を求める際に利用する。
  cofactor は入力された行列から余因子行列を求める関数になる。
  そして inverseMatrix が逆行列を求める関数になる。ここでは、cofactor 関数を呼び出し、逆行列の要素の元になる余因子行列を取得し、その逆行列を求める。そして、各符号を設定し、各要素に対して元の行列の行列式の逆数を係数に設定する。
 
　これに対して、次の行列を入れてみる。

  [[1,2,1],
   [0,-2,3],
   [2,1,5]]

  結果は次のとおり。

　[[(-2)/1,0/1,(-3)/1],
   [0/1,(-1)/1,(-3)/1],
   [(-1)/1,(-1)/1,(-4)/1]]

　これまでに作ってきた関数はいづれも線形代数の基本的な処理であるが。情報工学の技術では多くの場面において、これらが強烈なツールとなる。
　
　え、ライブラリがあるって？
　自分で作ったツールで出した結果の方が、絶対感動するだろうし。遊びでやっているんだから、車輪の再開発も大いに結構ではないか。

　次回は、これらを使った何か面白い計算をやってみる。