1つの応答曲面を遺伝的アルゴリズムで探査するPythonコード例

  • このエントリーをはてなブックマークに追加

機械学習で応答曲面を求めた後、その曲面(学習済モデル)の最小値や最大値を探査したい場合があります。ここではPythonによる実現の例として、とりあえずPyTorchで機械学習→PlatypusのNSGA-IIで探査といった方法を試してみます。

こんにちは。wat(@watlablog)です。ここではPyTorchとPlatypusで回帰と探査をやってみます

応答曲面法の概要と探査としての使い所

応答曲面法とは?

応答曲面法RSM : Response Surface Methotology)とは、シミュレーションや実験によってサンプリングされたデータを使って近似関数を作り、この関数を用いて計算を行う手法の事です。

以下は「Python機械学習!scikit-learnによる重回帰分析」で行った重回帰分析のメリットを説明する図ですが、応答曲面法のメリットも基本的にこれと同じです。

多数のデータが既にある場合で、設計変数間の未知の組み合わせにおける応答が知りたい場合に、まずは重回帰分析が有用な分析手法です。

重回帰分析の説明1

重回帰分析はデータ群に対して平面をフィットさせる手法です。このように平面をモデル(関数)で定義する事ができると、変数の連続的な変化による応答値の予測が可能となります。ちなみに設計変数が3つ以上の場合は超平面となり、もはや図で示す事が困難になります。

重回帰分析の説明2

応答曲面とは、重回帰分析における平面(または超平面)が曲面になったものと捉えて差し支えありません。

例えば以下の曲面は古典的機械学習の一種であるサポートベクターマシンを回帰モデル構築に利用した例です。データ群に対して曲面をフィッティングさせる事ができています。

※サポートベクターマシンについては以下が参考。
Pythonサポートベクターマシンで回帰分析!SVRの概要と実装
Python機械学習初心者用!サポートベクターマシンの概要と実装

応答曲面の例1

下図のような複雑な曲面もサポートベクターマシンによってフィットさせたものですが、このようにフィットされた曲面(つまりモデル)を使って応答を求める事で、この曲面は応答曲面として活用する事が可能です。

応答曲面の例2

曲面は平面に比べデータにフィットさせる自由度が増えるため、複雑な現象にもフィットする可能性を持ちます。もちろん重回帰分析と同様に設計変数が3つ以上の場合は超曲面で表現する事が可能です(図示は困難ですが)。

全ての応答を物理モデルで求めると時間がかかる

適切なモデルによるシミュレーションや確かな実験をする事で、物理現象に即した応答値を得る事ができます。しかし、より詳細な物理モデルを計算するにはスパコンを使った場合でも数日以上かかるケースがあります。また、実験はそもそも機器が高く段取り時間が多くかかったり、設計変数を正確かつ自由に振る事が困難であったりします。

単一の応答を得るだけならそれでも我慢できるかも知れませんが、最適化計算のように多目的で多数の計算結果を必要とするような検討には中々踏み込む事ができません。

応答曲面はそのような場合に活用可能です。
下図上段は任意の設計変数を物理モデルに設定して物理シミュレーションを行い、応答値を得る流れを示しています。

この物理モデル(あるいは実験)から少数の応答値をサンプリングし、そのサンプリングデータ点によくフィットするモデル(応答曲面)を使って応答値を得る方法が応答曲面法(下段)です。こうする事で、全ての応答を上段の物理シミュレーションで求めるよりは計算速度が圧倒的に速くなります。

応答曲面はあくまでデータから構成されるため、物理モデルと実験値を併用できる点もメリットの一つです。

物理シミュレーションと応答曲面法

但し、十分な設計空間を網羅していない場合や、サンプリング間隔が粗い場合等は、構築された応答曲面が予測した応答値の精度が著しく落ちてしまう事が注意点です。

応答曲面構築の手段

応答曲面は多数の手法で作成する事が可能です。簡単に列挙するだけでもステップワイズ回帰、フーリエ級数や多項式、様々な関数を使うパラメトリック曲面、ガウシアンプロセス…といった複数の方法があります。

色々ありますが、この記事ではニューラルネットワークによる応答曲面モデルの構築例を紹介します。

但しノーフリーランチ定理があるために応答曲面法も万能ではなく、どのような応答曲面構築手法を使ったとしても、出来上がったモデルの汎化性能が高いかどうかという観点は常に持ち続ける必要があります。

この辺の話は機械学習の勉強をすると嫌というほど出てきますので、是非当ブログの「AI」カテゴリをご覧ください。

設計変数が多いとひと目で最小値がわからない

続いて、応答曲面の最大値や最小値を求める方法(いわゆる最適化)の必要性を説明します。

例えば下図のように設計変数が2つ、応答である特性値が1つの場合は、設計変数をグリッドサーチして3Dプロットをする事でひと目で最大値と最小値を知る事が可能です。

設計変数が少ない場合

しかしながら、実際の工学の現場では設計変数が3つ以上ある場合、図示できないために簡単には最大値や最小値がわからない問題が多々あります。

設計変数が多い場合

