2012年2月29日水曜日

OpenBLAS のバグを回避

今日、OpenBLAS をコンパイルしていて気がついたが、
error: expected '=', ',', ';', 'asm' or '__attribute__' before
というエラーでコンパイルができないことがある。
このときには、common.h の最後のあたりで、

#ifdef __cplusplus
}
   
#endif  /* __cplusplus */
となっているところにセミコロンをつけて

#ifdef __cplusplus
};
   
#endif  /* __cplusplus */
としておいて、2回makeするとコンパイルできるようになる。
(なぜ2回 make しないといけないのかは、まだ調べていない)

GotoBLAS のときにも似たようなことが起きていたかもしれない。

今日の作業内容:SparsePOP チェック 3h + 証明 2h
今日のランチ:食堂 ライス+味噌汁+豚生姜焼き&から揚げ+春のバランス惣菜
明日の予測作業時間:5h


2012年2月28日火曜日

証明、行き詰り

先週からトライしている証明は、まだまだ分からず。
反例のようなものができたときに数値実験で確認しているが、反例ではないことが解ってしまう。

やはり、3x3以上の半正定値行列をいかに手軽に扱うか、がポイントになってくるようだ。

今日の作業内容:証明の続き 3h
今日のランチ:シッダルータ ベジタブルカレー
明日の予測作業時間:5h

2012年2月27日月曜日

なぜ主双対内点法の反復回数は20回程度なのか

主双対内点法の反復回数は、理論的にはO(sqrt(n) log (eplison)) となっている。
ところが、n >= 10000 で epsilon = 10^{-8} でも、SDPA などのソフトウェアの場合、たいていは 20 回程度、多くても 50 回程度で終わる。

なぜか。

理由は簡単で、理論的な解析が効く前に閾値に到達してしまうことにある。


たとえば、ステップ長がすべての反復で 0.5 以上できたとすると、一回の反復で実行可能性は2分の1以下になる。
そこで、(1/2)^k = 10^{-8} を k について解くと、閾値 10^{-8}を実行可能性が下回るまでの反復回数が分かるが、これが 26.5754 回なのである。

反復ではロングステップとショートステップがあって、ショートステップの方が理論的には少ない反復数であるとされているのに、実際に計算するとロングステップの方が速いのも、閾値からの影響が大きい。


コンピュータソフトウェアの場合、double など有限桁数で計算しているため、無限桁数で仮想している理論的な解析で高速なものが実際の計算でも高速とは限らない。



今日の作業内容:証明の続き 2h
今日のランチ:らく 焼き魚定食
明日の予測作業時間:4h

2012年2月24日金曜日

SDPA が testing から削除されてしまった

SDPA の mentors への dput がうまくできていない状態であるが、この段階で SDPA パッケージが testing から削除されてしまった。

つまり、testing にあったパッケージ (7.3.5.dfsg-1) では、MUMPS のバージョンアップに対応できておらず、build に失敗してしまうため、削除されたようである。

解決する手順としては、7.3.6.dfsg-2 を mentors に dput することである。
古いファイルがよくないのかもしれないと思い、SDPA を mentors から一度削除してみたが、それ以降も dput に失敗してしまう。
開発環境は Debian sid で apt-get dist-upgrade を施してあるのだが、これが 7.3.6.dfsg-1 をアップロードした後でどこかのパッケージにバグが混入しているのかもしれない。

ただ、mentors のアカウントの方に問題があるかもしれないので、mentors サーバーのサポートにもアカウントを調べてもらえるようにメールでお願いをしておいた。
やはり、エラーメッセージが出てこないので、どのあたりを調べたらいいのか、を調べる必要がある。

今日の作業内容:Debian 1h + 証明の続き 5h
今日のランチ:ちゅらさん ソーキのトマト煮
明日の予測作業時間:5h

2012年2月23日木曜日

証明難しい

