この記事は 数値計算 Advent Calendar 2018 5日目の記事です。
はじめに
数値解析アドベントカレンダーは大学時代を思い出させてくれるので懐かしい気持ちになります。そんな大学からの経歴についてと、本日お話することは下記の記事からの延長になるので、まずは一読をオススメします。
数値解析への想い
さて、本題からズレますが、数値解析アドベントカレンダーが始まってすぐに、私の研究テーマ（実際は手伝い）だったBFGS法のアドベントカレンダーも紹介されていて、テンションがあがります。
参考文献にも担当教授がいたりと、この分野の数値解析分野ならまだまだ負けない気がします💪
ちなみに、BFGS はメモリ効率を高めるための Limited-memory BFGS というアルゴリズムの方がより実践では使われることが多いです。7年以上前に実装した、Limited-memory BFGSのC言語のソースコードを公開しているので、もしよければお使いください。
Logic of Fast Inverse Square Root
前回の記事で Goで高速逆平方根計算アルゴリズム - Fast Inverse Square Root - についてコードを紹介しましたが、今回は少しコードを読み解いていければと思います。
func FastInverseSqrt(x float32) float32 { const threehalfs = float32(1.5) x2 := x * float32(0.5) y := x i := *(*int32)(unsafe.Pointer(&y)) i = 0x5f3759df - i>>1 y = *(*float32)(unsafe.Pointer(&i)) y = y * (threehalfs - (x2 * y * y)) return y }
数値解析の演算回数や誤差
const threehalfs = float32(1.5)
数値解析で研究を行うときは演算回数や情報落ち、桁落ちなどの誤差に気をつけなければなりません。
今回は const threehalfs = float32(1.5) として、 前もって3/2 = 1.5 として除算をするかしないかをコンパイラ依存（定数でなくても文脈から除算を前もってコンパイル時に解釈してくれるコンパイラもあります）をしないようにしたり、わかりやすいようにマジックナンバーではなくて定数名をつけたりして保守性を高めます。研究だと他の人にソースコードを渡すことが多々あるので、ここはリーダブルなコードを心がけましょう。
また、誤差伝播も気をつけなければなりません。
ポインターのキャスティング
x2 := x * float32(0.5) y := x i := *(*int32)(unsafe.Pointer(&y))
32ビット同士の整数と単精度浮動少数点数のポインターキャスティングを行っています。単精度の浮動少数点数の定義は
IEEE 754 での binary32 の定義は次の通りである。
符号ビット: 1ビット 指数部の幅: 8ビット 仮数部の幅: 23ビット
です。この、それぞれの符号、指数部、仮数部の情報を活用したまま32ビットの整数型にしています。
ある意味マジックナンバー 0x5f3759df
i = 0x5f3759df - i>>1 y = *(*float32)(unsafe.Pointer(&i))
Fast inverse square root - Wikipedia
ある意味マジックナンバーについては、First approximation of the resultに記載があります。
y = 1⁄√x を求解するために、logを取ってテイラー展開を行いそれを一次の精度で打ち切ります。これによって得られた近似式に対して、得られた数というのが 0x5f3759df になるのです。
ニュートン法の適用
y = y * (threehalfs - (x2 * y * y))
ここで、ニュートン法の反復を一度行っています。ここで、最大発生する計算誤差についても0.175%程度しか起きないため、そこそこの精度を保ったまま、一度の反復で逆平方根の近似値を得ることができます。これ以上の精度が欲しい場合は同じ計算を再度行うことによって、二回目の反復を実行します。
おわりに
数値解析のアルゴリズムは研究の精度を向上させることやより理解するためにも深掘りは重要になるので、なんとなくではなく、理論を知るのを怠らないようにしていければと思います。
