この記事は 数理最適化Advent Calendar 2023 の8日目の記事です。
はじめに
数理最適化Advent Calendar 2023の 1日目の記事 では、書籍「Pythonではじめる数理最適化」を共著で書かせていただいた @iwanaga-jiro さんより、書籍7章の内容を題材にした数理最適化問題の行列表現の実装について紹介されていました。
7章の「推薦商品のための興味のスコアリング」の題材は学生時代に @iwanaga-jiro さんと一緒に議論させていただいたテーマ1で、自身がデータ解析の仕事をはじめるきっかけとなりました。このモデリングには、観測データを絶対とするのでなくデータに期待する特性を柔軟に制約条件として組み込むことができるという面白さがあります。一方、書籍を手に取ってもらった知人からは行列表現が難しい...という声もききました。
そこで本記事では、Pythonのモデリング言語CVXPYを利用することで行列表現を用いず、なるべく数式そのままで凸2次計画問題とよばれるタイプの最適化問題2を実装する例を紹介します。
最後に、CVXPYのドキュメントを改めて読んで知ったRustとJuliaで実装されている凸最適化ソルバー Clarabel についても簡単に触れます。
CVXPYの特徴
CVXPYの公式ドキュメント では、いわゆる数理最適化の標準形3とよばれる形に変換することなく自然に数式の形で記述できることがCVXPYの特徴としてあげられています。
CVXPY is an open source Python-embedded modeling language for convex optimization problems. It lets you express your problem in a natural way that follows the math, rather than in the restrictive standard form required by solvers.
また、モデリング言語であるCVXPYを介してCVXOPTも含む複数のソルバーを呼び出すことができることも利点だとおもいます。
CVXPY relies on the open source solvers Clarabel, OSQP, SCS, and ECOS. Additional solvers are supported, but must be installed separately.
You can add solver names as “extras”; pip will then install the necessary additional Python packages.
pip install cvxpy [CBC,CVXOPT,GLOP,GLPK,GUROBI,MOSEK,PDLP,SCIP,XPRESS]
CVXPYを利用した実装
CVXPYのサンプルコードついては公式ドキュメントや日本語の解説記事もありますが、凸2次計画問題で目的関数でも制約条件でも行列表現を使っていない素朴な実装例は見当たらなかったので4、なるべく簡単な形で実装していきます。
本記事でも 1日目の記事 にならって、次の3つの決定変数$x_0,x_1,x_2$をもつ簡単な最小化問題を扱います5。（Pythonでの実装にあわせて添字は0はじまりにしています）
扱う最小化問題
-
決定変数の範囲
$ 0 \leq x_0 \leq 1 $
$ 0 \leq x_1 \leq 1 $
$ 0 \leq x_2 \leq 1 $ -
制約1（単調増加）
$ x_0 \leq x_1 $
$ x_1 \leq x_2 $ -
制約2（上に凸 → 増分が減少）
$ x_1 - x_0 \geq x_2 - x_1 $ -
目的関数（誤差最小化）
$ 100 \times (x_0 - 0.1)^2 + 50 \times (x_1 - 0.15)^2 + 10 \times (x_2 - 0.25)^2 $
CVXPYを利用すると「行列表現に持ち込むための二次計画問題の作成」や「二次計画問題の行列表現」は不要で直接実装をはじめることができます。
行列表現を利用しない凸2次計画問題の実装例
まず準備として、CVXPYのインストールと決定変数の定義をします。 現時点(2023/12/8)ではCVXPY1.4系最新なので以下のようにインストールします。
! pip install cvxpy==1.4
import cvxpy as cp
# 決定変数の定義
x = cp.Variable(3)
続いて、制約条件と目的変数をそれぞれ実装していきます。
-
決定変数の範囲
$ 0 \leq x_0 \leq 1 $
$ 0 \leq x_1 \leq 1 $
$ 0 \leq x_2 \leq 1 $
# 変数の範囲
constraints = []
constraints += [0 <= x[0], x[0] <= 1]
constraints += [0 <= x[1], x[1] <= 1]
constraints += [0 <= x[2], x[2] <= 1]
-
制約1（単調増加）
$ x_0 \leq x_1 $
$ x_1 \leq x_2 $
# 制約1（単調増加）
constraints += [x[0] <= x[1]]
constraints += [x[1] <= x[2]]
-
制約2（上に凸 → 増分が減少）
$ x_1 - x_0 \geq x_2 - x_1 $
# 制約2（上に凸、増分が減少）
constraints += [x[1] - x[0] >= x[2] - x[1]]
-
目的関数（誤差最小化）
$ 100 \times (x_0 - 0.1)^2 + 50 \times (x_1 - 0.15)^2 + 10 \times (x_2 - 0.25)^2 $
# 目的関数（誤差最小化）
objective = 100*(x[0] - 0.1)**2 + 50*(x[1] - 0.15)**2 + 10*(x[2] - 0.25)**2
実装した目的関数と制約条件を引数として解く問題を定義します。
# 最小化はcp.Minimize()。cp.Problem()の引数として目的関数と制約条件を与えて問題を定義する
prob = cp.Problem(cp.Minimize(objective), constraints)
solve()で問題を解きます。引数を空にすると目的関数値だけが返却されますが、引数として verbose=1 とすると詳細な実行結果が返ってきます。
# verbose=1で詳細を確認
prob.solve(verbose=1)
Google colabで実行した結果です。デフォルトの設定でOSQPというソルバーが利用されているようです。実行時間 (run time) は 3.63e-04s 、目的関数値 (optimal objective) は 0.0132 と出力されています。
===============================================================================
CVXPY
v1.4.0
===============================================================================
(CVXPY) Dec 07 12:47:13 PM: Your problem has 3 variables, 9 constraints, and 0 parameters.
(CVXPY) Dec 07 12:47:13 PM: It is compliant with the following grammars: DCP, DQCP
(CVXPY) Dec 07 12:47:13 PM: (If you need to solve this problem multiple times, but with different data, consider using parameters.)
(CVXPY) Dec 07 12:47:13 PM: CVXPY will first compile your problem; then, it will invoke a numerical solver to obtain a solution.
(CVXPY) Dec 07 12:47:13 PM: Your problem is compiled with the CPP canonicalization backend.
-------------------------------------------------------------------------------
Compilation
-------------------------------------------------------------------------------
(CVXPY) Dec 07 12:47:13 PM: Compiling problem (target solver=OSQP).
(CVXPY) Dec 07 12:47:13 PM: Reduction chain: CvxAttr2Constr -> Qp2SymbolicQp -> QpMatrixStuffing -> OSQP
(CVXPY) Dec 07 12:47:13 PM: Applying reduction CvxAttr2Constr
(CVXPY) Dec 07 12:47:13 PM: Applying reduction Qp2SymbolicQp
(CVXPY) Dec 07 12:47:13 PM: Applying reduction QpMatrixStuffing
(CVXPY) Dec 07 12:47:13 PM: Applying reduction OSQP
(CVXPY) Dec 07 12:47:13 PM: Finished problem compilation (took 6.354e-02 seconds).
-------------------------------------------------------------------------------
Numerical solver
-------------------------------------------------------------------------------
(CVXPY) Dec 07 12:47:13 PM: Invoking solver OSQP to obtain a solution.
-----------------------------------------------------------------
OSQP v0.6.2 - Operator Splitting QP Solver
(c) Bartolomeo Stellato, Goran Banjac
University of Oxford - Stanford University 2021
-----------------------------------------------------------------
problem: variables n = 6, constraints m = 12
nnz(P) + nnz(A) = 22
settings: linear system solver = qdldl,
eps_abs = 1.0e-05, eps_rel = 1.0e-05,
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 = 10000
check_termination: on (interval 25),
scaling: on, scaled_termination: off
warm start: on, polish: on, time_limit: off
iter objective pri res dua res rho time
1 0.0000e+00 2.50e-01 3.00e+01 1.00e-01 1.19e-04s
150 1.3154e-02 7.86e-06 8.48e-07 1.11e+00 2.78e-04s
plsh 1.3158e-02 0.00e+00 6.21e-17 -------- 3.63e-04s
status: solved
solution polish: successful
number of iterations: 150
optimal objective: 0.0132
run time: 3.63e-04s
optimal rho estimate: 1.97e+00
-------------------------------------------------------------------------------
Summary
-------------------------------------------------------------------------------
(CVXPY) Dec 07 12:47:13 PM: Problem status: optimal
(CVXPY) Dec 07 12:47:13 PM: Optimal value: 1.316e-02
(CVXPY) Dec 07 12:47:13 PM: Compilation took 6.354e-02 seconds
(CVXPY) Dec 07 12:47:13 PM: Solver (including time spent in interface) took 4.165e-03 seconds
0.013157894736842092
$x_0, x_1, x_2$の最適解は x.value とすると、arrayとして得られます。
x.value
> array([0.09736842, 0.16052632, 0.22368421])
上の例では簡単さのため冗長に記述しましたが、決定変数のベクトル x に対して一括で制約条件を設定したり、総和の関数 cp.sum を利用するとより簡潔な形で実装できます。
constraints = []
# 変数に対して一括で制約条件を指定
constraints += [0 <= x] + [x <= 1]
for i in range(2):
constraints += [x[i] <= x[i+1]]
for i in range(1):
constraints += [x[i+1] - x[i] >= x[i+2] - x[i+1]]
w = [100, 50, 10]
p = [0.1, 0.15, 0.25]
# cp.sum でリストの総和
objective = cp.sum([w[i]*(x[i] - p[i])**2 for i in range(3)])
prob = cp.Problem(cp.Minimize(objective), constraints)
prob.solve()
sum の他にも、自然に数式の形で記述するための sqrt, pow, log_sum_exp などの便利な関数が API Documantation に整理されています。
CVXPYでの行列表現を利用した凸2次計画問題の実装例
ここまでは行列表現を利用しないで実装するということを目的に進めてきました。しかしCVXPYでは行列表現が利用できないかというとそうではなく、公式ドキュメントのExamplesでも多くの例で行列表現が利用されています。
例えば、上記の簡単な例は1日目の記事にならって行列表現すると以下のように実装できます。
import cvxpy as cp
import numpy as np
x = cp.Variable(3)
G = np.array([[-1, 0, 0]
,[ 1, 0, 0]
,[ 0,-1, 0]
,[ 0, 1, 0]
,[ 0, 0,-1]
,[ 0, 0, 1]
,[ 1,-1, 0]
,[ 0, 1,-1]
,[ 1,-2, 1]])
h = np.array([0
,1
,0
,1
,0
,1
,0
,0
,0])
P = np.array([[2*100 ,0 ,0 ]
,[0 ,2*50 ,0 ]
,[0 ,0 ,2*10]])
q = np.array([-2 * 100 * 0.1
,-2 * 50 * 0.15
,-2 * 10 * 0.25])
objective = cp.Minimize((1/2)*cp.quad_form(x, P) + q.T @ x)
constraints = [G @ x <= h]
prob = cp.Problem(objective,constraints)
prob.solve()
こちらは標準形ではあるものの、目的関数が $Minimize~(\frac{1}{2})x^TPx + q^Tx$、制約条件が $Gx \leq h $ と、数式のまま自然に実装することができています。
余談：Rust/Julia実装の凸最適化ソルバーClarabel
本記事を執筆するにあたって改めてCVXPYのドキュメントを読んでいたところ、 CVXPY1.4 からRustとJuliaで実装されている Clarabel というソルバーが利用可能になっているようです。
もともとCVXPYが OSQP , SCS と並んでデフォルトでインストールしていた ECOS がエッジケースにて性能が不安定となる点を解消しているそうで、CVXPY1.5ではいくつかの問題タイプでECOSでなくClarabelがデフォルトで利用されることが計画されているそうです。
ECOS deprecation
CVXPY has used ECOS as the default solver for many years; however, it has known issues with performance and numerical stability in edge cases. Recently, a new solver, Clarabel, that improves the algorithm and implementation of ECOS has been under development.
In 1.5, CVXPY plans to start using Clarabel instead of ECOS by default for some categories of problems. In 1.6, we plan to no longer install ECOS as a CVXPY dependency. We have no plans to remove support for calling ECOS as a solver. As part of this transition, in 1.4 CVXPY will raise a warning whenever ECOS is called by default. We encourage you to try and use Clarabel instead, but if you’re dependent on ECOS’s exact behavior please explicitly specify it as a solver.
上述の問題に solver=cp.CLARABEL とソルバーを指定して実行すると、確かにClarabelが利用できるようでした。
prob.solve(solver=cp.CLARABEL, verbose=True)
===============================================================================
CVXPY
v1.4.0
===============================================================================
(CVXPY) Dec 07 12:47:13 PM: Your problem has 3 variables, 1 constraints, and 0 parameters.
(CVXPY) Dec 07 12:47:13 PM: It is compliant with the following grammars: DCP, DQCP
(CVXPY) Dec 07 12:47:13 PM: (If you need to solve this problem multiple times, but with different data, consider using parameters.)
(CVXPY) Dec 07 12:47:13 PM: CVXPY will first compile your problem; then, it will invoke a numerical solver to obtain a solution.
(CVXPY) Dec 07 12:47:13 PM: Your problem is compiled with the CPP canonicalization backend.
-------------------------------------------------------------------------------
Compilation
-------------------------------------------------------------------------------
(CVXPY) Dec 07 12:47:13 PM: Compiling problem (target solver=CLARABEL).
(CVXPY) Dec 07 12:47:13 PM: Reduction chain: Dcp2Cone -> CvxAttr2Constr -> ConeMatrixStuffing -> CLARABEL
(CVXPY) Dec 07 12:47:13 PM: Applying reduction Dcp2Cone
(CVXPY) Dec 07 12:47:13 PM: Applying reduction CvxAttr2Constr
(CVXPY) Dec 07 12:47:13 PM: Applying reduction ConeMatrixStuffing
(CVXPY) Dec 07 12:47:13 PM: Applying reduction CLARABEL
(CVXPY) Dec 07 12:47:13 PM: Finished problem compilation (took 3.399e-02 seconds).
-------------------------------------------------------------------------------
Numerical solver
-------------------------------------------------------------------------------
(CVXPY) Dec 07 12:47:13 PM: Invoking solver CLARABEL to obtain a solution.
-------------------------------------------------------------------------------
Summary
-------------------------------------------------------------------------------
(CVXPY) Dec 07 12:47:13 PM: Problem status: optimal
(CVXPY) Dec 07 12:47:13 PM: Optimal value: -2.737e+00
(CVXPY) Dec 07 12:47:13 PM: Compilation took 3.399e-02 seconds
(CVXPY) Dec 07 12:47:13 PM: Solver (including time spent in interface) took 5.739e-04 seconds
array([0.09736842, 0.16052632, 0.22368421])
おわりに
本記事では、書籍「Pythonではじめる数理最適化」で登場する簡単な凸2次計画問題の例を題材に、CVXPYを利用して行列表現を利用せず実装する方法を紹介しました。紹介した実装はGoogle Colab上にも公開しています。
またこの機会に、書籍のGitHub上のサポートページに7章の「推薦商品のための興味のスコアリング」についてCVXPYで実装したものをアップロードしましたので、別の例での実装に興味をもっていただけた方は是非ご確認ください。
参考記事
- CVXPYの公式ドキュメント
- Clarabelの公式ドキュメント
- 凸2次計画問題の標準形や応用例など
- CVXOPTを利用した簡単な凸2次計画問題の実装の解説
- 日本語のCVXPYに関する先行記事