Haskellで数値計算 線形代数編

BLASLAPACK

数値計算では、連立線形方程式に帰着させられる問題はとても多い。また、固有値問題は分野を越えて応用範囲が広い。そのため、これらの問題を扱うための線形代数ライブラリは古くから活発に開発されてきた。中でも代表的なものは基礎的な行列、ベクトル演算を行うBLASと、その上で線形方程式系や最小2乗問題、固有値問題などを解くLAPACKで、最適化されており高速に動作するので、多くの言語のライブラリでも背景で利用されている。

HaskellでもhmatrixパッケージがBLAS/LAPACKバインディングを提供してくれている。BLAS/LAPACKの関数の呼び出しで処理出来る問題の場合、自分の実装でこの速度を超えるのは難しいので、出来るだけ利用して行きたい。

目次

hmatrixの使い方

インストール

Haskellを使うのが初めての場合、まずはビルドツールstackから環境を構築するのがおすすめ。stack自体のインストールはこちらを参照。
Install/upgrade - The Haskell Tool Stack

hmatrixパッケージをインストールする。色々なパッケージを使うようになっても、プロジェクトごとに依存関係をstackが解決してくれる。

$ stack update
$ stack install hmatrix

cabalを直接使う場合は以下の通り。環境によってはsudoが必要。

$ cabal update
$ cabal install hmatrix

シェルでghciを起動して、:mでモジュールを読み込める。またファイルに書いたコードは、:lでロードできる。

$ stack ghci
> :m Numeric.LinearAlgebra
> :set prompt "> " -- promptの文字列を変えている ~/.ghciに書いてもよい
> :l [file]

この記事ではファイルに書いたコードをghciから:lでロードし、> に続くコードを実行して試していくことを想定している。main関数を持つファイルはghcコンパイルして実行可能バイナリにもできる。

MatrixやVectorの作成と基本的な演算

主な関数を表にまとめてみた。詳細はパッケージのドキュメントを参照。
hmatrix: Numeric Linear Algebra

モジュールのimport。

import Numeric.LinearAlgebra
import Numeric.LinearAlgebra.Data

-- ghci
> :m Numeric.LinearAlgebra
> :m Numeric.LinearAlgebra.Data

lts-12.0からPreludeの(<>)という同名の関数と競合するようになってしまったので、Numeric.LinearAlgebraの方をqualified importするか、Preludeの方をhidingする必要がある。この記事では簡単のため、Preludeの方をhidingする。古いltsを使っている場合は、この行を消せば動くはず。

-- 名前空間を分ける
-- Numeric.LinearAlgebraの関数はL.付きで呼ぶ必要がある
import qualified Numeric.LinearAlgebra as L

-- Preludeの方を隠す
import Prelude hiding ((<>))

リストからの生成やリストへの変換

関数 説明
(><) Int -> Int -> [a] -> Matrix a 行数と列数を指定する
toLists Matrix a -> [ [a] ]
fromLists [ [a] ] -> Matrix a
vector [a] -> Vector a
toList Vector a -> [a]
flatten Matrix a -> Vector a 行ベクトルを1列につなげる
reshape Int -> Vector a -> Matrix a 列数を指定する
toRows/Columns Matrix a -> [Vector a]
fromRows/Columns [Vector a] -> Matrix a
rand Int -> Int -> IO (Matrix Double) [0,1)一様分布擬似乱数
randn Int -> Int -> IO (Matrix Double) [0,1)ガウス分布擬似乱数

基本的な2項演算

演算 左辺値 演算子 右辺値
 \bf{u} \pm \bf{v} Vector + - Vector
 \bf{u} \cdot \bf{v} Vector <.> Vector
 \bf{u} \times \bf{v} Vector cross Vector
 A \pm B Matrix + - Matrix
 A \bf{u} Matrix #> Vector
 {\bf u} A Vector <# Matrix
 A B Matrix <> Matrix
 A \otimes B Matrix kronecker Matrix

要素ごとの操作など

操作 関数 説明
要素ごとの関数適用 cmap (e -> b) -> c e -> c b
indexでのアクセス ! 行列の場合は行ベクトルに
行列のスライス ?? 行と列の位置などを指定する
m ?? (Pos (idxs [3,4]), Pos (idxs [1,0]))
3行と4行、1列と0列を取る
m ?? (Take 3, Drop 2)
3行分を取り、2列分を捨てる
最大値
最大値のindex
maxElement
maxIndex
最小値
最小値のindex
minElement
minIndex
昇順にソート
その位置の値のソート後のindex
sortVector
sortIndex
L1/L2/L \inftyノルム (Vector)
作用素ノルム (Matrix)
norm_1
norm_2
norm_Inf
 \mid\mid {\bf x} \mid\mid_p = \left( \sum_i \mid x_i \mid^p \right)^{1/p}
 \mid\mid A \mid\mid = \sup_{\mid\mid {\bf x} \mid\mid_p = 1} \mid\mid A {\bf x} \mid\mid_p
この行列で写された任意の単位ベクトルのノルムの最大値
フロベニウスノルム norm_Frob  \mid\mid A \mid\mid_F = \sqrt{\sum_{i} \sum_{j} \mid a_{ij} \mid^2}
対角成分をVector takeDiag

特殊な行列の生成など

行列 関数 説明
 I_n ident n n次単位行列
 a_{ij} = \delta_{ij} x_i diag x 対角成分から行列をつくる
 ^t A tr 転置行列
 \det{A} det 行列式
 A^{-1} inv 逆行列
 (A + ^t A) / 2 sym 対称行列をつくる
 ^t AA mTm 転置行列との積

3次元の回転行列を使って座標を回転させてみる。

import Prelude hiding ((<>))
import Numeric.LinearAlgebra

rotateX :: Double -> Matrix Double
rotateX t = (3><3) [1.0,   0.0,    0.0,
                    0.0, cos t, -sin t,
                    0.0, sin t,  cos t]

rotateY :: Double -> Matrix Double
rotateY t = (3><3) [ cos t, 0.0, sin t,
                       0.0, 1.0,   0.0,
                    -sin t, 0.0, cos t]

rotateZ :: Double -> Matrix Double
rotateZ t = (3><3) [cos t, -sin t, 0.0,
                    sin t,  cos t, 0.0,
                      0.0,    0.0, 1.0]

> (rotateX pi <> rotateY pi <> rotateZ pi) #> vector [1.0,1.0,1.0]
[1.0000000000000004,1.0,0.9999999999999998]

連立線形方程式

連立線形方程式は行列を使って A \bf{x} = \bf{b}のように書けるので、左から A^{-1}を掛ければ解くことが出来るが、とても小さな行列以外では、クラメルの公式で逆行列を計算するのはコストが高い。そこで、 Aを扱いやすい行列の積に分解したり、要素の反復計算で漸近するなど色々な工夫が行われている。

hmatrixで用意されている主な線形ソルバは以下の通り。基本的に行列を特殊な行列の積に分解することで解いてくれる。引数は \bf{x} \bf{b}を行列化した AX = Bの形であたえられる。行列分解を1度計算しておけば、同じ Aに対する複数の \bf{x} \bf{b}の組み合わせを効率的に解くことができる。

関数 対象 説明
<\> 優決定系/劣決定系を含む 汎用線形ソルバ
式が多すぎる系や足りない系でもノルムを最小にするような解を計算してくれる
linearSolve 正則行列 正則でない場合を考慮してMaybeを返す
中身はluSolve
luSolve 正則行列  AをLU分解 (luPacked) して渡す
1000  \times 1000程度までの小さな正則行列なら基本的にこれを使えば問題ない
ldlSolve エルミート行列
(実対称行列を含む)
 AをLDL分解 (ldlPacked) して渡す
修正コレスキー分解
制約があるが速い
cholSolve 正定値のエルミート行列  Aをコレスキー分解 (chol) して渡す
制約があるが速い
triSolve 上三角行列または下三角行列 制約があるがとても速い
QR分解 (qr) した上三角行列も解ける
 QR{\bf x} = {\bf b}
 R{\bf x} = ^t Q {\bf b}
triDiagSolve 三重対角行列 制約があるがとても速い
cgSolve 正定値のエルミート行列
疎行列を想定したindex表現であるGMatrixであたえる
共役勾配法 (反復法の一種)
状況によってはとても速いが収束しないこともある
大規模な疎行列の時有利

