Haskellで実装する即席線形代数

プログラミングをしている時に行列積などのちょっとした線形代数の演算をしたいけど環境を整えるのが大変だしな… といった状況は稀にあるかと思います。本稿では即席線形代数と称して、そうした時にコピペで使える線形代数演算の実装例をまとめていきたいと思います。

もちろんHaskellでベクトルや行列を扱いたいのであれば

• hmatrix: BLAS/LAPACK の Haskell バインディング
• linear: Haskell による純粋な線形代数実装

といったライブラリの導入を検討するのが先決でしょう。ただこれらのライブラリを使用するには BLAS/LAPACK が先に必要だったり依存するライブラリが大きかったりして、何らかの理由でこれらを避けたい状況もあるでしょう。

Haskellでベクトルのようなデータ構造を扱うにはいくつか選択肢があり、本稿ではそれらのデータ構造による実装をそれぞれ掲載していきます。

リスト

まずはリストを使ってベクトルと行列を表現した場合の実装についてまとめていきます。Haskellのリストは連結リストの構造を持っているのでランダムな要素のアクセスには $O(n)$ の時間がかかるという欠点があります。一方でリストは関数型プログラミングにおける最も基本的なデータ型の一つであり、洗練されたコンビネータが豊富に用意されているという利点もあります。本稿ではこの利点を活かし、できるだけ簡潔な方法による実装を載せたいと思います。

型の定義

ベクトルは単なるリスト、行列は行ベクトルのリストとして定義します。

type Vector a = [a]
type Matrix a = [Vector a]

ベクトル

ベクトルは単なるリストなので値はリテラルや標準の関数を使って構築することができます。例えば定数ベクトルを作る場合は単純に replicate を使えば良いでしょう。

基本的な加減乗除は以下のように実装することができます（演算子は linear を参考）。

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

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

-- | ベクトルとベクトルの足し算
(^+^) :: Num a => Vector a -> Vector a -> Vector a
(^+^) = zipWith (+)

-- | ベクトルとベクトルの引き算
(^-^) :: Num a => Vector a -> Vector a -> Vector a
(^-^) = zipWith (-)

ベクトル同士の積としては要素同士の積 zipWith (*) も考えられますが、線形代数としては内積と外積が出てくることが多いでしょう。

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

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

ベクトルのノルムも以下のように実装することができます。ノルムによって必要な制約が変わってくるのは面白いですね。

-- | ベクトルの絶対値ノルム
norm1V :: Num a => Vector a -> a
norm1V = sum . map abs

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

行列

行列も単なる二重リストなのでリテラルや標準の関数を使って構築することができます。特殊な行列を作成するための関数だけ用意しておきましょう。

-- | 単位行列
identity :: Num a => Int -> Matrix a
identity n = map (\x -> map (\y -> if x == y then 1 else 0) [1..n]) [1..n]

また小数の行列はそのまま出力すると非常に見にくいので整形して表示する関数も定義しておきましょう。

import Data.List (intercalate)
import Text.Printf

-- | 行列を整形して表示する
displayM :: PrintfArg a
         => Int  -- 数値の表示幅
         -> Int  -- 有効数字
         -> Matrix a
         -> IO ()
displayM w p = putStrLn . intercalate "\n" . map (concatMap (printf "%*.*f" w p))

基本的な加減乗除は以下のように実装することができます。

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

-- | 各成分ごとの演算
elementWise :: (a -> b -> c) -> Matrix a -> Matrix b -> Matrix c
elementWise = zipWith . zipWith

-- | 行列の和
(!+!) :: Num a => Matrix a -> Matrix a -> Matrix a
(!+!) = elementWise (+)

-- | 行列の差
(!-!) :: Num a => Matrix a -> Matrix a -> Matrix a
(!-!) = elementWise (-)

-- | アダマール積
(!!*!!) :: Num a => Matrix a -> Matrix a -> Matrix a
(!!*!!) = elementWise (*)

行列とベクトルあるいは行列同士の積は以下のように実装できます。

import Data.List (transpose)

-- | 行列とベクトルの積
(!*) :: Num a => Matrix a -> Vector a -> Vector a
(!*) m v = map (dot v) m

