CVXPYによる行列表現を利用しない凸2次計画問題の実装

この記事は [数理最適化Advent Calendar 2023](https://qiita.com/advent-calendar/2023/mathematical-optimization) の8日目の記事です。

# はじめに
数理最適化Advent Calendar 2023の [1日目の記事](https://qiita.com/iwanaga-jiro/items/045a04abd47390904e27) では、書籍「Pythonではじめる数理最適化」を共著で書かせていただいた @iwanaga-jiro さんより、書籍7章の内容を題材にした数理最適化問題の行列表現の実装について紹介されていました。

7章の「推薦商品のための興味のスコアリング」の題材は学生時代に @iwanaga-jiro さんと一緒に議論させていただいたテーマ[^1]で、自身がデータ解析の仕事をはじめるきっかけとなりました。このモデリングには、観測データを絶対とするのでなくデータに期待する特性を柔軟に制約条件として組み込むことができるという面白さがあります。一方、書籍を手に取ってもらった知人からは行列表現が難しい...という声もききました。

そこで本記事では、Pythonのモデリング言語CVXPYを利用することで行列表現を用いず、なるべく数式そのままで凸2次計画問題とよばれるタイプの最適化問題[^2]を実装する例を紹介します。

最後に、CVXPYのドキュメントを改めて読んで知ったRustとJuliaで実装されている凸最適化ソルバー [Clarabel](https://oxfordcontrol.github.io/ClarabelDocs/stable/) についても簡単に触れます。 

# CVXPYの特徴
[CVXPYの公式ドキュメント](https://www.cvxpy.org/) では、いわゆる数理最適化の標準形[^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日目の記事](https://qiita.com/iwanaga-jiro/items/045a04abd47390904e27) にならって、次の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系最新なので以下のようにインストールします。
```python
! 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 $

```python
# 変数の範囲
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 $
```python
# 制約1（単調増加）
constraints += [x[0] <= x[1]] 
constraints += [x[1] <= x[2]] 
```
- **制約2（上に凸 → 増分が減少）**
$ x_1 - x_0 \geq x_2 - x_1 $
```python
# 制約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 $
```python
# 目的関数（誤差最小化）
objective = 100*(x[0] - 0.1)**2 + 50*(x[1] - 0.15)**2 + 10*(x[2] - 0.25)**2
```

実装した目的関数と制約条件を引数として解く問題を定義します。
```python
# 最小化はcp.Minimize()。cp.Problem()の引数として目的関数と制約条件を与えて問題を定義する
prob = cp.Problem(cp.Minimize(objective), constraints)
```

solve()で問題を解きます。引数を空にすると目的関数値だけが返却されますが、引数として `verbose=1` とすると詳細な実行結果が返ってきます。
```python
# 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` を利用するとより簡潔な形で実装できます。
```python
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](https://www.cvxpy.org/api_reference/cvxpy.html
) に整理されています。

## CVXPYでの行列表現を利用した凸2次計画問題の実装例
ここまでは行列表現を利用しないで実装するということを目的に進めてきました。しかしCVXPYでは行列表現が利用できないかというとそうではなく、公式ドキュメントの[Examples](https://www.cvxpy.org/examples/index.html)でも多くの例で行列表現が利用されています。

例えば、上記の簡単な例は[1日目の記事](https://qiita.com/iwanaga-jiro/items/045a04abd47390904e27)にならって行列表現すると以下のように実装できます。

```python
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](https://www.cvxpy.org/updates/index.html) からRustとJuliaで実装されている [Clarabel](https://oxfordcontrol.github.io/ClarabelDocs/stable/) というソルバーが利用可能になっているようです。
  
もともとCVXPYが [OSQP](https://osqp.org/) ,  [SCS](https://www.cvxgrp.org/scs/) と並んでデフォルトでインストールしていた [ECOS](https://web.stanford.edu/~boyd/papers/ecos.html) がエッジケースにて性能が不安定となる点を解消しているそうで、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が利用できるようでした。

```python
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上にも公開しています。

https://colab.research.google.com/drive/1KYI011hmRFS8g1xB1XmvfqIrExH8BTjO?authuser=1#scrollTo=S141uzdHtE3U

またこの機会に、書籍のGitHub上のサポートページに7章の「推薦商品のための興味のスコアリング」についてCVXPYで実装したものをアップロードしましたので、別の例での実装に興味をもっていただけた方は是非ご確認ください。

https://github.com/ohmsha/PyOptBook/blob/main/7.recommendation/recommendation_cvxpy.ipynb

# 参考記事
- CVXPYの公式ドキュメント
    - [CVXPY](https://www.cvxpy.org/)
- Clarabelの公式ドキュメント
    - [Clarabel jl/rs](https://oxfordcontrol.github.io/ClarabelDocs/stable/)
- 凸2次計画問題の標準形や応用例など
    - [学習用テキスト 非線形計画法 (1) ２次計画問題](http://www.me.titech.ac.jp/~mizu_lab/text/PDF-NLP/NLP1-QP-problem.pdf)
- CVXOPTを利用した簡単な凸2次計画問題の実装の解説
    - [書籍「Pythonではじめる数理最適化」の改訂版のお知らせ-数理最適化問題の行列表現の実装に入門-](https://qiita.com/iwanaga-jiro/items/045a04abd47390904e27)
- 日本語のCVXPYに関する先行記事
    - [【python cvxpy】モデリング言語cvxpyを利用して数理最適化問題を解く](https://note.com/fz5050/n/n4e3050b54ce3#dce6acdd-50f3-485a-bb7b-da5ad161a9cc)
    - [凸最適化ソルバーCVXPYの紹介](https://qiita.com/birdwatcher/items/885a96898304fddcbabc)
    - [Python による数理最適化モデリングツール CVXPY の初歩](https://mirucacule.hatenablog.com/entry/2019/10/22/163106)
    - [Python製凸最適化モデリングツールCVXPYの使い方](https://myenigma.hatenablog.com/entry/2016/11/06/114603)


[^1]:書籍「Pythonではじめる数理最適化」の推薦図書もあげられている[これ](https://www.sciencedirect.com/science/article/abs/pii/S1567422319300547)や[これ](https://www.sciencedirect.com/science/article/abs/pii/S002002551632134X)など
[^2]:応用例としては二乗誤差を最小化する回帰や分散を最小化するポートフォリオ選択問題などがあります

[^3]:凸2次計画問題の標準形については [学習用テキスト 非線形計画法 (1) ２次計画問題](http://www.me.titech.ac.jp/~mizu_lab/text/PDF-NLP/NLP1-QP-problem.pdf) が丁寧
[^4]:これも関連していてよいよ！という記事があれば教えてください
[^5]:@yappy0625 さんが同じ問題に対して [5日目の記事](https://yamaimo.hatenablog.jp/entry/2023/12/05/220000) で [PICOS](https://picos-api.gitlab.io/picos/examples.html) というモデリング言語を利用した凸2次計画問題の実装について紹介されていました