分析ノート

minimize関数の使い方の記事2本目です。

前回の記事で基本的な使い方を紹介しましたが、それは探索範囲に特に制限を設けず、各変数について実数値全体から探索するものでした。

しかし、現実の問題では特定の範囲に絞って探索したいってこともよくあります。一番よくあるのは特定の変数は正の値に絞るってケースでしょうか。また、各値がそれぞれの範囲内に絞られる、長方形や直方体系の制限だけでなく、円の内側のような指定も可能です。順番に見ていきましょう。

まずは、bounds 引数を使った境界制約です。
これは、単純に各変数の探索範囲にそれぞれ下限上限を指定できます。下限だけ指定したいとか上限だけ指定したい、といった場合はnp.inf で無限大を利用してください。

次のように、bounds引数に各変数の下限上限の配列を渡します。境界を指定しなければ最小値がない関数(ただの2変数の和)で実験してみましょう。

import numpy as np
from scipy.optimize import minimize


# 目的関数の定義
def sample_function(x):
    return x[0]+x[1]


# 初期値
x0 = [0, 0]

# 最適化の実行
result = minimize(sample_function, x0, method='L-BFGS-B', bounds=[[-2, 5], [3, np.inf]])

# 結果の表示
print("最適化の結果:", result)
print("最小値をとる点 x:", result.x)
print("最小値 f(x):", result.fun)
"""
最適化の結果:   message: CONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL
  success: True
   status: 0
      fun: 1.0
        x: [-2.000e+00  3.000e+00]
      nit: 2
      jac: [ 1.000e+00  1.000e+00]
     nfev: 9
     njev: 3
 hess_inv: <2x2 LbfgsInvHessProduct with dtype=float64>
最小値をとる点 x: [-2.  3.]
最小値 f(x): 1.0
"""

想定通り、それぞれの変数が最小値だったときに和も最小値でしたね。

methodは例として’L-BFGS-B’を使いましたが、boundsを指定できるmethodは数種類に限られます。もし何かしらエラーが発生したらドキュメントを確認して対応しているか確かめましょう。これは次の方法でも同じです。

次は、不等式を使って領域を絞り込む方法です。
これは制約を課す関数を定義し、その関数が「正の値」をとる範囲で探索します。

例えば、中心が(5, 3) 半径が2 の円の内側だけ探索するようにしてみましょう。
constraints引数に渡す方法はちょっとクセがあります。

# 目的関数の定義
def sample_function(x):
    return x[0]+x[1]


# 不等式制約
def constraint(x):
    return 4 - (x[0]-5)**2 - (x[1]-3)**2

inequality_constraint = {'type': 'ineq', 'fun': constraint}

# 初期値
x0 = [5, 3]

# 最適化の実行
result = minimize(sample_function, x0, method='SLSQP', constraints=inequality_constraint)

# 結果の表示
print("最適化の結果:", result)
print("最小値をとる点 x:", result.x)
print("最小値 f(x):", result.fun)
"""
最適化の結果:  message: Optimization terminated successfully
 success: True
  status: 0
     fun: 5.171572875250672
       x: [ 3.586e+00  1.586e+00]
     nit: 6
     jac: [ 1.000e+00  1.000e+00]
    nfev: 18
    njev: 6
最小値をとる点 x: [3.58578644 1.58578644]
最小値 f(x): 5.171572875250672
"""

method では今までと違って、 SLSQP を使いました。というのも、COBYLA, COBYQA, SLSQP  でしかサポートされてないのです。この辺りも気をつけて使う必要がありますね。

最近よく使っている、scipyの最適化関数の一つであるminimizeについて、まだ記事を書いてなかったので紹介します。

公式ドキュメントはこちらです。
参考: minimize — SciPy v1.14.1 Manual

これはタイトルの通りで、数値を返す関数を渡すとその関数が最小値をとる引数を探してくれるものです。ちなみに、最大値になる引数を探すメソッドはないので最大値を探したかったら、その関数に-1をかけて符号を反転させた関数を用意してください。

サクッと一つやってみましょう。正解がわかりやすいよう、2次関数でも例にとりましょうか。次の関数を使います。