-- | ベクトルと行列の積
(*!) :: Num a => Vector a -> Matrix a -> Vector a
(*!) v = map (dot v) . transpose

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

転置行列の計算には標準で定義されている transpose 関数をそのまま使うことができます。

線形代数においてこれらの積は次元が合っている場合にのみ定義されますが、以上の実装では次元が異なる場合にも実行されてしまいます。この問題は型に次元の情報を持たせることにより解決することができますが、それは別のデータ型で実装する際に実現することにしましょう。

最後に行列のトレースやノルムを定義しておきましょう。

-- | 行列の対角成分
diag :: Matrix a -> Vector a
diag m = zipWith (!!) m [0..]

-- | 行列のトレース
trace :: Num a => Matrix a -> a
trace = sum . diag

-- | 行列の絶対値ノルム
norm1M :: Num a => Matrix a -> a
norm1M = norm1V . concat

-- | 行列のフロベニウスノルム
norm2M :: Floating a => Matrix a -> a
norm2M = norm2V . concat

行列を行ベクトルを並べたベクトルに変換する関数として concat を用いることができます。

Vector

Vector は vector ライブラリで提供されているデータ型です。実体は配列でありランダムアクセスを $O(1)$ で行うことができます。Haskellには array という配列を扱うライブラリもありますが、vector は添字を Int に限定することで効率的な実装を伴った豊富なコンビネータを提供しているのが特徴です。さらに Vector 型に要素数の情報を持たせた型を提供する vector-sized というライブラリがあります。ここでは型で次元を管理して本来なら定義されない演算が実行されるのを防ぎながら、豊富コンビネータを活用することでリストとほとんど変わらない簡潔な方法で実装できることを見ていきたいと思います。

型の定義

ベクトルとしては vector-sized が提供する Vector 型をそのまま用います。行列は行ベクトルのベクトルとして定義します。

import Data.Vector.Sized (Vector)

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

vector-sized には Data.Vector.Unboxed.Sized というモジュールもあり、こちらのモジュールが提供する Vector 型の方が格納する値がボックス化されず効率が良いです。格納する値が Unbox のインスタンスである必要がありますが、線形代数の文脈であれば数値型を扱うことが多いので問題ないでしょう。Unbox な値を持つ Vector もまた Unbox のインスタンスになるので、こちらのモジュールを使う場合でも値が Unboxed のインスタンスであるという型クラス制約をつければ以下の実装でも基本的に動作すると思います。

ベクトル

ベクトルを作るには fromList を使ってリストから変換するのが基本的でしょう。fromList の返り値は Maybe に包まれているのでリストの長さが明らかな場合は fromJust で取り出すことも多いでしょう。定数ベクトルを作るなら replicate を使ったり、generate を使って [0..n-1] を元に生成するといった方法も便利です。

基本的な加減乗除はリストと同様に以下のように実装することができます。

import qualified Data.Vector.Sized as V

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

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

-- | ベクトルとベクトルの足し算
(^+^) :: Num a => Vector n a -> Vector n a -> Vector n a
(^+^) = V.zipWith (+)

-- | ベクトルとベクトルの引き算
(^-^) :: Num a => Vector n a -> Vector n a -> Vector n a
(^-^) = V.zipWith (-)

内積と外積もリストと同様に以下のように実装することができます。

-- | 内積
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

ベクトルのノルムも以下のように実装することができます。

{-# LANGUAGE DataKinds #-}

import GHC.TypeLits

-- | ベクトルの絶対値ノルム
norm1V :: Num a => Vector n a -> a
norm1V = V.sum . V.map abs

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

行列

行列をリストのリテラルから生成するための関数 fromList と、単位行列を作成するための関数を実装しておきます。単位行列の実装は行列のサイズに関する情報を型から取り出す事ができるのでリストのときよりむしろ簡潔になっていますね。

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

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

行列を整形して表示する関数も用意しておきます。

import Text.Printf

-- | 行列を整形して表示する
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)) ""

基本的な加減乗除は以下のように実装することができます。

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

