kashiの日記

2017/03/10(金)kv-0.4.41

kvライブラリを0.4.41にアップデートしました。

今回の変更は、キャストし忘れでode solverの内部使用型をddとしたときにコンパイル出来なくなっていた修正だけです。

また、ドキュメントにkvライブラリの応用例のページを追加しました。まだ2つだけしか載せていませんが、手元には載せたい応用例がたくさんあるので今後増やしていきたいと思います。また、このページでは数式を表示するプラグインのMathJaxを使ってみました。何か不具合があればお知らせください。

コメント（0件）

2017/01/27(金)kv-0.4.40

精度保証

kvライブラリを0.4.40にアップデートしました。

今回は、T大のT先生の依頼によるODE solver仕様変更が含まれています。その他、psaのドキュメントを最新にしたり、kraw_approxの1変数版を追加したりしました。

ODE solverの仕様変更は以下のようなものです。常微分方程式の計算を行なうodelong_maffineなどのodelong系の関数では、計算しながら何かをさせるためのcallback関数を指定する機能があります (正確に言うと、()演算子を持つcallback用のオブジェクト)。これの戻り値は今まではvoidだったのですが、それをboolに変更してその値をチェックするようにし、もしfalseなら計算を中止するようにしました。これにより、終了条件として時刻でなく、「値が何かの条件を満たしたら終了」のようなことが出来るようになりました。これは便利になっているのですが、もし今までのversionのkvを利用していてcallback関数を継承によりカスタマイズする機能を使っていた場合、コンパイルエラーになってしまうという申し訳ないことになってしまいました。修正は簡単で、operator()を定義する関数の戻り値をboolに変えて末尾にreturn true;を補うだけです。

こういうことはなるべく起きないように仕様を考えるときは気をつけていたのですが、今回はメリットが仕様変更のデメリットを上回ると考えてこうすることに決心しました。恐らくこの機能を実際に使っているのはT大のT先生だけだと思いますが、もし他に使っている方がいらっしゃれば大変申し訳ありません。

また、githubの方のソースコードも最新になっています。結局ローカルの開発体制をgitに移行するには至らず、version up時にまとめてgithubを更新するような何とも無意味な状態になっています。新しいものに馴染むのは難しい…

コメント（0件）

2016/12/30(金)kv-0.4.39

精度保証

kvライブラリを0.4.39にアップデートしました。

今回は特にバグフィックスは無く、主に使い勝手に関する機能追加です

大域最適化用のoptimize.hppにminimize/minimize_value/maximize/maximize_valueを追加。1変数用の関数も追加。
dka.hppのdka関数の収束条件を見直し。
高階微分用のhighderiv.hppに、高階微分の関数オブジェクトを生成するHighDerivと、偏高階微分の関数オブジェクトを生成するPartialDerivを追加した。
allsol.hppに1変数関数用のallsol関数を追加した。

大域最適化はあまり説明もなくひっそりと実装してあったのですが、本格的に利用したくなる場面があって、それには少し使いづらいので、いろいろと関数を追加しました。説明は、「大域最適化」に詳しく書きました。

dka関数の修正は、ルジャンドル多項式のゼロ点を求めようとすると収束条件が厳しすぎてループから出てこない問題点が見つかったのがきっかけです。高次のルジャンドル多項式は係数が非常に大きく、そのため残差に入る丸め誤差が大きくなるので収束しづらくなっていたようなので、係数の大きさを見て収束条件を加減するように変更してみました。

HighDerivおよびPartialDerivは、高階(偏)微分の関数オブジェクトを生成するもので、数値積分の誤差評価に出てくる高階微分の区間評価において、optimizeと組み合わせて過大評価を抑えることを意図して作ったものです。

いずれも誰かが使っているのを見て必要性を感じて作ったもので、ユーザのありがたさを感じています。皆様、もっともっと使って下さい！

コメント（0件）

2016/12/28(水)中点はどうやって計算するべきか

精度保証

今回は非常に単純な話で、IEEE754のdoubleで表現された数a,bが与えられたとき、その中点(a+b)/2はどう計算するのが最善でしょうか、という問題です。もちろん誤差は入りますが、出来ればnearestにしたい、すなわち、真の(a+b)/2をIEEE754のルールで丸めたものを計算することを目標にします。IEEE754の内部表現は2進数ですから、2で割る操作は指数部を1減らすだけなので普通は誤差は入りません。しかし、極端な場合を考えると案外難しいのです。

(a+b)*0.5

まずは普通に足して2で割る方法。「2で割る」は「0.5を掛ける」にした方がコンパイラに優しいでしょう。ごく普通の方法で安定度抜群ですが、例えば

a=2¹⁰²³
b=2¹⁰²³

の場合、a+bがオーバーフローしてしまう問題があります。

a0.5+b0.5

オーバーフローを避けて先に2で割る方法です。計算量は少し増えてしまいますが。こちらは、アンダーフロー領域で問題になります。例えば、

a=2^-1074
b=2^-1074

のとき、a*0.5、b*0.5はともにアンダーフローで0になってしまい、結果も0になってしまいます。

a+(b-a)*0.5