$$x^2-6x+5 = (x-3)^2 -4 $$

平方完成から分かる通り、$x=3$で最小ですね。

import numpy as np
from scipy.optimize import minimize

# 目的関数の定義
def sample_function(x):
    return x**2 - 6*x + 5

# 初期値
x0 = 0

# 最適化の実行
result = minimize(sample_function, x0, method='L-BFGS-B')

# 結果の表示
print("最適化の結果:", result)
print("最小値をとる点 x:", result.x)
print("最小値 f(x):", result.fun)
"""
最適化の結果:   message: CONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL
  success: True
   status: 0
      fun: -3.9999999999999982
        x: [ 3.000e+00]
      nit: 2
      jac: [ 1.776e-07]
     nfev: 6
     njev: 3
 hess_inv: <1x1 LbfgsInvHessProduct with dtype=float64>
最小値をとる点 x: [3.00000004]
最小値 f(x): -3.9999999999999982
"""

上のサンプルコードのように、最初の引数で最適化したい関数、2個目の引数で初期値、そしてオプションで利用するアルゴリズムなどの各種設定を行います。

ここで注意しないといけないのは、渡した関数の第1引数だけが最適化の対象ということです。なので、元の関数が複数の引数をとる多変数関数の場合、最適化したい引数等を一個の配列にまとめて渡す関数でラップしてあげる必要があります。また、最適化対象外の値を固定する引数については、argsで固定します。

例えば、$f(x, y, a, b) = x^2 + ax + y^2 + by$ みたいな関数があって、a, bは固定したとき、これを最小にする $x, y$を求めたい場合次のようにします。

# 元の関数（オリジナルの目的関数）
def original_function(x, y, a, b):
    return x**2 + a*x + y**2 + b*y


# 最適化用にラップする関数
# 元の関数のx, y を xという配列で渡してx[0], x[1]として内部で使っている
def optimization_target_function(x, *params):
    return original_function(x[0], x[1], *params)


# 最適化の実行
x0 = [0, 0]
result = minimize(optimization_target_function, x0, args=(4, -6), method='L-BFGS-B')

# 結果の表示
print("最適化の結果:", result)
print("最小値をとる点 x:", result.x[0])
print("最小値をとる点 y:", result.x[1])
print("最小値 f(x, y):", result.fun)

"""
最適化の結果:   message: CONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL
  success: True
   status: 0
      fun: -12.999999999999925
        x: [-2.000e+00  3.000e+00]
      nit: 2
      jac: [-1.776e-07 -7.105e-07]
     nfev: 9
     njev: 3
 hess_inv: <2x2 LbfgsInvHessProduct with dtype=float64>
最小値をとる点 x: -2.0000000804663682
最小値をとる点 y: 2.9999997404723753
最小値 f(x, y): -12.999999999999925
"""

最適化したい引数が何変数なのか？ってのを一見どこでも指定してないように見えて不安になるのですが、これは初期値の配列の要素数から自動的に判断してくれています。

結果は result.x に入っているのでここから自分で取得します。

本当に基本的な使い方は以上になります。

少し発展的な内容なのですが、最適化のアルゴリズムは何種類もある中から選べますが、中には導関数を利用するものもあります。このminimizeに導関数を渡さない場合は、数値微分でいい感じにやってくれるのですが、明示的に導関数を指定することも可能です。

やり方は簡単で、元の関数と同じ形で引数を受け取るように導関数を定義して、jac引数に渡すだけです。1変数の場合は簡単すぎるので2変数の例を出しますが、2変数の場合は第1引数による微分と第2引数による微分の2つがあるので、その結果を配列で返す関数を定義してそれを渡します。

例えば、$f(x,y)=x^2+y^2+4x+6y+13$ みたいなのを考えてみましょう。$x, y$での微分はそれぞれ、$2x+4$, $2y+6$ですね。

# 目的関数
def objective_function(vars):
    """
    2変数関数 f(x, y) = x^2 + y^2 + 4x + 6y + 13
    vars: [x, y]
    """
    x, y = vars
    return x**2 + y**2 + 4*x + 6*y + 13

