2017年12月4日月曜日

論文読み:Efficient Convex Optimization for Linear MPC

今回も Optimization Online のところから、論文をチェックしてみた。
今回は、

Efficient Convex Optimization for Linear MPC
http://www.optimization-online.org/DB_HTML/2017/11/6353.html

基本的には、線形モデル予測制御の問題が、目的関数が凸2次、制約が線形制約になるので、それをどのように効率的に解くか、という話。立ち位置としては、サーベイに近いような感じだと思う。
解法としては、内点法と active-set methods が書かれていて、内点法は比較的スタンダードなものに線形制約に現れる帯構造を使って時間短縮を図るというもの。
あと、active-set methods は warm-start のほうでメリットが大きい、とのこと。


個人的には、active-set methods をあまり使ったことがなかったので、もうちょっと調べてみようと思う。手元にはたくさんの対角ブロックになるような SDP があるので、それをうまく解けるといいなぁ、と思っている。



2017年12月2日土曜日

Adobe は、なぜ学習しないのか?

Adobe Acrobat Reader DC の注釈で書いていると気になるのだが、漢字変換の学習機能が効かない。

例えば、「この節では」と Windows で何回も入力していても、「この説では」と変換される。Adobe Acrobat Reader DC でも注釈で「この節では」と何回入力しても、学習されないようだ。

Adobe だと、なぜ学習しないのか、不思議だ。
ひょっとしたら変換辞書が別にあるのか?


ちなみに、普段はAdobeではない別の PDF リーダーを使っていて、そっちでは快適に漢字変換ができている。


PowerPoint の数式が思い

プレゼンテーションのファイルを作るのに、以前は Beamer を使っていたこともあるけど、PowerPoint + IguanaTeX の自由度の高さに気がついてからは Beamer には戻れなくなってしまった。
(個人的に思うところでは、Beamer を使い慣れている人には、PowerPoint + IguanaTeX をあまり勧めない。まぁ、Beamer と PowerPoint + IguanaTeX を比較しても大きな差はなく、例えていえば Beamer が普通の黒電話ぐらいで PowerPoint + IguanaTeX がスマートフォンぐらいという差の程度でしかない。無理して乗り換える必要もない。)

今回は、IguanaTeX を使わずにPowerPoint の数式だけでプレゼンテーションファイルを作ったらどうなるか、ということも試してみたくてやってみることにしてみた。
そもそも IguanaTeX だと TeX 環境をインストールしないといけないので、TeX がインストールされてないところでは編集ができない。そこで PowerPoint の数式だけでできるかを知っておくことも有益だろうと思ってやってみた。

結果として、PowerPoint の数式は、思っていた以上にスマート。入力に多少の癖はあるけど、TeX の数式を入力する知識がだいぶ役に立つので、Beamer を使っている人にもPowerPoint の数式は使いやすいと思う。

逆に困っているところは、以下の2つ。
(1) 行列を書こうとすると Unicode が入り込んでしまって、行形式での編集が極めて困難。
(2) 数式が多くなると、PowerPoint の処理が極端に遅くなる。例えば、"PowerPoint" という文字列を入力するだけでも 10 秒以上かかる。


まぁ、自分のところでは試していないけど、Office365 では LaTeX で数式を書けるようになるということなので、このあたりが解決するのは時間の問題かもしれない。
やはり LaTeX の数式はキーボードだけで書けるというメリットが大きいので、できれば Google スライドとか Office Online とかでも使えるようになるといいなぁ、とは思う。

2017年11月22日水曜日

IguanaTeX のバグを回避する

IguanaTeX は非常によくできているので、よく利用させてもらっているが、たまにバグがある。
例えば、一時ファイルを置く場所を Main Setting に入力しても、これが環境変数として ghostscript に渡らないため、TeX2imgc.exe などを使おうとしても途中の ghostscript で Unable to open the initial device, quittingとなって終了してしまう。

これを回避するには、TeX2imgc.exe を呼び出す bat ファイルを以下のように作って、TeX2imgc.exe の代わりに Main Setting で設定する。

%%%%%%%%%%%
set TEMP=[一時ファイルの場所]
set TMP=[一時ファイルの場所]
"[TeX2imgc.exeへのパス]\TeX2imgc.exe" %*
%%%%%%%%%%%
つまり、環境変数 TEMP, TMP を設定して、あとは TeX2imgc.exe にすべての引数を渡して実行する。

2017年11月17日金曜日

待ち行列の研究って、純粋数学になっていくよね、と思ったりすること

病院の患者などを研究している論文があって、その中に

1.患者の待ち時間
2.患者が病院にいる総時間
3.ベッドの使用率

をシミュレーションで解析している論文があった。

自分としては、「これって、待ち行列で解析すればいいのでは?」と思ったわけだが、論文ではシミュレーションソフトを使って解析していた。

そこで「待ち行列で解析」という考えを一度横において考えてみると、シミュレーションソフトを使って待ち行列を解析して、それで実用的な結果が得られる、というのは、確かにありかもしれない。

数理最適化についても、コンピュータの性能が上がればシミュレーションで十分になる範囲も格段に向上するだろうし、現段階では応用数学に分類されることが多い数理最適化も将来的には純粋数学に仲間入りすることもあり得るのかも、と思ったりもする。




2017年11月10日金曜日

Julia で返り値を複数の中から一つだけ取り出したいとき

たとえば、対称行列 A の固有値分解は

D, V = eig(A)

とすると得られるが、ここで V は不要だから D だけ欲しい、というときも結構ある。
その時は
D = (eig(A))[1]
とすると最初の引数だけが得られる。

もちろん、2番目がほしい時は
V = (eig(A))[2]
とすればよい。

ちなみに、配列の中で順番を入れ替えることもできて
V, D  = (eig(A))[[2 1]]
としたりすればよい。
ただし、[[2 1]]は2重カッコのように見えるが、(eig(A))[2 1] のように1重カッコだとエラーになって処理されない。

unimodular な SOCP

いま手元で作っている SOCP が良く良くデータを見てみたら、
||U*x|| <= a'*x + b
という形式で、行列の U のところが unimodular になっていた。

unimodular な場合っていうのは、何かに使えたりするんだろうか?


2017年11月8日水曜日

論文読み:From feasibility to improvement to proof: three phases of solving mixed-integer programs

今回は、Optimization Online とは別に、雑誌に掲載された論文をチェックしてみた。

From feasibility to improvement to proof: three phases of solving mixed-integer programs
http://www.tandfonline.com/doi/full/10.1080/10556788.2017.1392519

整数計画問題を厳密に解く時には分枝限定法が一般に用いられるが、この論文は以下のような良く見られる状況に着目している。
1.最初の実行可能解を得るまでの時間は、最適解を得るまでよりも短い
2.最適解を得るまでの時間は、その解が最適であると保証するまでの時間よりも短い

