本当は怖いHPC

HPC屋の趣味＆実益ブログ

数値微分（Juliaで）

Juliaで数値微分を書いてみます。

数値微分の基本

数値微分は、直感的には微小幅 h について関数の値を計算して差分を取れば、それが勾配になります。しかし、数学的にもう少しちゃんと根拠が示されています。

まず、微分を求めたい１変数関数を f、点を x_0とします。そして、微小幅 h > 0 について、関数 f(x + h) を x の周りでテイラー展開します。

{ \displaystyle f(x + h) = f(x) + h f^{(1)}(x) + \frac{h^{2}}{2} f^{(2)}(x) + \frac{h^{3}}{6} f^{(3)}(x) + \dots }

ただし、 \dots は、 h^{4} 以上の項です。同様に、関数 f(x - h) をxの周りでテイラー展開します。

{ \displaystyle f(x - h) = f(x) - h f^{(1)}(x) + \frac{h^{2}}{2} f^{(2)}(x) - \frac{h^{3}}{6} f^{(3)}(x) + \dots }

ここで上２つの式を引き算すると、都合の良いことに h^{2} の項が消えます。

{\displaystyle f(x + h) - f(x - h) = 2 h f^{(1)}(x) + \frac{h^{3}}{6} f^{(3)}(x) + \dots }

h^{3} 以上の項を無視することによって、

{\displaystyle f^{(1)}(x) = \frac{ f(x + h) - f(x - h) }{2h} }

という結果が得られました。これは中心差分の公式と呼ばれます。上の式変形からわかるように、 h^{1}の項まで考慮しているので（計算の途中で消えたけどｗ）、誤差のオーダーが O(h^{2})（二次精度）となります。

Juliaで実装

上式をJuliaで実装してみます。Float32をベースとして使って、 h としてJuliaの組み込み定数であるeps(Float16)を使ってみます。

コード自体は非常に単純です。

function diff(f)
  h = eps(Float16)
  return (x -> (f(x+h) - f(x-h))/(2h))
end

受け取った関数 f の1次の導関数を返します。試しに、いくつか解析的に自明な関数の計算をしてみます。

julia> f = x -> x  # 恒等関数：微分係数は常に1
(anonymous function)

julia> df = diff(f)
(anonymous function)

julia> df(0)
1.0f0

julia> df(1000)
1.0f0

julia> g = x -> x^2 + 3x + 1   # g(x) = x^2 + 3x + 1
(anonymous function)

julia> dg = diff(g)  # 導関数は g'(x) = 2x + 3
(anonymous function)

julia> dg(0)
3.0f0

julia> dg(5)
13.0f0

julia> h = x -> sin(x)  # h(x) = sin(x)
(anonymous function)

julia> dh = diff(h)  # 導関数は h'(x) = cos(x)
(anonymous function)

julia> dh(0)   # cos(0) = 1
0.9999998f0

julia> dh(pi/2)    # cos(π / 2) = 0
0.0

なんとなく良さそうです。

Juliaの16進浮動小数点リテラル

julia

Julia言語のチュートリアルを見ていて、16進浮動小数点リテラルなるものが存在するのを知りました[1]。書式としては、

  • 先頭に0xをつける（通常の16進リテラルと同じ）
  • 指数部の記号として p を使う
  • ただし、指数の底は2

ということで、例を見てみると

julia> 0x1.8p3
12.0

julia> (1.0 + (1/16*8)) * 2^3
12.0

ちなみに、これはJuliaだけの機能ではなくて、Javaにも存在するらしい[2]。

ふむふむなるほど。科学技術計算では使う機会がないこともなさそう…

[1] Integers and Floating-Point Numbers — Julia Language 0.5.1-pre documentation [2] Hexadecimal Floating-Point Literals (Joseph D. Darcy's Oracle Weblog)

最近の気になったニュース 2017/03/06

プロセッサ関係

Ryzen続報

AMDのRyzen CPUの続報がいろいろ出ています。

RyzenをHPCで使おうとしたベンチマーク - Qiita

【レビュー】AMD「Ryzen 7 1800X」はIntelの牙城を崩せるか? ～各種ベンチマークを実施 - PC Watch

【後藤弘茂のWeekly海外ニュース】AMD「Ryzen 7」の半導体チップの姿 - PC Watch

NVIDIA GeForce GTX 1080 Ti発表

NVIDIAがゲーム向けのコンシューマーGPUであるGTX 1080Tiを発表しました。Pascal世代のP102、GDDR5Xメモリ11GB、メモリ帯域が484GB/s、ピーク演算速度は単精度10.36TFlops。とうとう大台の10TFlopsに乗ってきました。

NVIDIA CEO Jen-Hsun Huang Unveils GeForce GTX 1080 Ti, World’s Fastest Gaming GPU | NVIDIA Blog

機械学習／AI関係

Line Announces Smart Assistant Clova, A Potential Alexa Competitor

メッセージングアプリのLINE株式会社が、Clovaという製品を発表しました。AmazonのAlexaと似ている据え置き型の対話型アシスタントです。

夏頃に、スマートフォン上のアプリと、音声端末である「WAVE」が日本と韓国で提供開始、ディスプレイ付きの「FACE」が冬ころの提供ということです。

Amazon Alexaは日本ではまだ提供されていませんし、GoogleのアプリもPixel以外のAndroid端末での提供は始まったばかりです。日本語の認識がどの程度なのか、またAmazonはGoogleに対抗できるのか興味深いです。

その他

AWS S3 の大規模ダウン

Amazon AWS S3サービス（Simple Storage Service）が4時間に渡ってダウンし、S3に依存している多くのサイトが影響を受けたそうです。直接の原因はヒューマンエラーで、しばらくトラブルを大規模なトラブルを経験していなかったことが逆に災いしてサービスの再起動に時間がかかったということです。また、その他の特徴としてはS3のダッシュボード（システムの状態を確認できるモニター）自身が影響をうけて正常動作しなくなったという点です。

毎度のことですが、現代のWebにおけるAmazon AWSの存在の大きさに気付かされます。

米国東部（バージニア北部、US-EAST-1）リージョンで発生した Amazon S3 サービス障害について

AWS、S3の大惨事の原因を公開―ヒューマンエラーが発端だった | TechCrunch Japan

【広告】