2015年12月17日木曜日

疫病の広がり方と研究資産

金曜日に参加した制御関係のワークショップで、疫病の広がり方についてエージェントベースのシミュレーションを用いている研究があり、非常に勉強になった。

ところで、この話を聞いて思い出したが、疫病の広がりを抑えるために、人的ネットワークのどこに医療予算をどのように配分するべきか、ということを半正定値計画問題で計算している、という論文があった。
この結果をシミュレーションのパラメータなどに利用したら、有効活用できるかと思われる。


ところで、疫病の場合は広がりを抑えることが計算目的となるが、同じような計算ステップを用いれば、研究分野でどの研究を普及させるためには誰に研究予算をつけるべきか、ということも計算できる。
たとえば、ドイツの「インダストリー4.0」などはキーワードとしてよく聞くようになっているが、研究をブランドとして確立していくためにも数理最適化のアプローチを使えるかと考えられる。
ネットワークデータをどのように入手すべきか、という課題があるが、計算手順としては開発可能なアプローチかと考えられる。

2015年12月7日月曜日

Julia の git 版をコンパイルしてみる

Debian で Julia を apt-get でインストールすると 0.3.2 というバージョンがインストールされ、少し古いバージョンである。

ということで、最新版を素っ飛ばして、git 版をインストールしてみる。
Julia はバージョンごとにファイルが整理されているので、複数のバージョンが共存しやすい。

まずは、
$ git clone https://github.com/JuliaLang/julia.git
としてソースを clone する。
そのあとは、

$ cd julia
$ make

として、数時間放っておくだけでコンパイルが終わる。
(apt-get で Julia をインストールしてあるので、コンパイルに必要なファイルなどは既にそろっているのではないかと推測している。)

実行するには
$ ./julia
とすればよく、Pkg.add などでパッケージもインストールできる。
(ただし、Version 0.4 とは別にインストールされるので、パッケージを一通りインストールしなおす必要がある。)

ちなみに、Version を表示してみると、
$ ./julia --version
julia version 0.5.0-dev
と表示される。

今日の作業内容:ワークショップ準備 2h + 樹木園計算 2h
明日の予測作業時間:3h

2015年12月2日水曜日

離散凸を勉強中

手元で最適化しようとしている関数が離散凸関数であるかどうかを判断するために、離散凸について、もう少し勉強中。

連続関数なら、ヘッセ行列が半正定値なら凸関数になることが分かるけど、離散凸では様子が違う、ということが下の資料に書いてある。
http://www3.grips.ac.jp/~tsuchiya/PSD/slides/murota140114.pdf

つまり、離散凸の場合の実行可能解である各点でヘッセ行列が半正定値でも離散凸関数とは限らないようだ。

このあたりをふまえて、もう少し調べてみることにする。

今日の作業内容:ワークショップ準備 1h, 離散凸を下調べ 2h
明日の予測作業時間:1h

2015年11月30日月曜日

離散凸解析、面白い

最近の計算していた内容に離散凸解析が使えるかも、ということで、教えてもらいながら勉強してみている。
同じ凸でありながら、連続の凸関数と離散凸関数では、それなりに様子が違っていて新鮮である。

たとえば、離散凸関数は、うまく拡張すると連続凸関数にできることは知っていたが、逆に連続凸関数だからと言って定義域を整数点だけにしても離散凸関数になっているとは限らないようだ。
(もう、てっきりそうなっているものだと思っていたが。)

まずは数値計算をしてみて、今の計算したい内容に使えるかどうかを見てみて、使えそうだったら、今の計算している関数が離散凸関数ってことだと判断することにしてみている。

今日の作業内容:ワークショップ準備1h+離散凸2h
明日の予測作業時間:2h

2015年11月28日土曜日

SOCP緩和とSDP緩和の間が必要そう

今日の数値実験だと、
SDP 緩和だと誤差が1%ぐらいだが、計算時間は10分。
SOCP 緩和だと誤差が4%ぐらいだが、計算時間は5秒。
となっていて、やはり SDP 緩和の誤差の小ささが際立つが、もう少し計算時間を短縮したいところ。

この間のあたりで、うまく使えそうなものがあるのかどうか、どうなのだろうか?

今日の作業内容:数値実験をこまごまと 2h
明日の予測作業時間:2h




2015年11月25日水曜日

Matlab での計算時間比較

いま、n 次元行列として A があって、V を 1 から n の添え字の部分集合としたとき、V の要素を行と列にもつ行列を A から作るには

A(V,V)

とすればよい。

ところで、この対角成分を取り出すとした時に、
v = diag(A(V,V))
とするよりも
diagA = diag(A); v = diagA(V);
とした方が圧倒的に高速である。
(このあたり、自動的に Matlab が変換してくれるかと思っていたが、どうやらそうでもないらしい)



今日の作業では、こういった高速化をいくつか入れることで、昨日1時間かかっていた部分を3秒程度まで短縮することができている。

今日の作業内容:プログラムの高速化 3h
明日の予測作業時間:2h




2015年11月24日火曜日

今日から従来モードに戻ります

今日は、数値実験のやり直しを進めておいた。
既存手法と、いま実装してみているヒューリスティクスを比較しやすいように、既存手法の結果を表としてまとめる作業を行っている。
いくつかのソースに分かれているので、複数ステップ必要となっているのが分かりにくいが、コツコツと進める感じになっている。

今日の作業内容:数値実験やり直し 2h
明日の予測作業時間:3h

2015年11月20日金曜日

樹木園計算、面白い

前から計算している樹木園計算が、なかなかに面白い状況となってきている。
どういう状況かというと、

1.既存手法を論文のアルゴリズム通りに自分で実装すると、計算が途中で破綻する。
2.しかしながら、既存手法を実装したソフトウェアがあって、それだと計算もできて、それなりの解が出る。(このソフトウェアは、ソースを公開していない。)
3.その一方で、自分で実装したものでも、データによっては計算が途中で破綻しないときもあり、その場合には既存ソフトウェアよりも30%程度高速に似たような解が作れる。


これから推測するに、「論文には載っていないヒューリスティクスを使っている」と考えられる。