したがって、分枝限定法を3段階に分けていて、
 1.最初の実行可能解を見つけるまで
 2.最適解を見つけるまで
 3.最適解が最適であると保証できるまで
これらの3段階でどのように推移していくか、というヒューリスティクスを議論している。
数値実験を見る限り、それなりにうまく動いているヒューリスティクスのようだ。

整数計画問題を近似的に解く方法として様々なヒューリスティクスが提案されているが、このような分枝限定法のヒューリスティクスを用いて他のヒューリスティクスを改良したら面白いかもしれない。


2017年11月7日火曜日

論文読み:Minotaur: A Mixed-Integer Nonlinear Optimization Toolkit

Optimization Online のところから、もうひとつ論文読み。

Minotaur: A Mixed-Integer Nonlinear Optimization Toolkit
http://www.optimization-online.org/DB_HTML/2017/10/6275.html

基本的には、タイトルにある通り、混合整数非線形最適化問題を解くソフトウェアを作りました、という内容。
アルゴリズムも分枝限定法を使っていて、内部的には CLP, Gurobi, CPLEX, BQPD, qpOASES, FilterSQP, IPOPT を呼んでいる。
インターフェースは C++ なので、実装するのは結構大変。

できれば、Julia にある JuMP.jl から呼び出せるようになってくれると便利だと思う。(JuMP.jlから呼び出せるライブラリは、まだまだ拡大の余地が大きいので。)

2017年11月6日月曜日

経営工学会でフラストレーション溜まりまくり

先週にパシフィコ横浜で行われた経営工学会の秋季大会に参加してきて、いろいろと発表を聞いてきた結果、フラストレーションが溜まりまくった、という話。

いろんな発表があるわけだけど、「この研究は、こういった数学を使ったら、もっと解析できるかも」というのが多くて、数理最適化としては研究の素になりそうな内容がゴロゴロと転がっていた。結果として「あれもやってみたいし、これもやってみたい。どれからやってみればいいのか分からん」というところになった。

やはり、一度に複数はできないので、とりあえず一つずつ検討してみたほうが良さそうにも思ったりしている。

2017年11月5日日曜日

論文読み:Enriching Solutions to Combinatorial Problems via Solution Engineering

11月の optimization online の更新が来ていたので、いくつかを論文読みをしてみた。
そのうちの一つが

Enriching Solutions to Combinatorial Problems via Solution Engineering
http://www.optimization-online.org/DB_HTML/2017/10/6252.html

簡単にいうと、組合せ最適化問題では同じ最適値を達成するような最適解が大量にある場合があって、そのようなときに大量の最適解から特徴的な最適解をどのように取り出すか、という話。

ざっくりと手法をまとめると、多様性を表すような目的関数を作って、別の最適化問題を解いている。
また、既に得られている最適解を切除平面で外すことで別の最適解を得たりしている。

個人的に思うに、「特徴的な解」に含まれる「特徴的な」の部分をどのように数学的に表すかは様々な方法があるから、割り当て問題などでの現実的なサイズになるとどうなるか、というあたりが結構大変そう、と思ったりする。

2017年10月24日火曜日

Chubanov の方法を凸結合したら面白そうね、って話

この前の RAMP で線形計画問題などを含む Chubanov の方法についての話があって、とても勉強になった。Chubanov の方法が流行っている、というのは既に聞いていたが、自分では勉強していなかったので、いい機会でもあった。

個人的な印象としては、
「漠然と想像していたよりも内点法に似通っている」
というところ。
基本的には、錐の中心へとスケーリングする、という計算方法がKarmarkar の内点法に近い印象になっているのではないかと思う。

対称錐の場合、homogeneous の性質が効くため、錐の中心とすべての実行可能内点は同じとみなすことができるので、Karmarkar の内点法で生成される内点の点列以外にも内点の点列が作れる、ということである。

単純に考えると、Karmarkar の内点法の生成する点列と Chubanov の内点法の生成する点列で凸結合を取れば、それら2つの方法を統合した計算手法が作れるんだろうなぁ、って思う。
そうすると、最初の反復では計算効率の良い Karmarkar 法を使って、解の近くに近づいたら Chubanov 法に切り替えて固有値計算をうまくやる、とかできるかとは思う。

2017年10月5日木曜日

Julia で Mosek 8 がこけるので、Mosek 7 を使えるようにする

Julia を Linux 上で動かしていて、ここから Mosek を呼びたいのだが、なぜか Segmentation fault が起こってしまい失敗する。
たぶん libc のバージョンとかがあっていないんだろうなぁ、と勝手に推測しているわけだが、Matlab から Mosek を呼んでも Segmentation fault が起こる。

そこで、Mosek.jl を build しなおして、Mosek 7 を利用するようにする。
まずは、julia で Mosek.jl をインストールする。
julia> Pkg.add("Mosek")
そのあとで、
$ cd ~/.julia/v0.6/Mosek/deps
に移動して、build.jl の最初にある
mskvmajor = "8"
mskvminor_minimum = 1

mskvmajor = "7"
mskvminor_minimum = 0
に修正する。
本来 Mosek.jl は Mosekのバージョンが 8.1 以降にしか対応していないので、そのバージョンチェックを 7 でも通過できるように修正していることになる。
(このあたりの不具合は、どっかで出てくるかもしれない。)

そのあとで、
$ cd ~/.julia/v0.6/Mosek/deps
$ MOSEKBINDIR=/home/username/mosek/7/tools/platform/linux64x86/bin julia -e 'include("build.jl")'
としてパッケージを build しなおす。このとき MOSEKBINDIR は Mosek のインストールされたディレクトリにあわせて適宜修正する。

これで何はともあれ Mosek 7 を Julia から呼べるようになる。



2017年10月3日火曜日

モデリング言語特集を読んでみて

今回の OPTIMA (つまり OPTIMA 103) はモデリング言語特集ということで、
JuMP (Julia), PuLP (Python), CVXPY (Python), ompr (R), ZIMPL
について書いてあって、それぞれの比較があって面白かった。

プログラミング言語としては、やはり Julia は R や Python よりも後発なだけあって使い勝手なども向上しているが、モデリング言語として考えると他との連携もあるので単純に比較するのは難しいので、結局最後は好みの問題とかこれまでの経験とかで決まるんだろうなぁ、と思ったりもする。

あと、それぞれのサンプルがあるが、サンプルとして数独が使われていることが多くて、研究者の中では数独って手ごろな問題としても面白いんだろうなぁ、って思ったりもする。

2017年10月2日月曜日

論文読み:Julia で Feasiblity Pump

Julia で Feasibility Pump ベースのヒューリスティクスを実装した、という内容があったので、チェックしてみた。