このようなケースは、モデルの全貌がわからないために暗闇の中で勾配情報のみを頼りにする方法や、進化的アルゴリズムを使う方法で最大値や最小値を探査する手法がよく使われます。

この記事では後者の1手法であるNSGA-IIを使った方法で探査を行う方法を紹介します。

それでは次の節から、

  • ニューラルネットワークで応答曲面を作る
  • NSGA-IIで最小値を探査する

2つの方法について、Pythonコードでトライしてみた結果を紹介します。

PyTorchで応答曲面を作るPythonコード

探査するための応答曲面(学習済モデル)を作るために、PyTorchの多層線形ネットワーク回帰コードを紹介します。

参考記事

線形ネットワークによる非線形関数の回帰方法

PyTorchは今ではディープラーニングのフレームワークとして有名になってきました。本記事は応答曲面の作成にこのライブラリを使います。

PyTorchによる線形回帰の方法は以下の記事で紹介していますので、参考にしてください。

PyTorchで色々な非線形関数を回帰してみたらすごかった

モデルの保存

PyTorchの学習済モデルはプログラムを終了させるとメモリから解放されてしまう儚い存在です。ファイルに保存していつでも再利用できるようにしておくためにcloudpickleというライブラリを使っています。

Python標準で付いているpickleだとできません。

cloudpickleについては注意点もあるため、「PyTorchモデルをcloudpickleで保存・読み込みする方法」をご覧ください。

学習してモデルを作るためのコード

サンプルの関数

今回は最適化ベンチマークとして使われる関数をトレーニングで回帰しますが、関数は以下のQiita記事を参考にしました(すごく網羅的にまとめられていて大変参考になります!)。

Qiita:最適化アルゴリズムを評価するベンチマーク関数まとめ

この中のFive-Well Potential Functionという次式の多峰性関数を使いました。

\[ \begin{eqnarray} f(x, y)=(1 &-& \frac{1}{1+0.05(x^{2}+(y-10)^{2})} \\ &-& \frac{1}{1+0.05((x-10)^{2}+y^{2})} \\ &-& \frac{1.5}{1+0.03((x+10)^{2}+y^{2})}\\ &-& \frac{2}{1+0.05((x-5)^{2}+(y+10)^{2})}\\ &-& \frac{1}{1+0.1((x+5)^{2}+(y+10)^{2})})(1 + 0.0001(x^{2}+y^{2})^{1.2}) \end{eqnarray} \]

すさまじい式です…。

全コード

以下にコピペ動作する全コードを示します。ここでは学習の過程を動画で記録するために「Pythonで複数画像からGIFを作る時に便利な処理まとめ」で紹介した方法でGif画像を作っていますが、動画が必要無い人はこの処理を呼び出さないようにする事で処理時間が速くなります。

出来上がるmovie.gifはこちら。ちょっと見難いですが、データ点に対してうねうねと動いているカラフルな面が学習によってデータにフィットしようとしている応答曲面です。

最小値を表現するために頑張っている様子が確認できますね。

学習の様子(動画)

学習済モデルをファイル読み込みして確認するコード

モデルをいきなり探査する前に、せっかく可視化できる次元でトライしているので、先ほど学習した応答曲面モデルを目視で確認してみましょう。

以下がファイル読み込みからプロットするまでのコードです。学習部分を省くとこんだけスリムなコードになるんですね(matplotlibのこだわり部分が長いだけ)。

下図のように表示されればOK。学習の度合いによって形は崩れると思いますが、今回は1万回のイタレーションで多峰が表現されたのでよしとしましょう。

学習したモデル(ファイル読み込みで確認)

Platypusの遺伝的アルゴリズムで応答曲面を探査するPythonコード

参考記事

Platypusを使った最適化でモデルを探査します。Platypusや最適化については以下の記事を参考にしてください。

Platypusで多目的最適化からパレートフロントを求める方法

応答曲面の最小値を探査するコード

先ほど紹介した記事とほぼ同じ構成であるため詳細解説は省きます。注意点は、PyTorchのテンソル型を使ったり応答曲面を作った時のデータshapeに揃えたりする部分です。

先ほど保存したmodel-A.ptを読み込んでNSGA-IIによる最小値探査を行います。

以下が結果です。おしい点がいくつかあるようですが、急峻な峰の探査は結構難しいのかも知れませんね。

最小値探査

特定の設計変数のペアを抽出したい時は、algorithm.resultをforループで回している部分でsolution.objectives[0]の値を基準にすれば良さそうです。

一度配列に格納、ソートすることで目的関数が最小値になる設計変数を抽出する事が可能です。

まとめ

本記事ではこの分野の初心者がPythonを使って応答曲面を作ったり、作った応答曲面を探査したりするコードを紹介しました。

もっと複雑なデータでトライしようと思いましたが、まず第一歩としては目視で確認できる次元の範囲でやってみました。

まだまだ実問題に適用するには課題がいくつかありそうです…。

PyTorchとPlatypusを連携させたいがための記事を書いてみました!
Twitterでも関連情報をつぶやいているので、wat(@watlablog)のフォローお待ちしています!

SNSでもご購読できます。

コメントを残す

*