2017年7月31日月曜日

平均が簡単なようでいて難しい

簡単に計算できるかな、というか解析的に求まるのかな、と思っていた平均が実はそれほど簡単でもなかった、というのがあって、以下のような問題が思っていた以上に難しい。


問題:「一辺が長さ 1 の正方形にランダムに2点を置いた場合の2点間の距離の平均はいくらか?」

それほど厳密な解がほしいわけでもないので、サクッとモンテカルロシミュレーションを Julia で作ってみたら、だいたい 0.5214 という解になった。
この 0.5214 って数字、なんだろう?


ちなみに、作ったプログラムは
sum0 = 0.0
T = 10000*10000
for i=1:T
    x1 = rand()
    x2 = rand()
    y1 = rand()
    y2 = rand()
    dist0 = sqrt((x1-x2)^2 + (y1-y2)^2)
    sum0 = sum0 + dist0
    if rem(i, 10000000) == 0
        println("iter = $i, average = $(sum0/i)")
    end
end

で 40 秒ぐらいで計算してみている。
julia> @time include("randdist.jl")
iter = 10000000, average = 0.521513779364118
iter = 20000000, average = 0.521471850234503
iter = 30000000, average = 0.5214132474790024
iter = 40000000, average = 0.5214231135462779
iter = 50000000, average = 0.521431221770643
iter = 60000000, average = 0.5214123606643466
iter = 70000000, average = 0.5213983220702314
iter = 80000000, average = 0.521385077693363
iter = 90000000, average = 0.5213964364889536
iter = 100000000, average = 0.5213935109322017
 41.304951 seconds (600.00 M allocations: 10.431 GB, 1.90% gc time)

ちなみに、1行で計算すると4倍ぐらい速かった。

julia> @time begin T = 10000*10000; x = rand(T,4); ave = sum(sqrt((x[:,1] - x[:,2]).^2 + (x[:,3]-x[:,4]).^2))/T; println("T = $T, ave = $ave"); end
T = 100000000, ave = 0.5213702192902951
  8.407286 seconds (45 allocations: 10.431 GB, 9.60% gc time)




2017年7月26日水曜日

新しい最適化の概念

数理最適化の分野では first-order method のような手法が「大規模なデータ」を扱えるということで多くの研究がされてきたが、昨日チェックしたもう1本の論文では、「first-order method が扱えるような小さなデータではなくて、ちゃんとした規模のデータを扱いたい」という話が載っていた。

そのために、今までの数理最適化の概念とは全く異なる最適化の概念を導入していて、簡単にいうと
「最適値や最適解を求めない」
というもの。

「なるほど、そういう考え方があったか」という気分で、目から鱗がパラパラ。

2017年7月25日火曜日

Home Health Care の論文をチェック

今回は、

OR problems related to Home Health Care: A review of relevant routing and scheduling problems,
Operations Research for Health Care (2017), http://dx.doi.org/10.1016/j.orhc.2017.06.001

という論文をチェックしてみた。
簡単にいうと、「自宅療養をしているような患者さんが複数いる場合に看護師がどういう順番で患者さんを見て回ったらいいか?」というような内容で研究している論文のサーベイだった。

結構面白いテーマかと思って読んだわけだけど、それぞれの論文でケースバイケースでアルゴリズムを作っている感じだった。ある意味で行き当たりばったりなアプローチが多く、うまく行く理由がたまたまなのか、ちゃんと科学的理由があるのか、検証がなされていないようだ。

一番致命的なのは、他の研究者がデータにアクセスできない(たぶん、公開を前提としてのデータ収集が行われていないためだろうな、とは思う)ことで、論文の結果を再現できないことだと思う。

このサーベイで指摘されているように、ベンチマークのようなデータセットがあって、それで試すことで本当に有益な手法なのかが分かるようになるのかもしれない。

2017年7月21日金曜日

クロスドッキングセンターの出庫エリアにおけるシュート・ドック割り当てとその解法の提案

論文を読んでみたので、思ったあたりをメモ書きしておく。

1.問題を mixed integer linear programming (MILP) で定式化する
2.MILP が小さい場合は Gurobi で解く
3.大きくなる場合は、遺伝的アルゴリズム+局所探索で解く
4.数値実験では、MILP が小さい場合は Gurobi で解いた解と同程度の解を遺伝的アルゴリズム+局所探索で得ることができる


読んでみた限りだと、遺伝的アルゴリズムと局所探索の実装は、問題にあまり特化せずに比較的スタンダードな実装が行われているようだ。
数理最適化の視点だと、ガッチリと問題の性質を使いまくった複雑なアルゴリズムを作って、その定理とか証明したりなっちゃうかもしれないけど、実用上では厳密解を求める必要がないから、こういうアプローチはありだなぁ、と思ったりする。


2017年7月19日水曜日

内点法の謎

内点法のソフトウェアで大きな謎なのが

「理論的には近傍を作ったほうがいいのに、数値実験では近傍を作らないほうがパフォーマンスがよい」

ということ。
近傍を作らなかった場合の解析というのは、まだ改良の余地があるのではないか?と思ったりする。

2017年7月16日日曜日

MUMPS がバージョンアップしていた

いままで気がついてなかったけど、SDPA の Debian パッケージがコンパイルできなくなっている、という連絡が来ていて、その理由として SDPA が利用している MUMPS の Debian パッケージの バージョンが上がっていることがあるようだ。

いままで何年もずっと 4.10 のままだったけど、5.1.1 になったようだ。

ちょうどいい機会でもあるので、SDPA も OpenBLAS に切り替えたりしつつ MUMPS の調整をするのがよさそうだ。
とはいえ、MUMPS の Debian パッケージも sid にあがってきたばかりのようなので、少し下調べしてみようと思う。


2017年7月4日火曜日

Facial Reduction 難しいなぁ

SDP に内点がない場合、Slater の制約想定を満たさないので、理論的には内点法ではうまく解けない。
そういったときに内点がある問題に SDP を修正する枠組みとして、Facial Reduction がある。
ということで、それを手元にある SDP に適用したりしてみているが、これがSDPソルバーで解いてみると、あんまり上手くいかない。簡単にいうと、精度が向上しない。

結構重要なこととして、
「SDP ソルバーで十分な精度で解けるかどうか、と、内点があるかないかというのは、それほど関係がない」
ということがある。
やはり、double 型が 64bit という影響は強くて、64bit 用の Facial Reduction というのを作ればいいんだろうなぁ、と思ったりしている。