Julia 言語で汎用一次元数値積分関数を書きました
二重指数関数型数値積分公式 を使った、非適応型の数値積分プログラムです。
高橋・森によって提案された変数変換を施すことによって、被積分関数を積分区間の両端で急激に減衰する新しい関数に変換します。このような関数は基本的に元の関数よりも数値積分が容易で、主に台形則などを使って高速かつ精度よく近似解を得ることができます。
被積分関数が積分区間全体にわたってなめらかに変化する場合に、特に少ない評価点数で高精度な値を返します。また、もう一つの特徴として端点での特異性に非常に強い点が挙げられます。
Tanh-sinh quadrature
積分区間や被積分関数の振る舞いによって多くの変形がありますが、最も汎用的かつ基本的な三種を実装しました。tanh-sinh quadrature はその一つで、任意の被積分関数の [-1, 1] の区間での数値積分を求めます。
\begin{align} \int_{-1}^{1} f(x) dx = \int_{-\infty}^{\infty} f\left(x(t)\right) \frac{dx}{dt} dt = \int_{-\infty}^{\infty} f\left(x(t)\right)w(t) dt \end{align}
についての被積分関数
を、
についての新しい被積分関数
へ変換し、それに伴い積分区間をうつしています。ここで、
および
は以下のような関数です。
\begin{eqnarray} x\left(t\right) = \sinh \left( \frac{π}{2} \sinh t \right) \\ w\left(t\right) = \frac{π}{2} \cosh t \cosh\left(\frac{π}{2} \sinh t \right) \end{eqnarray}
この変数変換が肝なので、以下の積分を例にとり被積分関数が実際にどのように変化するかを見てみましょう。
\begin{align} \int_{-1}^{1} \frac{dx}{\left(2-x\right)\left(1-x\right)^{1/4}\left(1+x\right)^{3/4}} \end{align}
この関数は積分区間の両端 で急速に発散しており、数値積分を難しくしています。しかし、変数変換を施したあとの関数
は正負の両方向で急速に減衰する関数に変換されていることがわかります。
元の関数の区間両端での発散が解消された代わりに両無限区間 [-∞, ∞] での積分になっています。しかし、変換後の被積分関数は正負の両方向へ向かうにつれ急激に減衰する関数となっているため、台形則の適用により無限和に近似し が十分小さくなった地点で打ち切ることができます。
\begin{align} \int_{-\infty}^{\infty} f\left(x(t)\right)w(t) dt \approx h\sum_{k = -\infty}^{\infty} f(x_k)w_k \end{align}
ここで および
は台形則における刻み幅
を用いて次のように表せます。
\begin{eqnarray} x_k = x(kh) \\ w_k = w(kh) \end{eqnarray}
および
は被積分関数に依存しないため、刻み幅
を決めておけば事前に
,
のテーブルを計算しておくことができます。これも速度が求められる場面では重要な特性といえるでしょう。また事前に計算しない場合でも、tanh-sinh quadrature については
が奇関数、
が偶関数なので、
に対する対称性 (
) を利用して計算量を減らすことが可能です。
さらに別の変数変換を適用することで任意の有限区間 [a, b] での積分をすべてカバーします。
\begin{align} \int_{a}^{b} f(x) dx = \frac{b - a}{2} \int_{-1}^{1} f\left(x(u)\right) du = \frac{b - a}{2} \int_{-\infty}^{\infty} f\left(x\left(u(t)\right)\right)w(t) dt \end{align}
\begin{align} where\ \ x(u) = \frac{b + a}{2} + \frac{b - a}{2}u \end{align}
DoubleExponentialFormulas.jl はほかに exp-sinh quadrature, sinh-sinh quadrature を組み合わせて任意の有限区間 [a, b]、片無限区間 [a, ∞] [-∞, b]、両無限区間 [-∞, ∞] に対応しています。
使い方
関数 quadde は Float64 の精度で数値積分を計算します。
I, E = quadde(f::Function, a::Real, b::Real, c::Real...;
atol::Real=zero(Float64),
rtol::Real=atol>0 ? zero(Float64) : sqrt(eps(Float64)))
I は得られた積分値、E は推定誤差です。たとえば、f(x) = 1/(1 + x²) の区間 [-1, 1] での積分は次のように計算します。
using DoubleExponentialFormulas using LinearAlgebra: norm f(x) = 1/(1 + x^2) I, E = quadde(f, -1.0, 1.0) I ≈ π/2 # true E ≤ sqrt(eps(I))*norm(I) # true
atol と rtol はそれぞれ目標とする絶対誤差および相対誤差で、E がこれらのうち大きいほうを下回るまで精度を上げながら反復を繰り返します。ただし、E はあくまで推定値であり I と真の値の正確な誤差ではありません。しかし、少なくとも I が収束した場合には E <= max(atol, rtol*norm(I)) が真となることが期待できます。言い換えれば、その条件が満たされていない場合、I は既定の反復回数のうちに収束しなかったことを意味しており、信頼できない値であると考えられます。
基本的には積分値の大きさに応じた精度を要求できる rtol を使うのが便利です。しかし、予想される積分値 I の絶対値が非常に小さい場合には rtol*norm(I) がアンダーフローを起こす恐れがあります。こうなると収束条件を満たせず反復回数が大きくなってしまうので、これを避けるために同時に atol にも小さい値を設定しておくと計算量を節約できる場合があります。atol が与えられた場合は、rtol は明示的に与えられない限り使われないので両方を明示的に指定しましょう。
区間には Inf を指定することができます。
# Computes ∫ 1/(1+x^2) dx in [0, ∞) I, E = quadde(x -> 1/(1 + x^2), 0.0, Inf) I ≈ π/2 # true # Computes ∫ 1/(1+x^2) dx in (-∞, ∞) I, E = quadde(x -> 1/(1 + x^2), -Inf, Inf) I ≈ π # true
また、a, b 以降にも区間が与えられた場合には積分区間を分割して計算した和を返します。
\begin{align} I = \int_{a}^{b} f(x) dx + \int_{b}^{c} f(x) dx + \cdots \end{align}
これは、積分区間に不連続点・発散する点を含む場合に有用です。例えば、関数 は
で発散しますが、区間を分割し
を端点にすることでうまく計算できるようになります。
# Computes ∫ 1/sqrt(|x|) dx in (-1, 1) # The integrand has a singular point at x = 0 I, E = quadde(x -> 1/sqrt(abs(x)), -1.0, 0.0, 1.0) I ≈ 4 # true
得意・不得意な関数の傾向
基本的には被積分関数 が積分区間全域にわたって緩やかに変化するような関数は効率よく計算するようです。逆に積分区間の一部だけ激しく変化するような関数の場合は、最も変化の激しい部分に引っ張られて全体の評価点数を増やすので効率が良くないです。これは非適応型の宿命といえます。具体的な例を挙げると、積分区間にスパイクが発生するような関数は評価回数が多くなってしまいます。
また、上でも述べられたように端点の特異性に関してはほかのアルゴリズムと比べても強いため、ほかのアルゴリズムでは対応できない場合に試してみるのはよいかもしれません。
QuadGK との比較
感想など
書く前は難しそうだと思っていたので、思ったよりもちゃんと動くな、というのが正直な感想です。いままで数値積分関数はブラックボックスとして扱っていたので勉強になりました。数値積分関数の説明には往々にして注意書きがいろいろあったり、「必ずしも正しい値を返すとは限りません」みたいな歯切れの悪い文言が並んでいたりするものなんですけれども、なんとなく理由がわかったような気がします。上のスパイクが発生するような関数なんかは二重指数関数型数値積分法に限らず数値積分関数にとっては悪意の塊みたいな関数なんだな、と。
汎用性のために一部速度や精度を犠牲にしている部分もあります。また、二重指数関数型積分法には被積分関数によっていろいろな変種が考案されているので、汎用の関数を書いておいてなんですが問題にあわせて最適化したものをつかうのが一番ですね。それでもほとんどの被積分関数に関して7～8桁（デフォルト）の精度での計算には問題なく動くようです。13桁かそれ以上を目指す場合には被積分関数によっては厳しいかもしれないです。
今のところはかなり素朴な実装 (v0.0.5) になっていて特にひねりもないので、何か思いついたら改良していこうと思います。
Julia の Language server を Vim でうごかす その２
以前、こういう記事を書きました。
しかし、準備が煩雑、必要なパッケージをグローバルにインストールしなければならないなど不満が多くありました。なので、これらの不満を解決するべく必要な設定やソースコードを一つの Vim plugin にまとめました。
必要な準備は julia コマンドにパスを通すこと、vim-lsp をインストールすることぐらいです。詳しくは vim-lsp-julia のREADME かドキュメントを読んでください。
注意としてはおそらく最初の起動時にはサーバーの準備ができるまでかなり時間がかかります。私のノートPCだと 20~30 分かかったような…？ほとんどは SymbolServer.jl がキャッシュを作るのにかかっていると思われます。なのでインストールしているパッケージの数にもよると思います。将来的にはよくなるかもしれません。
g:lsp_julia_depot_path 変数にパスを設定することでコンパイルキャッシュやログを外に切り出せるようにしたんですけど、肝心の SymbolServer.jl のキャッシュ (作成に時間がかかる) をまだ外に出せないので中途半端な感じです。これも将来的には解決できそうです。
また Vim のカラースキームを書いた
雪山かなんかの写真をみてなんかこんな感じ、と作りはじめました。
GitHub - machakann/vim-colorscheme-snowtrek: A light colorscheme for vim text editor
かなり明るめです。いままでは、カラースキームを作る場合 'background' オプションが light の場合と dark の場合の両面をつくっていたんですけど、今回はあっさり放棄しました。明るい色使いだけです。
わりと気に入ってて、特に vim の :help を読むのが気分いいです。個人的に。
Julia 言語の丸め関数
Julia 言語の丸め関数についてまとめました。
丸め関数
floor, ceil, trunc, round があります。
trunc(x)
原点 0 へ向かう方向に最も近い整数へと丸めます。
trunc(-1.1) == -1.0 trunc(-0.9) == -0.0 trunc( 0.9) == 0.0 trunc( 1.1) == 1.0
ceil(x)
小数点以下切り上げです。正方向に最も近い整数へと丸めます。
ceil(-1.1) == -1.0 ceil(-0.9) == -0.0 ceil( 0.9) == 1.0 ceil( 1.1) == 2.0
floor(x)
小数点以下切り捨てです。負方向に最も近い整数へと丸めます。
floor(-1.1) == -2.0 floor(-0.9) == -1.0 floor( 0.9) == 0.0 floor( 1.1) == 1.0
round(x)
いわゆる四捨五入ですが、デフォルトだと言葉通りの四捨五入にはなりません。次の "RoundingMode について" を参照してください。
floor(-1.1) == -1.0 floor(-0.9) == -1.0 floor( 0.9) == 1.0 floor( 1.1) == 1.0
RoundingMode について
round 関数は丸めたい浮動小数点数につづいて RoundingMode を指定することで挙動を変えられます。差異は小さいながら、時に重要なので違いを理解しておくことが重要です。
RoundNearest
round 関数のデフォルトの丸め規則です。半整数 ( ... -2.5, -1.5, -0.5, 0.5, 1.5, 2.5, ...) は 最も近い偶数 へ丸められるという、ちょっとトリッキーな動きです。
※2019/8/4 追記: Banker's rounding というらしいです。統計量の偏りを小さくすることを目的とした工夫みたいです。こちら の記事がわかりやすかった。
round(-1.5) == round(-1.5, RoundNearest) == -2.0 round(-0.5) == round(-0.5, RoundNearest) == -0.0 round( 0.5) == round( 0.5, RoundNearest) == 0.0 round( 1.5) == round( 1.5, RoundNearest) == 2.0
RoundNearestTiesAway
半整数は原点 0 から離れる方向に丸められます。絶対値を丸めているとも言えます。C/C++ 言語の round 関数と同じ挙動です。
round(-1.5, RoundNearestTiesAway) == -2.0 round(-0.5, RoundNearestTiesAway) == -1.0 round( 0.5, RoundNearestTiesAway) == 1.0 round( 1.5, RoundNearestTiesAway) == 2.0
RoundNearestTiesUp
半整数は正の方向へ切り上げられます。Java/JavaScript 言語の round 関数と同じ挙動です。
round(-1.5, RoundNearestTiesUp) == -1.0 round(-0.5, RoundNearestTiesUp) == -0.0 round( 0.5, RoundNearestTiesUp) == 1.0 round( 1.5, RoundNearestTiesUp) == 2.0
RoundToZero
trunc 関数と同じく原点 0 へ向かって丸めます。trunc(x) == round(x, RoundToZero) となります。
RoundUp
ceil 関数と同じく小数点以下切り上げです。ceil(x) == round(x, RoundUp) となります。
RoundDown
floor 関数と同じく小数点以下切り捨てです。floor(x) == round(x, RoundDown) となります。
RoundFromZero
ドキュメントによると BigFloat 型の数値を丸める場合にのみ使われるようです。よくわかりません。
個人的には、ほぼ整数だとわかっている数の Int 型への変換に RoundNearest を、ランダムな浮動小数点数を丸めたい場合に RoundNearestTiesUp を使っていることが多いような気がします。
# a, b は半整数 c = round(Int, a + b) # d ∈ [-3.0, -2.0, -1.0, 0.0, 1.0, 2.0, 3.0] d = round(cos(rand(Float64)*π)*3, RoundNearestTiesUp)
参考
肥大した $HOME/.julia を整理する
- Julia の REPL を立ち上げる
]キーを押して Pkg REPL モードに入るgcと入力してEnter
つまり、REPL を立ち上げて ]gc<Enter> 。簡単ですね！
必要なくなったバージョンのパッケージリポジトリのコピー (たぶん $HOME/.julia/packages/{package name}/ 以下にあるやつ) を削除するらしいです。Makie.jl, AbstractPlotting.jl, GLMakie.jl の master を追いかけていたら .julia フォルダのサイズが大変なことになってました。たまに自動で走らせたいけど頻度に悩む。
vim-sandwich で関数を消す
Vim の話です。
vim-sandwich というプラグインを以前書きまして、これを更新したよ、というのがこの記事の内容です。
これは何？
ざっくりといえば、vim-surround クローンであり、おおよそ対応する機能を持っています。つまり、テキストオブジェクトを括弧やクオーテーションなどで囲んだり、これらを削除・置換したりするためのプラグインです。
- 囲む
foo → (foo)
- 削除
(foo) → foo
- 置換
(foo) → [foo]
vim-surround よりも拡張性に主眼を置いており、複雑な動作をユーザーが定義できるようになっています。
何が変わったのか？
vim-sandwich の機能の一つに関数囲みの削除があります。(参考) 関数名まで含めて括弧を消す、というわけです。
func(foo) → foo
しかし、言語によっては関数名の記述はちょっと違ったりしますね。例を上げれば、Vim script を書くときには頻繁に s:func() という感じの関数名を書きます。このようなかんじで関数名のパターンもカスタムしたい、という要望がきたのでできるようにしました。
Allow user to customize function name to delete · Issue #71 · machakann/vim-sandwich · GitHub
以下のような変数を定義することでグローバルなルールを定義できます。foo::bar.baz() のようなものを関数として認識したい場合は次のような設定を vimrc などに書きます。
let g:sandwich#magicchar#f#patterns = [ \ { \ 'header' : '\<\h\k*::\h\k*\.\h\k*', \ 'bra' : '(', \ 'ket' : ')', \ 'footer' : '', \ }, \ ]
あるいは、b:sandwich_magicchar_f_patterns を使えばバッファローカルなルールを定義できます。これは主にファイルタイプに特有のパターンを記述するのに使われるでしょう。先の Vim script の例で言えば次のようになります。
augroup sandwich-filetype-vim autocmd Filetype vim let b:sandwich_magicchar_f_patterns = [ \ { \ 'header' : '\C\<\%(\h\|[sa]:\h\|g:[A-Z]\)\k*', \ 'bra' : '(', \ 'ket' : ')', \ 'footer' : '', \ }, \ ] augroup END
(ただし、上の設定は vim-sandwich が持つのでユーザーが自ら設定する必要はありません。)
リストの要素の辞書一つが、一種の関数のパターンに対応しており、四つのキー ('header', 'bra', 'ket', 'footer') を持っていなければいけません。'header' は関数名に一致する正規表現、'bra' は開き括弧、'ket' は閉じ括弧に一致する正規表現、'footer' は閉じ括弧に続くテキストに一致する正規表現を値に持ちます。
リストなので複数登録しておくことも可能です。これで関数囲みの削除が捗りますね！
vim-swap で関数の引数を並べ替える
Vim の話です。
vim-swap というプラグインを以前書きまして、これを更新したよ、というのがこの記事の内容です。
これは何？
もっとも代表的な用途は、関数の引数のようなコンマで区切られたテキストの順番を入れ替える、というものです。例えば、次のような三つの引数をとる関数があり Banana の上にカーソルがあるとします。
Basket(Apple, Banana, Cherry)
この時 g< と入力すると Banana と Apple を入れ替えます。
Basket(Banana, Apple, Cherry)
反対に、g> と入力すると Banana と Cherry を入れ替えます。
Basket(Apple, Cherry, Banana)
もう一つ機能があって、gs と入力すると対話的に入れ替えを行う swap mode を開始します。これは何が嬉しいのかと申しますと、ちょっと複雑な入れ替えを . で繰り返せるという点です。上の g< や g> も、もちろん . で繰り返し可能ですが、繰り返せるのは一つの要素の入れ替えのみです。gs は複数の要素の移動をまとめて繰り返せます。
例えば、次のバッファには Basket(Cherry, Banana, Apple) が三回現れますが、これをすべて Basket(Apple, Banana, Cherry) にするために g> や g< を使うと、. を繰り返し押す必要があります。
function! Basket(Cherry, Banana, Apple) abort return 0 endfunction call Basket(Cherry, Banana, Apple) call Basket(Cherry, Banana, Apple)
しかし、gs を使って入れ替えておくと、後は . 一回で変換できます。gif を見て何となく感じてください。
何が変わったのか？
Swap mode のコマンドを増やしました。gs と入力して swap mode に入り、以下のコマンドを使います。満足したら <Esc> キーで swap mode を抜けます。
Sort
s / S キーで要素をソートします。s は昇順、S は降順にソートします。
Reverse
r キーで要素を反転します。
Grouping
g キーでカーソル下の要素と右隣りの要素をグループ化して一つの要素として扱います。G でグループ化を解除します。
実装してみると、逆になぜ今までなかったのかと思うような機能ばかりでした。特にグルーピングがずっと欲しかった。部分的なソートや部分的な反転を実現したかったんですけど、力尽きたので切実に欲しくなってからにします。
vim-swap は最近 Github で 100 star を超えました🎉 ありがとうございます。