-- | 各成分ごとの演算
elementWise :: (a -> b -> c) -> Matrix m n a -> Matrix m n b -> Matrix m n c
elementWise = V.zipWith . V.zipWith

-- | 行列の和
(!+!) :: Num a => Matrix m n a -> Matrix m n a -> Matrix m n a
(!+!) = elementWise (+)

-- | 行列の差
(!-!) :: Num a => Matrix m n a -> Matrix m n a -> Matrix m n a
(!-!) = elementWise (-)

-- | アダマール積
(!!*!!) :: Num a => Matrix m n a -> Matrix m n a -> Matrix m n a
(!!*!!) = elementWise (*)

行列とベクトルあるいは行列同士の積は以下のように実装できます。これらの実装は型で次元が揃っていることを保証しているので先述のリストによる実装の場合に比べて安全です。

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

-- | 行列とベクトルの積
(!*) :: Num a => Matrix m n a -> Vector n a -> Vector m a
(!*) m v = V.map (dot v) m

-- | ベクトルと行列の積
(*!) :: (KnownNat n, Num a) => Vector m a -> Matrix m n a -> Vector n a
(*!) v = V.map (dot v) . transpose

-- | 行列積
(!*!) :: (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

行列の転置 transpose は sequenceA そのものです。その理由として Vector n a が Traversable であり Applicative のインスタンスが ZipList のように振る舞う ことがあります。

https://twitter.com/lotz84_/status/1536351912588869632

直感的には Vector n a は Finite n -> a というインデックスを取り要素を返す関数の型と同型になっており（つまり表現可能関手になっている）、行列の型 Matrix m n a（つまり Vector m (Vector n a) ）は Finite m -> Finite n -> a に同型です。sequenceA の型は

sequenceA :: (Traversable t, Applicative f) => t (f a) -> f (t a)

のようになっており、行列に適用すると引数の順番が Finite n -> Finite m -> a と入れ替わる形になります。すなわちこれは $ij$ 成分を $ji$ 成分に入れ替えているので、まさに転置の操作に対応していると考えることができるでしょう。一般的には関数の型は Traversable にはならず、"引数を入れ替える"概念はその双対である Distributive の distribute 関数が対応します。ですが Traversable の distribute は Distributive の sequenceA に対応する^[1]ので、今回の場合 sequenceA で実装しても変わりません。

最後に行列のトレースやノルムを定義しておきましょう。

{-# LANGUAGE NoStarIsType #-}

-- | 行列の対角成分
diag :: KnownNat n => Matrix n n a -> Vector n a
diag = (>>= id)

-- | 行列のトレース
trace :: (KnownNat n, Num a) => Matrix n n a -> a
trace = V.sum . diag

-- | 行列をベクトルに変換する
flatten :: Matrix m n a -> Vector (m * n) a
flatten = V.concatMap id

-- | 行列の絶対値ノルム
norm1M :: Num a => Matrix m n a -> a
norm1M = norm1V . flatten

-- | 行列のフロベニウスノルム
norm2M :: Floating a => Matrix m n a -> a
norm2M = norm2V . flatten

diag の実装はモナドの join そのものです。先ほど Vector の Applicative インスタンスは ZipList のように振る舞うと述べましたが、ZipList はモナドにはなりません。しかし Vector は ZipList とは異なり要素数が固定されているので対角要素を取るという操作が join として上手く振る舞うのでモナドにすることができます（さらに Vector は次数付きモナドになります）。

おわりに

本稿ではこれまでの記事中で実装した線形代数に関する演算をまとめてみました。種々の行列分解（LU分解、コレスキー分解、QR分解、SVD）や行列分解を用いて計算するのが効率の良い行列式・逆行列・ランクの計算、MVector を使った効率の良い計算などなどまだ書きたくて書けていないものは沢山ありますが、これ以上が必要であれば素直に hmatrix や linear を使うのが良いかもしれません。自分のチートシートとして使って行きたいと思っているので今後もアップデートしていくつもりです。

---

＼読んでいただきありがとうございました！／
この記事が面白かったら いいね♡ をいただけると嬉しいです☺️
バッジを贈っていただければ次の記事を書くため励みになります🙌

脚注

1. @viercc さんに証明してもらいました🙏✨