pyOpt の slsqp について

# 動機 下記のような最適化問題において、python のモジュールである pyOpt の slsqp アルゴリズムを使用してみた。 - 最適化問題 予測機を使って、仕入数から販売数を予測しながら、仕入の各店への配分を最適化する。 指標としては各店において、できる限り仕入数と販売数の差分が小さくなるように最適化する。 # 定式化 - 定数 総仕入数: $total\\_num$ 店舗数: $n$ 最低仕入数: $lower$ 各店舗のデータ: $df_i \in \mathbb{R}^k$, ( $k$ はデータの次元、$i$ は各店舗を表す添字 ) - 最適化定式 - 変数 : $x_i \in \mathbb{R}, i \in \{ i: 1 \leq i \leq n \} , ( 0 \leq x_i \leq 1.0 )$ - 補助定義 各店舗の仕入数: $y_i = total\\_num*x_i$ 販売数を出力する予測機 (機械学習などで学習して得た予測機): $pred : \mathbb{R}^k \times \mathbb{R} \to \mathbb{R} $ - 目的関数: 最小化 $ \sum_i { | y_i - pred( df_i,y_i ) |^2 )} $ - 制約: $ \sum_ix_i = 1$ ( $x_i$ が比率であることを意味する制約) $ y_i >= lower $ ( 仕入数は最低仕入数以上) # 上記問題の取り扱いにおける問題 予測機 pred は通常、数式的に与えられないため微分値を得るのが困難であり、最適化しずらい。 一般にこの種の問題の取り扱い方は二つある。 1. 数値微分を使う。 2. 関数の形を推論し、その推論した関数に対して最適化を行う。 1 の方が微分情報についてより直接的なので誤差は少ないと考えられる。ただし、ハイパーパラメータチューニングなどの場合、 そもそも関数の値を得るための計算時間が莫大になるために 2 の方法が選ばれる傾向にある。python の hyperopt モジュールは 2 の方法を取っている。 pyOpt のアプローチは 1 の「数値微分」を取り扱う方である。xgboost で学習した　pred の出力は比較的速いので有効と思われる。 数値微分をしているところについては本ブログの「pyOpt の微分情報について」を参照せよ。 # pyOpt のコード例 上記の定式化に対するコードを以下に掲載 ``` import numpy as np import pyOpt class FormulationForPuchaseOptimize: def __init__(self, model, df, purchase_total_num, lower, x_len ): self.model = model # 予測モデル. self.df = df # 各店舗のデータカラム. self.purchase_total_num = purchase_total_num # アイテムの全ての店舗の仕入数合計 self.lower = lower # 各店舗の最低仕入数. self.x_len = x_len slsqp = pyOpt.pySLSQP.SLSQP() # 1 にすると、multi process でやるときファイルエラーになる時がある. slsqp.setOption( "IPRINT", -1 ) self.slsqp = slsqp pass def create_formulation_prob( self ): # pyopt のインターフェイスに合わせるために下記のような関数を返すメソッドを用意した。 # (定式化に用いるパラメータを引数以外で関数に埋め込むため.) def get_func_value( x ): # x は 1 次元の np.array を想定 y = x * self.purchase_total_num # 割合が x の時の各店の仕入数. self.df[ "purchase_num" ] = pd.Series( y, index=self.df.index ) # out: 予測販売数 out = self.model.predict( self.df.values ) # ここが解析的にできない部分. f = mean_squared_error( out, y ) # 差分を目的関数とする. #f = mean_absolute_error( out, y ) # 差分を目的関数とする. # 等式、不等式制約を設定 g = np.zeros( 1 + self.x_len ) # 等式制約( == 0となるように記述する) g[ 0 ] = np.sum(x) - 1.0 # 合計=1, x が割合を表すようにするための制約. g[ 1: self.x_len+1 ] = self.lower - y # lower <= y, y は 1 以上という制約. fail = 0 return f, g, fail return get_func_value def prepare( self, seed ): # 変数初期化のための乱数シードの設定. np.random.seed( seed ) random.seed( seed ) opt_prob = pyOpt.Optimization( "purchase optimize", self.create_formulation_prob() ) # 変数の設定 x0 = np.random.rand( self.x_len ) # 0 ~ 1.0 の間の初期値をてきとうに選ぶ。 opt_prob.addVarGroup( 'x' , self.x_len, type='c', value=x0, lower=0.0, upper=1.0 ) # 販売予測値と仕入数の差分の最小化を目的関数とする。 opt_prob.addObj( "diff_pred_and_purchase" ) # sum(x) == 1 に相当する式の登録. opt_prob.addConGroup( "sum_x == 1", 1, "e" ) # lower <= y に相当する制約の登録 opt_prob.addConGroup( "lower_y", self.x_len, "i" ) return opt_prob def solve( self ,seed=1 ): # 定式化準備 opt_prob = self.prepare( seed ) [f_opt, x_opt, info] = self.slsqp( opt_prob, sens_type = "FD" ) if info['value'] != 0: print( "error! in %s", __name__ ) return None, None return f_opt, x_opt ``` create_formulation_prob は「 numpy.array である x が引数で渡された時に目的関数値 f と制約式の数値 g, 失敗したかどうかの数値 fail を返す関数」を返す関数である。 - なぜこんな面倒なことをしているかについての余談 slsqp のインターフェイスに合わせるためには引数は変えられない。引数以外からこの関数に対して、purchase_total_num のようなパラメータを埋め込みたい。python のように関数型言語的な機能がある場合、関数を返す関数で実装することになる(このように返された関数をクロージャとか呼んだりする)。 オブジェクト指向的なやり方ならば、そもそもインターフェイスとして slsqp クラスを継承して実装させるようにする。 普通の C 言語の場合、グローバル変数、スタティック変数、もしくは void* を引数に持たせて、内部で実装者がキャストして使うといったやり方になる。 このクラスのインスタンスを使用するコードは以下の通り ``` total_purchase = 30 lower = 1 with open( model_save_fpath, "rb" ) as f: model = pickle.load(f) df = pandas.read_csv( data_fpath, encoding="utf-8" ) optimizer = FormulationForPuchaseOptimize( model, df, total_purchase, lower, len(df) ) f_opt, x_opt = optimizer.solve() print(x_opt) ``` # 計算時間の目安 - 環境 下記の環境の docker 上で実行。 - 機種ID: MacBookPro14,2 - プロセッサ名: Intel Core i7 - プロセッサ速度: 3.5 GHz - プロセッサの個数: 1 - コアの総数: 2 - 二次キャッシュ（コア単位）: 256 KB - 三次キャッシュ: 4 MB - メモリ: 16 GB - 計測結果 　概ね 100 問に対して 15 秒程度 # pyOpt の実装 上記のコード例からもわかるように pyOpt では微分情報をどのように取っているのか明確ではない。pyOpt の実装を見ることでどのようにしているかがわかる。 - ソースコード ソースコードは svn コマンドにより次のようにして取ってくる。 ``` svn co http://svn.pyopt.org/trunk pyOpt ``` slsqp のメインとなるソースは pyOpt/pySLSQP にある。 pyOpt/pySLSQP/source ディレクトリにある fortran コードをコンパイルしたものを slsqp モジュールとし、pyOpt/pySLSQP.py から呼んでいるといった形式である。 - アルゴリズム 流れとしては初期点 x から目的関数値、制約式の数値、目的関数の微分値、制約式の微分値を用いて、 その周りの目的関数と制約式をそれぞれ二次関数、一次関数で近似した QP(quadratic programming) を解き、 変分 dx を決める。x = x + alpha* dx で x を更新。 といった形で、更新を進めていき、ある種の収束条件を満たしたら終了とする。 収束条件や alpha の詳細については以下の文献(ネットに落ちてる。)の chapter 2 による。 「Kraft D (1988) A software package for sequential quadratic programming. Tech. Rep. DFVLR-FB 88-28, DLR German Aerospace Center — Institute for Flight Mechanics, Koln, Germany.」 ただ、今だったら「Jorge Nocedal Stephen J. Wright, Numerical Optimization」(なぜかネットに落ちている。。)の chapter 18 にある 「Line search SQP method」を参考にしたほうが良いと思われる。 ちなみに nlopt の slsqp でも同様のソースを基にしている。 https://github.com/stevengj/nlopt/tree/master/src/algs/slsqp - QP について QP のアルゴリズム(lsq)は、「LAWSON & HANSON: SOLVING LEAST SQUARES PROBLEMS」(amazonで購入可能)という本に載ってるものそのまま。 この部分は40年くらい特に何もない。 # pyOpt の微分情報について 以下のようなコードで数値微分を行っている(説明の都合上、簡略化している)。デフォルトは sens_step = 1.0e-6 となっている。 pyOpt/pyOpt_gradient.py の class Gradient の getGrad メソッド ``` dh = sens_step xs = x k = 0 for i in mydvs: xh = copy.copy(xs) xh[i] += dh # ... 略 [fph,gph,fail] = opt_problem.obj_fun(xh, *args, **kwargs) if isinstance(fph,float): fph = [fph] for j in xrange(len(opt_problem._objectives.keys())): dfi[j,k] = (fph[j] - f[j])/dh for j in xrange(len(opt_problem._constraints.keys())): dgi[j,k] = (gph[j] - g[j])/dh k += 1 ``` # エラーが返されたのでデバッグ 店舗数が 30 以上といった場合に解けないということが生じていた. [f_opt, x_opt, info] = self.slsqp( opt_prob, sens_type = "FD" ) の info["text"] を見ると原因として、「Iteration limit exceeded」と出ている。 対処法としては以下の方法が考えられる。 1. 初期値を変えつつ何度もトライする。 2. 初期値を全て 1/total_num のように制約を満たすものにする。 3. slsqp.setOption( "MAXIT", 100 ) のようにイテレーション回数を増やす(デフォルトは 50) 4. slsqp.setOption( "MAXIT", 1.0e-4) のように判定誤差を甘くする。(デフォルトは1.0e-6 ) 5. sens_step(デフォルト 1.0e-6) を変える。(これは数値微分をとるときのステップサイズ) 6. 定式化を変える。 以下のように等式制約を変数にしてみた。 - 変数 : $x_i \in \mathbb{R}^{n-1} , ( 0 \leq x_i \leq 1.0 )$ - 補助変数(コード上では定義しない。) $ x_n = 1 - \sum_i(x_i) $ - 各店舗の仕入数: $ y_i = total\\_num*x_i$ - 目的関数: 最小化 $ \sum_i | y_i - pred( y_i ) |^2 $ - 制約: $\sum_{i\leq n-1}x_i \leq 1 $ ( $x_i$ が比率であることを意味する制約) $y_i \geq lower$ ( 仕入数は最低仕入数以上) - 結果 1,3,4 はどれも有効。現実としてはトライのたびに 3 のようにイテレーション回数を増やすといった実装が良いのではないか。 2 は試してないが解を偏らせるのではないかと心配。 4 も解を甘くしてしまうのは気になる所。 5 は試していない。predの感度に合わせてチューニングするのが良い。 6 は思ったよりうまくいかなかった。「LDL 分解する際に正定値でない」とかそんなエラーが生じる。