Juliaで2次計画法

この記事はJulia Advent Calender 2018 5日目の記事です。

2019-11-23 追記

JuMP v0.19.0以降のAPI変更を反映しました。

にて動作を確認しています。

Juliaで2次計画法

過去にPythonでOpenOpt, FuncDeisignerを使って2次計画法という記事を書いたのですが、今回はそのJuliaバージョンです。

2次計画法といえば金融分野のポートフォリオ最適化などをはじめ応用範囲が広い最適化問題の一つですが、線形計画法ほどは容易に解けず、オープンソースのソルバーでも線形計画法に比べると有名なソルバーが少ない印象です。
Juliaが様々な分野で広まっていくにつれてこうした問題を解く機会も増える（はず）なので、この記事が一助になれば幸いです。

今回もPythonのときと同様に以下のような問題を解くことを考えます。

\(
\begin{align*}
\min ~ &\frac{1}{2}\mathbf{x}^TH\mathbf{x} + \mathbf{f}^T x \\
\text{subject to} ~ ~
&\mathbf{l} \le \mathbf{x} \le \mathbf{u} \\
&A\mathbf{x} \le \mathbf{b} \\
&A_{eq}\mathbf{x} = \mathbf{b}_{eq}
\end{align*}
\)

例題としては

\(
\begin{align*}
\min ~ &2x_1^2 + x_2^2 +x_1x_2 + x_1 + x_2 \\
\text{subject to} ~ ~
&x_1 \ge 0 \\
&x_2 \ge 0 \\
&x_1 + x_2 = 1 \\
\end{align*}
\)

を実際に解いていきます。

パッケージとしてはJuMP及びOSQPを使います。JuMPは最適化問題のモデリングツールで、OSQPは特に二次計画法に特化したオープンソースのソルバーです。
OSQPはJuliaラッパーもあり、JuMPから容易に呼び出すことができます。（JuMPがOSQPに対応したのは2018年の2月から）

いつものごとく

(v1.2) pkg> add JuMP
(v1.2) pkg> add OSQP

でパッケージをインストールします。

スカラー変数の形式での記述

まず、例題の形をそのまま直訳（？）したものを書いてみます。

using JuMP, OSQP

model = Model(with_optimizer(OSQP.Optimizer))

# ソルバーの出力が煩わしいときは以下のset_silent()関数を用いる
# set_silent(model)

@variable(model, x₁ >= 0)
@variable(model, x₂ >= 0)

@constraint(model, x₁ + x₂ == 1)
@objective(model, Min, 2x₁^2 + x₂^2 + x₁*x₂ + x₁ + x₂)

print(model)  # 最適化問題の式を出力

optimize!(model)  # 計算結果やパラメータの表示

println("Objective value: ", JuMP.objective_value(model))
println("x₁ = ", JuMP.value(x₁))
println("x₂ = ", JuMP.value(x₂))

Juliaのマクロの強力さが如実に現れているかと思います。とても簡素に最適化問題を記述することができます。
このプログラムの実行結果は以下のようになります。

$ julia quadprog.jl
Min 2 x₁² + x₂² + x₁*x₂ + x₁ + x₂
Subject to
 x₁ + x₂ = 1.0
 x₁ ≥ 0.0
 x₂ ≥ 0.0
-----------------------------------------------------------------
           OSQP v0.6.0  -  Operator Splitting QP Solver
              (c) Bartolomeo Stellato,  Goran Banjac
        University of Oxford  -  Stanford University 2019
-----------------------------------------------------------------
problem:  variables n = 2, constraints m = 3
          nnz(P) + nnz(A) = 7
settings: linear system solver = qdldl,
          eps_abs = 1.0e-03, eps_rel = 1.0e-03,
          eps_prim_inf = 1.0e-04, eps_dual_inf = 1.0e-04,
          rho = 1.00e-01 (adaptive),
          sigma = 1.00e-06, alpha = 1.60, max_iter = 4000
          check_termination: on (interval 25),
          scaling: on, scaled_termination: off
          warm start: on, polish: off, time_limit: off