昨日の続きで SNL の定式化を見ているが、簡単そうに証明できそうなことがあったので、証明しようとしたら案外難しいことが分かった。
もともとのネットワークが連結でない場合は証明できないことが分かったが、そうでないときに成り立つのかどうかが難しい。
今日はこればかりだったが、もう少し考えてみる必要がありそうだ。


今日の作業内容:証明検討 5h
今日のランチ:角笛 鶏のささみ定食
明日の予測作業時間:5h

2012年2月22日水曜日

SNLのSDPをチェック

SNL の SDP 緩和の定式化を改めて見直しているが、これまでにあまり研究が進んでいない点があるので、そのあたりがどうなのかを式で変形してみている。

SFSDP だと本格的過ぎてトライアンドエラーをしづらいので、Matlab で 200 行ぐらいの基本的な実装をして、数値的にも確認できるようにしている。

とにかく、気になる点が4つほどあるので、これを一つ一つ片付けていきたい。

今日の作業内容:SDP緩和検討 4h
今日のランチ:つかさ サーモン照り焼き
明日の予測作業時間:5h

2012年2月21日火曜日

Debian の mentors にアップロードできず

SDPA の Debian パッケージを修正しており、mentors にアップロードしようと試みているが、なぜかうまくいっていない。
エラーメッセージもなく、「アップロードに成功」というメッセージが出るのであるが、実際にはアップロードができていない。
原因特定には、まだ時間がかかってしまいそうだ。

今日の作業:Debian パッケージ 2h
今日のランチ:らく 鶏の照り焼き定食
明日の予測作業時間:4h

2012年2月20日月曜日

POP の制約式の表現方法

この前の簡単に解けない問題については、制約式の表現方法を変えることで、sparsePOPで比較的にいい解が得られることが分かった。

つまり、x_i^2 = 1 という制約をそのまま SDP 緩和に入れるのではなく、
x_i^2 >= 0.95, -1<=x_i<=+1
のようにすると、比較的安定して解くことができる。
ここで、x_i^2 = 1 の解が -1,+1 の2種類しかないので、x_i^2>=1をx_i^2>=0.95 に少し緩めておいても、出てきた答えのx_i の符号を見れば -1 か +1 に修正しやすいことが要因であると思われる。

多項式制約問題は SDP を解くときよりも、どのように定式化するかに依存しやすく、その点に注意が必要だが、どのように定式化すべきかを体系的にまとめた内容はそれほど知られていないと思われるので、研究の一つになるかもしれない。

今日の作業内容:POPの制約チェック 2h + SDPA Debian パッケージ対応 2h
今日のランチ:味庵 鶏の唐辛子炒め
明日の予測作業時間:5h

2012年2月17日金曜日

簡単に見えて、どのソフトでも解けない問題

今やっている問題が解けないので、基本的な制約だけを残して定式化してみた。
それでもsparsePOP では解けなかった。
 NEOS に載っているほかのソフトでも試してみたが (GAMSを受け付けるものを試してみた)
BARON, LINDO, Gurobi, MOSEK, XpressMP, CONOPT, Ipopt, KNITRO, SNOPT, Bonnin, Conenne, PATHNLP
でも解けなかった。

問題は x \in R^5 が変数であって、
min c^T x : subject to e^T x = 1, x_i^2 = 1 (i=1,..,5)
というものである。
実行可能解としては、x_1=x_4=-1, x_2=x_3=x_5=1 があるので、実行不可能ではない。
しかし、NEOS のソフトの場合、たいていは Locally Infeasible あるいは Infeasible で終了してしまう。
sparsePOP の場合は、relaxOrder = 1 (通常の SDP 緩和)だと解くには解くが、x_1=x_4=0 などとなってしまい緩和が弱く、relaxOrder >= 2 だと SDP infeasible となる。

実際に使った GAMS のファイルは以下の通りである。
========


