CHOLMOD とユーザーマニュアル

CHOLMOD のパラメータを修正しようと試みたところ、ユーザーマニュアルを見ただけではよく分からず、結局ソースコードを解析して、修正するべきパラメータを特定した。
ただ、後になってみてユーザーマニュアルを見てみると、ぱっと見では分からないところに書かれていた。
CHOLMOD はいろいろと難しい。


今日の作業内容：SDPA-C 3h
今日のランチ：味庵 揚げ豆腐と野菜の炒め
明日の予測作業時間：5h

発表準備とSDPA-Cの実行

まだ直前ではないが、スケジュールが分かりづらいので、ある程度発表準備をしておいた。
だいたい８割ぐらい完成といったところ。
ざっと通しで時間を測ると20分目標のところ、18分ぐらいだったので、これで練習をしつつ、ゆっくりとしゃべれるようになればだいたい良さそうな感じである。

次に、SDPA-C を max-clique-problem で性能を確認してみたところ、
SDPA-C 6 が 2417 秒かかっていたところを 1500 秒で解けるところまで来ている。
ただ、まだ特にチューニングなどをしていないので、ここからチューニングが必要そうでもある。

今日の作業内容：発表準備 2h + SDPA-C 1h
今日のランチ：らく　鶏の照り焼き定食
明日の予測作業時間：5h

バグを取り除けた

ここまでのバグは、X のコレスキー分解のところに問題があった。
コレスキー分解は、通常は左上の要素から分解する必要があるが、SDPA-C の主問題側の分解は右下から行う必要があった。
これの方向をあわせて、順列などを再確認したところ、ある程度の問題は解けるようになった。
これからパフォーマンスチェックに入ることになる。

今日の作業内容：SDPA-C 3h
今日のランチ：つかさ　〆サバ
明日の予測作業時間：3h

バグ取れず２

コンパイルエラーや segmentation fault などのバグは取り切ったが、まだバグが残っているので、それを修正する必要がありそうだ。

現在の段階では、
(1) CHOLMOD チェック用のデモ
(2) SDPA-C プロトタイプ (Matlab 版)
(3) SDPA-C 本体
の順番で確認作業をしているが、(1) は正確にできているが、(2) の時点でバグが混入している様子である。
もう一度、(1)と(2)の間を確認しなおしたほうが良さそうである。

今日の作業内容：SDPA-C 4h
今日のランチ：信華園 レバー野菜炒め
明日の予測作業時間：4h

順列のバグが取れた

いろいろと調べた結果、順列についてのバグは取れた。
CHOLMOD の場合、順列についてのどう対応しているかがマニュアルなどに説明がないので、ここに書いておくことにする。

A x = b を解く場合、
matlab なら
[pp] = amd(A);
[dummy1, post] = etree(A);
pp = pp(post);
P = eye(n);
P = P(pp,:);
として順列 pp とそれに対応する行列 P が作れる。
この pp に対応するのが、CHOLMOD の cholmod_factor に入っている Perm である。

Ax = b については、
PAP^T = LL^T
となるようにコレスキー分解が入るので、
x = P^T*L^{-T}*L^{-1}*P*b
のように行われる。

したがって、CHOLMOD で Ax = b を解くには、

cholmod_dense* b2 = cholmod_solve(CHOLMOD_P, L, b, common);
cholmod_dense* b3 = cholmod_solve(CHOLMOD_LDLt, L, b2, common);
cholmod_dense* x = cholmod_solve(CHOLMOD_Pt, L, b3, common);
のようにして解くことになる。
（ここの CHOLMOD_P, CHOLMOD_Pt を書けるところが CHOLMOD のデモプログラムには明記されていない）

今日の作業内容：SDPA-C 4h
今日のランチ：やぶ A定食マグロの味噌焼き
明日の予測作業時間：4h

バグ取れず

順列のところのバグが簡単に取れるかと思いきや、案外手こずっている。
もういちど CHOLMOD の順列から、どういうふうに各行列に影響が出ているのかを確認しなおした方が速いのかもしれない。
まだまだ考える必要がありそうだ。

今日の作業内容：SDPA-C 3h + プレゼン準備 1h
今日のランチ：和泉庵 かつ丼とお蕎麦のセット
明日の予測作業時間：3h

バグに次ぐバグ

この前のステップ長がうまく計算できていないのは、CHOLMOD の順列を直すべきところがもう一か所あったためで、これを直したら問題なく実行できるようになった。
これで、clique がひとつになっているときは実行できることが分かった。
ただ、二つのcliqueに分解されるときに、まだバグが入っているようである。
次は、このバグを取ることを検討したい。

今日の作業内容：プレゼン資料作成 2h + バグ取り 3h
今日のランチ：花時計 生姜焼きと海老フライ
明日の予測作業時間：5h

cholmod_solve は自動的にはorderingを調整しない

CHOLMOD のデモプログラムを見るかぎりだと、cholmod_solve は右辺ベクトルと解ベクトルにordering を自動的に適用するように見えるが、調べてみるとそうなっていない。
このあたりは、ユーザマニュアルには書かれていないので、実際に実行して数値的に確認する以外に方法はない。
おそらく、ユーザマニュアルが実際のプログラムをきちんと反映できていないのであろう。
ユーザマニュアルを整備するのはプログラムを書くのと同じぐらい大変なので、このあたりの事情はよく分かる。

