2018/09/10(月)「精度保証付き数値計算の基礎」チュートリアル

今年の7月に

精度保証付き数値計算の基礎

という本を出しました。私が執筆したのは、第4章「数学関数の精度保証」と、第7章「常微分方程式の精度保証付き数値解法」です。

で、今回、

「精度保証付き数値計算の基礎」チュートリアル

で話させていただきました。そこで使ったスライドを上げておきます。内容は第4章の数学関数の精度保証の話で、発表時間は20分だったので、簡単な内容です。割と多くの方に来ていただいてとても嬉しかったです。

なお、今回のチュートリアルは6章までで、7章のODEと8章のPDEについては年末に別の機会を準備中です。

2018/06/17(日)kv-0.4.44

およそ半年ぶりに、kvライブラリを0.4.44にアップデートしました。

変更点は、
  • dd.hpp, rdd-hwround.hpp, rdd-nohwround.hppを修正し、主にddのオーバーフローに 関するバグを修正。
  • affine.hppのepsilon_reduceで、与えられたaffineベクトルが1次元だった場合に正しく 動かないバグを修正。
  • rkf78.hpp追加
  • その他細かい修正
です。

今回は、dd (double double) を使った区間演算に関する重大なバグフィックスを含んでいます。加減乗除の結果がオーバーフローするかしないかギリギリの値の場合の処理に問題があり、確率は低いものの、精度保証が正しく行われない可能性があるバグです。

以下、その詳細について説明します。ddでの区間演算は、例えば加算では真値と同じか真値よりも大きい値を計算するadd_up、真値と同じか真値よりも小さい値を計算するadd_downによって実現されています。例えばadd_downは、次のようなアルゴリズムによって実現されていました。簡単のため、python-likeに書いています。dd数(x1,x2)とdd数(y1,y2)の和を下向き丸めで計算するものです。
def dd_add_down(x1, x2, y1, y2) :
    z1, z2 = twosum(x1, y1)
    if z1 == -float("inf") :
        return z1, 0
    if z1 == float("inf") :
        return sys.float_info.max, sys.float_info.max * 2.**(-54)
    down()
    z2 = z2 + x2 + y2
    near()
    z1, z2 = twosum(z1, z2)
    return z1, z2
無誤差のtwosumはそのまま計算し、誤差の入る可能性のある「z2 + x2 + y2」の部分は下向き丸めで計算することで、下向き丸めの結果を得ています。これに加えて、赤字で示したオーバーフローの処理を行っています。ddライブラリではIEEE754と同様の無限大の取り扱いを出来るように設計しています。これは、端点に無限大を許す区間演算の実装には不可欠です。無限大は(inf,0)または(-inf,0)のような内部表現にしています。また、twosumは
def twosum(a, b) :
    x = a + b
    if (abs(a) > abs(b)):
        tmp = x - a
        y = b - tmp
    else:
        tmp = x - b
        y = a - tmp
    return x, y
のような実装になっていますが、a,bが無限大でなくa+bがオーバーフローした場合、(inf,-inf)を返すことが分かります。add_downのオーバーフロー処理は、
  • 負の方向へオーバフローした場合は、(-inf,0)を返す。
  • 正の方向へオーバーフローした場合は、下向き丸めなので無限大を返すのは誤りで、ddで表現できる正の最大数を返す。
という処理を行っていました。

自分はこれで抜けは無いと安易に考えていたのですが、研究室の修士の学生だった宮崎敏文さんにまずい場合があるのではないかと指摘され、確かに問題があることが分かりました。

問題は3つあります。
  1. 最初の問題は、「add_downで正の方向にオーバーフローするとddでの正の最大数を返すが、x1+y1は確かにddの正の最大数より大きいが、x2とy2が負の数で、全体の合計はddの正の最大数を下回ることがあるのではないか」というものです。宮崎さんの指摘はこの問題でした。
  2. 二つ目の問題は、「最初のtwosumの後は無限大の処理をしているが、二つ目のtwosumの後は何もしていない。最初のtwosumではオーバーフローしないが、二つ目のtwosumでオーバーフローする可能性があるのでは」というものです。
  3. 三つ目の問題は、「add_downの引数の片方がinf、もう片方が非infだったとき、その和はinfになるべきだが、ddの正の最大数が返ってくる」というものです。
これらの問題点を解決するため、次のように改修しました。
def dd_add_down(x1, x2, y1, y2) :
    if abs(x1) == float("inf") or abs(y1) == float("inf") :
        return x1 + y1, 0
    z1, z2 = twosum(x1, y1)
    if z1 == -float("inf") :
        return z1, 0
    if z1 == float("inf") :
        z1 = sys.float_info.max
        z2 = sys.float_info.max * 2.**(-54)
    down()
    z2 = z2 + x2 + y2
    near()
    z1, z2 = twosum(z1, z2)
    if z1 == -float("inf") :
        return z1, 0
    if z1 == float("inf") :
        z1 = sys.float_info.max
        z2 = sys.float_info.max * 2.**(-54)
    return z1, z2
