2016/05/09(月)kv-0.4.36

というわけで、予告通りkvライブラリを0.4.36にアップデートしました。

内容は前回の記事の予告通りで、maffine3を正式アルゴリズムに昇格させました。
  • maffine → 一応maffine0として残したが消滅候補
  • maffine2 → そのまま
  • maffine3 → maffineに改名
使い勝手は0.4.34以前と変わらず、単に性能がよくなっているはずです。

2016/05/08(日)新しいODE solverの性能評価

さて、前回の記事で書いた新しいアルゴリズムの性能を簡単な例で示します。例として選んだのはvan der Pol方程式です。
x''- μ(1-x2)x'+x = 0
この方程式はμ≥0の値によって計算のしやすさが変化し、μが1くらいまでなら計算しやすいnon-stiffな方程式ですが、μが大きくなると極めて計算しづらいstiffな方程式になります。ここでは、μ=1とμ=100でこの方程式の軌道を計算し、ステップ幅の変化を調べてみました。ステップ幅はある方法で局所誤差がmachine epsilon程度になるように計算していますが、stiffなODEの場合そこを変えて大きなステップ幅にしても精度保証の根幹となる不動点定理が成立しなくなってしまい小さなステップ幅への修正を強いられてしまいます。初期値は
x(0) = x'(0) = 1
とし、t=200まで計算しました。この軌道を相図(phase diagram)で示すと、次のような感じです。まずμ=1の場合:
dense1.png

周期解(limit cycle)に巻き付いている様子がよく分かります。次にμ=100の場合:
dense100.png

こちらは急激な変化を含んだ周期の長い軌道になります。

さて、これを、
  • ode-maffine (従来の標準的なアルゴリズム)
  • ode-maffine2 (従来の高速アルゴリズム。初期値に関する微分が出来ない欠点がある)
  • ode-maffine3 (新しいアルゴリズム)
  • awa (Lohnerによる有名な精度保証プログラム)
で計算して比較してみます。まずμ=1で:
step1.png

横軸は時刻tで、縦軸は上のグラフがx、下のグラフがODE Solverが選んだステップ幅です。non-stiffな方程式なので、どれでもあまり大きな違いがないことが分かります。計算時間は、
アルゴリズム計算時間(sec)
maffine2.276
maffine21.240
maffine31.798
awa3.815
でした。次に、μ=100の場合を示します。
step100.png

どの方法もμ=1の場合に比べると刻み幅が小さくなっています。変化の激しさに応じて刻み幅を頑張って適応させていることがよく分かります。刻み幅はawaが一番小さく、次にmaffine、maffine2とmaffine3はほぼ同じで比較的大きい刻み幅で計算できていることが分かります。結果として、計算時間にも次のように大きな差が出ました。
アルゴリズム計算時間(sec)
maffine53.820
maffine29.287
maffine313.478
awa176.78
maffine3はmaffineと完全に同じことが出来て性能が良さそうなので、maffineを残す意味はなさそうです。次期versionでは、
  • maffine → 消滅 (別名で残すかも)
  • maffine2 → 従来通り
  • maffine3 → maffineに改名
としたいと思います。これなら0.4.34以前と使い勝手は変わらずに単に性能が良くなったということになります。関数名の変更はユーザに混乱をもたらすので、現version (0.4.35) はなるべく短命にするつもりです。

2016/05/06(金)kv-0.4.35

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

まず、前回(0.4.34)に混ぜてしまったバグの修正です。1.00e-05みたいな数値が1.00e0-5のように表示されてしまうというバグを直しました。

また、常微分方程式のsolverに新しいアルゴリズム(maffine3)を追加しました。kv-0.4.30とstiffなODEで日記に書いたように、maffineはstiffなODEに弱く、maffine2はstiffに強いが初期値に関する精度保証された微分が得られないという状況でしたが、今回、「stiffに強くて初期値に関する微分もできる」maffine3を追加しました。今年のGWはほぼこの実装に費やしてしまいました。

