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

線形代数編

●hello, hmatrix!!

hmatrix はベクトルや行列を使った数値計算を高速に行うためのパッケージ
内部では LAPACK / BLAS を呼び出している
BLAS (Basic Linear Algebra Subprograms) はベクトルと行列の基本的な計算を行うパッケージ
LAPACK (Linear Algebra PACKage) は BLAS を使って連立 1 次方程式や固有値問題などを解くパッケージ
LAPACK / BLAS はフリーで利用でき、高速で信頼性の高いライブラリとして有名
hmatrix をインストールする前に、LAPACK と BLAS のインストールが必要になる
Ubuntu 18.04 (WSL) の場合、次のコマンドで簡単にインストールすることができる

$ sudo apt install liblapack-dev liblapack-doc

パッケージの依存関係で BLAS (libblas-dev) も同時にインストールされる
hmatrix のインストールは stack を使うと簡単

$ stack update
$ stack install hmatrix

モジュール名は Numeric.LinearAlgebra

簡単な実行例

Prelude> import Prelude hiding ((<>))
Prelude> :m + Numeric.LinearAlgebra
Prelude Numeric.LinearAlgebra> v1 = vector [1,2,3]
Prelude Numeric.LinearAlgebra> :t v1
v1 :: Vector R
Prelude Numeric.LinearAlgebra> v1
[1.0,2.0,3.0]
Prelude Numeric.LinearAlgebra> v2 = vector [4,5,6]
Prelude Numeric.LinearAlgebra> v2
[4.0,5.0,6.0]
Prelude Numeric.LinearAlgebra> v1 + v2
[5.0,7.0,9.0]
Prelude Numeric.LinearAlgebra> v2 - v1
[3.0,3.0,3.0]
Prelude Numeric.LinearAlgebra> v1 <.> v2
32.0

Prelude Numeric.LinearAlgebra> m1 = (2><2) [1,2,3,4] :: Matrix R
Prelude Numeric.LinearAlgebra> m1
(2><2)
 [ 1.0, 2.0
 , 3.0, 4.0 ]
Prelude Numeric.LinearAlgebra> :t m1
m1 :: Matrix R
Prelude Numeric.LinearAlgebra> m2 = (2><2) [5,6,7,8] :: Matrix R
Prelude Numeric.LinearAlgebra> m2
(2><2)
 [ 5.0, 6.0
 , 7.0, 8.0 ]
Prelude Numeric.LinearAlgebra> m1 + m2
(2><2)
 [  6.0,  8.0
 , 10.0, 12.0 ]
Prelude Numeric.LinearAlgebra> m2 - m1
(2><2)
 [ 4.0, 4.0
 , 4.0, 4.0 ]
Prelude Numeric.LinearAlgebra> m1 <> m2
(2><2)
 [ 19.0, 22.0
 , 43.0, 50.0 ]
Prelude Numeric.LinearAlgebra> inv m1
(2><2)
 [ -1.9999999999999998,                  1.0
 ,  1.4999999999999998, -0.49999999999999994 ]
Prelude Numeric.LinearAlgebra> (inv m1) <> m1
(2><2)
 [ 1.0000000000000002, 4.440892098500626e-16
 ,                0.0,    0.9999999999999998 ]

Vector a はベクトル、Matrix a は行列を表すデータ型
R は Double の別名 (他のデータ型は Numeric.LinearAlgebra.Data を参照)
<> は Prelude で定義されているので、import hiding で Prelude の定義を隠す
もしくは、import qualified で Numeric.LinearAlgebra に別名を付ける
参考 URL

hmatrix: Numeric Linear Algebra, (Hackage)
Haskellで数値計算線形代数編, (雨粒blog, lqtmirage さん)
数値計算パッケージの使い方 (PDF), (細田陽介さん)
LAPACK - Wikipedia

●ベクトルと行列の基礎知識

大きさ (長さ) と方向を持つ量のことを「ベクトル (vector)」という
大きさだけで方向を持たない量のことを「スカラー (scalar)」という
Haskell (hmatrix) ではデータ型 Vector a がベクトルを表す
ベクトルだけの演算であれば、行ベクトルと列ベクトルの区別は不要
行列とベクトルを演算する場合、行ベクトルと列ベクトルを考慮する必要がある

たとえば、行列 A とベクトル v の乗算 Au と uA を考える
Au は A #> u と記述するが、この場合 u を列ベクトルとして扱う
uA は u <# A と記述するが、この場合 u を行ベクトルとして扱う

大きさが 0 のベクトルを「ゼロベクトル (零ベクトル)」という
大きさ 1 のベクトルを「単位ベクトル」という

「行列 (matrix)」は要素 (element) を長方形または正方形に並べた「表 (table)」みたいなもの
要素は数、式、記号など
要素のことを成分ともいう
Haskell (hmatrix) ではデータ型 Matrix a が行列を表す
行列の横方向の並びを「行 (row)」といい、縦方向の並びを「列 (column)」という
すべての成分が 0 の行列を「零行列」という
行の数と列の数が等しい行列 (n 行 n 列の行列) を「正方行列」という
正方行列 A の左上から右下の対角成分が 1 で、それ以外の成分が 0 の行列を「単位行列」という

