CVXPY詳解：Pythonで始める数理最適化 🚀

CVXPYは、Pythonに組み込まれた、凸最適化問題のためのモデリング言語です。オープンソースで開発されており、誰でも無料で利用できます。スタンフォード大学の研究プロジェクトとして始まり、現在では多くの研究者やエンジニアによって開発が進められています。

「モデリング言語」と言われてもピンとこないかもしれませんが、要するに、数理最適化問題を、数学的な記述に近い自然な形でPythonコードとして表現できるツールだということです。通常、最適化問題をコンピューターで解くためには、問題を特定の「標準形」と呼ばれる形式に変換し、ソルバーと呼ばれる計算エンジンに入力する必要があります。CVXPYを使うと、この面倒な変換作業を自動で行ってくれるため、ユーザーは問題の定式化に集中できます。

数理最適化とは？

数理最適化（Mathematical Optimization）とは、与えられた制約条件の中で、ある目的関数を最小化（または最大化）する変数の値を見つけ出すことです。

• 目的関数 (Objective Function): 最小化または最大化したい対象（例：コスト、利益、誤差など）。
• 変数 (Variables): 最適化によって決定される値（例：生産量、投資配分、パラメータなど）。
• 制約条件 (Constraints): 変数が満たすべき条件（例：予算上限、資源の限界、物理法則など）。

特に、目的関数と制約条件が特定の「凸性（convexity）」という性質を持つ問題を凸最適化問題と呼びます。凸最適化問題は、局所的な最適解が大域的な最適解と一致するという優れた性質を持ち、効率的に解を見つけるアルゴリズムが存在します。CVXPYは、この凸最適化問題を主なターゲットとしています。

CVXPYの利点 ✨

• 自然な記述: 数学的な表現に近い形で問題を記述できるため、直感的で分かりやすい。
• 汎用性: 線形計画(LP)、二次計画(QP)、二次錐計画(SOCP)、半正定値計画(SDP)など、様々な種類の凸最適化問題に対応。幾何計画(GP)や準凸計画(Quasiconvex programming)も扱えます。
• ソルバー連携: Clarabel, OSQP, SCSといったオープンソースソルバーが標準で組み込まれているほか、GUROBI, MOSEK, CPLEXなどの高性能な商用ソルバーや他のオープンソースソルバー（CBC, GLPKなど）とも連携可能（別途インストールが必要）。
• Pythonエコシステム: NumPyやSciPyといった他のPythonライブラリとシームレスに連携できる。
• パラメータ化: 問題の一部をパラメータとして定義し、その値を変えながら繰り返し効率的に問題を解くことができる（Disciplined Parametrized Programming – DPP）。
• 活発なコミュニティ: 継続的に開発が進められており、ドキュメントや例題も充実。DiscordやGitHub Discussionsでコミュニティと交流できます。

CVXPYを利用するには、まずPython環境にインストールする必要があります。Python 3.9以上が必要です。

前提条件

CVXPYは以下のライブラリに依存しています。通常、CVXPYをインストールする際に自動的にインストールされますが、念のため確認しておくと良いでしょう。

• Python >= 3.9
• Clarabel >= 0.6.0 (バージョン1.6以降)
• OSQP >= 0.6.2
• SCS >= 3.0
• NumPy >= 1.21.6
• SciPy >= 1.11.0

依存関係の衝突を避けるため、`virtualenv`や`conda`を使って独立した仮想環境を作成することを強く推奨します。

pipを使ったインストール

最も一般的な方法は`pip`を使うことです。ターミナルまたはコマンドプロンプトで以下のコマンドを実行します。

pip install cvxpy

これにより、CVXPY本体とデフォルトのソルバー（Clarabel, OSQP, SCS）がインストールされます。

特定の追加ソルバー（例：GLPK, CBC, CVXOPT）も一緒にインストールしたい場合は、以下のように指定します。

pip install cvxpy[GLPK,CBC,CVXOPT]

利用可能なソルバーのリストや詳細は公式ドキュメントを参照してください。商用ソルバー（GUROBI, MOSEKなど）を使用する場合は、別途ライセンスを取得し、各ソルバーの指示に従ってインストールする必要があります。

デフォルトのソルバーを含めずに、CVXPYのコア機能だけをインストールしたい場合は、`cvxpy-base`をインストールします。

pip install cvxpy-base

condaを使ったインストール

Anacondaを利用している場合は、`conda`コマンドでもインストールできます。`conda-forge`チャンネルからインストールするのが推奨されています。