実装したアルゴリズムのアイデアは極めて単純です。常微分方程式の初期値問題を精度よく解くためには、いわゆるWrapping Effectを防ぐために解の初期値に関する微分を利用することがほぼ必須です。それを計算するのに必要なのが「初期値に関する変分方程式」というものです。すなわち、
x'(t)= f(x(t),t)
x(t0) = x0
に対して、x*(t)をこの式の真の解としたy(t)に関する行列微分方程式
y'(t) = fx(x*(t), t)y(t)
y(t0) = I (単位行列)
を初期値に関する変分方程式といいます。ある時刻t1におけるyの値y(t1)は、x(t0)にx(t1)を対応させる写像のヤコビ行列と一致します。

このyを計算するのに、従来は上の2つの式を連立させて、
x'(t)= f(x(t),t)
y'(t) = fx(x(t), t)y(t)
x(t0) = x0
y(t0) = I (単位行列)
のようなxとyの組を未知関数とするn+n×n次元の新しい初期値問題を作成し、これを解くという方法を使っていました。その方が実装が楽だったという事情もあります。

それを変更して、まずxだけの式を解いて真の解x*を区間関数として求め、それを代入してyの式を解く、素直な2段階法が、新しいアルゴリズムです。stiffでない普通の方程式だとどちらを使っても大差ありませんが、stiffな場合だと後者の方が刻み幅が大きくなり、結果として高速になるようです。

アルゴリズムの詳細や性能評価はまた改めて書く予定です。問題が無ければ、将来のバージョンではこの新しいmaffine3をmaffineとし、旧maffineは抹消してしまうことも考えています。

2016/05/01(日)emscriptenで遊んでみた

emscriptenという、c/c++のコードをLLVMのバイトコードを経由してjavascriptに変換するコンバータがよく話題になっています。以前、kvライブラリをemscriptenで変換して動かしたというメールを頂いたことがあって、ずっと気になっていました。そこで、ちょっとインストールして遊んでみました。

インストールしたばかりのubuntu 16.04で試しました。しかし、
sudo apt install emscripten
で入れてもうまく動かない模様です。clang/llvmのバージョンの問題っぽいがよく分かりません。

そこで、本家の Download and install — Emscripten 1.36.1 documentation に従ってインストールしてみました。そこではポータブル版というroot権限なしで入れる方法が紹介されていました。その方が周辺ツールのバージョンに左右されなくてインストールが楽なのかな。
sudo apt install cmake
Portable Emscripten SDK for Linux and OS X (emsdk-portable.tar.gz) をダウンロードして展開し、その中にcdして
./emsdk update
./emsdk install latest
./emsdk activate latest
します。最後の作業で/home/kashi/.emscriptenが作られました。これでインストールは終了。使う前に
source emsdk_env.sh
でPATHの設定が必要です。

適当なc++のプログラムを作って、
emcc test.cc
node a.out.js
で動くのを確認しました。

次に、kvのプログラムを動かしてみました。問題はboostで、emscriptenでは外部ライブラリを使ったプログラムはそのままでは動きません。外部ライブラリもソースから変換する必要があります。kvの場合は幸いなことにboostはheader onlyの範囲でしか使っていませんので、大丈夫なはずです。自分は、
ln -s /usr/include/boost ホームのどこか
とリンクを張ってしまって、-Iにそれを指定してごまかしました。これが正しいかどうかは分かりません。また、javascriptの中ではfesetroundによる丸めの向きの変更は効かないので、そのままでは精度保証出来ません。kvは-DKV_NOHWROUNDを付けると丸めの向きの変更を行わずエミュレーションで方向付き丸めを行なう機能を持っており、それをつけるとちゃんと精度保証付きで計算してくれるようです。
emcc -I(kvの場所) -I(boostの場所) -DKV_NOHWROUND hoge.cc
みたいな感じです。kvのサンプルはmpfrを使ったもの以外は大体動きそうです。実行速度はちゃんと測ってませんがネイティブのC++の30倍程度でしょうか。htmlに変換してfirefoxに食わせた方がnode.jsより速いようにも見えました。