ただ、このヒューリスティクスがなんなのか不明。
不明ではあるけど、うまく実行できている。
このヒューリスティクス、他のところにも使えるのではないか、と考えると、どんなふうに行っているのか、より一層気になるところである。

なかなかに面白い。


2015年11月4日水曜日

SparsePOP を Windows でコンパイルした時のエラー回避

SparsePOP を Windows 上の Matlab でコンパイルしようとすると、エラーが出てコンパイルができないことがある。

この時には、SparsePOP にある global.h で __func__ を __FUNCTION__ に置換するとコンパイルが通るようになる。


2015年11月2日月曜日

整数制約付き線形計画問題は、線形計画問題に定式化できる

いわゆる Mixed Integer Programming Problem、より正確には整数制約付きの線形計画問題は、線形計画問題に定式化できる。
したがって、指数オーダーで増加するの計算時間と計算メモリに耐えられれば、理論上はシンプレックス法で解くことができる。

定式化するには、
(1) copositive programming により定式化する
(2) copositive programming を sum of squares により定式化する
(3) sum of squares を非常に巨大な LP により定式化する
とすればよい。

2015年10月23日金曜日

RAMP のセッションオーガナイザーを終えて

先週の浜松で行われた RAMP で、無事にセッションオーガナイザーのお仕事が終わりました。

セッションオーガナイザーのお仕事は、その名前が示す通り、セッションをオーガナイズすることですが、そのお仕事を通じて、いろんな方のお世話になりました。

まずは、実行委員長とプログラム委員長のおかげで、セッションをオーガナイズすることに集中できました。

また、講演者の方には、準備段階から原稿の締め切りなど、ほぼきちんと守っていただいて、とても助かりました。おかげさまで、シンポジウムまでに原稿を20回ぐらい読み込むことができました。もちろん、シンポジウムでの発表も、「バリバリの非線形最適化」にふさわしい、どれも面白いものばかりで、勉強になることばかりでした。

帰りの新幹線では、サンテグジュペリの言葉を思い出したりしました。

自分がオーガナイザーをできるのは今回1回切りだと思いますが(いろんな方がオーガナイズすると発表の幅もできると思うので)、もう1回機会があったら是非再度やってみたい、と強く思える貴重な機会でした。


2015年9月29日火曜日

簡単そうに見えて、実は難しい、最適化問題

最小自乗問題
$$ \min \frac{1}{2} || Ax - b ||^2 \quad \mbox{s.t.} \quad x \in R^n $$
は、A の行が一次独立なら簡単に解ける問題。

次に、
$$ \min \frac{1}{2} || x - x^0 ||^2 \quad \mbox{s.t.} \quad x \in R^n,\ x \ge 0$$
と各変数が非負に制限されるときも簡単に解けて、解としては
$$ x^* = \max\{x,0\} $$
である。

ところが、この2つが組み合わさると
$$ \min \frac{1}{2} || Ax - b ||^2 \quad \mbox{s.t.} \quad x \in R^n,\ x \ge 0$$
とたんに、簡単には解けなくなる。


のように内点法ベースなら解けるのであるが、もっと簡単に解けそうにも見えて、そういった解法が見当たらない。
この「簡単そうなのに、より適した解法が案外見当たらない」というところが、面白い。

2015年9月28日月曜日

LPの今までとは違う解法

最近の論文をチェックしていたところ、LP のシンプレックス法や内点法とは異なる解法についての論文があった。
ぱらぱらっとめくったところでは、最適解が得られることなどは証明されている様子だが、どれだけの反復回数で収束するのか、のあたりはまだ検討されていない様子。ただ、「反復回数」の扱いも内点法と同じような基準でいいのかどうかもよく分からないので、単純な比較はできないかも知れない。

流れとしては、シンプレックス法よりも内点法に近い流れなので、中心パスなどとどう関係しているのか、などを考えてみると面白いかと思う。


2015年9月24日木曜日

unit simplex への射影

unit simplex、つまり
$$ \sum_{i=1}^n x_i= 1, x \ge 0 $$
への射影の方法が以下の論文に掲載されている。

A projected gradient method for optimization over density matrices
D.S. Gonçalves, M.A. Gomes-Ruggiero & C. Lavor
http://www.tandfonline.com/doi/abs/10.1080/10556788.2015.1082105?journalCode=goms20

Dykstra などの交互射影法とは異なり、有限回の計算で解が計算される点が面白い。

ところで、この方法を変形して、
$$ \sum_{i=1}^n x_i= 1 $$
にまずは射影してから計算すると効率が上がったりするのだろうか?
そのあたりを考えてみるのも面白そうではある。


2015年9月14日月曜日

著者が亡くなっても掲載される論文

Mathematical Programming Computation の Volume 7, Number 3 の Preface にいろいろと事情が掲載されている。

非常に意義深いことだな、と思う。

2015年9月1日火曜日

Matlab ソースを PDF に一括変換する

Matlab のソースを Matlab がインストールされていないパソコンで閲覧したいときには一度 PDF に変換したりすると便利だが、 Matlab の機能の publish だと HTML に変換して PDF に変換するのを一つ一つ行う必要があって、これがなかなかに面倒。

ということで、以下のようなシェルスクリプトで一括変換する。

for i in `find . | grep m$`
do
    LANG=C a2ps -1 --pro=color $i -o - | ps2pdf - $i.pdf
done

ここで、 LANG=C を入れているのは、a2ps で日本語が文字化けするので日付などが日本語だと文字化けしてしまうため。
また、--pro=color とすると、出力がカラーになるので、見やすくなる。

2015年8月26日水曜日

対角優位行列とゲルシュゴリンの定理

「対角優位行列は半正定値行列である」

という命題が成り立つ方法を考えていたが、ゲルシュゴリンの定理を使ったら、あっというまに示せることがわかった。

やっぱり、いろんな定理を知っている、ということも大事なものである。

2015年8月11日火曜日

Julia で多倍長計算

Julia でも多倍長計算できる、ということでチェックしてみた。

基本的には BigFloat を使えば簡単に行うことができて、with_bigfloat_precision を使えば特定の範囲だけ精度を変更できる。

with_bigfloat_precision(30) do
    x = BigFloat(0.1)
    i = 0
    while true
        i += 1
        x = x/2
        y = e^(-x^2)
        print(i)
        print(",")
        print(x)
        print(",")
        println(y)
        if y==1
            break
        end
    end

