QAP relaxation のつづき

今日はまとまった時間がとれなかったため、QAP を内点を持つように変形する方法をもう少し考えた。
またアイデアをもらったので、Matlab で簡単なプログラムを作ってチェックしているが、どこかに間違いがあるらしく、固有値にマイナスが含まれてしまう。
明日に続く展開になりそうだ。

今日の作業内容：QAP の変形
今日のランチ：味庵　白身魚のマーボー煮込み
明日の予測作業時間：5h

QAP relaxation

昨日のグラフ分割問題の制約緩和については、あのあともう少し考えてみたら簡単な数式で内点を表現できることが解かった。
SDPLIB の場合、内点がない問題は他に QAP の系統があるため、今日はこれを確認した。

問題の定式化については、いろいろと教えてもらったので、ほぼ大丈夫であるが、やはり内点を構成するのが難しそうである。グラフ分割と異なってストレートな制約緩和ができそうにないので、内点を先に構成してから、それを満たすような自然な制約緩和を入れるほうが良さそうである。

今日の作業内容：QAP relaxation
今日のランチ： 四川　汁無し坦々麺
明日の予測作業時間：2h

gpp relaxation

SDP 緩和のひとつにグラフ分割の問題がある。
SDPLIB にもグラフ分割の問題がいくつか含まれている。(gpp???.dat-sというもの)
グラフ分割から生成される SDP には、内点法から見ると大きな欠点があり、それは実行可能内点を持たないため、数値的に不安定になりやすい、ということである。

そこで、今日は実行可能内点をなくしている制約をすこし relaxation して実行可能内点を持たせてみることにした。これで、実際に SDPA-GMP でもきちんと解けている。
そこで、実行可能内点を数式で簡単に表現できるか？と思ってトライしてみたが、これが意外とストレートでないことがわかった。
これについては、もう少し調べてみようと思う。

今日の作業内容：gpp relaxation
今日のランチ：つかさ　サーモンハラス焼き
明日の予測作業時間：4h

MFPR と machine epsilon

mpackを使っての SDP を生成するプログラムを作成。
ただ、途中でMFPRの動きがすこし怪しい感じがある。

たとえば、1+e > 1 となる最小の e を求めると prec = 600 のときに 1.0e-200 程度であり、これが実質的な machine epsilon に相当するはずだが、mpfr_printf で表示させると 600 桁に近い情報を持っていることが表示される（足し算は600桁近い情報で行っている）。
つまり、四則演算に使っている桁数と比較演算に使っている桁数にずれがあるようなのである。
いずれにしても、SDP の生成の場合は、それを見越して prec を設定する。

SDPA-GMP で解くときを考えると、SDPA-GMP のパラメータを epsilon = 1.0e-20 とするなら、その約 30 倍の prec = 600 程度であれば桁数にゆとりを持って生成できる。epsilon = 1.0e-80 なら prec = 2400 である。もちろん、SDPA-GMP の パラメータの precision も調整しておく必要がある。 epsilon の6倍程度の桁数が必要である。epsilon = 1.0e-20 なら、precision = 120, epsilon = 1.0e-80 なら precision = 480 である。

SDPA-GMP の場合、出力ファイルに relative eps の行で machine epsilon を表示してくれるので、これの桁数が epsilon の桁数の倍以上の精度になるように precision を設定する必要がある。


今日の作業内容：MFPR, センサネットワーク
今日のランチ：　らく　チキンカツ
明日の予測作業時間：5h

MPFR と mpreal に苦戦

昨日に続いて MPFR でプログラムを組んでいるが、これの拡張版の mpreal クラスが扱うのがなかなかに難しい。
簡単にいうと、
printf("a = %d\n", a);
b = a;
printf("b = %d\n", b);
のようなことをすると、aとbで違った表示がされる。
MPFR が原因なのか、あるいは mpreal なのか解からないが、不思議な現象である。
このとき、a の桁数が打ち切られて b に入っているようでもあり、よく解からない。
とりあえずは、あとで valgrind をかけてメモリチェックをかけて、それから詳細を確認しようと思う。

それにしても、Matlab で１０秒でかけるコードを mpack で実装しようとすると５分はかかる。
mpack はマニュアルが未整備なので、関数を使うときにはソースコードを直接参照してから利用している。mpreal なども確認するためにソースコードの動きをざっと確認したりしており、これが結構時間がかかっている。

まだまだだいぶかかりそうだ。

今日の作業内容：MPFR プログラミング
今日のランチ：さか本　野菜天丼とおそばセット
明日の予測作業時間：4h

MPFR の sqrt

SDPA の関係で多倍長計算の MPFR を調べているが、sqrt の挙動でよくわからない点がある。

例えば、
  mpfr_set_default_prec(prec);
  mpf_class mtmp2 = 2.0;
  mtmp2 = sqrt(mtmp2);
  gmp_printf("sqrt(2)   in mpfr  %.80Fe\n",mtmp2.get_mpf_t());
  mpf_class residual = mtmp2 * mtmp2 - 2.0;
  gmp_printf("residual  in mpfr  %.80Fe\n",residual.get_mpf_t());
とした場合、ビット数 (prec) = 300 で
prec = 300
sqrt(2)   in mpfr  1.41421356237309504876000000000000000000000000000000000000000000000000000000000000e+00
residual  in mpfr  -1.07199229837017405496000000000000000000000000000000000000000000000000000000000000e-19
となるが 2000 にあげても
prec = 2000
sqrt(2)   in mpfr  1.41421356237309504876000000000000000000000000000000000000000000000000000000000000e+00
residual  in mpfr  -1.07199229837017405496000000000000000000000000000000000000000000000000000000000000e-19
と同じ結果が出てきている。
sqrt(2) の表示自体は桁数の表示の仕方に問題があるのかもしれないが、residual が変わらないのはおかしいような気がする。

これについては、明日続けて調べる必要がありそうだ。

今日の作業内容：MPFR チェック
今日のランチ：味庵 豚肉のキムチ炒め
明日の予測作業時間：5h

mpack のコンパイル

SDP の精度の計算用の問題を生成するのに、やはり精度の高い問題が必要と判断して、mpack を利用することにした。
ただ、かなりのコンパイル時間がかかるようで、Core2 でも５時間以上かかるようだ。
これについては、明日またチェックをする。

あと、昨日あったSDPAの細かいバグ（？）は、今日取ってもらった。
次のバージョンアップのときに公開しているソースに反映される。

SwissRoll のは久しぶりに計算をしてみたが、ある意味で致命的なミスに気がつく。
これは修正は難しいかも知れない。

今日の作業内容：mpack コンパイル
今日のランチ：　鶏こまち　鶏そぼろの蒸篭蒸し
明日の予測作業時間：３h