Haskell (hmatrix) では、n 次の単位行列を関数 ident n で生成する

単位行列は E や I と表記されることが多い

(2><3)        (3><2)      (3><3)         (3><3)
[ 1, 2, 3     [ 1, 2      [ 1, 2, 3      [ 1, 0, 0
, 4, 5, 6 ]   , 3, 4      , 4, 5, 6      , 0, 1, 0 
              , 5, 6 ]    , 7, 8, 9 ]    , 0, 0, 1 ]

 2 行 3 列   3 行 2 列     3 行 3 列       単位行列

\(\begin{array}{llll} \begin{pmatrix} 1 & 2 & 3 \\ 4 & 5 & 6 \end{pmatrix} & \begin{pmatrix} 1 & 2 \\ 3 & 4 \\ 5 & 6 \end{pmatrix} & \begin{pmatrix} 1 & 2 & 3 \\ 4 & 5 & 6 \\ 7 & 8 & 9 \end{pmatrix} & \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{pmatrix} \end{array}\)

行列 A と B において、行の数 m と列の数 n が等しいとき、A と B を同じ型であるという
同じ型の行列 A と B において、対応する成分がすべて等しいとき、A と B は等しい (A = B)
1 行だけの行列を「行ベクトル」、1 列だけの行列を「列ベクトル」という

●ベクトルの基本操作

ベクトルの生成には演算子 |> や関数 vector, fromList を使うと簡単

(|>) :: Foreign.Storable.Storable a => Int -> [a] -> Vector a
vector :: [R] -> Vector R
fromList :: Foreign.Storable.Storable a => [a] -> Vector a

n |> xs はリスト xs の先頭から n 個の要素を取り出し、それを Vector に格納して返す
このほかに、整数のベクトルを生成する関数 range と idxs がある

range :: Int -> Vector I
idxs :: [Int] -> Vector I

range n は 0 から n - 1 までの整数を格納した Vector を返す
I は CInt の別名 (Int32 と同じ)
ベクトルの要素数は関数 size で求める
要素のアクセスには演算子 ! または関数 atIndex を使う
関数 toList はベクトルをリストに変換する

toList :: Foreign.Storable.Storable a => Vector a -> [a]

Prelude Numeric.LinearAlgebra> v = vector [1,2]
Prelude Numeric.LinearAlgebra> v
[1.0,2.0]
Prelude Numeric.LinearAlgebra> w = fromList [3,4] :: Vector R
Prelude Numeric.LinearAlgebra> w
[3.0,4.0]
Prelude Numeric.LinearAlgebra> 5 |> [1..]
[1.0,2.0,3.0,4.0,5.0]
Prelude Numeric.LinearAlgebra> v ! 0
1.0
Prelude Numeric.LinearAlgebra> v ! 1
2.0
Prelude Numeric.LinearAlgebra> atIndex v 0
1.0
Prelude Numeric.LinearAlgebra> atIndex v 1
2.0
Prelude Numeric.LinearAlgebra> size w
2
Prelude Numeric.LinearAlgebra> toList w
[3.0,4.0]

●ベクトルの連結と分解

ベクトルの連結と分解する主な関数を以下に示す

vjoin :: Foreign.Storable.Storable t => [Vector t] -> Vector t
subVector :: Foreign.Storable.Storable t => Int -> Int -> Vector t -> Vector t
takesV :: Foreign.Storable.Storable t => [Int] -> Vector t -> [Vector t]

vjoin は引数のリストに格納されたベクトルをすべて連結したベクトルを返す
subVector の引数は 1 番目が始点 (s), 2 番目が要素数 (n) を表す
subVector は引数のベクトルの s 番目から n 個の要素を取り出し、それを Vector に格納して返す
takesV は第 1 引数のリストで取り出す部分ベクトルの要素数を指定する
takesV は引数のベクトルの先頭から順番に部分ベクトルを取り出し、それらをリストに格納して返す

Prelude Numeric.LinearAlgebra> v1 = vector [1,2,3]
Prelude Numeric.LinearAlgebra> v2 = vector [4,5,6]
Prelude Numeric.LinearAlgebra> v3 = vector [7,8,9]
Prelude Numeric.LinearAlgebra> v4 = vjoin [v1,v2,v3]
Prelude Numeric.LinearAlgebra> v4
[1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0]
Prelude Numeric.LinearAlgebra> subVector 1 4 v4
[2.0,3.0,4.0,5.0]
Prelude Numeric.LinearAlgebra> takesV [2,3,4] v4
[[1.0,2.0],[3.0,4.0,5.0],[6.0,7.0,8.0,9.0]]

●ベクトルの計算

スカラー k とベクトル v の要素の四則演算は k op v, v op k

k の型は整数ならば Num, 実数ならば Fractional であること
スカラーの指定には関数 scalar を使うと便利
なお、k が即値であれば Haskell が型変換する