end

このようにすると、30ビット分の精度で計算できるし、30を1000にすれば、もちろんもっと精度が出せる。

ちなみに、

X = [BigFloat(rand()) for i=1:5, j=1:5]

とすれば、5x5 の行列なども作れる。
足し算や inv など基本的な操作はすでに行えるが、eig や chol のあたりはまだ実行できない様子である。

行列周りが完全ではないけど、いろいろとできる範囲はそれなりに広くて、なかなかに面白そうな機能だ。


2015年7月31日金曜日

数値実験の続き

CPLEX で混合整数計画問題を解き始めて、はや1週間。
結局、CPUを 1500% 使用し続けたが上界と下界のギャップが8%程度までしか接近せず、n = 1000 の問題は断念。

もっと小さい問題、ということで n = 200 にして再度挑戦中。10% 程度までは比較的簡単だが、やはり本当の勝負はそこからの様子。

2015年7月24日金曜日

cplex が思っていた以上に時間がかかる

SOCP 付きの整数計画問題で厳密解がどうなっているかを知るために、あまり深く考えずに cplex にお任せしたら、想定以上に時間がかかっていることに気がついた。

変数の数は、0-1 変数が 1000, 連続変数も 1000 で、0-1 変数を 0<=x<=1 に緩和した SOCP は 0.1 秒もかからずに解けるので、それほど大きい問題ではないかと思い、そのまま cplex に丸投げしてある。

今のところ CPU 使用率が 1500% 程度を 24 時間ぐらいやってみて、上界と下界の差が 0.2% から縮まらない状態である。

分枝限定法での限定が上手くできないタイプの問題なのかもしれないが、もう少し定式化もいじってみようと思う。




2015年7月23日木曜日

Ipopt の Matlab インターフェースのコンパイルで No such file or directly と言われる

Ipopt を試すためにコンパイルしてみたが Matlab インターフェイスでコンパイルがうまく行かなかった。

より具体的には、
CoinIpopt/Ipopt/contrib/MatlabInterface/src
のディレクトリで
$ make
として ipopt.mexa64 を生成しようとすると、以下のエラーで停止する。

g++: error: matlabexception.o: No such file or directory
g++: error: matlabfunctionhandle.o: No such file or directory
g++: error: matlabjournal.o: No such file or directory
g++: error: iterate.o: No such file or directory
g++: error: ipoptoptions.o: No such file or directory
g++: error: options.o: No such file or directory
g++: error: sparsematrix.o: No such file or directory
g++: error: callbackfunctions.o: No such file or directory
g++: error: matlabinfo.o: No such file or directory
g++: error: matlabprogram.o: No such file or directory
g++: error: ipopt.o: No such file or directory

make のときのエラーを確認していったところ、.libs サブディレクトリにこれらのファイルがあることがわかったので、以下のように対応してみた。