Variables x1,x2,x3,x4,x5,objvar;
Positive Variables x1,x2,x3,x4,x5;
Equations e1,e2,e3,e4,e5,e6,e7;
* e1.. 0.2*x1 +0.7*x2 +0.1*x3 +0.9*x4 +0.4*x5 + objvar =E= 0;
e1.. objvar =E= 0.8*x1 +0.7*x2 +0.1*x3 +0.9*x4 +0.4*x5;
e2.. x1 + x2 + x3 + x4 + x5 =E= 1;
e3.. x1*x1 =E= 1;
e4.. x2*x2 =E= 1;
e5.. x3*x3 =E= 1;
e6.. x4*x4 =E= 1;
e7.. x5*x5 =E= 1;
Model m / all /;
m.limrow=0; m.limcol=0;
$if NOT '%gams.u1%' == '' $include '%gams.u1%'
Solve m using NLP minimizing objvar;

==
おそらくGAMSファイルは間違っていないはずで e2 の =E= 1 を =E= 5 にすると、最適解x_1=x_2=x_3=x_4=x_5 = 1 が得られる。
しかし、=E= 1 のままでは解くことができない。

もちろん、x_i^2 = 1 を x_i = {-1,1} にして Mixed Integer にすれば CPLEX などで簡単に解けるのであるが、本来解こうとしている問題にするために条件を追加する時点でうまくできなくなってしまう。

数理最適化の分野はいろいろなソフトが開発されているが、まだ完全とは言えないようで、だからこそ研究余地が大きく残っているともいえる。

今日の作業内容:SparsePOP, NEOS 確認 5h
今日のランチ:シッダルータ キーマカレー
明日の予測作業時間:2h

2012年2月16日木曜日

実行可能解が一つしかない

SDPA関係でオーストラリアから来ていた問い合わせが、ここ最近になってまた来るようになったので対応しているが、どうも数学の定式化のあたりがうまくできていないようだ。
制約条件を調べなおしたところ、実行可能解が一つしかないことが分かり、これが自動的に最適解にならざるを得ない。
おそらくモデル化で数式に変換できていないようである。

分野が違いすぎる関係で数式に落とす前の内容が分からないので、もう少し情報をもらう必要があるかもしれない。

今日の作業:SDPA 対応 3h + SDPARA チェック 2h
今日のランチ:信華園 トマト酸辣湯 \750
明日の予測作業時間:5h



2012年2月15日水曜日

sparsePOP のC++版を Ubuntu でコンパイル

sparsePOP の C++ 版を Ubuntu でコンパイルしたが、これが結構難しかった。
というのも、README に書いてある方法では、SDPA をコンパイルするときに同時にコンパイルされる MUMPS にリンクするようにしてあるのだが、MUMPS自体は apt-get でインストールしたものがあるので、これにリンクしようとしてライブラリ名がかみ合わずに失敗していた。

そもそも、MUMPS は single 版と並列計算用の parallel 版があって、single 版の場合には MPI のラッパーをかませることで基本的には対応している。
これが、Ubuntu だと libdmumps_seq.a にリンクすることになっており、標準的な libdmumps.a libpord.a は parallel 版のライブラリになっている。
つまり、libdmumps.a にリンクしてしまうと、MPI 関係のエラーが出て実行ができない。

このあたりは注意が必要なのだが、そもそも MUMPS を apt-get 由来のものにリンクするには configure.in を書き換えて autoreconf をかませる必要があって、ほとんどの人はこんなに面倒なことをしない。

今日の作業内容:sparsePOP のコンパイル 1h + SDPARA チェック 2h
今日のランチ:つかさ ヒラメのから揚げ
明日の予測作業時間:4h

2012年2月14日火曜日

SDPA のチェックポイント

SDPA では、途中までの解を出力しておいて、それを初期点として次の計算に利用できる。

たとえば、param.sdpa で maxIteration を 5 にすると
$ ./sdpa example.dat-s example.out
で example.out には 5 反復目までの解が入力される。
ここで、example.out の先頭行から xVec とある行までを削除したものを init.dat とすると、
$ ./sdpa example.dat-s example5.out init.dat
とすると、実質6反復目から計算できる。