FPBH.jl: A Feasibility Pump Based Heuristic for Multi-objective Mixed Integer Linear Programming in Julia
http://www.optimization-online.org/DB_HTML/2017/09/6195.html

手法としては、多目的関数最小化問題の Pareto をどうやっていくのか、というのがメインで、Feasibility Pump はサブぐらいの位置づけだった。
あと、Julia で整数計画問題を解くところは、CPLEX.jl を使っている。


2017年9月29日金曜日

Julia の now()

Julia の計算時間を測るには、@elapsed マクロが準備されているわけだけど、これだと begin-end でくくらないといけないので、Matlab のコードを tic, toc を移植したりするのが面倒だった。

調べてみたら、Julia には Julia で now() っていうナウい関数がちゃんとあった。例えば、
julia> start_time = now();
julia> X = rand(8,8);
julia> end_time = now()
julia> elapsed_time = end_time - start_time
10704 milliseconds
のようにできるので、Matlab の tic, toc を置き換えるのも比較的簡単に行うことができる。

(ちなみに、最初の now() だけは LLVM でのコンパイルがあるので時間が少しかかる)

2017年9月21日木曜日

論文読み:白血病治療に関する最適化

Optimization of Combined Leukemia Therapy by Finite-Dimensional Optimal Control Modeling

https://link.springer.com/article/10.1007/s10957-017-1161-9
Journal of Optimization Theory and Applications,
First Online: 15 September 2017

白血病治療に関する最適化の論文が出てきていた。
基本的には制御理論の計算手法を適用しているようで、細かいところは専門用語のために追い切れていないが、こういった使い方も最適化としてできる、という点が面白いと思う。

あと、感覚的にはこういった治療とかに関しての論文は、イスラエルの研究者も多いような気がする。どこかにそういったのをカウントしているデータとかあったりするのだろうか?


2017年9月20日水曜日

ソウルに行ってきた

先週のほとんどの時間を使って、ソウル滞在をしてきた。

日本から見ると、北朝鮮のミサイル問題などでソウルに行くのは危ない、と思うかもしれないが、実際にソウルに行ってみるとそれほど大騒ぎをしていなかった。日本のテレビのようにミサイルの発射画面を繰り返して写したりするわけでもないし、避難訓練とかの案内もまったくなかった。
むしろ、着実に生活水準を上げてきているソウルと、デフレ問題などで生活水準がジリ貧で下がってきている東京の差を、きちんと注視するべきだとは思う。なんだかんだ言っても、自分が頑張る場所は東京なわけだし。



ソウルでは、たまたまスケジュールが重なったドイツからの研究者の Max clique に関するSDP緩和の話を聞くことができて面白かった。自分が普段やってないことを聞くというのは、新しい視点ができていいものだと思う。

自分のほうは、SDP 緩和を SOCP で解く方法などを考えてみたりして、普段の環境とは違うと集中力なども違うので、いい転換になったかと思う。


2017年9月8日金曜日

DSDP という良くできたソフトウェア

数理最適化のコアとなる問題は LP であるが、LP に対する数値解法は SDP などに簡単に拡張できるので、LP に対する新しい論文が出てくると「SDP とかsymmetric optimization に使いました」っていう論文は良く出てくる。
面白いと思うのは、symmetric optimization に拡張したからといって、たいていの場合は内点法の時のように新しい知見が得られるわけではない、というところ。

これと同じように full NT direction についての論文も SDP に拡張されたりなど、いろいろと類似の論文がある。
これについては、
A New Full Nesterov–Todd Step Primal–Dual Path-Following Interior-Point Algorithm for Symmetric Optimization
という、簡単に拡張できるところは基本的にカバーしている内容があったりして、参考になる。(つまり、full NT direction でSDPの論文をさらに書くのは新規性の点で難しいってこと。)


ところで、SDP の探索方向には、NT, AHO, HKM といろいろとあるわけだけど、個人的には DSDP が使っている探索方向がベストな探索方向ではないかと思っている。
数理最適化問題の物理的構造を計算に取り込む、という点で非常にうまくできている。実装する時間が見つかれば SDPA/SDPARA に採用したいぐらいだ。これは10年後か20年後の研究に必ず役に立つと思う。





2017年9月4日月曜日

逆行列求めるって、難しいよね、っていう話。

例えば、sqrt(5)*I [ただし、Iは単位行列] の逆行列は、当然ながら (1/sqrt(5))*I である。
ところが、行列のサイズが n = 1000 ぐらいになると、数値計算ではうまく行かないことがある。余因子行列を用いる方法だと行列式が必要であるが、sqrt(5)*I の行列式が Inf になってしまうのである。Matlab なら


>> n = 1000; A = speye(n)*sqrt(5); det(A)
ans =
   Inf

となってしまう。

これが単純に逆行列ならば (1/sqrt(5))*I で処理すればいいので特に問題ないが、主双対内点法で解きたいような問題にこういった構造があったりすると、数値的不安定になったりして、話が非常にヤヤコシイ。


2017年8月30日水曜日

Julia のパッケージが壊れた

Julia でプログラムを組んでいるときに、あるソルバーの数値精度が悪い状況に遭遇。しかも古いバージョンではちゃんと解けていたので、これはJuliaのところで、そのソルバーのパッケージだけ古いのにすればうまくいくのでは?と思った。

その方法は、例えば、
Pkg.pin("JuMP", v"0.9.2")
などとする。(今回問題となっているソルバーは JuMP ではないけど、例として JuMP を使っている。)
ちなみに、バージョン情報は、
$ cd ~/.julia/v0.6/JuMP
$ git tag --format '%(refname) %(taggerdate)'
などとすれば分かる。

こういったことをいくつかのパッケージで行って遊んでいたら、なんとソルバーの関数を呼べなくなってしまった。
Pkg.free("JuMP")
で元のバージョンに戻してもダメになってしまった。

ちょっと謎だが、良く分からないので、
$ rm -rf ~/.julia/v0.6/
として全パッケージを削除して Pkg.add() でインストールしなおしてみた。



2017年8月10日木曜日

学術的「問い」とは

ある書類で、「学術的「問い」」について記述する項目が新設されていた。
今までにないタイプの質問ではあるけど、非常にいい項目ではないかと思う。

いい機会なので、「学術」について改めて調べてみると
https://kotobank.jp/word/%E5%AD%A6%E8%A1%93-460349
によると
#####
①学問。専門性の高いものをいうことが多い。 「 -論文」
②学問と芸術。また、学問と技術。学芸。
#####
とのこと。

あと、こういったものは「リサーチクエスチョン」とも呼ばれていて、それ自体に関する研究もたくさんある。

2017年8月3日木曜日

デスティニー

最近、研究をしていると
「あぁ、この研究をしているのって、デスティニーだよね」
と思うことがしばしばある。