自分で型変換することもできる
関数 fromIntegral は整数を Num に変換する
関数 realToFrac は実数を Fractional に変換する

ベクトルの要素に関数を適用する cmap もある (マッピング)

cmap :: (Container c e, Element b) => (e -> b) -> c e -> c b

ベクトル v, w の要素同士の四則演算は v op w
ベクトル v, w の内積は v <.> w または dot v w
ベクトル v の大きさ (ユークリッドノルム) |v| は関数 norm_2 v で求める

Prelude Numeric.LinearAlgebra> v
[1.0,2.0]
Prelude Numeric.LinearAlgebra> w
[3.0,4.0]
Prelude Numeric.LinearAlgebra> v + 10
[11.0,12.0]
Prelude Numeric.LinearAlgebr> 10 + v
[11.0,12.0]
Prelude Numeric.LinearAlgebra> v * 2
[2.0,4.0]
Prelude Numeric.LinearAlgebra> 2 * v
[2.0,4.0]
Prelude Numeric.LinearAlgebra> v - 10
[-9.0,-8.0]
Prelude Numeric.LinearAlgebra> 10 - v
[9.0,8.0]
Prelude Numeric.LinearAlgebra> v / 2
[0.5,1.0]
Prelude Numeric.LinearAlgebra> 1 / v
[1.0,0.5]
Prelude Numeric.LinearAlgebra> cmap (*2) v
[2.0,4.0]

Prelude Numeric.LinearAlgebra> v + w
[4.0,6.0]
Prelude Numeric.LinearAlgebra> w - v
[2.0,2.0]
Prelude Numeric.LinearAlgebra> v * w
[3.0,8.0]
Prelude Numeric.LinearAlgebra> v / w
[0.3333333333333333,0.5]

Prelude Numeric.LinearAlgebra> v <.> w
11.0
Prelude Numeric.LinearAlgebra> dot v w
11.0
Prelude Numeric.LinearAlgebra> norm_2 v
2.23606797749979
Prelude Numeric.LinearAlgebra> ev = v / scalar (norm_2 v)
Prelude Numeric.LinearAlgebra> ev
[0.4472135954999579,0.8944271909999159]
Prelude Numeric.LinearAlgebra> norm_2 ev
1.0

●行列の基本操作

行列の生成は関数 (><), fromLists, matrix などを使うと簡単

(><) :: Foreign.Storable.Storable a => Int -> Int -> [a] -> Matrix a
matrix :: Int -> [R] -> Matrix R
fromLists :: Element t => [[t]] -> Matrix t

m 行 n 列の行列を (><) で作成する場合、最初の引数が m で 2 番目の引数が n になる
このとき、(m><n) と指定することができる
matrix の第 1 引数は列の数 (行の要素数) を指定する
関数 size は配列の大きさをタプルで返す
関数 rows は行数を、関数 cols は列数を返す
要素のアクセスは演算子 ! を 2 回使う、または関数 atIndex を使う

行列 M の m 行 n 列の要素は M ! m ! n, または atIndex M (m, n)
M ! m とすると、m 行目の要素をベクトルに格納して返す

関数 row は引数のリストから行ベクトルを生成する
関数 col は引数のリストから列ベクトルを生成する
関数 toLists は行列をリストに変換する
rand は一様分布乱数 [0, 1.0) で、randn はガウス分布乱数で行列を作る

rand :: Int -> Int -> IO (Matrix Double)
randn :: Int -> Int -> IO (Matrix Double)

Prelude Numeric.LinearAlgebra> (><) 2 3 [1..6]
(2><3)
 [ 1.0, 2.0, 3.0
 , 4.0, 5.0, 6.0 ]
Prelude Numeric.LinearAlgebra> (2><3) [1..6]
(2><3)
 [ 1.0, 2.0, 3.0
 , 4.0, 5.0, 6.0 ]
Prelude Numeric.LinearAlgebra> fromLists [[1,2,3],[4,5,6]]
(2><3)
 [ 1.0, 2.0, 3.0
 , 4.0, 5.0, 6.0 ]
Prelude Numeric.LinearAlgebra> a = matrix 3 [1..6]
Prelude Numeric.LinearAlgebra> a
(2><3)
 [ 1.0, 2.0, 3.0
 , 4.0, 5.0, 6.0 ]
Prelude Numeric.LinearAlgebra> b = matrix 3 [4..9]
Prelude Numeric.LinearAlgebra> b
(2><3)
 [ 4.0, 5.0, 6.0
 , 7.0, 8.0, 9.0 ]
Prelude Numeric.LinearAlgebra> size a
(2,3)
Prelude Numeric.LinearAlgebra> rows a
2
Prelude Numeric.LinearAlgebra> cols a
3