# 勾配（偏微分値）
def gradient_function(vars):
    """
    勾配 ∇f(x, y) = [∂f/∂x, ∂f/∂y]
    vars: [x, y]
    """
    x, y = vars
    grad_x = 2 * x + 4
    grad_y = 2 * y + 6
    return np.array([grad_x, grad_y])  # 勾配ベクトルを返す

# 初期値
x0 = [0, 0]

# 最適化の実行 (勾配を指定)
result = minimize(objective_function, x0, jac=gradient_function, method='L-BFGS-B')

# 結果の表示
print("最適化の結果:", result)
print("最小値をとる点 [x, y]:", result.x)
print("最小値 f(x, y):", result.fun)

"""
最適化の結果:   message: CONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL
  success: True
   status: 0
      fun: 0.0
        x: [-2.000e+00 -3.000e+00]
      nit: 2
      jac: [ 0.000e+00  0.000e+00]
     nfev: 3
     njev: 3
 hess_inv: <2x2 LbfgsInvHessProduct with dtype=float64>
最小値をとる点 [x, y]: [-2. -3.]
最小値 f(x, y): 0.0
"""

関数が簡単で計算量も小さい場合は数値微分でも特に問題ないのですが、そうでない場合は明示的に導関数を渡すことでリソースを節約することもできます。

長くなってきたので今回の記事はここまでにしようと思います。境界条件とうの制約をつける方法とかを今後紹介したいですね。

最後にminimizeを使う時の注意点ですが、これ、極小値が複数あるような関数の場合、最小値ではなく極小値の一つを返してくることがあります。結果が初期値に依存してしまうのです。複雑な形の関数を最適化する場合は注意してください。

とある特殊な環境でPythonを書いていて、いくつかのライブラリがimportに失敗するという事態に遭遇しました。自分のローカルPC上であれば、pip freeze とかしてライブラリのバージョンを調べて原因を調査するのですが、その環境ではOSコマンドが打てず、同様の調査が不可能でした。importさえできれば {ライブラリ名}.__version__ みたいなプロパティから取得することもできたのですがimport自体が失敗するとあって調査に苦戦していました。

ところがどうやらimportを行わずにライブラリのメタデータにアクセスする方法がちゃんとあるようだったのでこの記事にまとめておきます。(正直、普通の環境であればpipで調べれれば済む話なので、ほとんどの人にとっては不要な知識だと思います。)

pkg_resources を使う方法

importせずにライブラリのバージョンを取得する方法の1個目は pkg_resources を使うものです。

次の例は、 arviz というライブラリのバージョンを調べたものです。

from pkg_resources import get_distribution


try:
    version = get_distribution("arviz").version
    print(f"version: {version}")
except pkg_resources.DistributionNotFound:
    print("ライブラリが見つかりません")

# version: 0.16.1

僕は上記の方法で一旦解決しました。

ただ、これはsetuptoolsに依存したライブラリなのですが、その公式ドキュメントを見ると、もう廃止されたからimportlib.metadataを使えと書いてあるのですよね。

ということで合わせてそちらを紹介します。

importlib を使う方法

importlib はPython3.8から標準ライブラリに含まれたライブラリです。

これはインストール済みのライブラリのメタデータにアクセスする機能を持っています。標準ライブラリになったわけなので、最近のPythonであればおそらくこちらを使うのが適切なんだと思います。

参考: importlib.metadata — パッケージメタデータへのアクセス — Python 3.13.0 ドキュメント

サンプルコードでもバージョンを取得していますね。それにならってやってみましょう。

import importlib.metadata


try:
    version = importlib.metadata.version("arviz")
    print(f"version: {version}")
except importlib.metadata.PackageNotFoundError:
    print("ライブラリが見つかりません")

# version: 0.16.1

同じような結果が得られました。

余談ですが、公式ドキュメントのサンプルコードでは、
from importlib.metadata import version
として versionという関数をインポートしています。僕も最初それに倣ってやったのですが、versionを結果の変数名で使いたいな、と思ったので上記のインポート方法にしました。
ただ、これはこれでイマイチな気もします。