conda install -c conda-forge cvxpy

condaでも、`cvxpy-base`のインストールが可能です。

conda install -c conda-forge cvxpy-base

インストール後、Pythonインタプリタで以下のコマンドを実行して、エラーが出なければ成功です。

import cvxpy as cp
print(cp.__version__)
print(cp.installed_solvers())

これにより、インストールされたCVXPYのバージョンと利用可能なソルバーのリストが表示されます。

CVXPYを使って最適化問題を解く基本的な流れを、簡単な線形計画問題(LP)を例に見ていきましょう。

基本的な構成要素

CVXPYで問題を記述するには、主に以下の要素を使います。

• 変数 (Variables): `cp.Variable()`で定義します。最適化によって決定される未知の値です。形状（スカラー、ベクトル、行列）や性質（整数、ブール値など）を指定できます。
• パラメータ (Parameters): `cp.Parameter()`で定義します。値が後から変更される可能性のある定数です。パラメータ化問題(DPP)で有用です。
• 式 (Expressions): 変数、パラメータ、定数（NumPy配列や数値）を演算子（+, -, *, @など）やCVXPYが提供する関数（`cp.sum`, `cp.abs`, `cp.norm`など）で組み合わせたものです。
• 制約条件 (Constraints): 式と等号（==）または不等号（<=, >=）で記述します。複数の制約条件はリストでまとめます。
• 目的関数 (Objective): 最小化または最大化したい式を`cp.Minimize()`または`cp.Maximize()`で指定します。
• 問題 (Problem): `cp.Problem()`で、目的関数と制約条件のリストをまとめて定義します。

例題：簡単な線形計画問題

以下の線形計画問題を考えてみましょう。

Maximize: 2x + 3y

Subject to:

• x + y ≤ 4
• 2x + y ≤ 5
• x ≥ 0
• y ≥ 0

CVXPYによる記述と求解

これをCVXPYで記述すると以下のようになります。

import cvxpy as cp
import numpy as np

# 1. 変数を定義
x = cp.Variable()
y = cp.Variable()

# 2. 目的関数を定義
objective = cp.Maximize(2*x + 3*y)

# 3. 制約条件を定義
constraints = [
    x + y <= 4,
    2*x + y <= 5,
    x >= 0,
    y >= 0
]

# 4. 問題を定義
problem = cp.Problem(objective, constraints)

# 5. 問題を解く
# solve()メソッドは最適値を返す。解の状態はproblem.statusに格納される。
optimal_value = problem.solve()

# 6. 結果を表示
print("解の状態:", problem.status)
# 'optimal'なら最適解が見つかったことを示す

if problem.status == cp.OPTIMAL or problem.status == cp.OPTIMAL_INACCURATE:
    print("最適値:", optimal_value)
    # 最適な変数の値は .value 属性に格納される
    print("最適解 x:", x.value)
    print("最適解 y:", y.value)
    # 制約に対応する双対変数の値は constraint.dual_value に格納される
    print("双対変数 (x + y <= 4):", constraints[0].dual_value)
    print("双対変数 (2x + y <= 5):", constraints[1].dual_value)
    print("双対変数 (x >= 0):", constraints[2].dual_value)
    print("双対変数 (y >= 0):", constraints[3].dual_value)
elif problem.status == cp.INFEASIBLE:
    print("問題は実行不可能です (Infeasible)")
elif problem.status == cp.UNBOUNDED:
    print("問題は非有界です (Unbounded)")
else:
    print("他のステータス:", problem.status)

このコードを実行すると、最適値と、それを達成するxとyの値が出力されます。`problem.solve()`を実行すると、CVXPYは問題を標準形に変換し、適切なソルバーを呼び出して解を計算し、結果を変数オブジェクトの`value`属性や制約オブジェクトの`dual_value`属性に格納します。

`problem.status`をチェックすることで、問題が正常に解けたか（`optimal`）、実行不可能なのか（`infeasible`）、非有界なのか（`unbounded`）などを確認できます。

ベクトルや行列を使う場合は、NumPy配列と組み合わせるのが一般的です。例えば、最小二乗問題を解く場合は以下のようになります。

import cvxpy as cp
import numpy as np

# 問題データの生成 (例)
m = 20
n = 15
np.random.seed(1)
A = np.random.randn(m, n)
b = np.random.randn(m)

# 変数の定義
x = cp.Variable(n)

# 目的関数 (最小二乗)
objective = cp.Minimize(cp.sum_squares(A @ x - b)) # '@' は行列乗算