iter   objective    pri res    dua res    rho        time
   1  -7.8808e-03   1.01e+00   2.00e+02   1.00e-01   9.50e-05s
  25   1.8745e+00   1.75e-04   3.32e-04   1.00e-01   1.31e-04s

status:               solved
number of iterations: 25
optimal objective:    1.8745
run time:             1.34e-04s
optimal rho estimate: 1.31e-01

Objective value: 1.8745175014907984
x₁ = 0.24994319406821397
x₂ = 0.7498813420106132

解析解が[0.25, 0.75]なのでだいたい合ってそうです。

ベクトルによる記述

今の例では変数をそれぞれx₁とx₂に分けて記述しましたが、一般的には変数はベクトルで表すのが普通です。
もちろんJuMPでもベクトルを用いてモデリングすることが可能です。

using JuMP, OSQP

model = Model(with_optimizer(OSQP.Optimizer))

# ソルバーの出力が煩わしいときは以下のset_silent()関数を用いる
# set_silent(model)

@variable(model, x[1:2] >= 0)

Aeq = [1, 1]
beq = 1
@constraint(model, Aeq'*x == beq)

H = [4 1
     1 2]
f = [1, 1]
@objective(model, Min, 0.5 * x'*H*x + f'*x)

print(model)

optimize!(model)

println("Objective value: ", JuMP.objective_value(model))
println("x = ", JuMP.value.(x))

上のコードの中ではベクトルfをf=[1,1]と定義した後に@objectiveの中でf'*xと転置して使っています。この部分を、先に転置した1x2行列f=[1 1]として定義してから@objectiveの中でf*xと使ってやるとエラーで動かなくなります。
私はこれにだいぶハマってしまいました…。ご注意ください。

コードの実行結果は以下です。
目的関数や制約条件の表記が上の例とはやや異なっているのが分かります。

$ julia quadprog_vec.jl
Min 2 x[1]² + x[1]*x[2] + x[2]² + x[1] + x[2]
Subject to
 x[1] + x[2] = 1.0
 x[1] ≥ 0.0
 x[2] ≥ 0.0
-----------------------------------------------------------------
           OSQP v0.6.0  -  Operator Splitting QP Solver
              (c) Bartolomeo Stellato,  Goran Banjac
        University of Oxford  -  Stanford University 2019
-----------------------------------------------------------------
problem:  variables n = 2, constraints m = 3
          nnz(P) + nnz(A) = 7
settings: linear system solver = qdldl,
          eps_abs = 1.0e-03, eps_rel = 1.0e-03,
          eps_prim_inf = 1.0e-04, eps_dual_inf = 1.0e-04,
          rho = 1.00e-01 (adaptive),
          sigma = 1.00e-06, alpha = 1.60, max_iter = 4000
          check_termination: on (interval 25),
          scaling: on, scaled_termination: off
          warm start: on, polish: off, time_limit: off

iter   objective    pri res    dua res    rho        time
   1  -7.8808e-03   1.01e+00   2.00e+02   1.00e-01   1.12e-04s
  25   1.8745e+00   1.75e-04   3.32e-04   1.00e-01   1.36e-04s

status:               solved
number of iterations: 25
optimal objective:    1.8745
run time:             1.41e-04s
optimal rho estimate: 1.31e-01

Objective value: 1.8745175014907984
x = [0.24994319406821397, 0.7498813420106132]

さいごに

前述の通り、Juliaはマクロが使えるのがとても強力なので、こうしたモデリング用途には大変適しているように感じます。
x₁やx₂のような下付き文字が使えるのも地味に可読性向上に寄与しているような気がします。

速度が注目されがちなJuliaですが、こうしたモデリングの見通しの良さも立派なアピールポイントですね。

参考