ただし、この場合は、theta などの一部の情報が引き継がれないので、これらがない状態での計算となり、ずれが生じる。
これらの情報を引き継ぐためには、ファイルに書き出して読み込み直す、という作業が必要となる。

今日の作業内容:SDPA チェックポイント 2h + 論文読み 1h
今日のランチ:角笛 さわらのから揚げ
明日の予測作業時間: 5h

2012年2月13日月曜日

unit disk graph

SNL関係の論文を読んでいたところ、unit disk graph という名前がついたグラフがあることが分かった。

Clark, Colbourn, Johnson, "Unit disk graphs", Discrete Mathematics, Vol 86 (1–3) pp 165–177, 1990

によると、5以上の奇数サイクルを持つグラフは unit disk graph であるという。
この性質はある最適化問題では重要な性質になっていることが知られている。

最適化からみると、予想外のところを結びつける内容なのかもしれない。

今日の作業内容:論文読み 2h
今日のランチ:四川 豚のもやし炒め
明日の予測作業時間:4h

2012年2月10日金曜日

A brief history of lift-and-project

整数計画問題を解く上で lift-and-project は重要な手法のひとつであるが、この手法の提案をした人が当時を振り返ってまとめた文章がある。
それが、

A brief history of lift-and-project
http://www.springerlink.com/content/b4h2n873j31t1448/
である。

これを読むと、legendary という形容詞がついている Egon Balas らと出会うところから始まり、Pittsburgh への異動など、研究の裏側でどのような生活があって、どのような出会いがあったか、などが分かって興味深い。

また、これを読むと、研究には運がつきもののようである。
一見うまくいっていなかった研究と結びついたりもするし、タイミングがよくなければダメであったりもする。

今日の作業内容:タンパク質データの読み方の勉強 2h + 論文探し 2h
今日のランチ: ちゅらさん しまなーチャンプル
明日の予測作業時間:4h

2012年2月9日木曜日

1F39

この前から調べていたたんぱく質の 1F39 だが、どうやら構造自体が難しいようだ。
pymol の cartoon だと alpha-helix, beta-sheet などの基本的な構造が表示されるが、1F39 の場合にはこれにカバーされていない領域が固まって存在している。
この部分が問題を難しくしているようだ。

今日の作業内容:たんぱく質データ確認 2h
今日のランチ:つかさ カンパチまかない丼
明日の予測作業時間:4h

2012年2月8日水曜日

たんぱく質データをチェック中

SNLの計算のために pymol で見たり、たんぱく質データをチェック中である。
pymol での操作は、リファレンスカード
http://s3.myopenarchive.org/season1_doc_org/94.pdf
も役に立つが、これは一通り操作を知っている人がコマンドを思い出すときのものであって、素人の自分には難しかった。
あと、特定の原子間の距離がどの程度あるのか pymol で表示させるのは難しく、
http://ameblo.jp/yare-yare-78/theme-10020442933.html
に載っている方法を利用するが、どうも特定の原子をマウスで指定するのが難しい。
(原子についている番号で指定できるようだが、3次元表示しているとどの原子が何番で、どう書けば指定できるのかがいまひとつ分かりにくい。)

結局のところ、pymolの場合は、マニュアルを見たりするよりも、ウェブで検索して試行錯誤をするのが一番早いことが分かった。



ところで、今扱っている13個のデータだが、調べてみると面白いことがあって、オングストロームに対してほぼ同じ距離で原子がある。
どういうことかというと、1オングストローム以内の原子ペアの数、2オングストローム以内の原子ペアの数 ... 7オングストローム以内の原子ペアの数を数え上げたときに、これらの原子ペアの総数に対して、それぞれのオングストローム内の原子ペアの数がほぼ同じなのである。