すなわち、まず引数に無限大が含まれている場合を先に処理し(3.の対策)、次にtwosum後に正方向にオーバーフローして正の最大数に置き換えた場合はreturnせずにそのまま計算を継続するようにし(2.の対策)、最後に2回目のtwosumに対しても同じようにオーバーフロー処理を行なう(2.の対策)ようにしました。

これによって、次のように問題が修正されました。

1.の問題は、例えば(x_1,x_2)=(2^{1023}-2^{970},-(2^{969}-2^{916})), (y_1,y_2)=(2^{1023},-2^{969})という入力で発生します。
#include <kv/interval.hpp>
#include <kv/dd.hpp>
#include <kv/rdd.hpp>

int main()
{
    std::cout.precision(32);
    kv::interval<kv::dd> x, y, z;

    x.lower().a1 = pow(2., 1023) - pow(2., 970);
    x.lower().a2 = -(pow(2., 969) - pow(2., 916));
    x.upper() = x.lower();
    std::cout << x << "\n";

    y.lower().a1 = pow(2., 1023);
    y.lower().a2 = -pow(2., 969);
    y.upper() = y.lower();
    std::cout << y << "\n";

    z = x + y;
    std::cout << z << "\n";
}
これを改修前のkvで計算すると
[8.9884656743115780417662938029053e+307,8.9884656743115780417662938029054e+307]
[8.9884656743115790396864485702651e+307,8.9884656743115790396864485702652e+307]
[1.797693134862315807937289714053e+308,inf]
となりますが、正しい値は1.797693134862315708145274237317049\cdots \times 10^{307}なので下端が真値より大きくなってしまっています。改修後は、
[8.9884656743115780417662938029053e+307,8.9884656743115780417662938029054e+307]
[8.9884656743115790396864485702651e+307,8.9884656743115790396864485702652e+307]
[1.797693134862315708145274237317e+308,inf]
と正しく計算されていることが分かります。

2.の問題は、例えば(x_1,x_2)=(2^{1023},2^{970}), (y_1,y_2)=(2^{1023}-2^{971},2^{969}-2^{916})という入力で発生します。
#include <kv/interval.hpp>
#include <kv/dd.hpp>
#include <kv/rdd.hpp>

int main()
{
    std::cout.precision(32);
    kv::interval<kv::dd> x, y, z;

    x.lower().a1 = pow(2., 1023);
    x.lower().a2 = pow(2., 970);
    x.upper() = x.lower();
    std::cout << x << "\n";

    y.lower().a1 = pow(2., 1023) - pow(2., 971);
    y.lower().a2 = pow(2., 969) - pow(2., 916);
    y.upper() = y.lower();
    std::cout << y << "\n";

    z = x + y;
    std::cout << z << "\n";
}
これを実行すると、改修前は
[8.988465674311580536566680721305e+307,8.9884656743115805365666807213051e+307]
[8.9884656743115780417662938029052e+307,8.9884656743115780417662938029053e+307]
[inf,inf]
のように[inf,inf]という不正な区間が発生していましたが、改修後は
[8.988465674311580536566680721305e+307,8.9884656743115805365666807213051e+307]
[8.9884656743115780417662938029052e+307,8.9884656743115780417662938029053e+307]
[1.797693134862315807937289714053e+308,inf]
のように[ddの最大数、inf]という正しい結果を返すようになりました。

3.の問題は、区間演算では発生しないのですが、add_downを明示的に呼ぶと再現できます。
#include <kv/interval.hpp>
#include <kv/dd.hpp>
#include <kv/rdd.hpp>

int main()
{
        std::cout.precision(32);
        kv::dd x, y, z;

        x = std::numeric_limits<kv::dd>::infinity();
        std::cout << x << "\n";
        y = 0;
        std::cout << y << "\n";
        z = kv::rop<kv::dd>::add_down(x, y);
        std::cout << z << "\n";
}
これを実行すると、
inf
0
1.797693134862315807937289714053e+308
のようにinf+0の下向き丸めが正の最大数になってしまいますが、改修でちゃんと
inf
0
inf
と正常になりました。

以前からこのバグには気付いていたのですが、多忙でなかなか改修できず、モヤモヤした気分が続いていました。確率が低いとはいえ、「精度保証が保証になっていない」バグを放置してはならないと考えていますので。しかし、double-double絡みのバグが目立つのは困りものです。針の穴を通すような繊細さが要求されると感じています。

2018/05/01(火)ubuntu 18.04 インストール (リンク集)

2018/05/01(火)ubuntu 18.04 インストール(10) リモートデスクトップ

ubuntu 18.04にリモートデスクトップを入れたときのメモです。まだ暫定版で、今後変わるかもしれません。

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

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

ubuntuは18.04になるとき、Xサーバを17.10で導入されたWaylandからXorgに戻したそうで、これはVNCやxrdpとの相性を考えてのことだそうです。また、unityからgnomeに戻った、更にパッケージのxrdpのバージョンが最新になった、ということもあり、xrdpが簡単に使えるようになったのではないかと期待していました。

しかし、いろいろ試したところあまり上手くは行かず、それでも何とかしてみた、というのが以下の記録です。