もちろん、destiny であって、よく似た単語の density ではない。

不思議なもので、何年もかかって自分のところに研究内容が来ていたりして、
「あぁ、これがご縁ってやつか」
とか
「あぁ、これを研究するために今まで勉強してたのか」
とか思ったりする。

なんとも不思議なものだ。

2017年8月1日火曜日

Julia が便利すぎる

今回計算しようと思ったのは、
「1辺が長さ1の正方形の中に n 個の点をランダムに置いた場合に、それらの巡回セールスマン問題の距離の平均はいくつ?」
ということ。

巡回セールスマン問題は NP-完全な問題なので厳密解を計算しようとしたら大変なので、モンテカルロシミュレーションを採用してヒューリスティクスを使って求めた場合の平均でいいことにした。ただ、巡回セールスマン問題のヒューリスティクスを実装したものが C や Python だと思ったほか見つからない(正確にいうと、かなりソースコードに自分で手を入れないといけない)。

ということで、Julia で探してみたら、巡回セールスマン問題のヒューリスティクスを実装してあるパッケージが既にあった。しかも、20行ぐらいを自分で書けば当面の目標は達成できる。

結局自分で書いたソースは以下のような感じで、Julia のパッケージを見つけてから10分ぐらいでコーディング完了。

##### ここから ####
using Distances
using TravelingSalesmanHeuristics

n_max = 20
for n = 1:n_max
    T = 10000
    sum0 = 0
    time1 = @elapsed begin
        for i=1:T
            pts = rand(2,n)
            distmat = pairwise(Euclidean(), pts)
            path, pathcost = solve_tsp(distmat)
            sum0 = sum0 + pathcost
        end
    end
    println("n = $n, average = $(sum0/T), time = $time1")
end
##### ここまで #####

で、実行結果は以下のような出力。
n = 1, average = 0.0, time = 0.135839147
n = 2, average = 1.0450777563744278, time = 0.182987882
n = 3, average = 1.5717746211246175, time = 0.240622039
n = 4, average = 1.8913379727886972, time = 0.309000624
n = 5, average = 2.125071816087441, time = 0.376932182
n = 6, average = 2.3193293493322282, time = 0.455320865
n = 7, average = 2.4779204158273442, time = 0.546624336
n = 8, average = 2.6150744388172296, time = 0.690397053
n = 9, average = 2.7494941691386257, time = 1.037431405
n = 10, average = 2.878494437536675, time = 1.229969995
n = 11, average = 2.984631760238664, time = 1.092407702
n = 12, average = 3.0957685529387713, time = 1.467921928
n = 13, average = 3.2022714958643905, time = 1.920219731
n = 14, average = 3.2997237925066214, time = 2.21402169
n = 15, average = 3.3949127954683385, time = 2.097072679
n = 16, average = 3.491072331414203, time = 2.212287591
n = 17, average = 3.5848683927044824, time = 2.481302509
n = 18, average = 3.6700061119169214, time = 2.758670548
n = 19, average = 3.757282880734828, time = 3.06751467
n = 20, average = 3.8405201350048874, time = 3.414271655

思っていた以上に Julia が便利。


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 というのを作ればいいんだろうなぁ、と思ったりしている。


2017年6月28日水曜日

TeX の数式を Word の数式に変換する

普段は TeX で数式を書いていると、たまに Word で数式を書かないといけないときが苦行である。
今日は、$$\sum_{i=1}^Z g_i x_i$$ という式を入れようとして、力尽きた。

で、いろいろと調べた結果、pandoc を使うと、ある程度は TeX を Word に変換できる。
例えば
$ pandoc -s sample.tex -o sample.docx
とすればよい。

もちろん、pandoc も完ぺきではなくて、数式の途中に色を使ってみたりとか、文章の途中で underline とか入れるとうまく処理ができなくなるようなので、pandoc では簡単な数式だけを処理させて、それを Word に貼りつけて色を調整したりとかするほうが楽ちんでもある。

2017年6月17日土曜日

TeX のソースコードをテキストに変換

TeX のソースコードをテキストに変換するスクリプトが
https://github.com/subbyte/tex2text/blob/master/tex2text.sh
にあった。

これを使うと、あとで英語の文法チェックのソフトなどにある程度一括でかけることができるので、だいぶ便利になる。

2017年6月16日金曜日

研究ネタになりそうなこと (SIAM Optimization 2017から)

先月参加してきた SIAM Conference on Optimization 2017 で、今後の内容になりそうな内容をメモとして残しておく。

ちなみに、セッション情報などは
http://meetings.siam.org/program.cfm?CONFCODE=OP17
で見ると分かる。
例えば、[MS1-2] はセッション MS1 の中の2番目の発表。

なお、ここにある情報は一部を取り除いてあって、完全版は別途あり。


== 2017/05/22 ==

[MS1-2] The Conic Optimizer in Mosek: A Status Report
Mosek が Ver 8 になって、40% 程度高速になったとのこと。
やはり、問題の再定式化などが要綱。
テスト問題として、以下の論文のものを用いていて、
Reduced-Complexity Semidefinite Relaxation of Optmial Power Flow Problems
http://ieeexplore.ieee.org/document/6692873/
データを作るプログラムがgithubで公開されている。
https://github.com/martinandersen/opfsdr

[MS1-3]
de Klerk の話で、First-order method で最もパフォーマンスが出ないのがどういう
状況なのかを解析するのに SDP を使うっていう話。
基本的には、内積は内積であればどの内積を使ってもよく、重要な情報は Gram 行列に集約される。
思うに、kernel trick のようなものなので、再生ヒルベルト空間の性質が聞いているのでは?
解析では、SDP を解くのではなく、特定のLagrange 乗数を使って KKT 条件を満たしている。
この Lagrange 乗数は、まだ幾何学的な解釈がないようなので、そのあたりも調べたら面白そう。
self-concordant にも似たように拡張していて、現在は内点法について拡張中。
たぶん、inexact Newton だったり、quasi Newton だったりで亜種を作れそう。
論文は Optimization Letter に To appear で、ドラフトとしては、
arXiv:1606.09365

[MS25-2]
SCIP-SDP の話。
SDPA も SDP ソルバーの一つに入れてもらっているけど、
最近の性能向上をさぼっていたためか、ちょっとダメな感じ。
もっと性能を上げるようにしないと。

[MS25-4]
Mixed-Integer Sum of Squares を Julia のパッケージで解く、という話。
Julia のパッケージとしては、MultivariatePolynomial.jl, PolyJuMP.jl,
SumofSquares.jl などがあって、これらを使うと、多項式最適化問題から
SDP 緩和を作るのがとても簡単にできる、ということ。
また、Steering a quadcopter through obstacles では、
Mixed-Integer Sum of Squares で
「障害物を避けながら進む」っていう計算ができていて、
これは計算結果をグラフィック表示で楽しむことができるので、
デモンストレーションなどにいいかもしれない。
(ただ、このデモンストレーションの計算については、論文などにはまだなっていない様子。)