つまり、7オングストローム内でお互いに存在する原子ペアの数を N として、aオングストローム以内の原子ペア数を N(a) とすると、N(a)/N (ただし、a=1,...,7) が13個のたんぱく質でほとんど同じ値になる。
これは不思議な現象でもあるので、PDBから必要なデータを取り出すルーチンでそのようにしているのかを確認したほうがよさそうだ。

今日の作業内容:論文読み 1h + たんぱく質データ計算 1h
今日のランチ:さか本 とろろめし丼膳
明日の予測作業時間:4h

2012年2月7日火曜日

昨日の論文の続き

今日は時間が比較的細切れだったので、昨日の論文読みの続きをしておいた。
一通り論文は読み終わったが、どうも Lemma 1 の r \ge 2 l sqrt(d) が今一つ正しくない、あるいはかなり大きく見積もっているようにも見える。
このあたりは、もう少し検討してみたい。

今日の作業内容:論文読み 2h
今日のランチ:味庵 鶏の黒胡椒炒め
明日の予測作業時間:2h

2012年2月6日月曜日

SNLの論文をチェック

今日は細切れの時間を利用して、SNLの新しめの論文をチェックし始めた。
読んでいるのは、
http://arxiv.org/abs/1010.2262
であって、アブストラクトからすると理論的な側面が強いようだ。
まだ前半部分だけなので、どのあたりが今後活用できそうかを検討しながら読んでおきたい。

今日の作業内容:論文読み 2h
今日のランチ:らく 豚生姜焼き
明日の予測作業時間:1h

2012年2月3日金曜日

論文斜め読み

今日は時間が比較的細切れとなったこともあって、5本ぐらいの論文を斜め読みした。
先週までに、ある程度基本的な論文を時間をかけて読んであるので、斜め読みにしても重要なキーワードなどで突っかかることがなくなりつつある。
今日調べた論文は、基本的な論文の拡張であって、自分の検討したい拡張とは別の方向性であったが、このあたりの知識も少しずつ増やしておくと、あとでアウトプットをするときに役に立ちそうである。

今日の作業内容:論文斜め読み 2h
今日のランチ:ちゅらさん タコライス
明日の予測作業時間:2h

2012年2月2日木曜日

Git の本

SDPA のソースコードの管理は Git で行っているが、Git についてまとめられている本がある。

[Gitによるバージョン管理]のAmazonへのリンク

Git は慣れてしまうと使い勝手がいいが慣れるまでが大変なので、この本が大いに参考になりそうだ。

今日の作業内容:Debian インストールの調整 1h + 論文読み 1h
今日のランチ:シッダルータ チキンカレー
明日の予測作業時間:3h

2012年2月1日水曜日

Debian testing のインストールを確認

今日は論文読みと並行して Debian testing のインストールを確認してみた。
基本的には
http://cdimage.debian.org/cdimage/daily-builds/daily/arch-latest/i386/iso-cd/debian-testing-i386-businesscard.iso
をダウンロードして VirtualBox にインストールしている。

すでに日本語でインストールできる状態になっており、kernel も 3.1 を選択できるようになっている。
パッケージの選択の時に何も選ばないとブラウザも入っていない最小構成を簡単に作ることができる。

今回はできるだけ小さい構成ということで lxde をインストールして、あとは ibus-anthy, TeX 関係、そして SDPA 関係をインストールしている。
このくらいでざっと 1000 個のパッケージとなり、仮想ディスクが 4GB ほど消費される。

VirtualBox のシームレスモードは VMware のユニティに相当するが、自分のところではユニティよりも高速に動くようだ。
シームレスモードで利用していれば、ブラウザはホストOSのものが使えるので、ブラウザをインストールする必要もない。

ところで、TeX に日本語が通せるか確認してみたが、これは今までの知識で対応できる範囲だった。
ただ、yatex だけが emacs で起動できず、debian の bug として登録しておいた。


今日の作業内容:Debian インストール 5h
今日のランチ:つかさ 鯛のから揚げ
明日の予測作業時間:3h