Prelude Numeric.LinearAlgebra> a ! 0
[1.0,2.0,3.0]
Prelude Numeric.LinearAlgebra> a ! 0 ! 0
1.0
Prelude Numeric.LinearAlgebra> a ! 1
[4.0,5.0,6.0]
Prelude Numeric.LinearAlgebra> a ! 1 ! 2
6.0
Prelude Numeric.LinearAlgebra> atIndex a (0,0)
1.0
Prelude Numeric.LinearAlgebra> atIndex a (1,2)
6.0
Prelude Numeric.LinearAlgebra> row [1,2,3]
(1><3)
 [ 1.0, 2.0, 3.0 ]
Prelude Numeric.LinearAlgebra> col [1,2,3]
(3><1)
 [ 1.0
 , 2.0
 , 3.0 ]
Prelude Numeric.LinearAlgebra> toLists a
[[1.0,2.0,3.0],[4.0,5.0,6.0]]
Prelude Numeric.LinearAlgebra> rand 3 3
(3><3)
 [  0.1819017982026105,    0.8054357985059898,  0.6154939842482535
 , 0.46874250353721086, 3.9034717734453606e-2, 0.27816141270015454
 ,  0.9679905604422049,   0.36524810612446074, 0.20691627413356503 ]

●行列の連結と分解

行列を連結する主な関数 (演算子) を以下に示す

(|||) :: Element t => Matrix t -> Matrix t -> Matrix t
(===) :: Element t => Matrix t -> Matrix t -> Matrix t
fromBlocks :: Element t => [[Matrix t]] -> Matrix t

A ||| B は A の右横に B を連結する
A === B は A の下に B を連結する
formBlocks は内側のリストに格納された行列を ||| で連結し、その行列を === で連結する

Prelude Numeric.LinearAlgebra> a = matrix 3 [1..6]
Prelude Numeric.LinearAlgebra> b = matrix 3 [7..12]
Prelude Numeric.LinearAlgebra> a
(2><3)
 [ 1.0, 2.0, 3.0
 , 4.0, 5.0, 6.0 ]
Prelude Numeric.LinearAlgebra> b
(2><3)
 [  7.0,  8.0,  9.0
 , 10.0, 11.0, 12.0 ]
Prelude Numeric.LinearAlgebra> a ||| b
(2><6)
 [ 1.0, 2.0, 3.0,  7.0,  8.0,  9.0
 , 4.0, 5.0, 6.0, 10.0, 11.0, 12.0 ]
Prelude Numeric.LinearAlgebra> a === b
(4><3)
 [  1.0,  2.0,  3.0
 ,  4.0,  5.0,  6.0
 ,  7.0,  8.0,  9.0
 , 10.0, 11.0, 12.0 ]
Prelude Numeric.LinearAlgebra> fromBlocks [[a, b], [matrix 3 [1..9], matrix 3 [11..19]]]
(5><6)
 [ 1.0, 2.0, 3.0,  7.0,  8.0,  9.0
 , 4.0, 5.0, 6.0, 10.0, 11.0, 12.0
 , 1.0, 2.0, 3.0, 11.0, 12.0, 13.0
 , 4.0, 5.0, 6.0, 14.0, 15.0, 16.0
 , 7.0, 8.0, 9.0, 17.0, 18.0, 19.0 ]

行列を分解する主な関数 (演算子) を以下に示す

takeRows :: Element t => Int -> Matrix t -> Matrix t
takeColumns :: Element t => Int -> Matrix t -> Matrix t
dropRows :: Element t => Int -> Matrix t -> Matrix t
dropColumns :: Element t => Int -> Matrix t -> Matrix t
subMatrix :: Element a => (Int, Int) -> (Int, Int) -> Matrix a -> Matrix a
(??) :: Element t => Matrix t -> (Extractor, Extractor) -> Matrix t

takeRows n A, A から n 行取り出す
takeColumns n A, A から n 列取り出す
dropRows n A, A から n 行取り除く
dropColumns n A, A から n 列取り除く
subMatrix (i, j) (m, n) A, 始点 (i, j) から m 行 n 列の行列を取り出す
(??) A (Extractor, Extractor), 取り出す行と列を Extractor で指定する
Extractor の指定方法

All
Range Int Int Int, 引数は start, step, end を表す
Pos (Vector I), 取り出す行 (列) を Vector I で指定する
PosCyc (Vector I)
Take Int, 先頭から n 行 (列) 取り出す
TakeLast Int, 末尾から n 行 (列) 取り出す
Drop Int, 先頭から n 行 (列) 取り除く
DropLast Int, 末尾から n 行 (列) 取り除く

Prelude Numeric.LinearAlgebra> a = matrix 5 [1..20]
Prelude Numeric.LinearAlgebra> a
(4><5)
 [  1.0,  2.0,  3.0,  4.0,  5.0
 ,  6.0,  7.0,  8.0,  9.0, 10.0
 , 11.0, 12.0, 13.0, 14.0, 15.0
 , 16.0, 17.0, 18.0, 19.0, 20.0 ]
Prelude Numeric.LinearAlgebra> takeRows 3 a
(3><5)
 [  1.0,  2.0,  3.0,  4.0,  5.0
 ,  6.0,  7.0,  8.0,  9.0, 10.0
 , 11.0, 12.0, 13.0, 14.0, 15.0 ]
