\( \def\vector#1{\boldsymbol{#1}} \) \( \newcommand{\argmax}{\mathop{\rm argmax}\limits} \) \( \newcommand{\argmin}{\mathop{\rm argmin}\limits} \)

Scala による数値演算アルゴリズム

2017年9月27日

序文

学生の頃に Numerical Recipes in C という多くの実数演算アルゴリズムの載った書籍があった。英語版は 3rd Ed. が出ているようだが日本語版は出ていない。アルゴリズム辞典や数値演算ライブラリを作る気はないのだが系統としてはあの本の方向で行ければ良いと思っている。

基本的に数式の解法や実装コードのフラグメント置き場である。いまさら実数演算という気がしなくもないが、実績のある数値演算ライブラリが GPU による並列処理 (ベクトル化) に対応していなかったり、独自の減衰関数を実装するような局面が最近は増えてきているように思われる。

実数演算のテクニック

誤差を削減し精度を上げるためのテクニック。

実数演算では誤差を回避するために定数を正確な数値リテラルで記述しなければならないことが多い。例えば \(\sqrt{2\pi}\) の近似は 2.50662827463100050 だが、math.sqrt(2*math.Pi) と書けば 2.5066282746310002 となり計算精度は低下する。
実数演算は単純な四則演算であっても評価順序によって結果が変わる。例えば (1e-308*1e10)/1e10 は 1.0E-308 だが (1e-308/1e10)*1e10 は 9.999987484956E-309 である。除算は割り切れないときに仮数部を全て使い切るため誤差を出しやすい。したがって (オーバーフローに注意して) 乗算を先に行い除算を後に行う方針は誤差を軽減することができる。

実装のための言語機能

データ型

Scala は Java と同等のデータ型を利用できる。整数は 8bit の Byte、16bit の Short、32bit の Int、64bit の Long でいずれも符号付き整数である。科学技術計算向けの浮動小数点は 32bit の Float と 64bit の Double を利用できる。

また Scala の特徴として多倍長整数の BigInt を利用できる。これは Java の BigInteger のラッパーで特別な何かを import しなくても利用できる。また四則演算の演算子オーバーロードが定義されているため数値データとして使用できる。

関数

Scala は第一級関数の言語機能を持つ。

def add(a:Int, b:Int):Int = a + b
add(1, 2)    // ==> 3

Java 8 でクロージャーが導入されたが Scala の関数もオブジェクトとして実装されている。上記の add() は以下のように変数に代入することが出来る。

val add = { (a:Int, b:Int) => a + b }
add(1, 2)    // ==> 3

ループと再帰呼び出し

関数型言語でのループは一般的な手続き型/オブジェクト指向言語での繰り返し処理とは根本のイメージが少し異なる。有限個の集合に対して iterative に処理を行うことが多い。モナドっぽく集合 \(\vector{A}\) から \(\vector{B}\) への射

for(i <- 0 until 10) ...

ループカウンターを使ったループが必要

var i = 0
while(i < 10){
  ...
  i += 1
}

@tailrec
def _loop(i:Int):Unit = if(i < 10){
  ...
  _loop(i + 1)
}
_loop(0)

部分適用とカリー化

確率密度関数のようにモデルを固定するためのパラメータとモデルに適用する確率変数とを明確に分けたい場合に部分適用を利用することができる。以下の例では正規分布のパラメータ \(\mu=0, \sigma^2=1\) を部分適用し、引数に確率変数 \(x\) を一つ取る標準正規分布の関数を作成している。

def norm(mu:Double, sigma2:Double)(x:Double):Double = {
  math.exp(- (x - mu) * (x - mu) / 2 * sigma2) / math.sqrt(2 * math.Pi * sigma2)
}
val stdnorm = norm(0, 1) _
stdnorm(0.0)                    // => 0.3989422804014327

Scala での部分適用の使い方は難しくない。関数定義の引数部分を 2 つ以上にすることと、部分適用する際に後方の引数が省略されていることを示すためのアンダースコア _ を付けることだけで良い。

複数の引数を一つづつの引数リストに分解することはカリー化 (currying) と呼ばれる。

トランポリンによる StackOverflow の回避

数値演算アルゴリズムでは再帰的に関数を適用することが多い。

def factorial(n:Int):BigInt = if(n < 2) 1 else n * factorial(n - 1)
factorial(100000)           // StackOverflow

import scala.util.control.TailCalls._
def factorial(n:Int):TailRec[BigInt] = if(n < 2) done(1) else {
  tailcall{ factorial(n - 1) }.map{ x => n * x }
}
factorial(100000).result    // 28242294079603478...

オーバーフローの判定

\(x \times y \gt {\rm max}\) より \(x \gt \frac{\rm max}{y}\) を判定する。また Scala では x.isInfinite で x がオーバーフローを起こしたかを判定することもできる。