[CP8-4]
サプライチェーンのロバスト版を考えることで SOCP が出てくる、ということ。
また、KKT 条件などから、均衡状態が示されている、ということなので、
SOCP の応用として覚えておくと面白そう。
というか、あっちの応用に組み合わせたら面白い。

== 2017/05/23 ==

[IP-3] Tom Luo の通信に関係する最適化問題の話。
interference の計算に WMMSE という座標降下法の一種を用いていて、
WMMSE はデファクトスタンダードとして用いられている、とのこと。
このときの目的関数はなかなかに複雑になっているので、
あとで細かいところを確認したら面白いかもしれない。

[MS39-1] min c^T x s.t. x in P \cap S で、
P は Ax <= b で S は closed set。
Disjunctive cut を outer-product-free という考えを用いて作られている。
テスト問題には GLOBALLib からとっている。
このGLOBALlib には多項式最適化のベンチマークも含まれるようなので、
いろいろとチェックしたら面白そうだ。

[MS41-2]
arXiv:1702.00526 https://arxiv.org/abs/1702.00526
問題は、min{f(x) | Qx = z, x in Conv(X), z in Z}
Subgradient, Proximal Point Method, Augmented Lagrangian Method,
Bundle Method との関係を調べることで、新しいアルゴリズムを作る。
Proximal SDM (Simplical decomposition method) では、
proximal point method の探索方向の近似を計算していて
column generation と同じ感じで探索領域が広がっていくようだ。
nonlinear SDM-GS (Gauss-Seidel) -ALM の方法を提案していて、
収束などを示しているようだ。
column geration が関係しているから、ほかのところにも使えるかも。

[MS41-3]
内点法ベースで Boyd らが作っているのが Forces というソフトウェア。
Augmented Lagrangian Filter という手法の大域的収束性を示している。
メモリが少ないような組み込み環境でも問題が解ける、ということなので、
に適したアルゴリズムとして考えれば、これは化けるかもしれない。
この話では、空冷システムのようなものの最適化計算を組み込み環境で計算していた。
http://ieeexplore.ieee.org/document/7902123/

[MS40-4]
Monte Carlo ので Adaptive Sampling Method を計算している。

[MS56-1]
Stochastic Mixed Integer Programming の計算を PIPS-SBB という
分子限定法の並列計算で行っている。
大まかに各ノードに割りあてする Coarse-grained と、もっと細かい制御で
並列計算する Fine-grained の2段階割り当てを行う。
ただ、200台ぐらいでパフォーマンスに限度が来る様子。
思うに、何台で最大パフォーマンスになるか、っていう予測があれば役に立ちそう。

[MS56-3]
並列計算の評価基準をどうするべきか、っていう話。
計算時間が遅くても、スケーラビリティーだと高く表示されてしまう危険性がある。

[MS57-3]
Birgin の発表。
ARC (adaptive regularization with cubics)
TRACE (a non-standard trust-region method that uses cubic descent)

f(x+s) <= f(x) - alpha*||s||^3
と最後が3乗になっているのがポイント。
これにより、2乗のときよりも全体の反復回数は抑えられる。
その一方で各反復の計算は2次計画問題を解く必要があって、コストが高い。
思うに、ニュートン法でも反復回数がかかるような問題に有効なのでは?
つまり、first-order < ニュートン法 < cubics のような感じ。

今後の研究では、min f(x) s.t. l <= x <= u のように
上限下限など制約を入れていったときにどうなるか、ということ。

The use of quadratic regularization with a cubic descent condition
for unconstrained optimization
(To appear SIAM J. Optim)
https://www.ime.usp.br/~egbirgin/publications/bmregu.pdf
https://www.ime.usp.br/~egbirgin/sources/bmregu/

[MS57-4]
問題の制約が線形制約な問題に対して、
Gradient Projection Algorithm (GPA) と
Linear Constrained Optimizer (LCO) を繰り返していく
PASA (Polyhedral Active Set Algorithm) というのを作った。
また、 Strong second-order assumption という仮定を導入して、
よりいろんなことを証明している。
Polyhedron への射影は PPPOJ という Hager/Zhang の方法を使う。
また、Dual Active Set Algorithm とも関係があり、
計算の途中で preconditined CG などを使ったりする。
あとは、Grippo のような Non-monotone なども使う。

テスト問題は、Cutest から持ってきている。
結果の比較は IPOPT (デフォルトパラメータ使用)と比較している。
あと、BIGBANK というテスト問題集があり、これには条件数が 10^8 などの
テスト問題が含まれている。
今回の結果は、いかに掲載されている、ということ。
An active set algorithm for nonlinear optimization with polyhedral constraints
https://link.springer.com/article/10.1007/s11425-016-0300-6

=== 2017-05-24===

[IP-4] Renegar の subgradient の話。
Plenary にもかかわらず、概要とか短めに切り上げて、思いっきりコアな数式に突入するあたりが、
うっとりするほど素敵。
通常の subgradient とは違った内容の subgradient なわけだけど、
自己双対錐の性質を使って単位ベクトルに移動してから計算するあたりが面白い。
話の内容としては、Optimization Online に乗っていたのがベース。

[CP12-1]
graph-valued version of total variation っていう話。
全てのラグランジュ乗数が定数になる場合には、実は制約なし最適化問題を解けば十分、
というのがとても面白い性質だった。
ちなみに、出てきている問題は2次計画問題に帰着していたけど、これが凸ではないので、
そう簡単には解けないはず。

[CP15-5]
カラテオドリの一般化の話。
このあたりは、理論的な枠組みと実際の数値計算がずれるかもしれないので、
もしずれるのであれば、それをどうやって修正するか、とか考えたら実用に使えるかも。

[CP15-6]
Conic Programming の実行可能でないとかいろいろと分岐するのを
ADMM  の計算過程に基づいて判別している。
このあたりは、SDPA の判別基準とも似たようなことになっているので、
その二つの論文を突き合わせて比較すると、もう少し詳細に踏み込めるはず。
(たぶん、SDPA のパラメータの omega についての考察とか。)
そういった情報をうまく取り入れていくと、面白いことがわかるかも。



2017年5月16日火曜日

IguanaTex の新しいバージョンと新しいバグ

IguanaTex のバージョンが新しくなっていて、なんとベクター画像を扱えるようになっていた。
いままでは基本的に png で画像ができて、それを貼っていたわけだけど、今度はベクターになるので、拡大や縮小したときにもスムーズにできる。