これもよく見る式です。IEEE754では関係ありませんが、内部演算が10進数の場合にこうやって計算することが推奨されているのを見かけたことがあります。しかし、

a=2¹⁰²³
b=-2¹⁰²³

でオーバーフローしてしまいます。

a+b0.5-a0.5

上の式でオーバーフローを避けたものです。今度こそ良さそうですが、

a=5*2^-1074
b=7*2^-1074

としてみましょう。このとき、IEEE754の偶数丸めルールの関係で、

a*0.5=2*2^-1074
b*0.5=4*2^-1074

となります。よって計算結果は7*2^-1074となってしまい、真の6*2^-1074になりません。

a+b0.5-a0.5 (ただし上向き丸めで計算)

INTLAB version 9で採用している計算方法です。a*0.5とb*0.5に入る丸めの方向が逆向きになることが無いので、先の例だとうまく計算できます。しかし残念ながら必ずnearestにはならないようで、Mさんがちゃんと計算できない例を発見しました。

a=(2⁵³-1)*2^-1073
b=7*2^-1074

とすると、

(nearest)=2^-1021+2^-1073
(計算値)=2^-1021+2^-1072

と値がずれてしまいます。a*0.5は無誤差ですが、b*0.5と加算で両方上向き丸めになった結果、真値(2^-1021+2^-1073+2^-1075)を丸めたものとずれてしまいました。

結論?

kvライブラリの区間の中点を計算するmid関数では、aとbの絶対値が両方共1より大きいときはa*0.5+b*0.5で、それ以外のときは(a+b)*0.5で計算しています。このようにifで場合分けをしていいなら話は簡単なのですが、matlabの場合はベクトル単位で一括処理しないとパフォーマンスが出ないのでなるべくifを使いたくないという事情があります。if文禁止で完璧な中点の計算アルゴリズムはまだ見つかっていません。誰か考えて！

コメント（0件）

2016/12/10(土)kv-0.4.38

精度保証

kvライブラリを0.4.38に更新しました。

今回は、区間演算にmidrad、max、minを追加しました。また、dka.hppの中にあるvdka (精度保証付きDKA法) で、解が重複する場合の処理をきちんと行なうようにしました。

midradとmax/minの追加は、どちらも同僚のS先生の要請によるものです。少し詳しく説明します。

midrad

midradは、区間 I に対して、

midrad(I, m, r);

とすると中心の近似mと、半径rを同時に計算します。これは、中心m、半径rの区間で I を包含する、すなわち、上端下端型の区間から中心半径型の区間への変換に使われることを意図しています。単純に

m = mid(I);
r = rad(I);

とすれば良さそうなものですが、中心mに丸め誤差が入るので、これだと I を包含しない可能性があるのです。例えば、ε=2^-52として、I=[1,1+3ε]とします。このとき、midとradを別々に計算すると、

m = 1+2ε
r = 1.5ε

となりますが、[m-r,m+r]は I を含まないことが分かります。midradを使うと、ちゃんと

m = 1+2ε
r = 2ε

になります。ミスが起こりやすいポイントなので気をつけて下さい。

max, min

max, minは区間同士のmax, minです。max(I1, I2)は、x1∈I1、x2∈I2としたときのmax(x1,x2)の取り得る範囲の区間を返します。minも同様です。計算は、

max([a,b],[c,d]) = [max(a,c),max(b,d)]
min([a,b],[c,d]) = min(a,c),max(b.d)]

のように行います。

この関数は、後述するようにVisual C++と少し相性が悪いこともあって実装していなかったのですが、実装しないと思わぬバグに遭遇する可能性があると同僚のS先生に指摘され、実装することにしました。例えば、max, minを実装していないversionのkvで、

#include <kv/interval.hpp>
#include <kv/rdouble.hpp>

typedef kv::interval<double> itv;

using namespace std;

int main()
{
    cout << max(itv(2., 4.), itv(3., 5.)) << endl;
    cout << min(itv(3., 5.), itv(2., 4.)) << endl;
}

のようなプログラムを動かすと、予想に反して

[2,4]
[3,5]

という結果が帰って来ます。これは、std::max, std::minが呼び出されてしまっており、恐らくmaxは

template <class T> 
T max(const T& a, const T& b)
{
    return (a < b) ? b : a;
}

のような実装になっていて、[2,4]<[3,5]は偽なので[2,4]が返されてしまっていると推測されます。よく考えれば分かるとは言え、これは事故を誘発しやすいので、ちゃんと区間用のmaxとminを実装しました。kv-0.4.38で上のプログラムを動かすと、ちゃんと

[3,5]
[2,4]

となります。

なお、Visual C++はwindows.hというヘッダファイルを持っており、これをincludeするとマクロでmaxとminを定義してしまうという大変極悪な仕様になっています。これがincludeされていると、max,minを定義している箇所がマクロ置換されてしまいinterval.hppがコンパイル出来ません。windows.hとkvをどうしても同時に使いたければ、

#define NOMINMAX
#include <windows.h>

のようにwindows.hをincludeする前にNOMINMAXを定義して下さい。

コメント（0件）