とりあえず、ここのところは修正したが、なぜかステップ長が本来の長さよりも短く算出されてしまう。
6次元の行列程度の問題なのに80反復ぐらいかかってしまう。
このバグを取り除くのが次の仕事になりそうだ。

今日の作業内容：SDPA-C 5h
今日のランチ：サイゼリヤ みぞれハンバーグ＆牛肉コロッケ
明日の予測作業時間：4h

example1, example2 は解けるようになった

ここまでのバグ取りで example1, example2 は解けるようになった。
昨日の問題点は、連立方程式を解く順番を変えたところ、問題なく解けるようになった。
今日のところでは、A が対称行列のときに x'*A*y を計算するときに、Aが対称なので上半分で２倍すればOKかと勘違いしていたが、これは x と y が異なるときだと正しい値にならない。
これを修正したところ、 example2 までは解けるようになった。

この次は mcq.1.dat-s という最大クリークの問題であるが、これは行列の Permutation がかかっており、この段階でバグが起きているようである。
なにはともあれ、それなりに前進したと考えられる。

今日の作業内容：SDPA-C デバッグ 4h
今日のランチ：らく 鶏の照り焼き定食（トマトソース）
明日の予測作業時間：5h

間違っていたかもしれない

プログラムのデバッグの最中で気がついたが、対称行列 X について
inv(chol(X,'lower')) と chol(inv(X),'lower')
では値が異なる。
これを検討しなおす必要がありそうだ。

今日の作業内容：SDPA-C デバッグ 3h
今日のランチ：味庵 バンバンジー
明日の予測作業時間：5h

emacs でのカーソルの色を白にする

emacs での設定で、目に優しいように背景を暗くしている関係でカーソルを白にしていたが、どうやら C-x 5 2 で新しいフレーム（ウィンドウ）を開いたときに黒になってしまう。
これを解決するには、~/.emacs.d/init.el に

(setq default-frame-alist
      (append (list '(cursor-color . "white")) default-frame-alist))
とすればよい。
普通に
 (set-cursor-color "white")
だけだと、最初のフレームは白になるものの C-x 5 2 で開いたフレームには受け継がれない。


SDPA-C のほうは、メインループを一通り書き終えて、数値があっているかをプロトタイプとして組んでいる Matlab と整合しているかを確認している。solveInfo がなぜか pdINF になってしまうので、ここに到達するまでにバグが入ってしまっているようだ。

今日の作業内容：SDPA-C 2h, SNL 1h
今日のランチ：シッダルータ キーマカレー
明日の予測作業時間：4h

Lancoz 法の確認

ステップ長の計算は Lancoz 法に依っているが、Dualは、cholmod のルーチンを使いまわせば、余計な一時配列を準備しなくても実行できそうな目処が立った。
逆に、primal 側は、各 clique ごとになるから、こちらは一時配列を確保しないといけないと考えている。


今日の作業内容：SDPA-C 1h, プレゼン準備 4h
今日のランチ：つかさ　サーモンハラス
明日の予測作業時間：5h

Lanczos 法の復習など

今日は、Lanczos 法の復習などに時間をあてた。
このあたりは、実装自体は比較的難しくないが、計算のためのワーキングスペースが必要であって、それをどのようにやりくりするか、がポイントとなる。

今日の作業内容：SDPA-C 1h + プレゼン資料作成 2h
今日のランチ：信華園 汐焼きそば
明日の予測作業時間：5h

探索方向まで計算

今日のところで、Schur 方程式を解いて探索方向(Predictor)の計算までは進んだ。
ただ、初期点は対角行列であるため、完全には確認できておらず、まだ油断は大敵である。
ここから次は探索ステップの計算になるので、Lanczos 法の実装に関係することになる。

今日の作業内容：SDPA-C 3h + プレゼン準備 1h
今日のランチ：角笛 チキントマト煮
明日の予測作業時間：3h

Schur 補完方程式のところを組む

まだデバッグしながらであるが、Schur 補完方程式の両辺までは計算できるようになった。
やはり、Schur 補完行列あたりは、SDPA7 のルーチンを使いまわせる部分があるので、それなりに楽でもある。
これで、Schur 方程式を解いて、どう進むかを確認していくことになりそうだ。

今日の作業内容：SDPA-C 4h + プレゼン準備 1h
今日のランチ：らく 鶏の照り焼き定食
明日の予測作業時間：3h

bool で返り値がない場合は false になりやすい

C++ で bool を返り値としたときに、うっかりと return で true か false を返さなかった場合には、不定のあたいとなるため、false になりやすい。
今日は、このバグで手間取った。

あと、今日でコレスキー分解のルーチンをとりあえず実装した。
前に別途確認したのであっているとは思うが、最初の反復は対角行列だけなので、これは2回目の反復以降でもう一度デバッグする必要がありそうだ。

今日の作業内容：SDPA-C 3h + プレゼン準備 2h
今日のランチ：シッダルータ ベジタブルカレー
明日の予測作業時間：5h