Prelude Numeric.LinearAlgebra> dropRows 1 a
(3><5)
 [  6.0,  7.0,  8.0,  9.0, 10.0
 , 11.0, 12.0, 13.0, 14.0, 15.0
 , 16.0, 17.0, 18.0, 19.0, 20.0 ]
Prelude Numeric.LinearAlgebra> takeColumns 3 a
(4><3)
 [  1.0,  2.0,  3.0
 ,  6.0,  7.0,  8.0
 , 11.0, 12.0, 13.0
 , 16.0, 17.0, 18.0 ]
Prelude Numeric.LinearAlgebra> dropColumns 1 a
(4><4)
 [  2.0,  3.0,  4.0,  5.0
 ,  7.0,  8.0,  9.0, 10.0
 , 12.0, 13.0, 14.0, 15.0
 , 17.0, 18.0, 19.0, 20.0 ]
Prelude Numeric.LinearAlgebra> subMatrix (1,1) (2,2) a
(2><2)
 [  7.0,  8.0
 , 12.0, 13.0 ]
Prelude Numeric.LinearAlgebra> a ?? (Pos (idxs [3, 1]), All)
(2><5)
 [ 16.0, 17.0, 18.0, 19.0, 20.0
 ,  6.0,  7.0,  8.0,  9.0, 10.0 ]
Prelude Numeric.LinearAlgebra> a ?? (Pos (idxs [3, 1]), Range 4 (-2) 0)
(2><3)
 [ 20.0, 18.0, 16.0
 , 10.0,  8.0,  6.0 ]
Prelude Numeric.LinearAlgebra> a ?? (Take 3, Drop 2)
(3><3)
 [  3.0,  4.0,  5.0
 ,  8.0,  9.0, 10.0
 , 13.0, 14.0, 15.0 ]
Prelude Numeric.LinearAlgebra> a ?? (TakeLast 2, DropLast 2)
(2><3)
 [ 11.0, 12.0, 13.0
 , 16.0, 17.0, 18.0 ]

●行列の計算

スカラー k と行列 M の要素の四則演算はベクトルと同じく k op M, M op k
行列の要素に関数を適用する cmap もある (マッピング)

cmap :: (Container c e, Element b) => (e -> b) -> c e -> c b

行列 M, N の要素同士の四則演算は M op N

Prelude Numeric.LinearAlgebra> a
(2><3)
 [ 1.0, 2.0, 3.0
 , 4.0, 5.0, 6.0 ]
Prelude Numeric.LinearAlgebra> b
(2><3)
 [ 4.0, 5.0, 6.0
 , 7.0, 8.0, 9.0 ]
Prelude Numeric.LinearAlgebra> a * 2
(2><3)
 [ 2.0,  4.0,  6.0
 , 8.0, 10.0, 12.0 ]
Prelude Numeric.LinearAlgebra> b / 2
(2><3)
 [ 2.0, 2.5, 3.0
 , 3.5, 4.0, 4.5 ]
Prelude Numeric.LinearAlgebra> cmap (*10) a
(2><3)
 [ 10.0, 20.0, 30.0
 , 40.0, 50.0, 60.0 ]

Prelude Numeric.LinearAlgebra> a + b
(2><3)
 [  5.0,  7.0,  9.0
 , 11.0, 13.0, 15.0 ]
Prelude Numeric.LinearAlgebra> b - a
(2><3)
 [ 3.0, 3.0, 3.0
 , 3.0, 3.0, 3.0 ]
Prelude Numeric.LinearAlgebra> a * b
(2><3)
 [  4.0, 10.0, 18.0
 , 28.0, 40.0, 54.0 ]
Prelude Numeric.LinearAlgebra> a / b
(2><3)
 [               0.25,   0.4,                0.5
 , 0.5714285714285714, 0.625, 0.6666666666666666 ]

行列の積は <>, 逆行列 A^-1 は関数 inv A を使う

逆行列 inv の引数は正方行列であること

ベクトル v と行列 M の積は演算子 #> と <# を使う

M #> v は v を列ベクトルとして扱う
v <# M は v を行ベクトルとして扱う

正方行列 A に逆行列 A^-1 が存在するとき、A は「正則行列 (regular matrix)」であるという
A が正則行列であれば、連立方程式 Ax = b の解 x は一意に定まる
実際には解 x を inv(A) * b で求めるよりも、linearSolve などの関数 (ソルバー) を使ったほうがよい

linearSolve :: Field t => Matrix t -> Matrix t -> Maybe (Matrix t)

Prelude Numeric.LinearAlgebra> a
(2><3)
 [ 1.0, 2.0, 3.0
 , 4.0, 5.0, 6.0 ]
Prelude Numeric.LinearAlgebra> c = matrix 2 [1..6]
Prelude Numeric.LinearAlgebra> c
(3><2)
 [ 1.0, 2.0
 , 3.0, 4.0
 , 5.0, 6.0 ]