単なる遊びであってあまり意味はないですが、楽しいです。

2016/04/30(土)ubuntu 16.04インストール(11) リモートデスクトップ(xrdp)

ubuntu 16.04にリモートデスクトップを入れたときのメモです。大量にいろいろ突っ込むのであまりお勧めしませんが一応記録を残しておきます。

リモートのubuntuに何かさせたいとき、普通は単にsshで、GUIなアプリを使いたければssh -Xで済みますが、稀にデスクトップ全体を転送したいことがあります。それを実現するものとしては、windowsのリモートデスクトップとVNCが有名ですが、素のwindowsで使えるし速度も速いのでリモートデスクトップを愛用しています。

ところが、unityを採用している最近のubuntuは、リモートデスクトップやVNCのサーバになれません。代わりに標準で「デスクトップの共有」という機能があるのですが(プロトコルはVNC)、これは本体にログインしている状態でしか使えず、ログアウトしてしまうとリモートからログイン出来ないというとても不便なものです。

以下は、ubuntu 16.04に無理やりリモートデスクトップをやらせたときの記録です。
sudo apt install xrdp
これで、サーバそのものは簡単に入ります。しかし、全く動作せず、灰色の画面が出るだけです。動作させるには、リモートから接続した場合はunityを捨てて他のデスクトップを使うしかありません。gnome2の操作性で軽く、見た目を重視しているということで最近人気のある、「MATE」を入れてみました。「メイト」ではなく「マテ」と読むらしいです。
sudo apt install mate-desktop-environment
sudo apt install mate-desktop-environment-extras
リモートでない本体の方でも使えるようになります。いったんログアウトするとデスクトップを選択するためのメニューが出るようになります。

さて、これだけではダメで、/etc/xrdp/startwm.sh を書き換えて、最終行の
. /etc/X11/Xsession
を消して
mate-session
にします。これで一応動いたようですが、日本語キーボードにならず、またかな漢字変換も出来ません。

まず、キーボードを何とかします。これは検索するとたくさん見つかって、/etc/xrdp で
sudo wget http://www.mail-archive.com/xrdp-devel@lists.sourceforge.net/msg00263/km-e0010411.ini
sudo mv km-e0010411.ini km-0411.ini
sudo ln -s km-0411.ini km-e0010411.ini
sudo ln -s km-0411.ini km-e0200411.ini
sudo ln -s km-0411.ini km-e0210411.ini

sudo systemctl restart xrdp
としたら直りました。

かな漢字変換は、fcitxがダメでibusにすれば大丈夫という情報がいくつか見つかりました。
sudo apt install ibus-mozc (すでに入っていた)
/etc/xrdp/startwm.shのmate-session起動前に
export GTK_IM_MODULE=ibus
export QT_IM_MODULE=ibus
export XMODIFIERS="@im=ibus"
ibus-daemon -d
を挿入したらうまく行きました。

ここまであれこれ入れると環境依存で不具合が出そうな気もしますが、参考になれば幸いです。

追記(2017年5月29日)

この方法を使ってNvidiaのGPUを積んだマシンをリモートで管理しようと思ったところ、なぜかNvidiaドライバを有効にするとxrdp経由で入ったときうまく動かず、firefoxが起動しなかったりかな漢字変換が起動しなかったりしていました。xorgドライバなら大丈夫なのですが、それだとcudaが使えず困ります。全く原因が分からずにしばらく放置していたのですが、偶然

nvidiaドライバを入れたubuntuにVNC接続する際に一部のプログラムが起動しない(segmentation fault)場合の解決策

を見かけて解決しました。
cd /usr/lib/x86_64-linux-gnu/
sudo ln -s /usr/lib/nvidia-375/libGLX_indirect.so.0 .
でうまくいきました。nvidia-375の数字の部分は今後のアップデートで変わっていく可能性があるので、そのときはまたやり直す必要があると思います。
OK キャンセル 確認 その他