2015年3月31日火曜日

Julia で sparse Cholesky 分解

Julia の場合、以下のような感じで sparse Cholesky factorization ができる。

julia> A = sprand(100,100,0.1);
julia> B = A*A' + 100*speye(100);
# ここまでで B に sparse で positive definite な行列を生成

julia> C = cholfact(B);
# これで CHOLMOD.Factor が生成される

julia> L = sparse(C);
# これで普通の sparse 行列として L を取り出せる。
# ただし、AMD による行と列の順番変更がある。
# この情報は C.Perm に入っているが、CHOLMOD の情報のままのためインデックスが 0 スタートであり、Julia の1スタートのために C.Perm + 1 など1を足す必要あり。

julia> vecnorm(L*L' - B[C.Perm+1, C.Perm+1])
1.428931331049563e-13
# これで数値誤差が評価できて、数値誤差を除けば正しく計算できていることが分かる。


2015年3月24日火曜日

ダメなときは、ダメなことが続くものだ

最近は、いろいろとダメなことが続いているので、
「あれもうまくいかず、これもうまくいかず」
という状態になっている。

まぁ、いろいろとやっていると、ダメなことが続く時期というのはやっぱりあるものなので、あまりジタバタとせずに我慢どころでもある。

とりあえずは身近にあるところからコツコツと進めることにして、
「失敗リスト」
と銘打ったリストを作ることで、どんな失敗が続いているのか、ということを見てみることにする。


2015年3月20日金曜日

SDPA-GMP を 7.1.3 に更新

最近の SDPA-GMP に関しての問い合わせがあった関係で、実際に SDPA-GMP をコンパイルしてみたところ、いくつかエラーが出てコンパイルが上手くできなかったので、それを修正したものをアップロードしておいた。

論文の場合は一度採択されてしまえば、それで終わりであって、メンテナンスされることはごく限られた論文以外ほとんどないが、ソフトウェアは再利用価値が高いので、メンテナンスなども必要である。

SDPA-GMP は他に類似ソフトウェアがないので、まだまだ使われていくことになりそうだ。

2015年3月10日火曜日

CVXOPT, CHOMPACK の Windows 64bit バイナリの生成方法

Python のライブラリである CVXOPT, CHOMPACK のビルドの方法は、手順さえ分かってしまえば簡単に行えるので、手順を載せておく。
ここでは、CHOMPACK をインストールする。なお、CHOMPACK は CVXOPT を使っているので、CVXOPT をあらかじめインストールしておく必要がある。(CVXOPT の簡単なインストール方法は下のほうに書いておく。)

1. WinPython をダウンロードして展開する
これで、python 関係だけでなく、mingw32 などコンパイルに必要なものも一通りインストールされる
2. https://github.com/cvxopt/chompack/archive/master.zip
からソースをダウンロードして zip ファイルを展開すると、chompack-master というフォルダができる
3. http://sourceforge.net/projects/openblas/files/ の最新バージョンの
OpenBLAS-v(バージョン名)-Win64-int32.zip をダウンロードして zip ファイルを展開すると、lib フォルダの中に libopenblas.a があるので、これを chompack-master フォルダにコピーする
4. chompack-master フォルダにある setup.py を編集して
BLAS_LIB_DIR = ['./']
BLAS_LIB = ['openblas']
LAPACK_LIB = ['openblas']
とする
5. WinPython のフォルダにある「WinPython Command Prompt.exe」 を実行してコマンドプロンプトを開く
このコマンドプロンプトでは、mingw32 などのパスが自動的に設定されている
6. コマンドプロンプトのcd コマンドで chompack-master に移動する
7. ビルドする
 python setup.py build --compiler=mingw32
8. ビルドされたバイナリファイルをwheel ファイルにまとめる
これにより、 pip での install, uninstall が簡単に行える
 python setup.py bdist_wheel
9. dist/ フォルダに出来上がるので、これを pip でインストールする
cd dist
pip install chompack-2.2.1-cp34-none-win_amd64.whl
これでインストール終了である。

pip でインストールしてあるので、アンインストールするときも
pip unintall chompack
でアンインストールできる。

なお、CVXOPT については、WinPython のものだと BLAS との連携の関係でうまく実行できないときがあるので(OpenBLAS でない Reference BLAS を使うと、うまく実行できないケースがある)、
pip unintall cvxopt
で一度削除したうえで
http://www.lfd.uci.edu/~gohlke/pythonlibs/
から
cvxopt‑1.1.7+openblas‑cp34‑none‑win_amd64.whl
をダウンロードして
pip install cvxopt‑1.1.7+openblas‑cp34‑none‑win_amd64.whl
としてインストールする。
(pip install cvxopt だけだと pip の標準のレポジトリからインストールされてしまうので、ダウンロードしてファイル名で指定する。)


2015年3月5日木曜日

CHOMPACK を Python3 で試してみる

sparse Cholesky に関係した CHOLMOD を Python から使えるようにしている CHOMPACK を使ってみているが、どうやら CHOMPACK は Python 2 用に書かれていて、そのままではPython 3 では動かない。
いくつか修正すべきところをリストアップしたので、作者に報告をしておいた。そのうちに反映されるかもしれない。

今回の CHOMPACK は mingw32 を使って Windows 用にコンパイルをして Windows 上で確認をしているが、OpenBLAS などもバイナリをダウンロードしてくれば自分でコンパイルする必要もなくバイナリを作れるし、wheel 形式にもすぐできるので pip によるパッケージ管理にも対応しており、非常に簡単になっている。

CHOMPACK は以前にチェックしたときよりも機能が豊富になっており、たとえば chordal graph による変換などを一括して行うプログラムなども実装されている。



2015年3月4日水曜日

Julia と Python で比較

Python は数値計算などでメジャーな言語になりつつあるが、Operations Research の研究者が作った数値計算用言語の Julia とで、sensor network localization (SNL) の SDP 緩和のプログラムを作ることで比較してみた。

結果:
Python --> 動作するプログラムを作れる
Julia --> あまりにも面倒で断念

Python がベターなところ:

  1. IDE の Spyder を使うと、その中でだいたいの作業ができる
  2. matplotlib があるので簡単に結果のグラフを作れる
  3. Numpy, Scipy, CVXOPT のデータ構造は、錐計画問題に使うにはそれぞれが不完全であるが、データ変換で連携させれば錐計画問題については一通りのことができる
  4. プロファイリングも簡単
  5. ドキュメントが整備されていて、Spyder の object inspector で簡単に見ることができる


Julia がベターなところ:

  1. よく使う行列演算がMatlabと同じ感じで使える (inv などが標準で使用可能)
  2. 行列の操作も Matlab に似ている ( A = [1 2; 3 4] で 2x2 の行列を作れる)
  3. パッケージ管理で、 Pkg.add("パッケージ名") と簡単に追加できる (git を内部的に利用している。)
ちなみに、Julia を今回断念した理由としては、

  1. グラフ描画の Gadfly が使いづらい。
    plot を複数回に分けて書けないので、グラフに必要なデータをすべて保持してから、一度の plot でグラフを完成させないとならない。
    また、サンプルにあるものは実行できるが、ドキュメントが足りないので、サンプル以外の内容を実行しようとするとソースを見ないと分からない。
    Python の matplotlib を呼び出すこともできるが、だったら Python のほうが便利。
  2. cd, pwd, ls (または dir) などでディレクトリの移動、ファイル名を取得するのが面倒。
  3. 標準的なIDE とされている Juno が使いづらい。
    Vista にも対応、とホームページには書いてあるが、実際には Vista では使用不可だったりする。
    また、Ctrl+Return でエディタに書いてある式を評価してくれるが、結果がエディタ部分に表示されたり、コンソール部分に表示されたり、両方確認する必要がある。むしろ、コンソール部分だけに集約してくれた方が表示されているかどうかを確認しなくて済む分、分かりやすい。
  4. ドキュメントが足りないので、ほとんどのことを Google などで検索する必要がある。
  5. プロファイリングでは時間での計測ではなく、関数呼び出し回数の計測。
といったあたり。
やはり、 Julia と比較するとmatlplotlib, Spyder, ipython が Python を使ううえで非常に強力で、このあたりを Julia に移植と Julia の使いやすさがアップするのでは?と考えてもいる。
今の Julia は Python のライブラリがないとうまく処理できない部分が多いので、Julia を使うためには Python も覚えたほうが効果的である。


2015年3月2日月曜日

論文チェック: From Predictive to Prescriptive Analytics

今朝のメールで来ていた Optimization Online への投稿の中から

From Predictive to Prescriptive Analytics
http://www.optimization-online.org/DB_HTML/2015/02/4773.html

についてチェックをしておいた。
「Machine Learning や OR/MS の知識で、prescribe な最適決定をする」ということが目的なようである。

基本的なフレームワークは、stochastic optimization にある Sample Average Approximation に基づいている。
これらをベースにして、よりいろんなことができるように拡張を考えている。
ちなみに、 Robust 最適化を取り込んだものが
http://www.nathankallus.com/DataDrivenRobOptv1.pdf
http://arxiv.org/pdf/1408.4445v1.pdf
で見ることができる。

この論文では、データの解析手法よりも、どこからデータを集めてきているか、というところが面白くて、IMDb などの映画に関するデータなども数値計算に取り込んでいる。