で、そのベクター画像の機能だけど、変なバグがあるようで、カラーにするとおかしなことが起こる。
たとえば、
\textcolor{red}{ここは赤色}
とかしても、なぜかカラーにならない。
ところが、Vector ではなくて Bitmap にするとちゃんと赤い文字になる。
そのあとでVectorとして画像を生成させると、今度はベクター画像で赤い文字になる。
なかなかに面白い現象だ。

いずれにしても、IguanaTex は beamer の100倍程度の便利さなので、いつまでもなくならないでほしいところだ。

Stochastic subgradient descent method for large-scale robust chance-constrained support vector machines

新しい論文
Stochastic subgradient descent method for large-scale robust chance-constrained support vector machines
https://link.springer.com/article/10.1007/s11590-016-1026-4
について読んでみた。

使っている手法としては、よくある感じのものだった。
ただ、理論的な解析がほとんどないので、どうしたんだろうと思っていたら、解いている最適化問題が非凸最適化問題だった。
したがって、first-order method によくあるような O(1/epsilon) などの収束レートは出てこない。

数値実験を見てみると、非凸最適化の定式化でも、そこそこ最適値を得られるようで、局所的最適解にはちゃんと到達できているのかもしれない。
そのあたりを詳しく解析してみると面白そう。

2017年5月2日火曜日

A doubly inexact interior proximal bundle method for convex optimization

Optimization Online に面白そうな内容が載っていたので、チェック。

A doubly inexact interior proximal bundle method for convex optimization
http://www.optimization-online.org/DB_HTML/2017/04/5960.html

基本的なアイデアとしては、inexact  Level-Set method と内点法を組み合わせたところ。
この中では、最適解に収束する、ということが証明されているが、こういった inexact の解法についてはどれくらいの収束レートになるか、ということを扱っている本があるので、それと合わせると収束レートが分かって面白いかもしれない。
気になるのは、内点法のように超一時収束できるのか、あるいはinexact が効いてあまり収束レートは良くないのか、というところ。

あと、主問題の解を構成していたり、Bregman distance を使っていたりするところは、他にも使ってみたい。

2017年5月1日月曜日

Model Building in Mathematical Programming

Model Building in Mathematical Programming
http://as.wiley.com/WileyCDA/WileyTitle/productCd-1118443330.html
をぱらぱらとチェックしてみた。

タイトルにある通り、数理最適化問題にどのように定式化するか、というのがメインだけど、いろんな定式化が載っていて為になる。

たとえば、x = {0} \cup [a, b] \cup {c} のように x の範囲がとびとびになってしまう場合に、どうやると整数計画問題にできるか、など書いてあったりもする。

実際の問題を数理最適化問題にするときには、いろんな定式化があって、どの定式化を選ぶかで計算時間に大きな影響があるので、こういうのを読んでおくっていうのは、いいことなんだろうなぁ、と思ったりする。

というか、日本語訳を出したらいいんじゃないだろうか、とも思う。
できれば最新の第5版で。(いろいろと追加されているし。)


2017年4月26日水曜日

Julia の Windows 用エディタとしては Notepad++ がベスト

Julia の Windows 用エディタとして、Atom や Visual Studio Code, IJupyter などを試してきたが、これらでは今一つであった。
最大の問題点は、動作が遅いことで、特に Atom はショートカットをクリックしてから30秒以上しないと Julia のソースファイル編集に入れない状態だった。

結局のところ、Notepad++ と REPL を使うのが現状ではベストのようだ。

便利なところとしては、

1.Julia 用の表示にできる

https://github.com/JuliaLang/julia/blob/master/contrib/Julia_Notepad%2B%2B.xml
このファイルをダウンロードして、Notepad++ のフォルダに入れて、[言語]からユーザー定義として「インポート」すれば、Julia 用にできる。

2.編集中のファイルのあるフォルダで Julia の REPL を起動できる。
実行コマンドの登録は、
http://so-zou.jp/software/tech/tool/editor/notepad-plus/commands/
に載っている方法でできるが、Julia の REPL を起動するときには、

cmd /c "cd /d $(CURRENT_DIRECTORY) & c:\Users\ユーザー名\AppData\Local\Julia-0.5.0\bin\julia.exe"

を登録すればよい。ただし、ユーザー名と Julia のバージョンで若干の修正が必要である。
このショートカットを、例えば Alt+F5 に割り当てれば、Alt+F5 で REPL が起動する。
REPL の中で cd で移動するのは大変なので、これが地味に便利である。

3.自動保存ができる。
Notepad++ のプラグインで Autosave を入れると、オプションで設定すれば Notepad++ からフォーカスが外れたタイミングで自動保存をできる。
つまり、Notepad++ で example1.jl を編集しているときに、Alt+Tab で REPL のウィンドウに移動したときに自動的に example1.jl が保存される。
そこで、include("example1.jl"); を実行すればよい。
(何回も実行するのであれば、ctrl+p で前のコマンドを REPL は呼び出せるので、Alt+Tab, ctrl+p, Enter で実行できる。)

4.速い。
Notepad++ を起動してから Julia のファイルの編集に入るまでに1秒以内でできる。自分のところでは、Atom の30倍以上高速である。


とりあえず、Notepad++ を使い込んでみようと思う。



2017年4月25日火曜日

scip には、合わないデータだったみたい

昨日までだいぶかかってインストールした SCIP に早速データを投入して数値実験。
結果としては、CPLEX のほうが性能が良かった。
このあたりは入力データ依存なので、たまたま SCIP に合わないデータなのかもしれない。

分かった点としては、
1.primal bound (実行可能点)が5時間ぐらいたっても1つも得られない。したがって duality gap も出てきていない。
   CPLEX の場合は、数分で実行可能点をとりあえずは発見している。ついでに、dual bound も CPLEX の値のほうがいい。
2.メモリ使用量が劇的に増える。
   分枝の方法があってないのかもしれないけど、とりあえず5時間ぐらいで300万ノードぐらいが生成されていて、100GB のメモリを使用。

このままだと「コンピュータに負担はかかるけど、有益な情報が得られない」ということになりそうなので、計算をストップすることにしてみた。

やはり、入力データにどれだけ性能が出せるかはやってみないと分からないものだ。

2017年4月24日月曜日

SCIP を Linux 上の Matlab から呼び出す、ってお話

SCIP を Linux 上 の Matlab から呼び出すのは、あまり推奨されていないのかもしれないが、コンパイルを頑張って、とりあえずはサンプルプログラムまでは実行できるようになった。(ちなみに、Matlab から呼び出すのでなく単純にコマンドラインから使うだけなら、それほど険しい道のりではない。)

ここにはインストールの記録を残しておく。なお、Ubuntu 系の Linux を想定しているので、ディレクトリ名などは適宜変更が必要かもしれないし、apt-get で必要なファイルなどをインストールする必要もあるかもしれない。
また、インストールするディレクトリは、とりあえず $HOME/lib としている。

