#Julia言語 ほとんどマクローリン展開なのですが、DS1～DS6の値をじっと見れば、高次の係数ほど大きくずらしてあって、ぴったりマクローリン展開にはなっていないことが分かります。
この微妙な点も結構面白いはず。
Juliaのリポジトリへのリンク↓
github.com/JuliaLang/juli…
この微妙な点も結構面白いはず。
Juliaのリポジトリへのリンク↓
github.com/JuliaLang/juli…
https://twitter.com/Atsushi_twi/status/1362760013941383172
#Julia言語 返答や引用ツイートを見ると、これをテイラー・マクローリン展開そのものだと誤解している人が多いな。
マクローリン展開ならDS1, DS2が
-1/6=-0.16666…
1/120=0.0083333…
と6,3がずっと続くことになるのですが、そうなっていない。
マクローリン展開そのものだと精度的に損。
マクローリン展開ならDS1, DS2が
-1/6=-0.16666…
1/120=0.0083333…
と6,3がずっと続くことになるのですが、そうなっていない。
マクローリン展開そのものだと精度的に損。
https://twitter.com/atsushi_twi/status/1362747637879894016
#Julia言語 #数楽
浮動小数点数の mod 2π での剰余を求めるのは非常に大変。
Juliaではそのためのコードが別ファイルになっていて、非常に複雑。
rem_pio2_kernel(x::Union{Float32, Float64}) はxをπ/2で割った剰余とxを2πで割った剰余がその何倍かを計算。
github.com/JuliaLang/juli…
浮動小数点数の mod 2π での剰余を求めるのは非常に大変。
Juliaではそのためのコードが別ファイルになっていて、非常に複雑。
rem_pio2_kernel(x::Union{Float32, Float64}) はxをπ/2で割った剰余とxを2πで割った剰余がその何倍かを計算。
github.com/JuliaLang/juli…
https://twitter.com/mkashi/status/1363083689182130184
#Julia言語 mod 2π の浮動小数点数での計算は大変なので sin(x) の実装は大変。
sin(πx) ならば圧倒的に楽。
Juliaのtrig.jlではsinpi(x)も定義されています。
sin(2πx) を計算したい場合にはJuliaでは sinpi(2x) と書いた方がお得です。sin(2πpx) なら sinpi(2p*x)
github.com/JuliaLang/juli…
sin(πx) ならば圧倒的に楽。
Juliaのtrig.jlではsinpi(x)も定義されています。
sin(2πx) を計算したい場合にはJuliaでは sinpi(2x) と書いた方がお得です。sin(2πpx) なら sinpi(2p*x)
github.com/JuliaLang/juli…
#Julia言語
Julia界では「初等及び特殊函数をpure Juliaで実装しておくことには自動微分でメリットがある」とよく言われています(した)。
現実のメリット→既存のライブラリをJulia版で置き換えると計算速度が上がることが多いのと、Julia界は数学強者が多いので、だったらJuliaで書いちゃえと。
Julia界では「初等及び特殊函数をpure Juliaで実装しておくことには自動微分でメリットがある」とよく言われています(した)。
現実のメリット→既存のライブラリをJulia版で置き換えると計算速度が上がることが多いのと、Julia界は数学強者が多いので、だったらJuliaで書いちゃえと。
#Julia言語
Juliaに限らず、特殊函数や高速フーリエ変換などのコードをどのように書くべきか(生成するべきか！)については、
JuliaCon 2019 | Keynote: Professor Steven G. Johnson
を視聴するとよいです。これ、何度も勧めているのですが、めちゃくちゃ面白いです。
Juliaに限らず、特殊函数や高速フーリエ変換などのコードをどのように書くべきか(生成するべきか！)については、
JuliaCon 2019 | Keynote: Professor Steven G. Johnson
を視聴するとよいです。これ、何度も勧めているのですが、めちゃくちゃ面白いです。
#Julia言語 Juliaの周辺の情報をあさると数学の勉強になってしまうことがよくあります。
高速計算の話題ではすぐにSIMDやコンパイラによる最適化などの話題にしたがる人が多いのですが、上で紹介したSteven G. Johnsonさんは全然違います。
高速計算の話題ではすぐにSIMDやコンパイラによる最適化などの話題にしたがる人が多いのですが、上で紹介したSteven G. Johnsonさんは全然違います。
#Julia言語 JuliaConでの講演なのでJuliaも話題に。
pure Juliaで書き直したerfinv函数はSciPyで使われているFortranライブラリよりも2～3倍速いらしい。
これを「Juliaのコンパイラの最適化がすごいのか？」のように解釈しちゃダメです。Johnsonさんはもっと普遍的な方法の解説をしています。
pure Juliaで書き直したerfinv函数はSciPyで使われているFortranライブラリよりも2～3倍速いらしい。
これを「Juliaのコンパイラの最適化がすごいのか？」のように解釈しちゃダメです。Johnsonさんはもっと普遍的な方法の解説をしています。
#Julia言語 erfinv函数は正規分布を扱うときには必須の基本特殊函数です。そういう基本特殊函数の普及しているライブラリが採用しているアルゴリズムは十分に最適化されていないっぽい。
そういうことがJuliaで特殊函数を実装し直すことによって判明してしまったわけです。
そういうことがJuliaで特殊函数を実装し直すことによって判明してしまったわけです。
#Julia言語 SpecialFunctions.jlのerfinvのコード
こちらでは多項式近似ではなく、有理函数近似が使われています。
Horner法で分子分母の多項式を効率よく計算している。
この方法はJulia以外の言語でもほとんどそのまま使えます。
github.com/JuliaMath/Spec…
こちらでは多項式近似ではなく、有理函数近似が使われています。
Horner法で分子分母の多項式を効率よく計算している。
この方法はJulia以外の言語でもほとんどそのまま使えます。
github.com/JuliaMath/Spec…
#Julia言語 erfinv関連 g++との比較
https://twitter.com/genkuroki/status/1312208871776153605
#Julia言語 SGJさんはMITでの講義で指数積分函数E₁(z)の実装を宿題で出しています。その模範解答(Julia v0.6版)が
nbviewer.jupyter.org/github/steveng…
(↑これも非常に面白い！)
にあり、Julia v1でも動くようにしたものが
nbviewer.jupyter.org/gist/genkuroki…
にある。あるFortranライブラリの5～6倍の速さ！
nbviewer.jupyter.org/github/steveng…
(↑これも非常に面白い！)
にあり、Julia v1でも動くようにしたものが
nbviewer.jupyter.org/gist/genkuroki…
にある。あるFortranライブラリの5～6倍の速さ！
既存の基本特殊函数ライブラリの最適化が不十分であることは、 #Julia言語 について調べた人が気付いてしまう驚くべき事実の一つだと思います。
誰かCでも使えるような(ゆえにほとんどの環境で使えるような)、爆速の基本特殊函数ライブラリを書くべきだと思います。
誰かCでも使えるような(ゆえにほとんどの環境で使えるような)、爆速の基本特殊函数ライブラリを書くべきだと思います。
例えばベイズ統計関連の数値計算(例えばMCMC法)では確率密度函数の対数の計算が膨大な回数繰り返されます。その計算では基本特殊函数が使われます。
仮にその部分が数倍高速化がすればMCMC法の計算も体感できるくらい速くなると思う。
基本特殊函数ライブラリの最適化は重要です。
仮にその部分が数倍高速化がすればMCMC法の計算も体感できるくらい速くなると思う。
基本特殊函数ライブラリの最適化は重要です。
ただし、CやFortranで書かれた基本特殊函数ライブラリが十分に最適化されていなかったことには理由があると思う。
MITでの宿題の答え nbviewer.jupyter.org/github/steveng… を見ると、特殊函数のコードの最適化を気分よく行うためには、コードの自動生成(マクロ)、視覚化、その他諸々の道具が必須なことが分かる。
MITでの宿題の答え nbviewer.jupyter.org/github/steveng… を見ると、特殊函数のコードの最適化を気分よく行うためには、コードの自動生成(マクロ)、視覚化、その他諸々の道具が必須なことが分かる。
#数楽 #Julia言語 #統計 google.com/search?q=minim…
minimax polynomial approximation of sine を検索
「ミニマックス」は「最悪の(最大の)誤差を最小化するようにパラメータを調節すること」を意味しています。
統計的推測のミニマックスゲームの例↓
minimax polynomial approximation of sine を検索
「ミニマックス」は「最悪の(最大の)誤差を最小化するようにパラメータを調節すること」を意味しています。
統計的推測のミニマックスゲームの例↓
https://twitter.com/genkuroki/status/1362386206781153288
最悪の場合の害を最小化すること(ミニマックス)は多くのゲームの基本。
例えば、相手が自分にとって最悪の手を指して来ても対応できるような手を探すことが将棋や囲碁の基本。
そういう考え方は最悪の誤差を最小化したい場面でも出て来る。
例えば、相手が自分にとって最悪の手を指して来ても対応できるような手を探すことが将棋や囲碁の基本。
そういう考え方は最悪の誤差を最小化したい場面でも出て来る。
#数楽 私も知らないので、自分でテキトーにやってみました。
* チェビシェフ多項式のミニマックス性を使って、sin((π/4)x)の十分高次のマクローリン展開を区間[-1, 1]で最良近似する13次の多項式を数値的に求める。
* そこからsin(x)の近似多項式を作る。
gist.github.com/genkuroki/2f82…
* チェビシェフ多項式のミニマックス性を使って、sin((π/4)x)の十分高次のマクローリン展開を区間[-1, 1]で最良近似する13次の多項式を数値的に求める。
* そこからsin(x)の近似多項式を作る。
gist.github.com/genkuroki/2f82…
https://twitter.com/tchaikovsky1026/status/1363390971614928903
#数楽 高次の多項式を低次の多項式で区間[-1,1]上のsupノルムの意味で最良近似するにはチェビシェフ多項式を使えます。
だから、原理的には、sin((π/4)x)を[-1,1]でものすごい精度で近似する多項式を高次のマクローリン展開で作って、低次の多項式で最良近似すれば欲しいものが得られます。
だから、原理的には、sin((π/4)x)を[-1,1]でものすごい精度で近似する多項式を高次のマクローリン展開で作って、低次の多項式で最良近似すれば欲しいものが得られます。
#数楽 チェビシェフ多項式を使った計算はBigFloatで遂行。
添付画像の上半分のグラフはそのようにして私が作った多項式函数の値とJuliaの例のsinの差の10^16倍のグラフです。ほぼ一致！
下半分のグラフは13次のマクローリン展開の場合です。x=π/4の近くで誤差が大きい
gist.github.com/genkuroki/2f82…
添付画像の上半分のグラフはそのようにして私が作った多項式函数の値とJuliaの例のsinの差の10^16倍のグラフです。ほぼ一致！
下半分のグラフは13次のマクローリン展開の場合です。x=π/4の近くで誤差が大きい
gist.github.com/genkuroki/2f82…
#数楽 #Julia言語
添付画像は上から順に
* BigFloat表示での私が求めた13次の近似多項式
* Juliaのtrig.jlにある係数
* Float64に変換した私が求めた13次の近似多項式
* Float64に変換した13次のマクローリン展開
Juliaのと私のでは結構値が違います。
gist.github.com/genkuroki/2f82…
添付画像は上から順に
* BigFloat表示での私が求めた13次の近似多項式
* Juliaのtrig.jlにある係数
* Float64に変換した私が求めた13次の近似多項式
* Float64に変換した13次のマクローリン展開
Juliaのと私のでは結構値が違います。
gist.github.com/genkuroki/2f82…
#数楽 #Julia言語
何も知らないので(チェビシェフ多項式は知っていた)、全部自分でテキトーにやってみた↓
gist.github.com/genkuroki/2f82…
sin(x)の13次の多項式による近似
結果的にJuliaのsinの値とほぼ同じ値を区間[0,π/4]上で計算できるようになりました。
自分でやってみると色々分かる。
何も知らないので(チェビシェフ多項式は知っていた)、全部自分でテキトーにやってみた↓
gist.github.com/genkuroki/2f82…
sin(x)の13次の多項式による近似
結果的にJuliaのsinの値とほぼ同じ値を区間[0,π/4]上で計算できるようになりました。
自分でやってみると色々分かる。
#数楽
#Julia言語 のtrig.jlにおけるsin(x)を近似する13次多項式と私が求めた13次多項式の係数がかなり違う理由：
* Juliaの側は最大相対誤差を最小化
* 私の側はsupノルムで測った最大絶対誤差を最小化
添付画像
① 絶対誤差
② 相対誤差
S(x)は私が求めたもので、sin_Float64はJuliaのやつ。
#Julia言語 のtrig.jlにおけるsin(x)を近似する13次多項式と私が求めた13次多項式の係数がかなり違う理由：
* Juliaの側は最大相対誤差を最小化
* 私の側はsupノルムで測った最大絶対誤差を最小化
添付画像
① 絶対誤差
② 相対誤差
S(x)は私が求めたもので、sin_Float64はJuliaのやつ。
https://twitter.com/genkuroki/status/1363421270981812225
#数楽 あまりにも違う答えが出て来たので悩みました。
最大相対誤差の最小化と最大絶対誤差(supノルムで測った距離)の最小化で、結果が違うのは当たり前でした。
でも、13次多項式によるsinの2つの近似の差はFloat64の精度では大したことないです。
最大相対誤差の最小化と最大絶対誤差(supノルムで測った距離)の最小化で、結果が違うのは当たり前でした。
でも、13次多項式によるsinの2つの近似の差はFloat64の精度では大したことないです。
#数楽 マクローリン展開をそのまま使うことの効率の悪さが印象的。同じ計算量で最悪の誤差が100倍以上悪化する。
https://twitter.com/genkuroki/status/1363420271638552579
• • •
Missing some Tweet in this thread? You can try to
force a refresh
