MATLABの最適化関数(fmincon)の高速化Tips
はじめに
この記事は MATLAB/Simulink Advent Calendar 2023 のシリーズ2の 3 日目の記事です .
はじめまして,hook(@hookedonmas)といいます.某機械メーカーで自動化に関する研究開発をしている社会人一年目です. どうぞよろしくお願いします. 今回は,2年前に書いて下書きとして眠っていたものを,せっかくの機会なので掘り出してきました.
私が調査した内容を備忘録としてまとめたもので,実際に全てを試したわけではありません.
そのため,間違えている,または効果のない内容が含まれている可能性があります.
それを踏まえて参考として読んでいただけると助かります.
高速化方法
ヤコビ・ヘッセ行列の計算
通常,最適化の中で用いられる目的関数・制約の勾配(ヤコビ行列)の推定に使用されるのは有限差分であるため,勾配の計算は正確ではありません.
多くの場合,解の更新の回数が増え,計算時間が長くなります.
そこで,ヤコビ行列を戻り値とする関数を用意することで,より高速に計算することができます(詳細).
(ただ,勾配の計算が複雑な場合は,有限差分の方が速いことがあります.)
また,信頼領域 Reflective 法と内点法の場合は,ヘッセ行列を与えることができます. (詳細: 信頼領域 Reflective 法の場合と内点法の場合) ただ,ヘッセ行列の場合は有限差分による推定だけではなく,bfgs(準ニュートン近似)による計算も可能なので、問題によってはこちらの方が適している場合があります.
ヤコビ行列(とヘッセ行列)をfminconに与える方法はいくつかあり,以下にそれらの方法と参考リンクを挙げます.
自分の手で微分して計算
自分の手でヤコビ行列(とヘッセ行列)を計算する関数を用意する方法です.
Symbolic Math Toolboxがあれば,次の「シンボリック微分によって計算」の方法が楽でかつ正確です.
しかし,シンボリック微分で扱えないケースや,手計算により導いた導関数の方が評価が高速であるケースがあります.
Symbolic Math Toolboxが手元になくて,手計算によって導関数を求める必要がある場合,Wolfram Alphaが便利です.
シンボリック微分によって計算
シンボリック微分は数式微分とも呼ばれる方法で,自動微分とは別の方法です.
簡単に説明すると,"x2"を具体的なxの値を決めずに変数のまま扱って,微分することで"2x"を得て,その後"x=1"を代入することで微分係数を計算する方法です.
言わば,人間が微分係数を求める行為をそのままコンピュータが行うような方法です.
有限差分とは違い,正確な微分係数を求めることができます.
シンボリック微分を用いることで,目的関数と制約を定義するだけで,その導関数まで計算してくれます(詳細).
また,導関数自体は事前に計算しておくことができるので,fmincon自体の計算を短くできます.
ただし、微分したい関数が大規模で複雑だと、微分に時間がかかるので、注意が必要です。
自動微分によって計算
自動微分はシンボリック微分とは異なる方法で,微分の連鎖律を利用することで,微分の計算を効率的かつ正確に行う方法です(詳細).
深層学習などで用いられる誤差逆伝播法は自動微分の一部(一つのタイプ)です.
問題ベースの最適化問題solveに書き換えることで,自動微分を使うことができます(詳細).
並列計算
Parallel Computing Toolboxがある場合,並列計算によって計算時間を短くできる可能性があります.
基本的には,for文の中身を並列に動作できる場合(インデックスの順番を変えても結果が同じ場合)に,for文を並列化する並列計算が主です.
目的関数と非線形制約関数の有限差分を自動で並列に計算する設定にすることができます(詳細).
ただヤコビ行列・ヘッセ行列を与えた場合には自動で並列計算をしてくれないので,自分で並列化コードを実装する必要があります[5].
また,並列化の準備(parpoolの実行)に数秒時間がかかる上,並列化のオーバーヘッドが存在するので,対象が大規模な問題で,かつ十分なコア数を持つコンピュータの場合にのみ有効です.
そうでない場合は,計算が遅くなってしまう可能性があります.
MEXファイルの利用
計算量の多い関数をMATLABで実装するよりも,C++などで記述してからC++ MEX APIを用いる(参考)か,既存のC++ライブラリを呼び出す(参考)ことで高速化できる可能性があります.特にその関数がMATLABから提供されていない場合や,効率的な計算量のアルゴリズムでない場合,この方法が有効です.
私自身これを実際に利用したことはありませんが,何度かこのような利用をされているケースを見たことがあります.
例えば,この質問の質問者の方は,目的関数と制約およびそれらの導関数の計算に必要な微分方程式の数値解を,CVODESというC言語で書かれたライブラリをC MEXインターフェースによって使えるようにして,得ているようです.
アルゴリズム [6]
デフォルトのアルゴリズム(内点法)ではなく、問題に合わせてアルゴリズムを選択する。 (詳細)
その他
目的関数と制約を同じ関数で定義する
デフォルトだと目的関数と制約は異なる関数で定義します.
しかし,目的関数と制約の中の一部で同じ処理を行う場合,異なる関数で定義しているために,2度同じ処理を行うことになり冗長です.
このような場合に,目的関数と制約を同じ関数で定義する方法があります(詳細).
初期値をより最適解に近づける(warm start)
初期値が最適解に近いと反復回数が少なくて済むので,高速化されます.
また,目的関数・制約によっては解が局所解に陥ってしまい,いい最適解が得られないため,初期値にはできるだけ最適解に近い値を与えた方が良いです.
if を使わない [4]
収束基準 [4]
完全に収束するまで待たなくていいなら,許容値を大きくする.
別の関数を試す [4]
参考
[1] Optimization Toolbox での並列計算とは
[3] Parallel Computing Toolbox を使用した、時間のかかる最適化問題の最小化
[4] How to design optimisation and predictive control to be more computationally efficient?
[5] UseParallel Option with fmincon
[6] アルゴリズムの選択
[7] 勾配を含める