Prelude Numeric.LinearAlgebra> d = a <> c
Prelude Numeric.LinearAlgebra> d
(2><2)
 [ 22.0, 28.0
 , 49.0, 64.0 ]
Prelude Numeric.LinearAlgebra> vector [1,2] <# a
[9.0,12.0,15.0]
Prelude Numeric.LinearAlgebra> a #> vector [1,2,3]
[14.0,32.0]

Prelude Numeric.LinearAlgebra> e = ident 2
Prelude Numeric.LinearAlgebra> e
(2><2)
 [ 1.0, 0.0
 , 0.0, 1.0 ]
Prelude Numeric.LinearAlgebra> d <> e
(2><2)
 [ 22.0, 28.0
 , 49.0, 64.0 ]
Prelude Numeric.LinearAlgebra> e <> d
(2><2)
 [ 22.0, 28.0
 , 49.0, 64.0 ]
Prelude Numeric.LinearAlgebra> d <> e == d
True
Prelude Numeric.LinearAlgebra> e <> d == d
True
Prelude Numeric.LinearAlgebra> e <> e
(2><2)
 [ 1.0, 0.0
 , 0.0, 1.0 ]
Prelude Numeric.LinearAlgebra> e <> e == e
True

Prelude Numeric.LinearAlgebra> a1 = matrix 2 [2,3,4,7]
Prelude Numeric.LinearAlgebra> a1
(2><2)
 [ 2.0, 3.0
 , 4.0, 7.0 ]
Prelude Numeric.LinearAlgebra> ra1 = inv a1
Prelude Numeric.LinearAlgebra> ra1
(2><2)
 [  3.5, -1.5
 , -2.0,  1.0 ]
Prelude Numeric.LinearAlgebra> a1 <> ra1
(2><2)
 [ 1.0, 0.0
 , 0.0, 1.0 ]
Prelude Numeric.LinearAlgebra> ra1 #> vector [4,6]
[5.0,-2.0]
Prelude Numeric.LinearAlgebra> ra1 <> matrix 1 [4,6]
(2><1)
 [  5.0
 , -2.0 ]
Prelude Numeric.LinearAlgebra> linearSolve a1 $ matrix 1 [4,6]
Just (2><1)
 [  5.0
 , -2.0 ]

det A, 行列 A の行列式 (determinant)
行列 A が正則行列であれば、det A は 0 にならない
rank A, 行列 A のランク (rank)
行列 A が n 次の正則行列であれば、rank A の値は n になる
行列 A のトレース (trace, 跡) は行列の対角線上にある要素の和のこと

n 次の正方行列 A のトレースは \(a_{11} + a_{22} + \cdots + a_{nn}\)
hmatrix にトレースを求める関数はなさそうだ (M.Hiroi の探し方が悪いのかもしれない)
関数 takeDiag で対角成分を取り出し、その和を関数 sumElements で求めればよい

takeDiag :: Element t => Matrix t -> Vector t
sumElements :: Container c e => c e -> e

Prelude Numeric.LinearAlgebra> a1
(2><2)
 [ 2.0, 3.0
 , 4.0, 7.0 ]
Prelude Numeric.LinearAlgebra> det a1
2.0
Prelude Numeric.LinearAlgebra> rank a1
2
Prelude Numeric.LinearAlgebra> b1 = matrix 3 [1..9]
Prelude Numeric.LinearAlgebra> b1
(3><3)
 [ 1.0, 2.0, 3.0
 , 4.0, 5.0, 6.0
 , 7.0, 8.0, 9.0 ]
Prelude Numeric.LinearAlgebra> det b1
6.661338147750939e-16
Prelude Numeric.LinearAlgebra> rank b1
2

Prelude Numeric.LinearAlgebra> trace m = sumElements $ takeDiag (m::Matrix R)
Prelude Numeric.LinearAlgebra> trace a1
9.0
Prelude Numeric.LinearAlgebra> trace b1
15.0

●ベクトルと行列の変換

行列をベクトルに変換する主な関数を以下に示す

toRows :: Element t => Matrix t -> [Vector t]
toColumns :: Element t => Matrix t -> [Vector t]
flatten :: Element t => Matrix t -> Vector t

toRows は行の要素をベクトルに格納し、それらをリストに格納して返す
toColumns は列の要素をベクトルに格納し、それらをリストに格納して返す
flatten は toRows の結果を平坦化したもの

Prelude Numeric.LinearAlgebra> a = matrix 3 [1..9]
Prelude Numeric.LinearAlgebra> a
(3><3)
 [ 1.0, 2.0, 3.0
 , 4.0, 5.0, 6.0
 , 7.0, 8.0, 9.0 ]
Prelude Numeric.LinearAlgebra> toRows a
[[1.0,2.0,3.0],[4.0,5.0,6.0],[7.0,8.0,9.0]]
Prelude Numeric.LinearAlgebra> toColumns a
[[1.0,4.0,7.0],[2.0,5.0,8.0],[3.0,6.0,9.0]]
Prelude Numeric.LinearAlgebra> flatten a
[1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0]

