2015-12-20

これは，Julia Advent Calendar 20日目の記事です．この記事では，JuliaOptについて取り上げてみようと思います．

JuliaOptとは

JuliaOptとは，Juliaで数理最適化を実行するためのパッケージを集めたプロジェクトです*1．
JuliaOpt: Optimization packages for the Julia language.
http://www.juliaopt.org/

JuliaOptは，Juliaの中でも最も精力的に開発が進められているプロジェクトの一つと言っても過言ではありません．
上記の公式ページは解説や関連する論文などが充実しています．
後ほど説明するモデリングツールであるJuMPやConvex.jlのサイトにもドキュメントが豊富です．
JuMP: https://jump.readthedocs.org/en/latest/
Convex.jl: http://convexjl.readthedocs.org/en/latest/

さらに，Jupyter Notebook形式で豊富な実行例が提供されています．
http://nbviewer.ipython.org/github/JuliaOpt/juliaopt-notebooks/tree/master/notebooks/

JuliaCon2015でもプレゼンテーションやワークショップが行われており，githubには各種の資料が，Youtubeには動画が上がっています．

スライド等の資料
https://github.com/JuliaOpt/JuliaCon-2015
JuliaOpt --Optimization packages in Julia

プレゼンテーション
"JuliaOpt Optimization related projects in Julia"

ワークショップ
"Solving optimization problems with JuliaOpt"

JuliaOptの概観

JuliaOptは複数のパッケージから構成されます．公式ページの中段付近にあるパッケージの関係図がわかりやすいのでご参照ください．
http://www.juliaopt.org/
これらのパッケージは，最適化問題を数式で定式化するモデリングツールと実際に最適化問題を解くソルバーに大別されます．
また，モデリングツールとソルバーをつなぐ抽象化のレイヤーにはMathProgBase.jlが配されています．

モデリングツール

モデリングツールは，最適化問題を数式で定式化します．JuMPとConvex.jlの2つが提供されています．

JuMP
線形計画法，2次計画法，非線形制約最適化問題のための代数モデリング言語．商用のモデリングツールと同程度の速度でモデルを生成するとともに，ソルバーのコールバックなど高度な機能も提供する．
Convex.jl
凸プログラミングのための代数モデリング言語．

それぞれのモデリングツールが対応している問題，解法は下表のとおりです．
JuMPが線形計画法，混合整数線形計画法，2次錐計画法，混合整数2次錐計画法，半正定値計画法，非線形計画法，混合整数非線形計画法に対応しているのに対して，Convex.jlは非線形計画法，混合整数非線形計画法には対応していません．

	線形計画法	混合整数線形計画法	2次錐計画法	混合整数2次錐計画法	半正定値計画法	(制約付き)非線形計画法	混合整数非線形計画法
JuMP	○	○	○	○	○	○	○
Convex.jl	○	○	○	○	○	--	--

ソルバー

ソルバーは，Gurobi, CPLEX, Mosek等の商用ツール，SCS, ECOS等のフリーツールなど多くのものが対応しています．
詳しくは，JuliaOptの公式ページを参照してください．

JuMPの使用例

JuMPの使用例について，例題を通して見ていきましょう．以下の実行環境はOS X 10.10，Julia 0.4.2です．
まずは簡単な例として，JuMPのページにも例題として載っている以下の問題を解いてみます．
$\begin{eqnarray}\max \, \, &5x + 3y \\ s.t. 0 \leq &x \leq 2, \, 0 \leq &y \leq 30, x + 5y &\leq 3\end{eqnarray}$

この問題は，JuMPを用いて次のように記述します．線形計画法を解くソルバーとして，ここではClp.jlをインストールしておきます．

julia> Pkg.add("JuMP")
julia> Pkg.add("Clp")
julia> using JuMP
julia> # モデルのインスタンスの生成
julia> m = Model()
Feasibility problem with:
 * 0 linear constraints
 * 0 variables
Solver set to Default
julia> # 変数の定義
julia> @defVar(m, 0 <= x <= 2)
x
julia> @defVar(m, 0 <= y <= 30)
y
julia> # 目的関数の設定
julia> @setObjective(m, Max, 5x + 3y)
5 x + 3 y
julia> # 制約条件の追加
julia> @addConstraint(m, x + 5y <= 3.0)
x + 5 y &#8804; 3
julia> # モデルの確認
julia> m
Maximization problem with:
 * 1 linear constraint
 * 2 variables
Solver set to Default
julia> # 最適化の実行
julia> status = solve(m)
INFO: Precompiling module Clp...
:Optimal
julia> println("Objective value: ", getObjectiveValue(m))
Objective value: 10.6
julia> println("x= ", getValue(x))
x= 2.0
julia> println("y= ", getValue(y))
y= 0.2

要点は以下のとおりです．

Model()でモデルのインスタンスを生成．
@defVarで変数を定義する．第一引数にはモデルのインスタンス，第二引数には変数の範囲を指定する．
@setObjectiveで目的関数を設定する．第一引数にはモデルのインスタンス，第二引数には目的関数の最大(Max) or 最小(Min)，第三引数には目的関数の式を指定する．
@addConstraintで制約条件を追加する．第一引数にはモデルのインスタンス，第二引数には制約条件の式を指定する．
solve関数で最適化を実行する．
最適化の実行後，目的関数の値はgetObjectiveValue関数，各変数の値はgetValue関数で取得する．

続いて，非線形最適化の例として，Rosenbeck関数 $(1-x)^2 + 100(y - x^2)^2$ を最小化する $x, y$ を求めてみましょう．非線形最適化を実行するIpopt.jlをインストールしておきます．

julia> Pkg.add("Ipopt")
julia> using JuMP
julia> m = Model()
Feasibility problem with:
 * 0 linear constraints
 * 0 variables
Solver set to Default
julia> @defVar(m, x, start = 0.0)
julia> @defVar(m, y, start = 0.0)
julia> @setNLObjective(m, Min, (1-x)^2 + 100(y-x^2)^2)
julia> solve(m)
julia> println("x = ", getValue(x), " y = ", getValue(y))
x = 0.9999999999999899 y = 0.9999999999999792

Convex.jlの使用例

Convex.jlの使用例として，次の問題を解いてみましょう．
$\begin{eqnarray} \min \, \, &\||Ax - b\||^2 \\ s.t. x_{i} &\geq 0 \, (i = 1, \, \dots, \, N) \end{eqnarray}$
以下では， $A$ は各成分が一様乱数からなる $4$ 行 $5$ 列の行列， $b$ が4つの一様乱数からなるベクトルとして，Convex.jlを用いて最適なベクトル $x$ を求めています．

julia> Pkg.add("convex")
julia> Pkg.add("SCS")
julia> using Convex, SCS
julia> set_default_solver(SCSSolver(verbose=0))
SCS.SCSSolver(Any[(:verbose,0)])
julia> srand(71)
julia> m = 4;  n = 5
julia> # 行列Aとベクトルbの設定
julia> A = randn(m, n); b = randn(m, 1)
julia> # 最適化する変数 (ベクトルx)
julia> x = Variable(n)
julia> # 最適化問題の定式化
julia> problem = minimize(sumsquares(A * x - b), x >= 0)
julia> # 最適化の実行
julia> solve!(problem)
julia> # ベクトルxの推定値
julia> x
Variable of
size: (5, 1)
sign: Convex.NoSign()
vexity: Convex.AffineVexity()
value: 5x1 Array{Float64,2}:
 -2.8087e-6 
  0.26101   
  1.99817   
 -1.21593e-6
  0.135552
julia> # 目的関数の値
julia> problem.optval
1.166663382501328