基本的に、
の記事に従ってやってみました。
sudo apt install xrdp
これで、サーバそのものは簡単に入ります。カーソル回りに不具合があるらしく、/etc/xrdp/xrdp.iniで、
new_cursors=true
を、
new_cursors=false
に書き換える
sudo systemctl restart xrdp
が必要です。

上のサイトによれば、~/.xsessionrc に、
export GNOME_SHELL_SESSION_MODE=ubuntu
export XDG_CURRENT_DESKTOP=ubuntu:GNOME
export XDG_DATA_DIRS=/usr/share/ubuntu:/usr/local/share:/usr/share:/var/lib/snapd/desktop
export XDG_CONFIG_DIRS=/etc/xdg/xdg-ubuntu:/etc/xdg
と書くとgnomeが使えるとのこと。しかし、試してみると、本体でログアウトしておかないと使えないことが分かりました(少なくともうちの環境では)。本体使用中に接続しようとしてもすぐに画面が消えてしまいます。クライアントAで使っている時にクライアントBから接続すると、Aの画面が閉じてBで続きができる、という、windowsに似た動作になります。また、どこかのクライアントで使用中のときは、本体でログイン出来なくなってしまいます。

これでは不便なので、16.04のときと同じく、「MATE」を使う作戦を試してみました。MATEは、「メイト」ではなく「マテ」と読み、gnome2の操作性で軽く、見た目を重視しているということで最近人気があるデスクトップ環境です。
sudo apt install ubuntu-mate-desktop
インストール中、display managerをgdm3とlightdmのどちらにするか聞かれましたが、とりあえずMATEの標準と思われるlightdmにしました (gdm3でどうなるかは未確認)。再起動し、本体の方でログイン時に「MATE」と「Ubuntu(デフォルト)」(gnome)が切り替えられ、どちらでも正常に使えることを確認しました。

まず、/etc/xrdp/startwm.sh を書き換えて、最後の2行のXsessionを起動している場所の前に
exec mate-session (これを挿入)

test -x /etc/X11/Xsession && exec /etc/X11/Xsession (これは元から)
exec /bin/sh /etc/X11/Xsession (これは元から)
のようにMATEの起動コードを挿入します。これで一応動きましたが、起動時に「Could not acquire name on session bus」と変なウインドウが出てしまいます。これは、mate-sessionの起動前に
unset DBUS_SESSION_BUS_ADDRESS
を挿入したら直りました。

16.04のときと違ってキーボードは正常に使えます(xrdpが新しいおかげ)。かな漢字変換は作動しないので、16.04のときに倣ってmate-session起動前に
export GTK_IM_MODULE=ibus
export QT_IM_MODULE=ibus
export XMODIFIERS="@im=ibus"
ibus-daemon -d
を挿入したらうまく動きました。

なお、xrdp接続するとthinclient_drivesというフォルダがホームにでき、接続が切れるとこれがアクセスできないフォルダになってlsの度に警告が出て大変鬱陶しいです。一応ホームで
sudo umount thinclient_drives
とすれば直りますが、面倒です。これを何とかしようと、/etc/xrdp/sesman.iniで、
FuseMountName=thinclient_drives
を、
FuseMountName=.thinclient_drives
に変更
sudo systemctl restart xrdp
として見えなくしました。

18.04でのxrdpに対する期待とは裏腹にかなり面倒なことになってしまいましたが、一つの例として記録を残しました。もう少し賢い方法がありそうな気もするので、少し様子見ですかね。

2018/05/01(火)ubuntu 18.04 インストール(9) その他

その他入れた細々としたもの。
sudo apt install openssh-server
これでsshログイン出来るようになります。
sudo apt install git
後は個人的に必要なもの。
sudo apt install lv
sudo apt install checkinstall
以前あった、bashでccとc++に対してファイル名補完が効かない問題は修正されました。

ホームディレクトリに作られる「ダウンロード」などのディレクトリが日本語だと何かと不便なので、英語表記に直します。
LANG=C xdg-user-dirs-gtk-update
として、「Don't ask me this again」をチェックして「Update Names」をクリックします。これでホームディレクトリが
examples.desktop  テンプレート  ドキュメント  ピクチャ      公開
ダウンロード      デスクトップ  ビデオ        ミュージック
から
Desktop    Downloads  Pictures  Templates  examples.desktop
Documents  Music      Public    Videos
に変わりました。元に戻すにはLANG=Cなしで単に「xdg-user-dirs-gtk-update」。

さらに、ふと目を離すとロックされてパスワードが要求されるのが嫌なので、設定→プライバシー で画面ロックをオンからオフに変更しました。

また、警告音がうるさいので、設定→サウンド→音響効果 をオフにしました。

これで大体日頃使っている環境が出来た気がします。Intlabやopenblas関連については後ほど。

追記(5/24)

画面上部に表示されるのが時刻のみで残念だったので、
gsettings set org.gnome.desktop.interface clock-show-date true
sudoなしで実行して、日付と時刻が両方表示されるようにしました。
gsettings reset org.gnome.desktop.interface clock-show-date
で元に戻せます。
OK キャンセル 確認 その他