1. SoPlex をコンパイル
cd $HOME/lib
wget http://soplex.zib.de/download.php?fname=soplex-3.0.0.tgz
tar xzf soplex-3.0.0.tgz
cd soplex-3.0.0
make SHARED=true

(SHARED=true は Matlab の mex にリンクするときに必須なので忘れないように。)

2. ZIMPL をコンパイル
cd $HOME/lib
wget http://zimpl.zib.de/download/zimpl-3.3.4.tgz
tar xzf zimpl-3.3.4.tgz
cd zimpl-3.4.4
make SHARED=true

(このあたりの wget は、ライセンスの関係でうまくいかないかもしれないので、ホームページからダウンロードしたほうがいいかもしれない。)

3. IPOPT をコンパイル
cd $HOME/lib
wget  https://www.coin-or.org/download/source/Ipopt/Ipopt-3.12.7.tgz
tar xzf Ipopt-3.12.7.tgz
cd Ipopt-3.12.7
./configure
make

(make install は特に必要なし。)

4. SCIP をコンパイル
cd $HOME/lib
wget http://scip.zib.de/download.php?fname=scip-4.0.0.tgz
tar xzf scip-4.0.0.tgz
cd scip-4.0.0
make SHARED=true

(このあたりのコンパイルで、ZIMPLE や SOPLEX, IPOPT がどこにあるか指定する必要あり。画面上のメッセージは分かりにくいけど、20回ぐらいトライ&エラーをすれば、なんとかなった。)

5. OPTI をコンパイル
cd $HOME/lib
wget https://github.com/jonathancurrie/OPTI/archive/master.zip
unzip master.zip
cd OPTI-master/Solvers/Source/scip

ここで、以下のようなスクリプトでコンパイル。

******************************
#!/bin/sh -x

SCIP_DIR=$HOME/lib/scip-4.0.0
SOPLEX_DIR=$HOME/lib/soplex-3.0.0
OPTI_DIR=$HOME/lib/OPTI-master
IPOPT_DIR=$HOME/lib/Ipopt-3.12.7
mex -largeArrayDims -output scip.mexa64 -I$SCIP_DIR/src -I$SOPLEX_DIR/src -I$OPTI_DIR/Solvers/Source/scip/Include \
    -I$IPOPT_DIR/Ipopt/src/Common -I$OPTI_DIR/Solvers/Source/opti\
    -L$SCIP_DIR/lib/shared \
    scipmex.cpp scipeventmex.cpp scipnlmex.cpp -lscip -lnlpi.cppad -lut
******************************

基本的には、リンクで見つからない関数があれば、そのつど nm コマンドで *.so を全解列挙して、grep で検索すればよい。

6. Matlab でパスを追加
addpath(genpath(getenv('HOME'),'/lib/OPTI-master'));

7. Matlab を以下のスクリプトで起動

******************************
export LD_LIBRARY_PATH=/usr/lib/gcc/x86_64-linux-gnu/6.3.0:$HOME/lib/
scip-4.0.0/lib/shared:/usr/lib/x86_64-linux-gnu:$HOME/lib/scip-4.0.0/
lib/shared:$LD_LIBRARY_PATH

LD_PRELOAD=/usr/lib/gcc/x86_64-linux-gnu/6.3.0/libgcc_s.so:/usr/lib/gcc/x86_64-linux-gnu/6.3.0/libstdc++.so:/usr/lib/gcc/x86_64-linux-gnu/6.3.0/libgfortran.so:$HOME/lib/scip-4.0.0/lib/shared/libscip.so:/usr/lib/x86_64-linux-gnu/libreadline.so:$HOME/lib/scip-4.0.0/lib/shared/libnlpi.cppad.so:$HOME/lib/scip-4.0.0/lib/shared/libscipsolver.so  /usr/local/bin/matlab -nodisplay $*
******************************

実行時にうまく行かない場合も、nm と grep で探し当てて、LD_LIBRARY_PATH と LD_PRELOAD に載せれば何とかなる。
(なぜか LD_PRELOAD だけでは、うまく行かない時があるので、LD_LIBRARY_PATH にもフォルダをあわせて書く方がよいようだ。)

7. サンプルソースの実行
あとは、$HOME/lib/OPTI-master/Examples/Global_NonlinearProgramming.m にある SCIP 用のサンプルを実行すればよい。
具体的には、Example 4 - White Box Quartic となっているところである。



途中途中であまりスリムになっていないところもあるが、これでサンプルプログラムを実行できる。
Matlab から $ help scip とすれば使い方も分かるようになっている。

たぶん、Mac でも同じような手順でコンパイルできるかと思うが、Mac は OS 自体が複雑なので、もっと大変なのかもしれない。


2017年4月20日木曜日

SCIP がインストールできずに、てこずる話

どうやら SCIP というのが性能がいいらしい、という話を今更ながらに聞いてきたので、今日インストールしようとし始めた。

扱いたいデータは Linux サーバーのMatlab 上にあるので、SCIP を Linux の Matlab から呼び出すようにしないといけない。
どうやら、YALMIP を使って OPTI 経由で SCIP を呼び出すのが流儀のようなので、インストールしてみる。

今日のところインストールしたのは、
SCIP
Soplex
ZIMPL
YALMIP
OPTI
Ipopt
まで。

で、OPTI は Windows 版はあるけど Linux 版がないので、scip のインターフェイスを直接 mex でコンパイルしようと模索。
すると、なんと libscip.a が -fPIC がついていなかった。

明日のうちに SCIP のコンパイルからやり直さないと。。。

できれば、SCIP が apt-get でインストールできると便利ではあるけど、apt-get でインストールできるようになるのは、現状では難しいと思う。

2017年4月12日水曜日

救急車再配置の実例

今朝のニュースで、リアルタイムでの救急車再配置を実際に行っている例が取り上げられていた。
前に調べた時には、日本では実例が見つからなかったので、大変参考になった。
実際にどのように行っているのか、情報がどこかにあるか調べてみようと思う。


2017年4月7日金曜日

IT Text 数理最適化

世の中には、非常に不思議だと思うこともあって、「IT Text 数理最適化」の Amazon のページにレビューが1件もなかったりする。

数理最適化の基礎を勉強するうえでとてもいい本だと思うし、自分もちょくちょく参考にしていたりするし、他の人から「数理最適化って、どの本で勉強したらいいですか?」って聞かれたら「この本がおススメ」って伝えていたりもするし。

なんでレビューが1件もないのか謎だ。不思議だ。



2017年3月31日金曜日

TeX で使われていない数式番号を見つけ出す

