hmatrix について(前書き)

高速な行列計算ライブラリーです。基本的な使い方はほかの方のqiita記事にも詳しく書いてありますので省略しますが、私が躓いたところに関して、身の丈に合った解説をします。

一言だけ初歩的解説しますと、import Numeric.LinearAlgebraで呼び出します。

大きく2点躓きました。1つは複素行列の定義、2つは外部ファイルからの行列読み込みです。そしてそれらの合成にも躓きました。

複素行列

hmatrixは複素行列の計算にも対応しいる。（というより、していないと閉じない計算が多い）
定義する場合は、後で複素行列であることを明示すればよい。

let a = (2><2)[1,2,3,4] :: Matrix (Complex Double)
let b = (2><2)[1:+2,2,3,4:+(-4)] :: Matrix (Complex Double)

よくある間違いとして（私がよくやる間違い）

let a = (2><2)[1..4] :: Matrix (Complex Double)
--複素数なので".."表現が使えない。
let b = (2><2)[1:+2,2,3,4:+(-4)]
--Matrix Double として扱おうとしてエラーになる。ちょま、

haskellはコンパイラが型推論してくれるが、Matrixに関しては自分で変換する必要がある。
実行列 a と複素行列 b の積は

let a = (2><2)[1,0,0,1]::Matrix Double
let b = (2><2)[-1,0,0,1]::Matrix (Complex Double)

let c = (complex a) <> b

とする。ちなみに、Complex Doubleに complex b などと書いても大丈夫なので、場合によっては複素行列になるといった場合にも対応できる(可読性は下がる)。

saveMatrix "FilePath" "%lf" a

列は半角スペース区切り行は改行区切りで"FilePath"に出力される。カンマでも区切りたいなら、

saveMatrix "FilePath" "%lf," a

のようにする。こうすると半角スペース、改行に加えカンマが毎回挿入される。同様に任意の区切りが可能。

ファイルから行列をロードする方法

let a = loadMatrix "FilePath"

でとることができる。しかし、a の型は IO (Matrix Double) のようなIO型になるため、

svd a
a + a

のようなことはできない。そこでIO型からMatrix Doubleを取り出すため、"<-"を用います。
(System.IO.UnsafeのunsafePerformIOという関数を用いる方法もあります。)

a <- loadMatrix "FilePath"
svd a

これにより、ファイルから読み込んだ行列用いて計算し行くことができる。
大きさは自動で判別されるが、列は半角スペース、行は改行で区切られている必要がある。(カンマ区切りはダメ)

saveMatrix, LoadMatrixは実行列Matrix Doubleにしか対応していない。複素行列を保存する場合、自前で作るか、実部、虚部を分けるかする必要がある。ここでは後者を紹介する。

toComplex, fromComplexを用いると比較的簡単に書くことができる。

let a = (2><2)[1:+2,2,3,4]::Matrix(Complex Double)
let (a_r,a_c) = fromComplex a
saveMatrix "real.txt" "%lf" a_r
saveMatrix "complex.txt" "%lf" a_c

2行で書くと

let a = (2><2)[1:+2,2,3,4]::Matrix(Complex Double)
saveMatrix "real.txt" "%lf" (fst$fromComplex a)
saveMatrix "complex.txt" "%lf" (snd$fromComplex a)

プチ解説するとfst, sndは2要素のタプルに対して、第1成分と第2成分を取るものです。

a_r <-loadMatrix "real.txt"
a_c <-loadMatrix "complex.txt"
let a = toComplex (a_r,a_c)

今後も、このテーマに関して新たな見地を得たら加筆修正しようと思います。忌憚なきご意見を頂ければ幸いです。(近いうちに複素行列を使う例を挙げようと思います。)

大きな修正は記録しておきます。(加筆はここには書きません)