エルミート行列は転置して複素共役を取ると自身と一致するような行列。虚部がすべて0なら実対称行列になる。正定値行列は固有値がすべて正の実数であるような行列。基本的に対象となる行列への制約が強いソルバほど速いので、解こうとしている行列が制約をみたすなら対応するソルバを使うと良い。

乱数で初期化した1000  \times 1000行列をLU分解を使って解いて時間を測ってみる。簡単のためGHCのStrict拡張で正格評価にしている。

{-# LANGUAGE Strict #-}

import Prelude hiding ((<>))
import System.CPUTime
import Numeric.LinearAlgebra

randomLU :: Int -> IO ()
randomLU n = do
  -- 乱数で初期化
  a <- rand n n
  b <- rand n 1
  -- 時間を計測
  start <- getCPUTime
  let x = luSolve (luPacked a) b
  end <- getCPUTme
  -- 誤差を計算
  print $ fromIntegral (end - start) / 1e12
  print $ norm_2 (a <> x - b) / norm_2 b

> randomLU 1000
8.3574e-2
3.1152329739293184e-13

おそらくLAPACKの関数の実行速度が出ている。LU分解ではどの言語のどのような実装でも並列化しなければこれ以上速くするのは難しい気がする。さらに速度を求める場合は、(収束性を考慮する必要が出てくるが) 共役勾配法やSOR法のような反復法を使ったほうが良いかもしれない。

参考文献

連立線形方程式の数値解法についてはこの辺りがわかりやすかった。

http://ri2t.kyushu-u.ac.jp/~watanabe/RESERCH/MANUSCRIPT/KOHO/GEPP/GEPP.pdf
http://nkl.cc.u-tokyo.ac.jp/16n/SolverDirect.pdf
http://nkl.cc.u-tokyo.ac.jp/16n/SolverIterative.pdf
Eigen: Linear algebra and decompositions
Eigen: Benchmark of dense decompositions
Eigen: Solving Sparse Linear Systems

固有値固有ベクトル

 A \bf{x} = \lambda \bf{x} ( \bf{x} \neq \bf{0}) を満たすような \lambda固有値といい、 \bf{x}固有ベクトルという。 A n \times n行列の時、複素数値も含めて、 n個の固有値と対応する n個の固有ベクトルが存在する (重解でなければ) 。これらをすべて求めるのが固有値問題

関数 対象 説明
eig 一般 固有値Vector固有ベクトルのMatrixをタプルで返す
eigenvalues 一般 固有値のみを求める場合
eigSH エルミート行列 制約があるが速い
eigenvaluesSH エルミート行列 固有値のみを求める場合
geigSH エルミート行列の一般化固有値問題  A{\bf x} = \lambda B {\bf x}
 A Bはエルミート行列、 Bは正定値

主成分分析を固有値問題として解いてみる。

import Prelude hiding ((<>))
import Numeric.LinearAlgebra
import Numeric.LinearAlgebra.Data

pca :: Matrix Double -> Vector Double
pca x = let n = fromIntegral $ cols x
            -- 各成分の平均
            avg = vector . map ((/n) . sumElements) $ toRows x
            -- 中心化
            ctr = fromColumns . map (+ (-avg)) $ toColumns x
            vcm = cmap (/n) $ ctr <> tr ctr
            -- 共分散行列の固有値と固有ベクトルを求める
            (s, w) = eig vcm
            -- 分散が最大になる方向を選択する
        in  cmap realPart $ toColumns w !! maxIndex s

randomPca :: Int -> Int -> IO ()
randomPca l m = do
  -- l次元の成分を持つ標本をm個つくる
  x <- rand l m
  let v = pca x
  -- 第1主成分ベクトル
  print v
  -- 第1主成分軸方向への射影
  print . map (<.> v) $ toColumns x

> randomPca 3 10

1000  \times 1000行列について、eig関数の処理時間も測ってみる。基本が非正格評価なので式ごとの処理時間の計測などは簡単ではなさそうなのだが、GHCのStrict拡張やseq関数などを使うと、妥当そうな結果が得られることもある。

randomEigen :: Int -> IO ()
randomEigen n = do
  a <- rand n n

  start <- getCPUTime
  let (s, w) = eig a
  end <- seq w getCPUTime

  print $ fromIntegral (end - start) / 1e12
  print $ norm_2 (complex a <> w - w <> diag s) / norm_2 s

> randomEigen 1000
2.90387
2.2744873540326555e-15

参考文献

固有値問題の数値解法についてはこの辺りがわかりやすかった。

http://nalab.mind.meiji.ac.jp/~mk/labo/text/eigenvalues.pdf
http://www.r-ccs.riken.jp/r-ccssite/wp-content/uploads/2014/03/sps14_enshu2-1.pdf
Eigen: Catalogue of dense decompositions

numpyやEigenとの比較

同様の扱いやすいライブラリであるPythonのnumpyやC++のEigenについてもサンプルコードを書いてみた。また、関数呼び出し部分のみ手元の環境での実行時間を測ってみた。100回計測して平均を取っている。

線形方程式系
1000  \times 1000
線形方程式系
200  \times 200
固有値問題
200  \times 200
説明
hmatrix 0.08489044s 0.00156016s 0.05794971s 並列化しないなら
repaより速い
numpy 0.0847463s 0.0014532s 0.0577511s 処理はできるだけ
numpyで完結したい
Eigen 0.0618023s 0.00073476s 0.0619551s 速度的にはBoost
よりEigenが良い

環境によって多少の前後はあるかもしれないが、すべて同等の速度だと見ていいと思う。numpyは後ろでLAPACKを使っていた気がするので当たり前の結果かもしれない。Eigenは選択できる模様。

Eigen: Using BLAS/LAPACK from Eigen

小さなコードでは書きやすさもあまり変わらないが、大きなアプリケーションになった時はどうだろうか。計測に使用したコードは以下の通り。

hmatrix

{-# LANGUAGE Strict #-}

import Prelude hiding ((<>))
import Numeric.LinearAlgebra
import System.CPUTime
import System.Environment (getArgs)

randomLU :: Int -> IO Double
randomLU n = do
  a <- rand n n
  b <- rand n 1

  start <- getCPUTime
  let x = luSolve (luPacked a) b
  end <- seq x getCPUTime

  let t = fromIntegral (end - start) / 1e12
  print t
  print $ norm_2 (a <> x - b) / norm_2 b
  return t

randomEigen :: Int -> IO Double
randomEigen n = do
  a <- rand n n

  start <- getCPUTime
  let (s, w) = eig a
  end <- seq w getCPUTime

  let t = fromIntegral (end - start) / 1e12
  print t
  print $ norm_2 (complex a <> w - w <> diag s) / norm_2 s
  return t

main :: IO ()
main = do
  [n,r] <- map read <$> getArgs
  t <- mapM randomLU $ replicate r n
--  t <- mapM randomEigen $ replicate r n
  print $ sum t / fromIntegral r

$ stack ghc -- -o lh -O2 la.hs
$ ./lh 1000 100

numpy

import numpy as np
import numpy.linalg as la
import sys
import time

def randomSLE(rank):
    a = np.random.rand(rank, rank)
    b = np.random.rand(rank)

    start = time.process_time()
    x = la.solve(a, b)
    end = time.process_time()

    t = end - start
    print("duration: {:.7f}".format(t) + "s")
    print(la.norm(np.dot(a, x) - b) / la.norm(b))
    return t

def randomEigen(rank):
    a = np.random.rand(rank, rank)

    start = time.process_time()
    s, w = la.eig(a)
    end = time.process_time()

    t = end - start
    print("duration: {:.7f}".format(t) + "s")

    d = np.matmul(a, w) - np.dot(w, s)
#    d = np.matmul(a, w) - np.matmul(w, np.dot(np.identity(rank), s))
    print(la.norm(d) / la.norm(s))
    return t

rank = int(sys.argv[1])
repeat = int(sys.argv[2])

t = 0
for i in range(0, repeat):
    t += randomSLE(rank)
#    t += randomEigen(rank)
print("average: {:.7f}".format(t / repeat) + "s")

$ python3 la.py 1000 100

Eigen

#include <chrono>
#include <iostream>
#include <Eigen/Dense>

#define NDEBUG
#define EIGEN_NO_DEBUG

using namespace Eigen;

double randomLU(const int rank) {
  MatrixXd a = MatrixXd::Random(rank, rank);
  MatrixXd b = MatrixXd::Random(rank, 1);

  const auto start = std::chrono::system_clock::now();
  MatrixXd x = a.lu().solve(b);
  const auto end = std::chrono::system_clock::now();

  const double elapsed = std::chrono::duration<double>(end - start).count();
  std::cout << "duration: " << elapsed << "s" << std::endl;

  MatrixXd d = a * x - b;
  std::cout << d.norm() / b.norm() << std::endl;

  return elapsed;
}

double randomEigen(const int rank) {
  MatrixXd a = MatrixXd::Random(rank, rank);

//  ComplexEigenSolver<MatrixXd> es;
  EigenSolver<MatrixXd> es;

  const auto start = std::chrono::system_clock::now();
  es.compute(a);
  MatrixXcd s = es.eigenvalues();
  MatrixXcd w = es.eigenvectors();
  const auto end = std::chrono::system_clock::now();

  const double elapsed = std::chrono::duration<double>(end - start).count();
  std::cout << "duration: " << elapsed << "s" << std::endl;
  MatrixXcd d = a * w - w * s.asDiagonal();
  std::cout << d.norm() / s.norm() << std::endl;

  return elapsed;
}

int main(int argc, char *argv[]) {
  const int rank = std::atoi(argv[1]);
  const int repeat = std::atoi(argv[2]);

  double t = 0.0;
  for (int i = 0; i < repeat; i++) {
    t += randomLU(rank);
//    t += randomEigen(rank);
  }
  std::cout << "average: " << t / repeat << "s" << std::endl;

  return 0;
}

$ clang++ -Ofast -march=native -std=c++14 la.cpp -o lc
$ ./lc 1000 100

Free Monadで言語内DSLをつくろう

抽象構文木の構築とその評価の分離

インタプリタコンパイラでは最初に、コードの文字列を構文解析して抽象構文木 (AST) を構築する。その具体的な表現形式は変化して行くことも多いが、基本的にはこの抽象構文木に対して様々な最適化がほどこされる。インタプリタの場合にはこれが評価され実行されるが、コンパイラの場合にはこれがアセンブリの命令に置き換えられ、スケジューリングされ、レジスタが割り当てられ、最終的にはバイナリが生成される。こうした構文解析構文木の構築と、最適化、その評価の分離は、各工程をそれぞれ独立に作成し、必要に応じて取り替えられるという意味で、とても便利だ。

言語内DSLの処理系をそこまで重くするのは本末転倒感があるが、少なくとも構文木を構築しその評価と分離することは、抽象化の手段として適切な場合もあると思う。Haskellでは、言語内DSLのコードとその評価の分離には色々な方法 (型クラスによる抽象など) が考えられるが、ここではFree Monadを取り上げたい。自分で定義した、DSLを表現するデータ型をFunctorにし、Free Monadに入れれば、Monad計算で構文木を構築できる。これを評価するインタプリタは独立に複数用意できるし、評価先の型も自由に変えられる。

目次

Free Monadで言語内DSLをつくる

構文木を表現するデータ型の準備

まずは、自分のDSL構文木を表現するデータ型を定義する。Freeに包むためにはFunctorである必要があるが、この場合その定義は自明なので、簡単のためGHCのDeriveFunctor拡張を使ってderiving Functorしている。(この導出では型コンストラクタの最後の型変数に対応する値を写すようなfmapをつくってくれる。) またこのあと使う予定のモジュールもimportしておく。Freeには色々な実装や、関連するライブラリ、派生形などがあるが、ここでは最初の素朴な物を使っている。(処理効率を考えて実装されたControl.Monad.Free.Churchや、ASTのデータ型からFreeのアクションを導出してくれるControl.Monad.Free.THなどを参照されたし。)

{-# LANGUAGE DeriveFunctor #-}

import Control.Monad.Free
import Control.Monad.State.Lazy
import Control.Monad.Writer.Lazy

data MyAST r = Open r
             | Elephant String r
             | Close r
             | End
  deriving Functor

値コンストラクタはノードに対応し、型変数rは残りの木を表している。Freeに包むと、FreeとMyASTが交互に入ることになる。

構文木構築のための、FunctorからのFreeの作成

次に、MyASTの値から対応するFreeのアクションを作成する。Freeの値コンストラクタを使ってもいいが、liftFを使うと少し簡単に書ける。これらのFreeをつなげれば、Free MyAST a型の構文木が構築されていく。

{-
data Free f a = Pure a | Free (f (Free f a))
instance Functor f => Monad (Free f) where
  Free x >>= f = Free $ fmap (>>= f) x
  Pure a >>= f = f a
  return = Pure
liftF = Free . fmap return
-}

-- (対応するFreeの値) = liftF (DSLのデータ型の値)
-- (open = Free $ Open (Pure ()))
open  = liftF $ Open ()
close = liftF $ Close ()
end   = liftF End

putElephant :: String -> Free MyAST ()
putElephant name = liftF $ Elephant name ()

扉を開けて象を入れて扉を閉めてプログラムを終わらせるような命令を用意してみる。Freeは構文木の構築手段を用意しているだけなので、実際に何が行われるかはインタプリタが決める。

インタプリタの作成

最後に構文木を評価するためのインタプリタをつくる。Free MyAST a型の構文木を渡すと、最初のノードにパターンマッチして、何か処理を行い (例えばIntegerやStringなどの値に評価してもいいし、Monad計算に評価してもいい)、残りの木に対する処理を再帰的に呼び出す。葉はPureで表される。

interpretPrint :: Free MyAST a -> IO ()
interpretPrint f = case f of
  Free (Open r)       -> putStrLn "open" >> interpretPrint r
  Free (Elephant e r) -> putStrLn ("elephant: " ++ e) >> interpretPrint r
  Free (Close r)      -> putStrLn "close" >> interpretPrint r
  Free End            -> putStrLn "end"
  Pure _              -> return ()

DSLのプログラムを書く

準備ができたので、幾つかDSLのプログラムを書いて、インタプリタで評価してみる。

program1 :: Free MyAST ()
program1 = do
  open
  putElephant "Sally"
  putElephant "Emet"
  end

program2 :: Free MyAST ()
program2 = do
  open
  putElephant "Emet"
  close
  putElephant "Sally"
  end

先ほど定義した、MyASTの値が入ったFreeをつなげることで、このプログラムを表現する構文木が作成される。

> interpretPrint program1
open
elephant: Sally
elephant: Emet
end

これをインタプリタで評価すれば結果が返る。このインタプリタでは標準出力にプログラムの構成を表示するだけ。

色々なインタプリタをつくる

操作のログを取る

構文木の構築とその評価が分離されているので、同じプログラムを評価する色々なインタプリタをつくることができる。interpretPrintはプログラムをIOにしていたが、interpretNamesは象の名前のリストに評価してくれる。

interpretNames :: Free MyAST () -> Writer [String] ()
interpretNames f = case f of
  Free (Open r)       -> interpretNames r
  Free (Elephant e r) -> tell [e] >> interpretNames r
  Free (Close r)      -> interpretNames r
  Free End            -> return ()
  Pure _              -> return ()

WriterにMonoidをtellすると、後ろでmappendしてくれる。Monoidは単位元memptyを持ち、結合法則を満たす2項演算mappendが使えることを要請する型クラスで、例えばリストのmemptyは[]でmappendは(++)。完成したWriterをexecWriterすれば、このMonoidの演算結果を返してくれる。

> execWriter $ interpretNames program2
["Emet","Sally"]

操作の解析をする

冷蔵庫の状態を表現する型を用意してStateと一緒に使うことで、冷蔵庫が適切に使われているか、最後にどんな状態になっているかを調べてくれるようなinterpretRefもつくってみよう。

data Refregirator = Opened | Closed | Broken deriving (Show, Eq)

interpretRef :: Free MyAST () -> State Refregirator ()
interpretRef f = case f of
  Free (Open r) -> do x <- get
                      case x of Broken -> return ()
                                _ -> put Opened >> interpretRef r
  Free (Elephant _ r) -> do x <- get
                            case x of Opened -> interpretRef r
                                      _ -> put Broken
  Free (Close r) -> do x <- get
                       case x of Broken -> return ()
                                 _ -> put Closed >> interpretRef r
  Free End -> return ()
  Pure _ -> return ()

閉まっているのに象を入れると壊れてしまうことにした。Stateはputやget、stateで明示的に状態を操作できるMonadで、execStateに完成したStateと初期状態をあたえれば、最後の状態を返してくれる。

> execState (interpretRef program1) Closed
Opened
> execState (interpretRef program2) Closed
Broken

DSLのコードはそのままに、IOから切り離して、インタプリタを取り替えて色々なテストや検証もできそうだ。

動的な制御と評価時の環境の利用

分岐する木の構築と評価時の動的な制御

ここまでは構文木と言いながら、リストに近いものしか扱ってこなかった。プログラムは1列に並んだ命令列のような構造を持っていた。それというのも言語内DSLではその言語の制御構文がそのまま使えるので、制御フローが静的に決定できるなら、構文木の構築時に分岐のない1つの命令列を選択できるためである。(分岐条件などの値は、インタプリタでの評価時に構文木の構築時点に送り込むこともできる。) 一方で、どうしてもインタプリタでの評価時に制御フロー自体を動的に決定したい時は、構文木を分岐する木として表現したいこともあるかもしれない。

data ArithAST a r = Unary (a -> a) r
                  | Binary (a -> a -> a) r r
                  | Branch (a -> Bool) r r r
                  | Val a
  deriving Functor

type FA a = Free (ArithAST a) ()

unary :: (a -> a) -> FA a
unary op = liftF $ Unary op ()

binary :: (a -> a -> a) -> FA a -> FA a
binary op l = Free $ Binary op l (Pure ())

branch :: (a -> Bool) -> FA a -> FA a -> FA a
branch cp t f = Free $ Branch cp t f (Pure ())

val :: a -> FA a
val x = liftF $ Val x

ここでは簡単のため、単項演算と2項演算と条件分岐を含む1つの式を表現する構文木をつくり、Branchの制御フローを変えたインタプリタを用意してみる。

type ASTF = Free (ArithAST Double) ()

interpretArith :: ASTF -> Double
interpretArith f = case f of
  Free (Unary f r)      -> f $ interpretArith r
  Free (Binary f r l)   -> f (interpretArith r) (interpretArith l)
  Free (Branch f l r c) -> if f (interpretArith c)
                             then interpretArith l else interpretArith r
  Free (Val x)          -> x
  Pure _                -> 0

interpretArithBr :: ASTF -> Double
interpretArithBr f = case f of
-- same as interpretArith ...
  Free (Branch f c l r) -> if f (interpretArith c)
                             then interpretArith l else interpretArith r
-- same as interpretArith ...

コードを書いてそれぞれのインタプリタで評価してみる。ノードに演算子があり、木が下向きに延びているので、計算の流れは下から上になる。

-- ((*3) (-2))
calc1 :: ASTF
calc1 = do
  unary (*3)
  val (-2)

-- (log (max calc1 ((**2.5) 0.5)))
calc2 :: ASTF
calc2 = do
  unary log
  binary max calc1
  unary (**2.5)
  val 0.5

-- if (>0) 1.0 then calc1 else calc2 -- interpretArith
-- if (>0) calc1 then calc2 else 1.0 -- interpretArithBr
calc3 :: ASTF
calc3 = do
  branch (>0) calc1 calc2
  val 1.0

> interpretArith calc3
0.9535708819095106
> interpretArithBr calc3
-4.75415180475803e-2

評価時の環境からの変数束縛

それは入れ子の無名関数でしかないが、do記法の変数束縛記法は直感的で便利だ。ASTを表現するデータ型 (ここではStAST) で関数のフィールド (これがfmapされるようにする) を使えば、インタプリタでの評価時の環境から構文木の (実際の) 構築時点に値を送り込み、Freeのdoの中で変数束縛できる。DSLの中で、例えば評価時にIOから取り出した値を利用したり、StateのようにStASTの文脈から取り出した値を利用したりできる。

data StAST s cs r = Read (s -> r)
                  | Const (cs -> r)
                  | Write s r
  deriving Functor

readS :: Free (StAST s cs) s
readS = liftF $ Read id
-- (readS = Free $ Read (\x -> Pure x))

constS :: Free (StAST s cs) cs
constS = liftF $ Const id

writeS :: s -> Free (StAST s cs) ()
writeS x = liftF $ Write x ()

StASTに変化する状態を表現するフィールドsと、変化しない定数を表現するフィールドcsを作り、それぞれを操作するためのFreeのアクションを用意する。

interpretStAST :: Free (StAST s cs) a -> s -> cs -> s
interpretStAST f x c = case f of
  Free (Read g)    -> interpretStAST (g x) x c
  Free (Const h)   -> interpretStAST (h c) x c
  Free (Write y r) -> interpretStAST r y c
  Pure _           -> x

execStateに相当するインタプリタを作成する。FreeのプログラムとStASTの初期状態をあたえると、最後の状態を返してくれる。プログラムを書いて評価してみよう。

sproc1 :: Free (StAST Double Double) ()
sproc1 = do
  x <- readS
  if x < 0.0 then constS >>= writeS
    else writeS $ sqrt x
  y <- readS
  writeS (x + y)

> interpretStAST sproc1 2.0 0.0
3.414213562373095
> interpretStAST sproc1 (-2.0) 1.0
-1.0

こうした方法を使えば、StASTの状態をDSLの中から確認して制御を変えたりもできる。

まとめと参考文献

Freeを使うと、DSL構文木構築とその評価を分離できる。同じDSLのコードに対して複数のインタプリタを使い分けられるので、例えばIOから切り離してテストや検証などにも利用できる。Freeの中では評価時の状態を制御に組み込めるし、どうしても必要なら制御フロー自体を評価時に選ぶこともできる。Freeの便利さや使い方については、より進んだ内容についてもたくさんの良い記事が書かれているので、参照されたい。

free: Monads for free
Haskell for all: Why free monads matter
Haskell for all: Purify code using free monads
Dave Laing - Free for DSLs, cofree for interpreters
Free and Freer Monads

Haskellの代数的データ型

多相性の実現

Haskellでは型クラスによってアドホック多相が実現でき、型変数によってパラメータ多相が実現できる。また、部分型付けは型変数に対して型クラス制約を課すことで近しいことが実現できる。個人的な印象では、Haskellの型システムは強い型付けという意味で安全で、静的な強い型付けにも関わらず簡単に多相性を実現できるという意味で高機能だと思う。多相再帰も簡単に書けて便利なことがあった。

アドホック多相は同じ名前の関数を引数の型ごとに使い分ける抽象化。例えばDoubleでもIntegerでも(+)が使える。また、パラメータ多相は型によらない共通の処理を行う抽象化。例えばmapはどの型のリストに対しても使える。部分型付けは (オブジェクト指向言語で) 型の派生関係を利用して、共通して持つ処理を実際の型ごとに使い分ける。Haskellでは完全に同等な機能はないが、型クラス制約を付けた型変数を利用することに近い。また、GHC存在量化拡張または一般化代数的データ型拡張を利用すれば、部分型付けを実現できる。
ポリモーフィズム - Wikipedia

目次

型変数を持つ代数的データ型

Haskellの代数的データ型では、要素を列挙したり、それぞれの要素に幾つかのフィールドを持たせたりできる。フィールドは再帰的にも定義できる。また、フィールドには名前を付けることもできる。

data (型名) = (値コンストラクタ名(型名と同じでも良い)) (フィールドの型名 ..) | ..
data Type0 = A | B | C
data Type1 = Recursive Integer String Type1 | Terminal
data Type2 = Type2 { filed1 :: Type0, field2 :: Integer -> Integer }

代数的データ型にも型変数を持たせることができて、例えばリストや木のような構造を定義して、中に入っている型に依らない抽象的な処理を記述することができる。

data (型コンストラクタ名) (型変数名 ..) = (値コンストラクタ名) (フィールドの型名(型変数名を含む) ..) | ..
data Tensor a = Tensor [Tensor a] | Elem a deriving (Eq)

instance Show a => Show (Tensor a) where
  show (Elem x) = show x
  show (Tensor []) = "[]"
  show (Tensor xs) = "[" ++ init (concatMap ((++ ",") . show) xs) ++ "]"

instance Functor Tensor where
  fmap f (Elem x) = Elem (f x)
  fmap f (Tensor xs) = Tensor $ map (fmap f) xs

型変数には自分で定義した代数的データ型でも関数でも何の型でも入れられるので、それをコンテナと解釈するか、何かの文脈と解釈するか、可能性は自由に広がっている。型変数を持つデータ型は、その型変数に入るHaskellの任意の型に対して、ラベルを付けているとも言える。そこで、こうしたラベル付きの型に対して、関数を適用してラベル付きのまま計算していく方法を考えてみたい。まず取り上げたいのが、そのデータ型の中身に関数を適用できるFunctorと、そのデータ型の中に入っている関数に、そのデータ型の中に入っている引数を適用していけるApplicativeという型クラス。(型クラスは、そのインスタンスとしたデータ型が指定の関数を持つことを要請してくれる。)

FunctorとApplicative Functor

データ型の中身に関数を適用して、別のものに写せることを要請する型クラスがFunctor。Functorのインスタンスはfmap関数を実装し、その関手としての意味からFunctor則を満たすことが要請されている。(fmap id == id (恒等変換で不変)、fmap f . fmap g == fmap (f . g) (関数合成の保存) Data.Functor)

class Functor f where
  fmap :: (a -> b) -> f a -> f b

cite #source

例としてTypeというデータ型を定義し、TypeをFunctorのインスタンスにし、fmapを実装してみる。

data Type a = Content a | Empty deriving (Show)

instance Functor Type where
  fmap f (Content x) = Content (f x)
  fmap _ Empty = Empty

> fmap (^2) (Content 7)
=> Content 49
> fmap (*2) Empty
=> Empty

Typeに入っていても、Haskellの普通の型 (ここではTypeの外の世界の型という意味。Typeも普通の型。) の値であるかのように、Typeの中身にfmapで普通の関数を適用できるようになった。いったい何が嬉しいのかわからないかもしれないが、例えばリストはFunctorでリストのfmapはmapである。そのデータ型の意味に応じて適切に定義されたfmapによって、データ型の中身を普通のHaskellの関数で写せることは役に立ちそうだ。

Functorでは、fmapに渡す関数はTypeの中から1つの引数しか取れないが、複数のTypeの中身を引数としてあたえたい時はどうしたら良いだろうか。Typeの中に関数を入れて、次に渡されたTypeの中身を部分適用した関数を、Typeに入れて返せるような仕組みを作れば良さそうである。こうした仕組みはApplicativeという型クラスが要請している。(ApplicativeのインスタンスはFunctorのインスタンスで、Applicative則を満たす必要がある。Control.Applicative)

class Functor f => Applicative f where
  pure :: a -> f a
  (<*>) :: f (a -> b) -> f a -> f b

cite #source

instance Applicative Type where
  pure = Content
  Content f <*> Content x = Content (f x)
  _ <*> _ = Empty

> pure (^) <*> Content 2 <*> Content 4
=> Content 16
> mod <$> Content 14 <*> Content 3
=> Content 2

pureは引数をTypeの中に入れ、<*>はTypeの中に入っている関数に次のTypeの中に入っている値を適用する。<$>はfmapの別名。Typeの中に部分適用された関数が入れば続けて適用を繰り返せる。Typeの中に入っていても、Haskellの普通の型の値のように、複数の引数を取る普通の関数に適用できるようになった。

FunctorやApplicative Functorはデータ型の中身を、そのデータ型が意味する文脈に応じて、普通の関数で操作するという意味合いが強いと思うが、Haskellの代数的データ型の表現力はそれに留まらず、文脈を付与されたデータ型に対して、その文脈において意味のある計算を自由に定義し、これを連ねていくことを記述できる。

MonadMonad変換子

MonadインスタンスはApplicativeのインスタンスであって、Monad則をみたす必要がある。(Control.Monad)

class Applicative m => Monad m where
  (>>=) :: forall a b. m a -> (a -> m b) -> m b
  
  (>>) :: forall a b. m a -> m b -> m b
  m >> k = m >>= \_ -> k
  
  return :: a -> m a
  return = pure

cite #source

forall a bはaとbが型変数であることを明示的に宣言しているがここでは省略できる。Monadが要請するのはbindと呼ばれる演算子>>=の実装で、Monadの中身を取り出し、あたえられた関数 (a -> m b) に適用し、Monadとして返す。>>=にあたえる関数はMonadアクションと呼ばれている。>>は左の値を利用しないbindで (文脈の処理は行われる)、returnはpureと同じもの (つまり値を中に入れるだけで、文脈は操作しない)。

これだけでは何のことかわからないかもしれないが、Monadは元々はプログラムを純粋な関数によって構成しながらも、文脈 (状態) を伴った計算をすっきり記述するために導入されている。これを実現しようとした時、まずは文脈を引数で渡して行くことが考えられるが、常に引数で持ち回すのは煩雑だし、文脈を明示的に意識せず表の計算や計算の構成に集中したいこともある。そこで文脈を型に入れて、背景で自動的に処理される2つ目の計算ラインに移してしまったのがMonadだ。

FunctorやApplicativeとの違いは、Monadでは>>=によって背景で自動的に文脈の処理が行われることや、その時、アクションが文脈に働きかけることができる点にある。適切なアクションによる値の操作と文脈の供給を用意しておけば、文脈の処理はMonadにまかせることができ、値の計算や、アクションの組み合わせという計算の構成に集中することができる。つまり、>>=を定義しMonadアクションを用意することは、Monadの文脈を処理する計算を構成するためのDSLを用意することに近しい。

Typeを拡張して、Emptyが出てこない限り、計算の点数を数えるような計算という意味をTypeに追加してみる。

data Type a = Content a Int | Empty deriving (Show)

instance Functor Type where
  fmap f (Content x n) = Content (f x) n
  fmap _ Empty = Empty

instance Applicative Type where
  pure x = Content x 0
  (Content f n) <*> (Content x m) = Content (f x) (n + m)
  _ <*> _ = Empty

instance Monad Type where
  Empty >>= f = Empty
  (Content x n) >>= f = case f x of
                          Empty -> Empty
                          (Content y m) -> Content y (n + m)

このMonadでの文脈の処理として、>>=で中身を取り出しアクションにあたえて、その後ろで点数を足し合わせていく。>>では点数の加算のみが行われる。また、1度でもEmptyが出て来たら、常にEmptyになる。

アクションを定義して、計算しながら掛け算の数を数えてみる。

square :: Num a => a -> Type a
square n = Content (n * n) 1

fact :: Integral a => a -> Type a
fact n = Content (f n) (fromIntegral n - 1)
  where f m | m < 2 = 1
            | otherwise = m * f (m - 1)

> square 2 >>= fact >>= square
=> Content 576 5

addCount :: Int -> Type ()
addCount n = Content () n

calc :: Type Integer
calc = do
  x <- square 2
  y <- fact x
  addCount 3
  square y

> calc
=> Content 576 8

do記法はMonad計算を直感的に書くための構文糖衣。各式はそのMonadの型で、>>で結ばれていると見なせる。また、Monadの中身は<-で取り出して変数に束縛できる。 (その場合は、>>=で値を取り出し、無名のアクションの引数として変数に束縛しているような状態。) もっとたくさんのアクションが用意され、calcのdo記法の中身のような計算が大規模になって複雑になっても、点数の計算と失敗の記録はTypeが後ろで自動的に行ってくれるので、このMonadを使う人は、ただ値の計算に集中することができる。

こんな風に、自分のデータ型にあたえたい計算の意味に応じて、適切に>>=とアクションを定義すれば、文脈に沿った処理が簡単に書けるようになる。関数合成演算子 (.) のように、Monad計算におけるアクションの合成を行う演算子 (f >=> g = \x -> f x >>= g) を考えると、Monad則は、>=>による計算の合成において、returnが左単位元かつ右単位元となることと、>=>が結合法則をみたすことを意味している。(つまり、returnは文脈を変更せず、連続する>=>はどの2項から計算を合成しても結果が変わらない必要がある。) このように考えると、Monadをつくることは、単位元を持ち結合法則を満たす新しい計算の体系をつくっているとも捉えられる。自然数の掛け算で類推するなら、自然数Monadアクションで、*が>=> 、1がreturnになる。

Haskellのライブラリには、状態のStateやReaderやWriter、入出力のIO、失敗可能性のMaybeやEither、DSL構文木処理のFreeなど便利なMonadが色々用意されていて、必要に応じてMonad変換子で複数を組み合わせて使うこともできる。

import Control.Monad.Trans.Maybe
import Control.Monad.Trans.Class (lift)
import Control.Monad (guard)
import Text.Read (readMaybe)

readPositiveInteger :: MaybeT IO Integer
readPositiveInteger = do
  l <- lift getLine
  x <- (MaybeT . return) (readMaybe l)
  guard (x > 0)
  lift $ putStrLn "thank you for correct value!"
  return x

> runMaybeT readPositiveInteger
hello
=> Nothing
> runMaybeT readPositiveInteger
1.0
=> Nothing
> runMaybeT readPositiveInteger
1
thank you for correct value!
=> Just 1

MaybeTの中にIOが入っているが、MaybeTのdoの中でも、liftを使えば内側のIOの関数を使うことができる。guardは零元を持つ2項演算を抽象化したMonadPlus (零元はmzero、演算はmplus) の関数で、引数がFalseの時、零元になる。MaybeTでの実装では、MaybeT (return Nothing) が入る。ちなみにMonadPlusのApplicative版にAlternativeがあり、また単位元を持つ2項演算を抽象化したMonoid (単位元はmempty、演算はmappend) などもある。特定のMonadに依らない、抽象化されたMonad計算を用意したい時にはこうした型クラスも使える。Monadが便利に使われている身近な例としては、リストの内包表記もリストMonadのdo記法だったりする。Monadを使うと純粋な関数のみでも表現力の高い言語内DSLを作ることができるので、そのうちそうした記事も書いてみたい。

実用的なMonad変換子の使い方についてはこちらの記事がわかりやすかった。
http://dev.stephendiehl.com/hask/#monad-transformers

MonadHaskellの文法によって特殊な地位をあたえられている訳ではなく (構文糖衣doはあるが) 、こうした記述能力は型変数を持つHaskellの代数的データ型と型クラスの能力によって自然に実現されている。

メモ

Haskellモナドの話には圏論の言葉がよく出てくるのですが、ざっくりとした理解しかできていません。自分の理解をまとめるためにメモを残します。間違っているかもしれないので、鵜呑みにはしないでください。圏論はなんとなく集合論をさらに抽象化したような印象です。

  • 圏は対象と、対象から対象への射を持った構造
  • 関手はある圏ともう1つの圏の間で、対象を対象に写し、射を射に写す
  • その際、恒等な射と、射の結合法則を保存しなければいけない

これは、HaskellのFunctorにあてはめると以下のようになります。

  • Haskellのすべての型の値を対象として、関数を射とすることで、圏にできる
  • 関手はFunctorの値コンストラクタとfmapで、これで写すことで、Haskellのすべての型の圏はFunctorの圏に対応させることができる
  • 関手で写す前後の、恒等な射と、射の結合法則の保存は、Functor則そのもの
fmap id ≡ id
fmap (p . q) ≡ (fmap p) . (fmap q)

HaskellMonadについてはもう少し複雑で、Haskellのすべての型の圏からMonadの圏への関手 (MonadHaskellの型なので、自分から自分への自己関手になる) が、自然変換を射とすることでまた圏を為し、この関手の圏において、モノイドの性質を持つような対象であるMonadの値コンストラクタ (と自然変換の組) がモナドです。

  • Haskellのすべての型の値を対象として、関数を射とすることで、圏にできる
  • Monadの値コンストラクタを関手とすることで、Haskellのすべての型の圏はMonadの圏に対応させることができる
  • この関手であるMonad (の値コンストラクタ, 以下では省略) に対して、自然変換を射として導入することで、圏にできる
  • 自然変換は構造を保存したままMonadMonadに写すような変換で、ここでは恒等関手をMonadに写すような変換returnと、Monadによって写したMonadを、Monadに写すような変換join
return :: a -> m a -- 自身に対応させる (何もしない) 関手を、mに写す関手に変換していると捉えられる
join :: m (m a) -> m a
return x >>= f ≡ f x
m >>= return ≡ m
(m >>= f) >>= g ≡ m >>= (\x -> f x >>= g)
  • joinはbindと以下のように変換できるので、bindをつなげていくことは、Monadの適用を繰り返していくことに等しくなる。
m >>= f ≡ join $ fmap f m

まとめるとMonadの値コンストラクタはHaskellのすべての型の圏からMonadの圏への関手で、この関手自体が自然変換returnとjoinを射として圏をなし、モノイドの性質を持っています。モノイドなのでどこから演算を始めても結果は変わりません。

意識していませんでしたが、ここでは、型クラスFunctorやMonadインスタンスである型の値が対象である圏をそれぞれFunctorの圏やMonadの圏と呼んでいます。値コンストラクタについても同様です。

自分の語彙の使用は正確ではないかもしれない…。もっと正確で詳細な内容が知りたい方はこちらなどをご参照ください。
Haskell/圏論 - Wikibooks

Haskellで数値計算

普段はJavaFortran数値計算をしているのだが、使える範囲でHaskellも使ってみたくなった。宣言的に書けるのでコードが簡潔になることや、副作用が厳格に管理されている上に、静的な型付けがしっかりしていてコンパイラが賢いので、バグが見つけやすいことなどがうれしい。ただ、自分で管理できるワークステーションならともかく、スパコンGHCを入れて多次元の流体計算などしたいとなると、どうだろうか。Fortranよりも計算資源を効率的に使えていると説明できないといけないかもしれない。一方で、MathematicaPythonでやっているような計算の一部は、Haskellでより柔軟に速くできることも多そうだ。ここでは、よく使われる簡単なアルゴリズムを書いてみた。

目次

準備と方針

処理系はGHC 8.0.2を使っている。また、行列計算のためにhmatrixパッケージを入れた。(ここではcabalを直接使っていますが、状況に応じてStackなどを使ってください。) 行列分解や連立線形方程式を解くアルゴリズムもいずれ紹介してみたい。

> cabal update
> cabal install hmatrix

ここでは、そこまで速度やメモリにシビアではない状況を考えている。速度やメモリを気にする場合は、リストの代わりにData.Vector.UnboxedのVector (参照ではなく値の配列で、要素を辿るようなライブラリ関数はFusionが単一のループにしてくれたりするらしい) を使ったりすると良い。使い勝手もリストとそれほど変わらない。
Numeric Haskell: A Vector Tutorial - HaskellWiki
vectorパッケージはhmatrixを入れると付いてくる気がするが、個別にインストールすることもできる。

> cabal update
> cabal install vector

また、アルゴリズムの選択は何より性能に影響をあたえると思うが、ここでは見通しを優先して逆行列を計算するライブラリ関数を素直に使っていたりするので、その辺りは必要に応じて適宜修正してほしい。hmatrixの逆行列関数はLU分解で逆行列を計算しているので、1000  \times 1000程度までの行列を数十回計算する程度なら問題なさそうではある。並列計算についてもそのうち紹介してみたい。

4段4次の古典的Runge-Kutta法

まずは連立常微分方程式 (ODE) を解いてみる。4段4次の古典的Runge-Kutta法は簡単で精度や速度もそれなりなので、よく使われている。ただ陽解法なので、問題によっては適さないこともある。その場合、hを小さくしたり可変にしたりしてみて、それでもだめなら、ステップごとに非線形方程式を解くことになるが、1段2次から3段6次程度までの陰的Runge-Kutta法なども使えるかもしれない。

import Data.List

kf :: Double -> Double -> [[Double] -> Double] -> [Double] -> Double -> [Double] -> [Double]
kf tc xc fs (t:xs) h ks = map ($ rps) fs
  where rps = st:sx
        st = t + tc*h
        sx = zipWith (\x k -> x + xc*h*k) xs ks

rungekutta_44 :: [[Double] -> Double] -> Double -> Double -> [Double] -> IO ()
rungekutta_44 dfs h lim ps@(t:xs) = do
  print ps
  if t > lim then return ()
  else rungekutta_44 dfs h lim ((t + h):ss)
    where k1s = map ($ ps) dfs
          k2s = kf 0.5 0.5 dfs ps h k1s
          k3s = kf 0.5 0.5 dfs ps h k2s
          k4s = kf 1.0 1.0 dfs ps h k3s
          ss = zipWith5 sx xs k1s k2s k3s k4s
          sx x k1 k2 k3 k4 = x + (h/6.0)*(k1 + 2.0*k2 + 2.0*k3 + k4)

高階微分の方程式を解く場合、新しい変数を導入して1階微分の方程式を連立させることになるので、これを自動化してくれる関数も書いてみた。

rungekutta_44_s :: ([Double] -> Double) -> Double -> Double -> [Double] -> IO ()
rungekutta_44_s df h lim ps = rungekutta_44 (dfn ++ [df]) h lim ps
  where dfn = map dndt [0..(length ps - 3)]
        dndt n = (!! (n + 2))

簡単な方程式でテストしてみる。方程式や初期値はリストであたえているので、数や順番は自分で気をつけないといけない。ライブラリにする場合は別の方法を考えた方がいいかもしれない。

> -- 調和振動
> rungekutta_44 [\[t,x,v] -> v,\[t,x,v] -> -x] 0.1 100.0 [0.0,1.0,0.0]
> -- 減衰振動
> rungekutta_44_s (\[t,x,v] -> -x - 0.1*v) 0.1 100.0 [0.0,1.0,0.0]

f:id:lqtmirage:20170529112439p:plainf:id:lqtmirage:20170529232044p:plain
もう少し面白いローレンツアトラクタなども描いてみる。

> let l_dxdt [t,x,y,z] = -10.0*x + 10.0*y
> let l_dydt [t,x,y,z] = -x*z + 28.0*x - y
> let l_dzdt [t,x,y,z] = x*y - 8.0/3.0 * z
> rungekutta_44 [l_dxdt,l_dydt,l_dzdt] 0.01 100.0 [0.0,1.0,1.0,1.0]

f:id:lqtmirage:20170529232335p:plain

Newton-Raphson法

人間誰しも連立非線形方程式を解きたくなることがあるはずだ。Newton-Raphson法は微分係数を使って解に近づいていく収束が速いアルゴリズムで、適切な関数に対して適切な初期値をあたえれば、近くにある解を1つ見つけてくれる。ただ、収束しないこともあるし、自動的にすべての解を見つけてくれるわけではない。
まずは1変数の場合。微分係数は差分で近似している。deltaを小さくし過ぎると数値誤差で結果が悪化することもあるので注意が必要。数値微分で精度が必要な場合は中心差分や高階差分なども検討されたい。

delta = 1.0e-9 :: Double

derivative :: (Double -> Double) -> Double -> Double
derivative f x = (f (x + delta) - f x) / delta

newton :: Int -> Double -> (Double -> Double) -> Double -> Maybe Double
newton lim err f x
  | lim <= 0 = Nothing
  | abs fx < err = Just x
  | otherwise = newton (lim - 1) err f sx
  where sx = x - (fx / dfdx)
        fx = f x
        dfdx = derivative f x

簡単な関数でテストしてみる。発散したらすぐに打ち切るべきか。

> newton 1000 1.0e-9 (\x -> x^2 - 4.0) 4.0
=> Just 2.0000000000000053
> newton 1000 1.0e-9 (\x -> x^2 - 4.0) 0.0
=> Nothing

多変数の場合に拡張する。微分係数をJacobi行列 (偏微分係数を並べた行列) に置き換えて、関数とパラメータをベクトル化しただけで、ほとんど同じに書ける。速度が必要な場合はリストをData.Vector.UnboxedのVectorやData.Vector.Unboxed.MutableのMVectorに置き換えて、ランダムアクセスやin placeな書き換えも使うと良いかもしれない。

import qualified Numeric.LinearAlgebra as L

delta = 1.0e-9 :: Double

derivative_p :: ([Double] -> Double) -> [Double] -> [Double]
derivative_p f xs = zipWith dfdx xs [0..]
  where dfdx x i = (f (replace xs i (x + delta)) - f xs) / delta
        replace xs i x = take i xs ++ [x] ++ drop (i + 1) xs

newton_m :: Int -> Double -> [[Double] -> Double] -> [Double] -> Maybe [Double]
newton_m lim err fs xs
  | lim <= 0 = Nothing
  | L.norm_2 fx < err = Just xs
  | otherwise = newton_m (lim - 1) err fs sx
  where sx = L.toList $ L.vector xs - (L.inv j L.#> fx)
        fx = L.vector $ map ($ xs) fs
        j = L.fromLists $ map (flip derivative_p xs) fs

適当な関数でテストしてみる。

> let vf0 [x,y,z,t] = 2.0*x^2 - y^3 + z - exp t
> let vf1 [x,y,z,t] = 3.0*x - 6.0*y^5 - 2.0*z*t
> let vf2 [x,y,z,t] = 2.0*x^2 - y - 2.0*z*t
> let vf3 [x,y,z,t] = x^4 + y*z - 2.0*t
> newton_m 1000 1.0e-9 [vf0,vf1,vf2,vf3] [1.0,2.0,3.0,4.0]
=> Just [1.066252215687592,0.7772738508955951,0.7861552762455988,0.9517927114657457]
> let a = [1.066252215687592,0.7772738508955951,0.7861552762455988,0.9517927114657457]
> map ($ a) [vf0,vf1,vf2,vf3]
=> [1.3322676295501878e-15,1.5543122344752192e-15,1.9984014443252818e-15,1.5543122344752192e-15]

Double程度の精度で解けていそうだ。適宜deltaやerrを調整する必要があるかもしれない。

Newton-Cotes型定積分

等間隔で関数値があたえられている時に積分値を計算する、Newton-Cotes型と呼ばれる一連のアルゴリズムがある。なかでも2次関数で補間するSimpson法はよく使われている気がする。4次関数で補間するBoole法もためしてみた。

newton_cotes :: [Double] -> Double -> (Double -> Double) -> Int -> Double -> Double -> Double
newton_cotes ces ce f n a b = ((b - a) / (r ^ n)) * ce * ig n a b
  where r = fromIntegral $ length ces - 1
        ig n a b
          | n == 1 = sum $ zipWith (*) (map f cs) ces
          | otherwise = sum $ zipWith (ig (n - 1)) cs (tail cs)
          where cs = map (\i -> a + i * ((b - a) / r)) [0.0..r]

simpson :: (Double -> Double) -> Int -> Double -> Double -> Double
simpson = newton_cotes [1.0,4.0,1.0] (1.0/3.0)

boole :: (Double -> Double) -> Int -> Double -> Double -> Double
boole = newton_cotes [7.0,32.0,12.0,32.0,7.0] (2.0/45.0)

分割位置のリスト生成について、数値誤差で長さが変わらないように冗長に見える方法をとってみた。簡単な関数でテストしてみる。

> simpson (\x -> x^4) 4 0.0 1.0
=> 0.20000203450520831
> boole (\x -> x^4) 2 0.0 1.0
=> 0.2

Boole法のほうが次数は高いが必ず精度が良いとは限らない。

Schemeの継続と環境

前回の記事で、call/cc は呼ばれた時点での継続を取り出してくれるということを書いたのだが、その時点での環境は保存されているのだろうか? 前回の記事に出て来た wayback 関数を少し書き換えて実験してみる。

(define cont #f)

(define (capture x)
  (call/cc (lambda (cc) (set! cont cc) x)))

(define (wayback2 n)
  (let* ((a 0) (x (capture n)))
    (display (list a x))
    (newline)
    (set! a 1)
    (if (> x 0) (cont (- x 1)) x)))

let* は前に束縛された変数を、後の変数束縛で参照したい時に使えるのだが、それゆえに実行順も決められる。ここでは、aに0が束縛された後、xに (capture n) が束縛される。capture はその時点での継続をcontに保存し、引数をそのまま返すような関数。その後、aとxを表示した後、aに1を代入する。これを実行すると、GaucheでもChibi-Schemeでも以下のような結果になった。

> (wayback2 3)
(0 3)
(1 2)
(1 1)
(1 0)
=> 0

継続 cont が呼び出されると、xへの引数の束縛以降の処理が行われるが、aへの0の束縛はそれ以前に終わっているので繰り返されない。しかし、もし継続が環境も保存しているなら、何度呼び出してもaは0に戻るはずだ。だが、結果はそうはならず、(set! a 1) の影響がその後の継続呼び出しにも反映されている。継続は環境を保存してはくれない。

評価順序と継続
ここに詳しく書いてあった。

Schemeと継続

目次

継続とcall/cc

継続とは

Schemeでは継続を手続きとして取り出して利用することができる。継続とはプログラムの (基本的にはトップレベルに戻るまでの) 残りの処理のことで、例えば、以下の [] の部分には、(square 2) をかけて1を足すという処理が残されている。

(+ 1 (* [] (square 2)))

call/ccは継続を取り出してくれる

call/cc (call-with-current-continuation) を使うと、この残りの処理を手続きとして取り出し、この手続きを引数として関数を呼び出してくれる。

> (+ 1 (* (call/cc (lambda (cc) (cc 3))) (square 2)))
=> 13

現在の継続が取り出され、これを引数として (lambda (cc) ...) が呼び出されている。ここではこのccが継続で、それは普通の手続きである。ccが呼ばれると、(call/cc ...) が呼ばれた時点のその全体がccの引数で置き換えられる。先ほどの [] の部分がccの引数で埋められるとも言える。

(lambda (cc) ...) も普通の関数で、普通に呼び出されているので、もしccが呼ばれないなら最後に評価された式の値が返される。

> (+ 1 (* (call/cc (lambda (cc) 0)) (square 2)))
=> 1

継続はファーストクラスのオブジェクト

継続はファーストクラスのオブジェクトなので、変数に束縛したり、関数の引数や返り値として渡すこともできる。

> (define cont #f)
> (+ 1 (* (call/cc (lambda (cc) (set! cont cc) 0)) (square 2)))
=> 1
> (cont 3)
=> 13

時をさかのぼる

さて、継続が面白いのは、(その定義通り) 継続が取り出された時点のプログラムの残りの処理が保存されていることである。

(define cont #f)

(define (wayback n)
  (let ((x (call/cc (lambda (cc) (set! cont cc) n))))
    (begin (display `("wayback" ,x))
           (newline)
           (if (> x 0) (cont (- x 1)) x))))

xへの変数束縛以降の継続がcontに保存されており、1度waybackを抜けた後も、contを使って何度でも呼び直すことができる。

> (wayback 2)
(wayback 2)
(wayback 1)
(wayback 0)
=> 0
> (cont 3)
(wayback 3)
(wayback 2)
(wayback 1)
(wayback 0)
=> 0

まるで時間をさかのぼって、別の未来への処理をやり直しているように見える。継続の保存期間には制限がないので、束縛されている限りは何度でも呼び出せる。

どこでも取り出せる継続

あまり迷う人もいないかもしれないが、クロージャは定義された位置での環境を参照するのに対して、call/ccは呼ばれた時点での継続を取り出す。そこで、こんな手続きを使うこともできる。

(define cont #f)

(define (capture x)
  (call/cc (lambda (cc) (set! cont cc) x)))

(define (func n)
  (cond ((= n 0) (capture 1))
        (else (* n (func (- n 1))))))

captureは呼ばれた時点での継続をcontに束縛し、xを返す。そのため、コード中の好きな式aを (capture a) に置き換えれば、結果を変えずにその位置での継続をcontに束縛できる。

> (func 4)
=> 24
> (cont 5)
=> 120

Schemeではcall/ccを使えば、どこでもその時点の継続を取り出すことができ、それを変数に束縛したり、関数に渡したり返したりして、遠く離れた場所からも利用できる。これはなんだか大変な力に思える。

継続の利用

それでは、具体的には継続は何に使えるのだろうか?ここでは簡単なサンプルコードを示して継続の使い方を紹介してみたいと思う。

大域脱出

継続のわかりやすい使用例は大域脱出だ。継続がなくても同様の処理は書くことができるが、深い位置からの脱出ではコードがより簡潔になることもあるかもしれない。

(define (escape x lst)
  (call/cc (lambda (return)
             (map (lambda (l)
                    (if (= l 0) (return '()) (/ x l)))
                  lst))))

returnにescape関数の最も外側の継続が保存されているので、returnが呼ばれた時点で (残りの処理を飛ばして) その引数が直ちにescape関数の返り値となる。それゆえに、returnはCのような言語のreturn文のように使うことができる。

> (escape 2 '(1 2 3))
=> (2 1 2/3)
> (escape 2 '(0 1 2))
=> ()

処理の一時中断と再開

継続の力がより発揮される例の1つは、好きな位置で一時中断して、後で停止した位置から再開するような処理ではないだろうか。まずは以下のようなコードを考えてみる。

(define resume #f)

(define (resumable n)
  (call/cc (lambda (return)
             (display `("resumable" ,n))
             (newline)
             (call/cc (lambda (cc)
                        (set! resume (lambda () (cc #t)))
                        (return n)))
             (resumable (+ n 1)))))

resumable関数は内側のcall/ccが呼ばれた時点の継続をresumeの中に保存し、returnでnを返している。なので、resumeを呼び出せば、内側の (call/cc ...) が#tで置き換えられ、(resumable (+ n 1)) が呼ばれる。

> (resumable 0)
(resumable 0)
=> 0
> (resume)
(resumable 1)
=> 1
> (resume)
(resumable 2)
=> 2
...

ただ、このresumableはnを返して止まってしまうので、返り値を使うことはできず、例えば以下のようにリストを作ることはできない。

> (import (srfi 1))
> (fold (lambda (i l) (cons (resume) l)) '() (iota 10))
(resumable 3)
=> 3

それというのも、(resume) は前回の resumable 途中の継続に飛んで行ってしまって、呼び出し元の継続には帰って来ないからである。返り値を使いたい場合は、(resume) の中で取り出された継続に帰って来てくれないといけない。そこで、以下のように改良してみる。

(define yield #f)

(define (generator proc)
  (define current proc)
  (define (next)
    (call/cc (lambda (return)
               (set! yield (lambda (x)
                             (call/cc (lambda (cc)
                                        (set! current cc)
                                        (return x)))))
               (current))))
  next)

(define (endless)
  (let loop ((i 0))
    (begin (display `("endless" ,i))
           (newline)
           (yield i)
           (loop (+ i 1)))))

generator内のnext関数は、yieldの呼び出し時の継続をcurrentに保存し、yieldの引数xをreturnに返すような関数をyieldにセットし、currentを呼び出す。そのため、currentの呼び出し先でyieldが呼ばれると、currentはyieldへ戻るための継続に更新され、nextはxを返す。こうして、nextは呼ばれる度に、yieldまでの処理を実行し、xを返し、一時中断する。

> (define gen (generator endless))
> (gen)
(endless 0)
=> 0
> (gen)
(endless 1)
=> 1
...
> (import (srfi 1))
> (fold (lambda (i l) (cons (gen) l)) '() (iota 10))
=> (11 10 9 8 7 6 5 4 3 2)

こうした技術を使うと、重い処理や無限に続く処理などを適当なところで止めておいて、必要に応じて再開することができる。また、もっと対称的に2つの処理がお互いの継続を呼び合うような方法を考えることもできる。これらは実行順や範囲を意識的に制御される平行処理のようなものを実現する。

(define (coroutine-1 cont)
  (let loop ((i 3))
    (if (>= i 0)
      (begin (display `("c-1" ,i))
             (newline)
             (set! cont (call/cc cont))
             (loop (- i 1))))))

(define (coroutine-2 cont)
  (let loop ((i 0))
    (if (< i 4)
      (begin (display `("c-2" ,i))
             (newline)
             (set! cont (call/cc cont))
             (loop (+ i 1))))))

それぞれの (set! cont (call/cc cont)) ではその時点での自身の継続を引数に相手の継続を呼び出し、帰ってきた相手の継続をcontに束縛している。

> (coroutine-1 coroutine-2)
(c-1 3)
(c-2 0)
(c-1 2)
(c-2 1)
(c-1 1)
(c-2 2)
(c-1 0)
(c-2 3)

まとめと参考文献

Schemeではファーストクラスの継続を簡単に取り出して利用することができる。継続を使うと、プログラムを任意の時点で中断したり、再開したりすることができる。そのため、大域脱出や例外処理、ジェネレータやコルーチン、探索時のバックトラッキングなど、複雑な制御を比較的簡単に実現することができる。

ここまで継続について学んだことをまとめてみた。あまり実用的な例を紹介することはできなかったが、もし自分のような勉強中の方の参考になったらさいわいだ。最後に、参考にしたサイトや本を紹介しておきたい。ここで紹介したコードの多くは、これらの参考文献をもとに書かれている。より進んだ内容や、より正確な内容が知りたい方は参照してほしい。

call-with-current-continuation
Scheme/継続の種類と利用例 - Wikibooks
使いたい人のための継続入門
なんでも継続
独習 Scheme 三週間 Teach Yourself Scheme in Fixnum Days (Chapter 13, 14)
もうひとつの Scheme 入門 (16, Appendix 3)
On Lisp (20-23)
R7RS 日本語訳 (報告書 P49-52)

Schemeによる記号処理入門 | 森北出版株式会社
O'Reilly Japan - プログラミングGauche