TeX で \label をつけて数式番号を振ってあるけど、その数式番号をどこでも使ってなかった、ということが、よくある。
とりあえずは、\label がついているものすべてに \ref が使用されているかどうかを見ればわかるので、Linux なら次のようなスクリプトで \label と \ref の diff を取ることができる。

(ただし、paper.tex は適宜変更のこと)
#!/bin/sh
grep -o -e "ref{[^}]*}" paper.tex | sort | uniq  | sed "s+ref{++" | sed "s+}++" > __ref.txt
grep -o -e "label{[^}]*}" paper.tex | sort | uniq  | sed "s+label{++" | sed "s+}++" > __label.txt

diff __ref.txt __label.txt
rm __ref.txt __label.txt

2017年3月3日金曜日

Optimality condition and complexity analysis for linearly-constrained optimization without differentiability on the boundary

Optimization Online で

Optimality condition and complexity analysis for linearly-constrained optimization without differentiability on the boundary
Gabriel Haeser, Hongcheng Liu, Yinyu Ye
http://www.optimization-online.org/DB_HTML/2017/02/5861.html

というのが掲載されていたので、チェックをしてみた。
基本的には、KKT 条件を epsilon 誤差だけ許容するような条件に拡張した場合にどうなるか、という点が議論されていた。

実際の数値計算では、厳密に(誤差0で)KKT条件を満たせることは、ほぼないので、epsilon-KKT 条件のようなものが計算されているはずで、その理論展開が普通の KKT 条件と似たようなところで展開されている、というのが面白い。
つまり、「KKT 条件が成り立たないはずの数理最適化問題が、KKT条件を想定して作られたソフトウェアでも、なぜか答えが求まっちゃう」ということが実際には良くあるわけだが(例えば、Slater 条件を満たさない SDP の解を SDPA が求めちゃったり)、こういう研究を見てみると、そういう現象の背景が推察できたりする。

2017年2月28日火曜日

Newton Sketch

SIAM Optimization に「Newton Sketch」というタイトルの論文が出ていたので、チェックをしてみたが、なかなかに面白かった。

普通に Newton 法を実行すると Hessian $$\nabla^2 f(x)$$ が重たいわけだけど、ある構造を満たすような正定値対象行列 $$S$$ を用いて $$S \nabla^2 f(x)^{1/2}$$ のように変更する。
ここで、$$SS$$ は確率的に選んでいて、しかも行列積などが高速に行える構造を持っているようなものを作っている。

このようにすると、stochastic gradient よりも速い収束が得られる、という話だった。

たしかに、stochastic gradient は反復回数が非常に多くて十分な収束が得られないので、こういったように改善する、というのは面白いアイデアだと思う。

論文では、self-concordant なども議論していて、SDP なども解けるとのこと。
ただ、数値実験では SDP が出ていないので、たぶん SDP の構造にはあまりあてはまっていない可能性もある。
そういった可能性もあるが、SDP でどうやって性能を引き出せるか、ということに注目すれば、新しいテーマにできるかもしれない。



2017年1月16日月曜日

みかん最適化分割問題

今朝の NHK で、「みかんの皮を切っていって、酉の形ができあがる」というスゴ技をやっている人がいました。

「あたらしいみかんのむきかた」という本(Amazon の「あたらしいみかんのむきかた」)にはみかんの皮のむき方が書いてあるようなので、実際の絵などは、そちらを見ると分かりやすいかと思います。

これがどうすごいか、というのを数学的な視点からいうと

1.みかんの皮に余りがない(つまり、酉の形になっている皮を戻すと球形に戻る)
2.酉の形は連結である
3.みかんのヘタの部分が酉の目の位置に来ている

の3点が挙げられます。
数理最適化で考えてみると、「3次元多様体から連続的に変化させて得られる2次元多様体で、もっともイメージに近いような変化を求める」というようなことになるかと思われるわけです。

この場合の計算方法としては、次のうちのどちらかなぁ、と思います。
(1)多様体から多様体への写像を求める
(2)皮を離散化して、座標ごとの連結情報を維持しながら、各座標の写像を計算する


うーん、こういったのを頭の中で計算できてしまうとは、世の中には色んな達人がいるものです。


2017年1月13日金曜日

Julia の型変換は、なかなかに難しい (SDPA を Julia で実装する:その3)

Julia では、function の引数に型指定できるので、間違った入力が来ても、引数チェックの時点ではじくことができる。

たとえば、
function func1(x::Int, y::Float64)
とすれば、x は Int, y は Float64 が指定できる。
具体的には、
func1(x::Int, y::Float64) = x+y
等と定義すれば、func1(3, 4.5) で計算できる。


ここで、Int の 3 と Float64 の 4.5 の足し算をすると、Float64 の 7.5 が帰ってくる。
そこで、y のほうに整数を入れても自動的に変換できるだろう、と思うと、実はそうではない。
func1(3,4) とすると、「型が合わない」というエラーが出る。

特に難しいのはベクトルの定義で、例えば
a = [-2; 3; 1]
とすると、Array{Float64} では型が合わないので、基本的には Array{Int}, Array{Float64} のそれぞれで関数を作るか、あるいは float で変換するか、というあたりである。

このあたりを、どうすると簡単に SDPA を使えるようにできるか、検討している。




2017年1月11日水曜日

SDPA を Julia で実装してみる(その2)

SDPA sparse format の読み込みの続きを行った。
Matlab と Julia では、文字列処理が異なるので、単純にコピーするわけにもいかず、このあたりが結構手間がかかる。

Julia では、replace で文字列の置換、split で文字列分解ができて、このあたりは強力な操作ができる。

あと、疎行列については Matlab であればメモリー領域を先に確保しておくことができるが、Julia では同じような機能がない様子。
今は小さいファイルしか読み込んでいないので特に問題ないが、大きい問題になると読み込みに時間がかかりそうだ。

2017年1月10日火曜日

SDPA を Julia で実装してみる(その1)

最近になって、SDPA を Julia で実装しなおしている。
主な理由としては、SDPA のメンテナンスコストを抑えるため。

やはり、SDPAのインストールなどを考えると、コンピュータに関する知識がある程度は必要だし、かといって、OpenBLAS などを毎回コンパイルするのも大変。
そこで、プログラミング言語としても、それなりにスピードがある Julia で実装を行ってみている。


ということで、まずは最初のところとして、SDPA sparse format を Julia に読み込むようにしている。
SDPA では、fscanf の文字処理が強力なので、それに任せているが、Julia だと fscanf を使えないので(libcを経由すれば使えるけど、それだとlibcがないところで使えない)、正規表現で読み込むようにしている。

とりあえずファイルからは読み込めたので、あとは疎行列形式に持たせるようにすれば、OK になるはず。
Matlab に近いインターフェイスにできるだけしたいけど、Matlab だとスピードが出ないところもあるので、そのあたりは適宜修正になりそう。