ベクトルを行列に変換する主な関数を以下に示す

reshape :: Foreign.Storable.Storable t => Int -> Vector t -> Matrix t
asRow :: Foreign.Storable.Storable a => Vector a -> Matrix a
asColumn :: Foreign.Storable.Storable a => Vector a -> Matrix a
fromRows :: Element t => [Vector t] -> Matrix t
fromColumns :: Element t => [Vector t] -> Matrix t

reshape の第 1 引数は列数 (行の要素数) を表す (関数 matrix と同じ)
asRow はベクトルを行ベクトルに、asColumn は列ベクトルに変換する
fromRows は引数のベクトルが行の成分になる
fromColumn は引数のベクトルが列の成分となる

Prelude Numeric.LinearAlgebra> reshape 3 $ vector [1..9]
(3><3)
 [ 1.0, 2.0, 3.0
 , 4.0, 5.0, 6.0
 , 7.0, 8.0, 9.0 ]
Prelude Numeric.LinearAlgebra> asRow $ vector [1,2,3]
(1><3)
 [ 1.0, 2.0, 3.0 ]
Prelude Numeric.LinearAlgebra> asColumn $ vector [1,2,3]
(3><1)
 [ 1.0
 , 2.0
 , 3.0 ]
Prelude Numeric.LinearAlgebra> fromRows [vector [1,2,3], vector [4,5,6]]
(2><3)
 [ 1.0, 2.0, 3.0
 , 4.0, 5.0, 6.0 ]
Prelude Numeric.LinearAlgebra> fromColumns [vector [1,2,3], vector [4,5,6]]
(3><2)
 [ 1.0, 4.0
 , 2.0, 5.0
 , 3.0, 6.0 ]

●ベクトルと行列の生成

関数 konst はすべての成分が e のベクトル (行列) を生成する

konst :: Konst e d c => e -> d -> c e

第 1 引数が成分の値
第 2 引数がベクトル (行列) の大きさ
Int だとベクトル、(Int, Int) だと行列を生成する

関数 build は添字を関数 f に渡してベクトル (行列) の成分を生成する

build :: Build d f c e => d -> f -> c e

第 1 引数がベクトル (行列) の大きさ (指定方法は konst と同じ)
第 2 引数に成分を生成する関数 f を渡す
関数 f の引数はベクトルだと \i -> ..., 行列だと \i j -> ...

関数 assoc は行列の成分を連想リストで指定する

assoc :: Container c e => IndexOf c -> e -> [(IndexOf c, e)] -> c e

第 1 引数がベクトル (行列) の大きさ (指定方法は konst と同じ)
第 2 引数が成分のデフォルト値
第 3 引数が成分を指定する連想リスト
IndexOf は添字を表すデータ型 (ベクトルだと Int, 行列だと (Int, Int) になる)

関数 accum は連想リストの値とベクトル (行列) の成分を関数に渡して評価し、その結果で置き換えたベクトル (行列) を返す

accum :: Container c e => c e -> (e -> e -> e) -> [(IndexOf c, e)] -> c e

第 1 引数がベクトル (行列)
第 2 引数が更新関数 \x y -> ...

第 1 引数 x が連想リストの値, 第 2 引数 y がベクトル (行列) の成分

第 3 引数が連想リスト

関数 linspace は区間 [s, e] を n - 1 等分した成分を持つベクトルを生成する

linspace :: (Fractional e, Container Vector e) => Int -> (e, e) -> Vector e

第 1 引数がベクトルの大きさ
第 2 引数が区間を表すタプル
ベクトルの大きさを n とすると、始点と終点と区間を n - 1 等分した n - 2 個の点がベクトルの成分となる

Prelude Numeric.LinearAlgebra> konst 0 10 :: Vector R
[0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0]
Prelude Numeric.LinearAlgebra> konst 1 10 :: Vector R
[1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0]
Prelude Numeric.LinearAlgebra> konst 2 10 :: Vector I
[2,2,2,2,2,2,2,2,2,2]
Prelude Numeric.LinearAlgebra> konst 0 (3, 3) :: Matrix R
(3><3)
 [ 0.0, 0.0, 0.0
 , 0.0, 0.0, 0.0
 , 0.0, 0.0, 0.0 ]
Prelude Numeric.LinearAlgebra> konst 1 (2, 3) :: Matrix R
(2><3)
 [ 1.0, 1.0, 1.0
 , 1.0, 1.0, 1.0 ]
Prelude Numeric.LinearAlgebra> konst 1 (2, 2) :: Matrix I
(2><2)
 [ 1, 1
 , 1, 1 ]

