はじめに

つまり、数値計算を自分で実装する場合の注意点です。

特に必要がなければ車輪の再発明をせず、実装済みのソフトウェアやライブラリやソースコードを探しましょう。
そうではない、以下のような人向けです。

精度や速度はそこまで必要としない
実装することでアルゴリズムを理解したい
実装済みのプログラムが見つからない、あるいは使用できない

なお、自分用のメモですので、内容についてはお察しください。

覚書

命名規則

変数名や関数名はとにかく数式と完全に一致させる。

プログラムでよく使われるものとは異なる場合があったり、大文字小文字が混ざっていたりするが、気にせずそのままとした方がいい。

なお、記号（ダッシュ記号'やキャレット記号^など）は用途や読みが異なる場合があるが、自分の宗教に合わせる。

必要であれば、あとで変更する。

添え字

数式上は、$0$始まりや$1$始まり、$n$終わりや$n-1$終わりなどある。
プログラムでは$0$始まり$n-1$終わりが一般的だが、上記のため必ずしもそうはならないことに注意する。

一般的にプログラムでは、ループ回数が明確になるような書き方をすることが多い。
そのため、ループする値自体は記述とずれがあることを意識しなければならない。

例えば次のような式

\sum_{k=1}^{n} f(k)

は、C言語のようなfor文では

for (int i = 0; i < n; i++) {
  value += f(i + 1)
}

Pythonでは

for i in range(n):
    value += f(i + 1)

と+ 1が必要であり、少し直観と反する書き方になる。
なお、ループ範囲を$+1$だけずらしてもよいが、Pythonでは

for i in range(1, n + 1):
    value += f(i)

となり、個人的には気持ち悪い書き方だと思う。

なお、MatlabやVBAなどはまた異なるので注意が必要。

割り算

割り算は値のスケールが大きく変化するので十分注意する。
処理に工夫をすることで精度が改善することはよくある。

整数の割り算

例えば、組み合わせの計算

{}_n C_k = \frac{ n! }{ (n - k)! k! }

は必ず整数となるが、分子・分母を別々に計算するのではなく、まとめて計算した方がよい。
つまり例えばPythonでは

value = 1
for i in range(k):
    value *= n - i
    value /= i + 1

としたほうがよい。
なお、C言語などのように、整数を整数で除算すると整数に切り捨てられる場合は、計算順序に十分に注意する。

既存の関数の使用

主に逆正接atanを計算する場合。
JavaScriptにはMath.atan2(y, x)が存在するので、こちらを積極的に使用する。

ログスケール

掛け算と割り算が大量に発生する場合は、ログスケールでの計算を実施してから元のスケールに戻すことも考える。
例えば、Lanczos近似によるガンマ関数$\Gamma(x)$の計算を用いたベータ関数$B(x, y)$の計算においては、$\log \Gamma(x)$の計算が容易なため、ログスケールでの計算が有効である。

B(x, y) = \frac{\Gamma(x)\Gamma(y)}{\Gamma(x + y)} = \exp ( \log \Gamma(x) + \log \Gamma(y) - \log \Gamma(x + y) )

定義済みのもの

言語によっては、頻出の数学的な関数や定数を定義している場合がある。
一通り確認しておくとよい。

例えばJavaScriptでは、

Math.LN2 // 2の自然対数
Math.SQRT1_2 // 1/2の平方根

等がある。

特殊な値

二重階乗の$0!!$など、定義や特定の条件で計算方法が変わる場合は、その値に注意する。

また、逆正規累積分布関数のように定義域が存在する場合は、域外の値について適切に対応する。

多次元変数の演算

大きさの意識

各次元の大きさを意識して実装すると間違いが少なくなる。

例えば、連続分布型HMMで各状態の出力をガウス分布とした場合は、そのガウス分布の分散は $状態数 \times 出力次元数 \times 出力次元数$ の三次元変数となる。
更新式では、更新値が必ずこの大きさになる。

簡略化

行列演算では、実装は数式に対して計算量を削減できる場合がある。

例えば、重み付き最小二乗法の解析解は

\hat \beta = ( X^{\top} W X )^{-1} X^{\top} W y

であるが、この$W$は対角行列であるため、$X^{\top} W$は行列積を計算する必要はなく、要素ごとの積で十分である。

計算順序

計算順序を変えるだけで、計算速度が改善する場合がある。

最小二乗法の解析解は

\hat \beta = ( X^{\top} X )^{-1} X^{\top} y

であるが、$X$の列数よりも$y$の列数の方が少ない場合、数式通りに左から順番に計算するよりも、先に$X^{\top} y$を計算する方が計算量が小さくなる。

逆行列

逆行列は真面目に計算しなくても問題ないことがある。

例えば

a = X^{-1} y

は逆行列を計算する必要はなく、ガウスの消去法などによって効率的に$a$を求めることができる。

数式と完全に一致しなくてよい場合

一部、数式と数値計算プログラムの処理とが一致していなくても問題ない場合がある。

ゼロ除算

除算での割る側の値、あるいは$-1$乗する値が$0$になる場合は、次のような対応を実施するとよい場合が多い。
行列についても逆行列の計算について、正則行列であることが約束されない場合も同様。大抵は対称行列かどうかで判断する。

微少量を足す
行列の場合は、対角成分に微少量を足す
$1$とする（割り算を実施しない）
エラーとする

ほとんどの場合、微少量を足すことで何とかなる。

単位ベクトルや離散確率分布のように、合計が$1$であることが要求される場合は、状況に応じて適切に対応する。

微分値

例えば機械学習のニューラルネットワーク界隈で有名なReLU関数は、数学では原点において微分不可能である。
しかし、ニューラルネットワークの誤差逆伝播時にはその微分値を、$0$あるいは$1$のどちらにしても問題なく動作する。

値の確認

数値計算で一番確認が難しいのは、計算結果が間違っていることである。
この部分はしっかりと確認する。

例えば数式の値や入出力を、高精度計算サイトや他の文献、あるいは既存のライブラリの出力と比較する。
楕円有理関数などのように別の関数を組み合わせてできる関数は、一つ一つ個別に比較するとなお良い。

値が比較確認できない場合は、値の性質を満たすことを確認する。

なお、精度は必要なだけ一致していればよい。
$0$とならなければならない部分が1.0e-16などとなっていることはよくある。

アルゴリズムの選択

実装が面倒な数式には大抵、複数の計算アルゴリズムが存在する。
また、近似値でいい場合は近似アルゴリズムも確認する。

アルゴリズムによって前提条件や計算速度、数値の収束速度が異なる場合が多いので、最終的には前提条件を満たすものの中から速度と精度のバランスを見て選ぶ。

有名なのは固有値計算で、Jacobi法、LU分解法、QR分解法などがあるが、非対象行列においてはQR分解法が比較的良い精度とされるらしい。

言語や実行環境によって速度と精度にばらつきが出るので注意する。