# 制約条件 (例: box constraints)
constraints = [0 <= x, x <= 1]

# 問題の定義と求解
problem = cp.Problem(objective, constraints)
problem.solve()

# 結果の表示
print("解の状態:", problem.status)
if problem.status == cp.OPTIMAL or problem.status == cp.OPTIMAL_INACCURATE:
  print("最適値:", problem.value)
  # print("最適解 x:", x.value) # 値はベクトルとして得られる
else:
  print("問題は解けませんでした。")

このように、NumPyの知識があると、より複雑な問題を簡潔に記述できます。

CVXPYは、線形計画問題(LP)や二次計画問題(QP)だけでなく、より高度な凸最適化問題にも対応しています。

対応している主な問題クラス

• 線形計画問題 (Linear Programming – LP): 目的関数、制約条件がすべて線形の問題。
• 二次計画問題 (Quadratic Programming – QP): 線形制約のもとで、凸な二次関数を最小化（または凹な二次関数を最大化）する問題。
• 二次錐計画問題 (Second-Order Cone Programming – SOCP): LPやQPを包含する、より一般的なクラスの問題。ノルム項を含む制約などを扱える。
• 半正定値計画問題 (Semidefinite Programming – SDP): 変数が対称行列で、行列の半正定値性を制約とする問題。制御理論や組み合わせ最適化の緩和問題などに応用される。
• 幾何計画問題 (Geometric Programming – GP): 特定の形式（posynomial）を持つ目的関数と制約条件を持つ問題。対数変換により凸計画問題に変換できる (Disciplined Geometric Programming – DGP)。
• 混合整数凸計画問題 (Mixed-Integer Convex Programming – MICP): 変数の一部が整数値を取ることを許す凸計画問題。CBC, GLPK_MI, CPLEX, GUROBI, MOSEK, SCIPなどの整数計画に対応したソルバーが必要。
• 準凸計画問題 (Quasiconvex Programming – QCP): 目的関数や制約関数が準凸性を持つ問題。特定の条件下で解くことができる (Disciplined Quasiconvex Programming – DQCP)。CVXPY 1.6以降、ClarabelがDQCPに対応している。

具体的な応用例

CVXPYは様々な分野で活用されています。

公式ドキュメントのExamplesセクションには、これらの分野を含む多数の具体的な応用例とコードが掲載されています。

CVXPYには、基本的な使い方以外にも、より効率的かつ柔軟に問題を扱うための機能が備わっています。

Disciplined Convex Programming (DCP)

CVXPYが問題を正しく凸計画問題として認識し、標準形に変換できるためには、ユーザーが問題を特定のルールに従って記述する必要があります。このルール体系がDisciplined Convex Programming (DCP)です。

DCPの核心は、CVXPYの式（Expression）が持つ2つの属性「符号 (Sign)」と「曲率 (Curvature)」に基づいています。

• 符号: 正 (positive)、負 (negative)、ゼロ (zero)、不明 (unknown)
• 曲率: 定数 (constant)、アフィン (affine)、凸 (convex)、凹 (concave)、不明 (unknown)

CVXPYは、基本的な要素（変数、定数）の符号と曲率、そしてそれらを組み合わせる関数や演算子のルール（例：「凸関数 + 凸関数 = 凸関数」、「非負の定数 * 凸関数 = 凸関数」など）に基づいて、構築された式の曲率を自動的に判定します。

目的関数や制約条件がDCPルールに従っていない場合（例えば、凸関数を最大化しようとしたり、非アフィンな等式制約を使ったりした場合）、CVXPYはエラー (`DCPError`) を発生させます。

パラメータ化問題 (Parameterized Problems – DPP)

問題に含まれる定数の一部を `cp.Parameter()` で定義すると、そのパラメータの値を変えながら同じ構造の問題を繰り返し高速に解くことができます。これはDisciplined Parametrized Programming (DPP)と呼ばれる機能です。

最初に `problem.solve()` を呼んだ際に、CVXPYは問題の解析と標準形への変換（コンパイル）を行います。パラメータの値を変更した後は、`problem.solve()` を再度呼び出すと、コンパイル済みの情報を使って、より高速に解を求めることができます。これは、例えばパラメータを変えながら感度分析を行ったり、リアルタイム制御で繰り返し最適化計算を行ったりする場合に非常に有効です。

import cvxpy as cp
import numpy as np
import matplotlib.pyplot as plt