Prelude Numeric.LinearAlgebra> build 10 id :: Vector R
[0.0,1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0]
Prelude Numeric.LinearAlgebra> build 10 (*2) :: Vector R
[0.0,2.0,4.0,6.0,8.0,10.0,12.0,14.0,16.0,18.0]
Prelude Numeric.LinearAlgebra> build 10 (**2) :: Vector R
[0.0,1.0,4.0,9.0,16.0,25.0,36.0,49.0,64.0,81.0]
Prelude Numeric.LinearAlgebra> build (3, 3) (+) :: Matrix R
(3><3)
 [ 0.0, 1.0, 2.0
 , 1.0, 2.0, 3.0
 , 2.0, 3.0, 4.0 ]
Prelude Numeric.LinearAlgebra> build (3, 3) (\i j -> i * 3 + j) :: Matrix R
(3><3)
 [ 0.0, 1.0, 2.0
 , 3.0, 4.0, 5.0
 , 6.0, 7.0, 8.0 ]
Prelude Numeric.LinearAlgebra> build (3, 3) (\i j -> j * 3 + i) :: Matrix R
(3><3)
 [ 0.0, 3.0, 6.0
 , 1.0, 4.0, 7.0
 , 2.0, 5.0, 8.0 ]

Prelude Numeric.LinearAlgebra> assoc 5 0 [(1,10), (3,20)] :: Vector R
[0.0,10.0,0.0,20.0,0.0]
Prelude Numeric.LinearAlgebra> assoc 5 1 [(1,10), (3,20)] :: Vector R
[1.0,10.0,1.0,20.0,1.0]
Prelude Numeric.LinearAlgebra> assoc 7 2 [(1,10), (3,20), (5,30)] :: Vector R
[2.0,10.0,2.0,20.0,2.0,30.0,2.0]
Prelude Numeric.LinearAlgebra> assoc (3, 3) 0 [((0,0),1), ((2,2),2)] :: Matrix R
(3><3)
 [ 1.0, 0.0, 0.0
 , 0.0, 0.0, 0.0
 , 0.0, 0.0, 2.0 ]
Prelude Numeric.LinearAlgebra> assoc (3, 3) 1 [((1,0),11), ((2,1),12)] :: Matrix R
(3><3)
 [  1.0,  1.0, 1.0
 , 11.0,  1.0, 1.0
 ,  1.0, 12.0, 1.0 ]
Prelude Numeric.LinearAlgebra> assoc (3, 3) 2 [((1,0),21), ((0,2),22), ((2,1),23)] :: Matrix R
(3><3)
 [  2.0,  2.0, 22.0
 , 21.0,  2.0,  2.0
 ,  2.0, 23.0,  2.0 ]

Prelude Numeric.LinearAlgebra> a = vector [1,2,3,4,5,6,7,8]
Prelude Numeric.LinearAlgebra> a
[1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0]
Prelude Numeric.LinearAlgebra> accum a (+) [(2,10), (5,20)]
[1.0,2.0,13.0,4.0,5.0,26.0,7.0,8.0]
Prelude Numeric.LinearAlgebra> accum a (-) [(1,10), (6,20)]
[1.0,8.0,3.0,4.0,5.0,6.0,13.0,8.0]
Prelude Numeric.LinearAlgebra> accum a const [(2,10), (4,20), (6, 30)]
[1.0,2.0,10.0,4.0,20.0,6.0,30.0,8.0]
Prelude Numeric.LinearAlgebra> b = matrix 3 [1..9]
Prelude Numeric.LinearAlgebra> b
(3><3)
 [ 1.0, 2.0, 3.0
 , 4.0, 5.0, 6.0
 , 7.0, 8.0, 9.0 ]
Prelude Numeric.LinearAlgebra> accum b (+) [((0,0),10), ((2,1),20), ((1,2), 30)]
(3><3)
 [ 11.0,  2.0,  3.0
 ,  4.0,  5.0, 36.0
 ,  7.0, 28.0,  9.0 ]
Prelude Numeric.LinearAlgebra> accum b (*) [((0,0),10), ((2,1),20), ((1,2), 30)]
(3><3)
 [ 10.0,   2.0,   3.0
 ,  4.0,   5.0, 180.0
 ,  7.0, 160.0,   9.0 ]
Prelude Numeric.LinearAlgebra> accum b const [((0,0),10), ((2,1),20), ((1,2), 30)]
(3><3)
 [ 10.0,  2.0,  3.0
 ,  4.0,  5.0, 30.0
 ,  7.0, 20.0,  9.0 ]

Prelude Numeric.LinearAlgebra> linspace 5 (1,5) :: Vector R
[1.0,2.0,3.0,4.0,5.0]
Prelude Numeric.LinearAlgebra> linspace 9 (1,5) :: Vector R
[1.0,1.5,2.0,2.5,3.0,3.5,4.0,4.5,5.0]
Prelude Numeric.LinearAlgebra> linspace 10 (0,9) :: Vector R
[0.0,1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0]
Prelude Numeric.LinearAlgebra> linspace 19 (0,9) :: Vector R
[0.0,0.5,1.0,1.5,2.0,2.5,3.0,3.5,4.0,4.5,5.0,5.5,6.0,6.5,7.0,7.5,8.0,8.5,9.0]

Functional Programming