$ cd CoinIpopt/Ipopt/contrib/MatlabInterface/src
$ cp .libs/*.o .
$ make

これで無事に ipopt.mexa64 が生成できる。

2015年7月3日金曜日

IguanaTeX で日本語が通るようになった!

http://karat5i.blogspot.jp/2015/04/iguanatex-powerpoint-tex.html
に書いてある通り、 IguanaTeX で日本語が通るようになった。
自分のところでも試してみたが、正常に動作している。

今までは日本語をうまく使えていなかったが、これで platex でコンパイルできるものをPowerPoint にも使えるようになり、とても便利になった。

素晴らしい。

学術雑誌で投稿ができない状態になっているものがある

数理最適化関係の学術雑誌で、投稿数が増えすぎて新規投稿を reject している雑誌がある。
別の雑誌に投稿するように、という意味での reject になっている。
ただ、ホームページにそういったことが書かれてはいないので、投稿するまではその雑誌がそういう状態になっていることを把握するのは難しそうだと思う。

それにしても、査読などをたまに引き受けたりもするが、やっぱり査読者を見つけてくるという Editor や Associate Editor のお仕事は相当大変なものだろうから、そのあたりはほんと感謝だな、と思ったりする。


2015年6月26日金曜日

Julia でループ展開しても、あまり速くない?

疎行列と密ベクトルの積を計算するような以下のコードを書いてみた。

==== Ax.jl =====

function Ax(n, density, seed)
    srand(seed)
    println("Generating A")
    time_A = @elapsed @time A = sprand(n,n,density)
    x = rand(n,1)
    println("Computing b = A*x")
    time_b = @elapsed @time b = A*x

    println("Computing b2 = A*x")
    time_b2 = @elapsed @time begin
        b2 = zeros(n,1)

        for j=1:n
            const xj = x[j]
            for index_i = A.colptr[j]:(A.colptr[j+1]-1)
                b2[A.rowval[index_i]] += A.nzval[index_i]*xj
            end
        end
    end

    println("norm(b-b2) = ", norm(b-b2))
    println("time (A) = ", time_A, " time (b) = ", time_b, " time(b2) = ", time_b2)

    return A, x, b, b2

end


==============

実行結果は以下の通り。

julia> include("Ax.jl"); (A, x, b,b2) = Ax(2000,0.8,1024); (A,x,b,b2) = Ax(2000,0.8,1024);
Generating A
elapsed time: 0.173287485 seconds (76853144 bytes allocated, 24.33% gc time)
Computing b = A*x
elapsed time: 0.011240257 seconds (16064 bytes allocated)
Computing b2 = A*x
elapsed time: 2.677362552 seconds (498951072 bytes allocated, 9.45% gc time)
norm(b-b2) = 0.0
time (A) = 0.173453157 time (b) = 0.011356873 time(b2) = 2.677506154
Generating A
elapsed time: 0.141352727 seconds (76853144 bytes allocated, 12.73% gc time)
Computing b = A*x
elapsed time: 0.010914271 seconds (16064 bytes allocated)
Computing b2 = A*x
elapsed time: 2.708701738 seconds (498951072 bytes allocated, 9.03% gc time)
norm(b-b2) = 0.0
time (A) = 0.141499582 time (b) = 0.011023018 time(b2) = 2.708858507


b の計算はもともとあるもので、b2 のほうがループを展開して計算している。
ちなみに、2回実行しているのは、Julia のコンパイルで時間がかかると遅くなるかもしれない、ということで2回実行してみた。

うーむ、やはり numpy のように、できるだけ行列やベクトルの演算に持ち込んだ方が
速いのだろうか?

そうそう、@elapsed, @time マクロは begin - end で挟むと、その間の時間を測定してくれる。
@elapsed マクロのほうは返り値で経過時間を返してくれるのでそれを代入しているが、@time マクロのほうが細かい情報が表示される。

2015年6月25日木曜日

Julia で簡単な Newton 法の実装

関数 f_1(x) = x[1]^4 - 2x[1]x[2]+ 2x[2]^2 - 4x[1] + 2x[3] についての Newton 法を実装してみると以下のようになり短いコードで書ける。

面白い点としては、

  1. C や Matlab なら 2*x[1]*x[2] と書かなければいけないコードが 2x[1]x[2] と * を省略して書ける
  2. UTF-8 でファイルを作っていれば、∇も普通の文字として扱える (全角空白が紛れ込むと、かなり厄介だけど)

実行するには
julia> include("newton.jl")
とする。
ちなみに、このコードなら短いので 15 分程度で書ける。

==== newton.jl ====

# -*- coding:utf-8 mode:Julia -*-

function f_1(x)
    return x[1]^4 - 2x[1]x[2]+ 2x[2]^2 - 4x[1] + 2x[3]
end
function ∇f_1(x)
  return [4x[1]^3-2x[2]-4; -2x[1]+4x[2]+2]
#    return [4*x[1]^3-2*x[2]-4; -2*x[1]+4*x[2]+2]
end
function ∇2f_1(x)
    return [12x[1]^2 -2; -2 4]
end


ε = 1.0e-8
k = 0
x = [1; 1]
norm_∇f = 0.0
while k < 100
    ∇f = ∇f_1(x)
    norm_∇f = norm(∇f)
    if norm_∇f < ε
        println("Convergence!!")
        break
    end
    ∇2f = ∇2f_1(x)
    @printf "k = %d : x = [%.2e, %.2e], norm(∇f) = %.2e\n" k x[1] x[2] norm_∇f
    x = x - ∇2f \ ∇f
    k = k+1
end
@printf "k = %d : x = [%.2e, %.2e], norm(∇f) = %.2e\n" k x[1] x[2] norm_∇f

==========






2015年6月23日火曜日

内点法で新しい結果

arXiv をチェックしていたら、内点法の新しい結果が出ていて、なんと infeasible interior-point method のオーダーが低くなっている。

内点法や線形計画問題の理論的な性質は 2000 年頃までの研究でほとんどなされており、それ以降はソフトウェア開発などが中心で理論的な前進はあまりない状況だった。
理論的な論文に関しても、「今までと同じオーダーの解析ができました」が主流であったが、今回の論文は違うようだ。

非常に面白そうなので、あとで読んでみることにしよう。




2015年6月20日土曜日

Julia で最小自乗法を組んでみる

簡単な最小自乗法なら、以下のような感じでプログラムを書ける。
Matlab にそっくりな感じ。

==== least1.jl =====
function leastsquare(m, n, seed)

    if m < n
        error("Too small m")
    end
    srand(seed)
    A = rand(m,n)
    b = rand(m)

    xhat = (A'*A) \ (A'*b)
    e = A*xhat - b
    norm_e = norm(e)
    @show A
    @show b
    @show e
    @printf "norm_e = %.3e\n" norm_e

end
===============

実行するには、
julia> include "least1.jl"
julia > leastsquare(5, 2, 1024)
といった感じ。




2015年6月8日月曜日

なんだか良く分からないバグが取れた

C++で組んでいたときに、

if (data[419] ! = data2[419]) {
   printf("diff = %d\n", data[419]-data2[419]);
}

としたら、なぜか diff = 0 と表示される、というバグが発生。
なぜなのか良く分からないが、data, data2 にデータを入れる前に全部の要素を0で初期化するとバグが取れた。


おそらくメモリリークしているのかと推測されるが、他のところがあまりに単純なコード過ぎてメモリリークを発見できず。
しかも、Linux だと発生せずに Windows 上だと発生する、というたちの悪いバグになっており、valgrind などでもうまく特定できず。

なんかよく分からないけど、とりあえずバグが取れたのでいいにしておこうかと思う。


2015年6月4日木曜日

TeX でå (a の上に丸がついたもの)を使うマクロ

TeX で a の上に丸がついた å を使うとき、本来であれば、
\r{a}
とすればOKだが、\r を別のマクロで使っていたりすると、これが困る。

この回避方法は、automake の texinfo.tex を見ていて分かったが、
\def\ringaccent#1{{\accent23 #1}}
とマクロを定義してから
\ringaccent{a}
とすると回避できる。


ドイツ語のウムラウトなども\accent24 などで定義されていたりするので、これをうまく使うと \o などを別のマクロに使ったりもできて便利でもある。

2015年6月2日火曜日

Julia で Rosenbrock 関数の最急降下法

練習用に実装をしてみた。
以下のファイルを rosen.jl として保存して、
include("rosen.jl"); test01(3);
とすれば、Rosenbrock の n=3 のときを最小化してくれる。
(ただし、4154 反復かかる。)

やはり、Julia で便利なのは
1. リスト内包ができる
2. UTF-8 の文字が変数名に利用できる(∇も利用できる)
3. @show マクロが使える
のあたりかと思う。
Rosenbrock の関数式を調べるところからはじめて、この文章を書くまでで、だいたい1時間といったところ。

そういえば、C 言語なら x == 0 || y == 0 のようなところの || の前に改行が入ってもOKだが、Julia の場合は NG になる様子。

== ここから rosen.jl ==

# -*- coding:utf-8 mode:Julia -*-

function Rosen(x)
    sum([ (1-x[i])^2 + 100*(x[i+1] - x[i]^2)^2 for i = 1:length(x)-1])
end

function ∇Rosen(x)
    N = length(x)
    [[-2*(1-x[i])+2*100*(x[i+1]-x[i]^2)*(-2*x[i]) for i=1:length(x)-1]; 0] + [0; 2*100*[x[i+1] - x[i]^2 for i=1:length(x)-1]]
end

function minimizer(n, f, ∇f)
    count_f   = 0
    count_∇f = 0
    η = 0.05

    x = zeros(n,1)
    current_f = f(x)
    count_f += 1
    max_iter = 10000
    iter = 1
    while iter <= max_iter
        gf = ∇f(x)
        norm_gf = norm(gf)
        count_∇f += 1
        α = 1.0
        next_f = 0.0
        next_x = 0
        while true
            next_x = x + α*(-gf/norm_gf)
            next_f = f(next_x)
            count_f += 1
            if next_f <= current_f - η*α*norm_gf || α < 1.0e-10
                break
            end
            α *= 0.5
        end
        @printf("iter: %d, f: %.4e, gf = %.2e, α = %.2e\n", iter, current_f, n\
orm_gf, α)
        if norm_gf < 1.0e-4
            println("Norm of gradient is small")
            break
        end
        if α < 1.0e-10
            println("Step length is too short")
            break
        end
        if abs(current_f - next_f)/max(abs(current_f),1.0) < 1.0e-10
            println("No improvement of f")
            break
        end
        current_f = next_f
        x = next_x
        iter += 1
    end
    x
end

function test01(n)
    minimizer(n, Rosen, ∇Rosen)
end



2015年5月22日金曜日

Julia の発展途上なところ

この前の簡単な内点法などをしてみて、以下のところが数理最適化に使うには Julia はまだまだ発展途上な印象を受ける。

(1) 遅い
C, C++ などと比較しても遜色ない速さ、とうたっているが、これはある意味で当然で、そもそも実行前にコンパイルしている(インタープリターではない)。

また、実際にプログラミングをしてみると
using PyPlot
など頻出する using の遅さが地味に結構痛い。

ちなみに、これについては、PyPlot などを ~/.juliarc でロードしておけば、それなりに快適にできる。Matlab が起動に時間がかかっている作業を Julia でも行っているようなもので、起動に時間をかけることで、そのあとの作業を効率化できる。


(2) Profiling がカウンターであって秒で計測できない。
たとえば 1/1000 秒ごと、と指定してもサンプルカウントをしているので、かなり荒い情報しか得られない。だいたいのボトルネックを把握するには有効であるが、大まかなチューニングしか行えない。
しかも @time などの時間測定のコマンドを組み入れた関数を @profile で測定すると、サンプルカウントがうまくいかないことがある。

(3) IDE の機能がいまひとつ。
Juno を使ってみたが、数理最適化に多い反復計算には向いていなかった。現状では、emacs とコマンドラインが、一番便利。
ただし、Julia の REPL は結構強力であるので、コマンドラインからの利用でもそれなりに行けたりする。



ここに挙げたうちで一番重要そうなのは (3) で、たぶん IDE に Python の spyder なみの環境ができれば Julia の利用可能性はぐっと高まるかと思う。



2015年5月21日木曜日

TeXインストーラ3での「mktexlsr の実行に失敗しました」を修正する

TeX インストーラ3で TeX を一括で更新した時に、なぜか
\usepackage{color}
などで
! LaTeX Error: File `color.sty' not found.
と怒られるときがある。
これは、TeX インストーラ3 の更新の時に mktexlsr の実行に失敗しているためのようで、abtexinst_log.txt にもエラーメッセージとして記録されている。

これを修正するには、コマンドプロンプトを「管理者権限」で開いて
$ mktexlsr
と実行する。
(「管理者権限」でなくても mktexlsr は実行できて成功するが、ファイルには保存されないので実質的に「管理者権限」が必要である。)


2015年5月14日木曜日

実用的な最適化とは?

シミュレーテッドアニーリングで max cut を解く、という論文が arXiv の

Solving maximum cut problems by simulated annealing

で読むことができる。
シンプルな実装(ソースコードも途中に書いてあるぐらい)になっているが、数値結果を見てみると、それなりに良い結果になっているようにも見える。


最適化の論文には難しい理論が書かれているマニアックなものが多くて、そういうのが多いのも面白いが、こういう simple is best みたいなことも研究の方向としてはあるべき一つなのかと思ったりもする。
この論文がどこの雑誌に採択されるか、気になるところだ。

2015年5月13日水曜日

Julia で戸惑うところ (Julia と Matlab の違い)

Julia は Matlab に似たように書けるので、ついつい同じように書いてハマることがあるので、そういったことについてのメモをしておく。

(1) 変数の deepcopy に気をつけろ!
たとえば、Matlab で
x = [3 4 5];
y = x;
y(2) = 7;
としても、x は変わらず x = [3 4 5];である。
しかし、Julia で
x = [3 4 5];
y = x;
y[2] = 7;
とすると、y につられて x も x = [3 7 5] に変化する。
要するに、Matlab では y = x で別の配列が準備されて配列の中身もコピーされるが、Julia の y = x は単純に「 x という配列に y という名前でもアクセスできる」ということで中身も同じである。
Matlab のようにするには、
y = deepcopy(x)
とする。


(2) integer と float に気をつけろ!
Matlab の場合は基本的に float (double) で計算されているが、Julia の場合は厳密に違う型として認識されている。
たとえば、
x = [3 4 5]
とすると出来上がる型は 1x3 Array{Int64,2}: である。
したがって、
x[1] += 2.3
とすると、float の値となる 3 + 2.3 = 5.3 を x[1] に代入しようとして
ERROR: InexactError()
で失敗する。

このときには、
x = float([3 4 5])
と float をかませておけば回避できる。


(3) max と maximum の違いに気をつけろ!
Matlab の max 関数は、Julia では max と maximum に機能の違いで分割されている。
max は要素を引数として並べることできる状態、maximum なら配列が引数となる。
たとえば、

max(2,3) と直接 2と3が引数となる場合は max, maximum([2,3]) と配列[2,3]が引数ならmaximum である。
逆に max([2,3]), maximum(2,3) とするとエラーで失敗する。




2015年5月11日月曜日

Julia の面白いところ(数理最適化のプログラミングから見て)

この前、内点法を実装してみて分かった、面白いところを列挙してみる。

1. UTF-8 の文字列を変数に使える。例えば
  α = β + γ
なんて式が可能。

2. A ? B : C が使える。(C っぽい)

3. リスト内包が使える。(python っぽい)
これは、SDP のときに便利で、例えば内積の A_k \bullet X というのを A_1,...,A_m で行いたければ
float([ sum(A[:,:,k].*X) for k=1:m ])
とすれば一行で作れる。
(float をかけているのは、リスト内包で上のリストを作るとなぜか Any のリストになってしまって、その後に足し算引き算などができなくなってしまうので、float に型を強制変更している。)

4. Cell にできる。(Matlab っぽい)
たとえば、
Z = cell(3,2)
Z[1,2] = [3 4; 5 6]
Z[3,1] = spzeros(7,8)
のように違うデータを入れ込むことができる。


3 と 4 は SDP のときには相当便利で、これがあれば SDPA-M 入力と SeDuMi 入力の両方のアドバンテージを一つのデータでこなすことができる。

2015年5月5日火曜日

Julia で SDP の内点法を実装してみた

次のようなソースファイルを sdpa01.jl として保存しておいて
julia > include("sdpa01.jl")
とすると内点法で SDP を解ける。例題は、SDPA にある example1.dat-s を解いている。


# input data
m = 3;
n = 2;
A = zeros(n,n,m)
A[:,:,1] = [10 4; 4 0]
A[:,:,2] = [0 0; 0 -8]
A[:,:,3] = [0 -8; -8 -2]
b = [48, -8, 20]
C = [11 0; 0 -23]

# parameters
ε = 1.0e-6
γ = 0.9
β = 0.1
λ = 1.0e+4
maxIter = 30

# initialization
iter = 0
X = λ*eye(n);
Y = λ*eye(n);
z = zeros(m,1);

M = zeros(m,m); # Schur complement
# main loop
while iter < maxIter
iter += 1
# compute residual
p = b - float([sum(A[:,:,k].*X) for k=1:m])
   # float is necessary here
D = C - Y - sum([ A[:,:,k]*z[k] for k=1:m])
max_p = maximum(abs(p));
max_D = maximum(abs(D))
μ = sum(X.*Y)/n
residuals = max(max_p, max_D, μ);
if  residuals <= ε
@printf "iter = %3d, r_p = %.2e, r_d = %.2e, μ = %.2e, α = %.2e\n" iter max_p max_D μ α
println("Converge !!")
break
end

Yinv = inv(Y)
# build Schur complement
for j=1:m
XAY = X*A[:,:,j]*Yinv
M[:,j] = float([sum(A[:,:,i].*XAY) for i=1:m])
end

# search direction
dz = M \ (p - float([ sum(A[:,:,k].* (β*μ*Yinv - X - X*D*Yinv)) for k=1:m]) )
dY = D - sum([ A[:,:,k]*dz[k] for k=1:m ]);
dX = β*μ*Yinv - X - X*dY*Yinv
dX = (dX + dX')/2;

# step length
Lx = chol(X,:L); Ly = chol(Y,:L)

(Dx, Vx) = eig(inv(Lx)*dX*inv(Lx)')
eigx = minimum(Dx)
αp = eigx < 0 ? -1.0/eigx : 100.0
(Dy, Vy) = eig(inv(Ly)*dY*inv(Ly)')
eigy = minimum(Dy)
αd = eigy < 0 ? -1.0/eigy : 100.0

α = min(γ*αp, γ*αd, 1.0)
# update
X += α*dX
Y += α*dY
z += α*dz
@printf "iter = %3d, r_p = %.2e, r_d = %.2e, μ = %.2e, α = %.2e\n" iter max_p max_D μ α
end

if iter <= maxIter
println("--Optimal Solution--")
println("X = ", X)
println("Y = ", Y)
println("z = ", z)
else
println("Failed to converge during ", iter, " iterations.")
end

2015年4月28日火曜日

エルデシュナンバーの平均値はどうなるんだろうか?

ふと思うと、エルデシュナンバーの平均値は、年を追うごとに変化しているはずで、平均値はどのような変動をするのか予測できたら、それはそれで面白そうでもある。

各研究者でいえば、エルデシュナンバーは単調減少する。
逆に現状では2より小さくなることはできないし、今から100年経過したら3より小さくするのは困難になるはずで、200年も経たないうちに4より小さくはできないはずだろうから、下限は単調増加でもある。
この2つのバランスがどっち側に傾くのか、そのあたりが面白そうでもある。

2015年4月24日金曜日

現時点での SDPARA のコンパイル方法

温故知新、ということで現時点での Debian 上のSDPARA のコンパイル方法をまとめておく。

(1) ほとんどは Install.txt にある情報でコンパイルできる
(2) MUMPS は、4.8 だとエラーが起きてコンパイルを通せなかったので、4.10.0 に変更。使用するバージョンに合わせて mumps/Makefile の build/lib/libdmumps.a も書き換える。
(3) openmpi のファイルが変更になっているので、以下のような修正を Makefile に行う。


SCALAPACK_LIBS  = /usr/lib/libscalapack-openmpi.a /usr/lib/libblacsCinit-openmpi
.a /usr/lib/libblacs-openmpi.a /usr/lib/libblacsCinit-openmpi.a
FORTRAN_LIBS = -lgfortran -lmpif77

特に、FORTRAN_LIBS の -lmpif77 を入れないと MUMPS の一部のファイルで「関数が見つからない」というエラーが起こる。

これで一通りコンパイルも通過して、実行可能となる。

Harmony Search Algorithm

メタヒューリスティクスのひとつに Harmony Search Algorithm というものがある。

A self-adaptive global best harmony search algorithm for continuous optimization problems
http://www.sciencedirect.com/science/article/pii/S0096300310001128

音楽における即興作曲などでどのように曲が変わるか、というところからアナロジーを得て作られているアルゴリズムのようだ。
遺伝的アルゴリズムと似たような雰囲気もあるが、遺伝的アルゴリズムよりも速く計算ができる、とのこと。



2015年4月20日月曜日

整数計画問題による定式化

JIMA の Vol. 66, No.1 に載っている論文のうち、

1.サプライチェーン途絶に対する最適リスク緩和方策
2.災害救助活動のためのロジスティクス・モデルに関する研究
3.人間ドックにおける最適スケジュール作成について

の3つにおいて整数計画問題による定式化が行われていて面白い。

特に、1,2は大震災などに対してどのようにことで数理最適化が役にたつのか、という意味でも勉強になるし、3のようなことも社会に使われる最適化、という意味合いがある。


2015年4月15日水曜日

ぷよぷよと確率最適化

以前に「ぷよぷよ」が NP 完全であることが示されている。
ふと思うに、確率最適化の理論を使ってぷよぷよを解いたら、対戦相手としてはどれくらい強くなれるのであろうか?
将棋や囲碁よりも計算時間の制限がきついので、それなりに難しいかも。
むしろ車の自動運転とかの考えを使ったほうが計算時間的には有利?


2015年4月14日火曜日

非対称錐での内点法

非対称錐についての内点法についての論文が

A homogeneous interior-point algorithm for nonsymmetric convex conic optimization
http://link.springer.com/article/10.1007%2Fs10107-014-0773-1

に出てきていることに気がついた。

ポイントとしては、

  1. 主問題側だけのバリア関数を使っている
  2. homogeneous な定式化を用いている
  3. 計算量は通常の内点法と同じオーダーになる
  4. Runge-Kutta の方法を用いて、探索方向を修正している
といったあたり。

パッと思いつく、改良できそうな方向としては、
  1. DNN に特化して計算方法を改良する
  2. homogeneous ではなくて、infeasible interior-point method にできるかどうか
  3. 並列計算で大きな問題を解く
のあたり。


2015年4月7日火曜日

東日本大震災とオペレーションズ・リサーチ

最近の OR/MS Today にインタビュー記事が掲載されていて、東日本大震災のあとにオペレーションズ・リサーチの研究者(本人がいうには数学者)が来ていて、避難計画の策定などにも関係している、ということが分かる。

http://viewer.zmags.com/publication/e9a496eb#/e9a496eb/38

このインタビューでは、東日本大震災だけでなくエボラ出血熱との関わりも述べられている。(OR/MS Today はアメリカのINFORMSなので、エボラ出血熱についての記述のほうが記事では多い。)

オペレーションズ・リサーチを通してどのような貢献ができるのか、という点でとても参考になる。まったく同じ貢献でなくても、このような貢献をするためにはどのように自分は進むべきか、ということを考える起点にもなりそうである。また、折を見て、読み返してみようと思う。

2015年4月3日金曜日

Renegar の first-order method

Nesterov の first-order method を、Renegar が SDP 用にしたものが arXiv で読むことができる。

Efficient First-Order Methods for Linear Programming and Semidefinite Programming
James Renegar
http://arxiv.org/pdf/1409.5832v2.pdf

SDP の定式化を変形することによって、Nesterov の手法を使えるようにする、というのがミソで、この変形を行うと制約式が線形制約だけになる。(問題のむずかしさが目的関数に集約される感じか。)


なかなか面白い手法ではあるが、その一方で
Practical first order methods for large scale semidefinite programming
Stephen Tu, Jingyan Wang
http://www.cs.berkeley.edu/~stephentu/writeups/first-order-sdp.pdf
を見るとソフトウェアへの実装は難しいようである。
どのあたりに難しさがあるのか、ということを確認すると、新しい研究に結びつくかもしれない。


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 などの映画に関するデータなども数値計算に取り込んでいる。

2015年2月27日金曜日

Stochastic optimization の基本的なキーワード

A Tutorial on Stochastic Programming
Alexander Shapiro and Andy Philpott

を参考にすると、 Stochastic Optimization 関係のキーワードとしては、

  • Two-stage optimization
  • 複数のシナリオを発生させてからの最適化
  • Sample Average Approximation (SAA)
  • chance constraints
    • VaR (Value at Risk)
    • C-VaR (Conditional Value at Risk)
  • Robust Optimization
といったあたりが重要。

2015年2月23日月曜日

Debian では -fPICが非推奨とのことらしい

Linux では、lib???.a ファイルは重要なファイルであるが、Debian ではこの lib???.a を作る際に -fPIC を使うのが非推奨、とのことらしい。
これによって、たとえば apt-get でインストールされる libABC.a というライブラリを使用して libXYZ.so などの shared library の作成ができない、ということが標準のようである。

したがって、たとえば、octave の *.mex を作成するには、*.mex の実態が shared library であることから、libABC.a などの利用に制限がかかるのが標準、ということらしく、apt-get でインストールすのとは別途 libABC.a をソースからコンパイルする必要がある。

2015年2月20日金曜日

素晴らしい対応の Python Extension の Unofficial Windows Binary

Windows での python の状況をチェックしていて、
Unofficial Windows Binary
http://www.lfd.uci.edu/~gohlke/pythonlibs
を利用し始めていたが、今までの CVXOPT バイナリにはバグが残っていて、自分のプログラムではエラーが頻発して使い物にならなかった。

このバグは Netlib の BLAS/LAPACK をリンクしているためであって、MKL とリンクするbinary はバグを取り除けるが CVXOPT のライセンスで配布ができないとのことであった。
そこで、OpenBLAS をリンクすればライセンスの問題は回避できるかも、というメールを送ってみた。

すると、わずか1時間後には CVXOPT が OpenBLAS とリンクされた binary がアップロードされていた。
http://www.lfd.uci.edu/~gohlke/pythonlibs/#cvxopt
しかも、複数の python バージョンで、である。

驚くべき素晴らしい対応すぎて、自分も見習いたいものである。


ちなみに、今回の binary で自分のプログラムは問題なく実行できるようになった。
Conic Programming の研究では、Matlab を使わなくても python で 90% 以上のことができているのでは?と思うが、そうなってきているのも、こうやって binary などをメンテナンスしてくれている方がいるからこそである。



2015年2月12日木曜日

Google の OR-Tools

Google にもオペレーションズリサーチ関係のソフトウェアパッケージがあるらしく、or-tools という名前になっている。

https://code.google.com/p/or-tools/

ページを見てみると、どうやら入っているのは

  • 制約プログラミング
  • トラック配送問題
  • GLPK, Gurobi, SCIP などのインターフェイス
  • ナップサック問題
  • グラフ問題(最短路計算、最大流問題など)
ということの様子。
インストールについては、現段階ではソースからのビルドのみのようで、まだまだ発展途上の段階のようでもあるが、そのうちに使い道が出てくるかもしれないので、少し注目しておくことにする。

2015年2月10日火曜日

Matlab から python への移行を検討し始める

最近、Matlab の周辺環境がいろいろと使いにくくなりつつあって、市販ソフトウェアとしての限界が出てきつつあるので、python への移行を検討し始めている。

とりあえず検討すべきこととしては、
1. SDPA などの主力ユーザは Windows なので、Windows でどれくらい簡単に利用できるかを把握する(python は Linux 中心で発展してきているので、Linux については特に問題ないかと推測している)
2. Matlab からどれくらい移行が大変かを確認する。文法的な差異はそれほど問題ではないが、BLAS/LAPACK や疎行列へのアクセスがどれくらいできるのか、という点が気になるところ。


2015年2月6日金曜日

永久機関と余弦定理

SIAM News の記事に、永久機関と余弦定理が関係していることが書かれている記事があった。

http://sinews.siam.org/DetailsPage/tabid/607/ArticleID/359/Perpetual-Motion-and-the-Theorem-of-Cosines.aspx

「永久機関が存在しない、ということで式を立てると、そこから余弦定理が導ける」といういもので、意外なところの結びつきが面白い。


2015年2月4日水曜日

Windows 上での Python モジュールのビルド

Windows 上で Python モジュールをどのようにビルドするのか、という内容が、

http://docs.python.jp/2/extending/windows.html

にまとめられている。
Matlab の Mex ファイルも、Python モジュールと同様に実体は DLL ファイルなので、とても参考になる。(Matlab の場合は、拡張子を mexw32, mexw64 としている)

SDPA-M の Windows 版は Linux 上でクロスコンパイルしているが、上記ページの「4.2. Unix と Windows の相違点」が勉強になる。


2015年2月2日月曜日

主双対内点法でまだ解析されていないこと

SDPA や SeDuMi は主双対内点法を実装しているが、これらのソフトウェアを実行していると実際の性能と理論的枠組みには大きなズレがあることもわかる。特に重要な点として、

1.主双対内点法は最適解の近くだと2次収束するが、実際のソフトウェアでは数値誤差の影響もあり、1次収束までの動作しか得られない。わかりやすくいうと、C言語の double で実装した場合、主双対内点法では 10^{-8} 程度の数値誤差が出てしまうが、2次収束を開始するのは 10^{-8} よりも小さいところである。(乱数で問題を生成すると問題の性質がよくなり、もう少し早いところから2次収束する傾向が見られやすい。)

2.ところが1次収束する他のアルゴリズムと比較すると、最初の数反復での収束の速さは速く、その後の1次収束するときの係数も小さい。

などがある。
特に2のところは不思議で、理論的な収束としては1反復で duality gap が $$ 1-\sqrt{n} $$ ぐらいの比でしか小さくならないはずであるが、実際に動かしてみるとは 0.1 ~ 0.5 ぐらい出ているのでは?と思える。



2015年1月20日火曜日

python の起動と Matlab の起動を時間で比較

Python を Anaconda 経由でインストールしてみたが、起動に時間がかかるのでどのくらいかかるのかを Matlab の起動と比較してみた。
結果としては、

Python
Spyder のロゴが出るまで 20秒
Spyder のWindow が出るまで +40秒
Spyder でコマンド入力できるようになるまで + 10秒
合計で1分10秒かかっている。

Matlab
Matlab のロゴが出るまで 2秒
Matlab の Window が出るまで + 30秒
Matlab でのコマンドが入力できるようになるまで +3秒
合計で 35秒かかっている。

Matlab は、起動にやけに時間がかかるときがたまにあって安定した時間ではないのだが、Python の Spyder はざっと Matlab の2倍程度の起動時間がかかるようだ。

2015年1月14日水曜日

maxNumCompThreads の代わり

Matlab で最大スレッド数を知るために maxNumCompThreads を利用しているが、実査に使うと

警告: maxNumCompThreads は将来のリリースで削除される予定です。この関数のインスタンスをコードから削除してください。
> In maxNumCompThreads at 26

という警告文が出る。もう3,4年は出ている警告なのでいつ削除されるのか、そもそもほんとに削除するつもりがあるのかどうかよくわからないが、警告文が出る。

この場合には、

>> numberOfThreads = feature('NumCores');

とすると、同等の機能を得ることができる。

2015年1月13日火曜日

ルービックキューブの最大必要手数は20である

今回の SIAM Review には、「ルービックキューブの最大必要手数は20である」という論文が掲載されている。

The Diameter of the Rubik's Cube Group Is Twenty,
T, Rokicki, H. Kociemba, M. Davidson, J. Dethridge,
SIAM Review, Vol. 56, No. 4, pp 645--670, 2014

単純に数学的な側面の議論だけでなく、ソフトウェア実装なども議論されている(gcc によるアセンブラコードなども登場)あたりが面白いところかと思う。



2015年1月8日木曜日

SDPAフォーマットへの定式化についての対応

SDPA 関係での問い合わせは、比較的コンスタントにあって、最近もいくつかあるが、やはり SDPA フォーマットへの定式化についてのものも多い。

このような場合は、YALMIP を使うと簡単に SDPA フォーマットに変更できるので、YALMIP をおススメする場合もある。
また、最終的には Matlab を使わずに C++ などだけで実装したい、というような場合もあって、そのような場合には、まずは YALMIP でsavesdpafileにより SDPA フォーマットを生成してから、その SDPA フォーマットと照らし合わせながら作業を進めたり、ということもある。


2015年1月5日月曜日

詳しくSDP の双対定理についての証明過程が載っている本


SDP の理論展開において、双対定理は重要であるが、いくつかの証明などを組み合わせて証明する必要がある。
この証明の過程をフローチャートでまとめてあるのが、

Gartner, Matousek,
Approximation Algorithms and Semidefinite Programming,
Springer, 2012

の Section 4.8 The Largest Eigenvalue の Figure 4.5 (P69) である。
この図を見ると、どのような証明が組み合わさって双対定理が示されているかが分かりやすくなっている。