# パラメータを定義
rho = cp.Parameter(nonneg=True) # nonneg=True で非負パラメータを指定

# 変数と定数
n = 10
x = cp.Variable(n)
a = np.random.randn(n)
b = np.random.randn()

# 問題定義 (パラメータ rho を含む)
objective = cp.Minimize(cp.sum_squares(x - a) + rho * cp.norm(x, 1))
constraints = [cp.sum(x) == b]
problem = cp.Problem(objective, constraints)

# パラメータを変えながら解く
rho_values = np.logspace(-2, 2, 50)
x_values = []

for val in rho_values:
    rho.value = val          # パラメータに値を設定
    problem.solve()          # 問題を解く (2回目以降は高速)
    if problem.status == cp.OPTIMAL or problem.status == cp.OPTIMAL_INACCURATE:
        x_values.append(x.value)
    else:
        # 解けなかった場合の処理 (例: None を追加)
        x_values.append(np.full(n, np.nan))

# 結果のプロット (例)
# plt.figure()
# for i in range(n):
#     plt.plot(rho_values, [xv[i] for xv in x_values])
# plt.xlabel('rho')
# plt.ylabel('x_i')
# plt.xscale('log')
# plt.title('Solution path')
# plt.show()
print("DPPによる計算が完了しました。")

DPPを使うためには、問題がDPPのルールセットに従っている必要があります。これもDCPと同様に、パラメータの符号や曲率に関するルールがあります。詳細はDPPのドキュメントを参照してください。

ソルバーの選択と設定

CVXPYは、問題の種類に応じて自動的に適切なデフォルトソルバーを選択します（例：QPならOSQP、SOCPならClarabel）。しかし、`solve()`メソッドの`solver`引数を使って、使用するソルバーを明示的に指定することも可能です。

# SCSソルバーを指定して解く
problem.solve(solver=cp.SCS)

# MOSEKソルバーを指定して解く (別途インストールが必要)
# problem.solve(solver=cp.MOSEK)

インストールされている利用可能なソルバーは `cp.installed_solvers()` で確認できます。

また、ソルバー固有のパラメータ（例えば、最大反復回数や許容誤差など）を設定することもできます。`solve()`メソッドの`kwargs`を通じて、辞書形式でソルバーオプションを渡します。

# SCSの最大反復回数を5000に設定
problem.solve(solver=cp.SCS, max_iters=5000)

# OSQPの許容誤差を設定し、詳細なログを出力
problem.solve(solver=cp.OSQP, eps_abs=1e-7, verbose=True)

`verbose=True` を指定すると、ソルバーの実行過程に関する詳細なログが出力され、デバッグに役立ちます。利用可能なオプションはソルバーごとに異なるため、各ソルバーのドキュメントやCVXPYのソルバー設定ドキュメントを確認してください。

その他の高度な機能

• 双対変数 (Dual Variables): `constraint.dual_value` 属性を通じて、制約条件に対応する双対変数の値を取得できます。これは感度分析などに利用できます。
• N次元配列のサポート (v1.6+): 変数やパラメータをN次元配列として定義し、より自然に多次元データを扱えるようになりました。
• 微分可能プログラミング: パラメータに関する解の微分係数を計算する機能（`problem.derivative()`）。PyTorchやTensorFlowと連携して、最適化問題を組み込んだニューラルネットワークの学習などが可能です。

CVXPYは、Pythonで凸最適化問題を扱うための非常に強力で使いやすいライブラリです。

• 数学的な記述に近い自然な構文で問題をモデル化できます。
• LP, QP, SOCP, SDP, GP, MICP, QCPなど、幅広い種類の最適化問題に対応しています。
• 複数のオープンソースおよび商用ソルバーと連携し、問題に応じて最適なソルバーを選択できます。
• DPPにより、パラメータ化された問題を効率的に繰り返し解くことができます。
• DCPというルール体系により、問題の凸性を保証しやすくなっています。
• NumPyやSciPyなどのPythonエコシステムとシームレスに連携します。

数理最適化は、データサイエンス、機械学習、金融工学、制御工学、オペレーションズ・リサーチなど、非常に多くの分野で不可欠なツールとなっています。CVXPYを使えば、これらの分野における複雑な問題を、より直感的かつ効率的に解決へ導くことができます。

この記事が、皆さんのCVXPYを用いた数理最適化への第一歩となれば幸いです。ぜひ公式サイトのドキュメントや例題を参考に、様々な問題に挑戦してみてください！ 💪