hmatrixはBLAS/LAPACKをベースにしたHaskellで行列を扱うライブラリです
調べるとlinearという線形代数のためのライブラリがありますが依存関係を見ているとlinearはlinearはOpenGL等のグラフィックの処理に使われることが多く数値計算にはhmatrixが使われているみたいです。
ベクトルと行列
hmatrixでのベクトルと行列の作り方と簡単な使い方を紹介します。
ghciでNumeric.LinearAlgebraをimportした前提で書いていきます。
> import Numeric.LinearAlgebra
Vector はベクトルを表す型でvector :: [R] -> Vector R を使って作るのが簡単です。ここでRはDoubleのエイリアスです
> vector [1..5]
[1.0,2.0,3.0,4.0,5.0]
> size (vector [1..5])
5
Matrixは行列を表す型です。サイズや行サイズだけを指定して生成する方法などいろんな作り方が用意されています。
> (2><3) [1..6]
(2><3)
[ 1.0, 2.0, 3.0
, 4.0, 5.0, 6.0 ]
> matrix 5 [1..15]
(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 ]
> size ((2><3) [1..6] :: Matrix R)
(2,3)
> ident 3
(3><3)
[ 1.0, 0.0, 0.0
, 0.0, 1.0, 0.0
, 0.0, 0.0, 1.0 ]
VectorやMatrixはNumのインスタンスになっているので四則演算を行うことができます。これらは基本的に成分同士の演算になりますが内積や行列積ももちろん用意されています。
> let v1 = vector [1,2]
> let v2 = vector [3,4]
> v1 + v2
[4.0,6.0]
> v1 * v2
[3.0,8.0]
> v1 <.> v2
11.0
> let m1 = (2><2) [1..4] :: Matrix R
> let m2 = (2><2) [5..8] :: Matrix R
> m1 + m2
(2><2)
[ 6.0, 8.0
, 10.0, 12.0 ]
> m1 * m2
(2><2)
[ 5.0, 12.0
, 21.0, 32.0 ]
-- 行列積
> m1 <> m2
(2><2)
[ 19.0, 22.0
, 43.0, 50.0 ]
-- 行列とベクトルの積
> v1 <# m1
[7.0,10.0]
> m2 #> v2
[39.0,53.0]
他にも様々な関数が用意されています。
> let v = vector [1,2]
> let m = (2><2) [1..4] :: Matrix R
-- フロベニウスノルム
> norm_1 v
3.0
-- ユークリッドノルム
> norm_2 v
2.23606797749979
-- 転置行列
> tr m
(2><2)
[ 1.0, 3.0
, 2.0, 4.0 ]
-- 逆行列
> inv m
(2><2)
[ -1.9999999999999996, 0.9999999999999998
, 1.4999999999999998, -0.4999999999999999 ]
-- 階数
> rank m
2
-- 行列式
> det m
-2.0
-- 固有値・固有ベクトル
> eig m
([(-0.3722813232690143) :+ 0.0,5.372281323269014 :+ 0.0],
(2><2)
[ (-0.8245648401323938) :+ 0.0, (-0.4159735579192842) :+ 0.0
, 0.5657674649689923 :+ 0.0, (-0.9093767091321241) :+ 0.0 ])
-- 最小二乗法
> m <\> v
[1.3877787807814457e-16,0.4999999999999998]
> m #> (m <\> v)
[0.9999999999999997,1.9999999999999996]
-- 特異値分解
> svd m
((2><2)
[ -0.40455358483375703, -0.9145142956773042
, -0.9145142956773045, 0.4045535848337568 ],
[5.464985704219043,0.3659661906262574],
(2><2)
[ -0.5760484367663209, 0.8174155604703631
, -0.8174155604703631, -0.5760484367663209 ])
-- LU分解
> lu m
((2><2)
[ 1.0, 0.0
, 0.3333333333333333, 1.0 ],
(2><2)
[ 3.0, 4.0
, 0.0, 0.6666666666666667 ],
(2><2)
[ 0.0, 1.0
, 1.0, 0.0 ],-1.0)
-- コレスキー分解
> chol (sym (tr m <> m))
(2><2)
[ 3.1622776601683795, 4.427188724235731
, 0.0, 0.6324555320336748 ]
などなど！
まだまだあるのでぜひHaddockの方も読んでみてください
ニューラルネットワーク
ディープではありません(笑)
作るものはとても単純な三層パーセプトロンで
- 中間層の活性化関数はReLU
\phi(u) = \max(0,u)
- 出力層の活性化関数はソフトマックス関数
\phi(u_k) = \frac{\exp^{u_k}}{\sum_i \exp^{u_i}}
- 出力層の誤差関数はクロスエントロピー
E(w) = -\sum_{n,k} d_{n,k}\log y_{n,k}
- 学習方法は誤差逆伝播法
こんな感じでXORの分類器を作ってみます。
Main.hs
import Numeric.LinearAlgebra
relu :: R -> R
relu = max 0
softmax :: Vector R -> R -> R
softmax xs x = exp x / sumElements (cmap exp xs)
tlp :: (Matrix R, Vector R) -- 入力層から中間層への重み行列と定数項
-> (Matrix R, Vector R) -- 中間層から出力層への重み行列と定数項
-> Vector R -- 入力
-> Vector R -- 出力
tlp (w1, b1) (w2, b2) input =
let hidden = cmap relu $ w1 #> input + b1
in softmax >>= cmap $ w2 #> hidden + b2
training :: R -- 学習係数
-> (Vector R, Vector R) -- 学習データ
-> ((Matrix R, Vector R)
,(Matrix R, Vector R)) -- 学習前の重み行列と定数項
-> ((Matrix R, Vector R)
,(Matrix R, Vector R)) -- 学習後の重み行列と定数項
training a (input, correct) ((w1, b1), (w2, b2)) =
let hidden = cmap relu $ w1 #> input + b1
output = softmax >>= cmap $ w2 #> hidden + b2
-- 出力層の学習
loss = output - correct
b2' = b2 - scalar a * loss
w2' = w2 - scalar a * (loss `outer` hidden)
-- 中間層の学習
h x = if x < 0 then 0 else 1
loss' = loss <# (w2 <> diag (cmap h hidden))
b1' = b1 - scalar a * loss'
w1' = w1 - scalar a * (loss' `outer` input)
in ((w1', b1'), (w2', b2'))
main :: IO ()
main = do
let trainingData = concat . replicate 1000 $
[ (vector [0, 0], vector [1, 0])
, (vector [0, 1], vector [0, 1])
, (vector [1, 0], vector [0, 1])
, (vector [1, 1], vector [1, 0])
]
hDim = 30
-- パラメータの初期化
w1 <- randn hDim 2
b1 <- flatten <$> randn 1 hDim
w2 <- randn 2 hDim
b2 <- flatten <$> randn 1 2
-- 学習
let (h,o) = foldr (training 0.2) ((w1,b1),(w2,b2)) trainingData
-- 結果
print $ tlp h o (vector [0, 0])
print $ tlp h o (vector [0, 1])
print $ tlp h o (vector [1, 0])
print $ tlp h o (vector [1, 1])
$ runhaskell Main.hs
[0.9960662495329804,3.933750467019611e-3]
[3.66632024396632e-3,0.9963336797560337]
[3.5199543634252773e-3,0.9964800456365748]
[0.9960076764061874,3.992323593812549e-3]
うーん、教科書通りに実装したのですがうまく動いてそうです
VectorやMatrixはFunctorのインスタンスにはなってませんがContainerという型クラスのメソッドcmapを使うことで一つ一つの要素に対して関数を適用することができます。これをつかって活性化関数を適用しています。
VectorやMatrixの値を定数倍する時はscalarという関数が便利です。
> scalar 2 * vector [1,2,3]
[2.0,4.0,6.0]
> scalar 2 * (2><2) [1..] :: Matrix R
(2><2)
[ 2.0, 4.0
, 6.0, 8.0 ]
outerはベクトル同士の要素を掛け合わせたような行列を作ります
> vector [1,2,3] `outer` vector [0.1, 1, 10]
(3><3)
[ 0.1, 1.0, 10.0
, 0.2, 2.0, 20.0
, 0.3, 3.0, 30.0 ]
randnはランダムな要素をもつ行列を作れる便利な関数です。
> randn 2 2
(2><2)
[ 0.9147819402345773, -0.49367886497544466
, 0.30457323650647944, -1.5975938635395308 ]
以上、hmatrixの使い方を中心に実装例をまじえて簡単に説明してみました。