MATLABで最適化の途中経過を可視化＆動画化する

はじめに

突然ですが、最適化っていきなり答えが導出されるので味気ないですよね？
もっとこう、頑張ってる感じがしたほうがよいと思いませんか。
マクロ組むとズルしてるとか言われる世界線だと特に

そんなわけでMATLABにおける可視化＆動画化の方法を紹介します。
バージョンは2010b以降、ツールボックスは一切不要です。

〇想定する読者
　MATLAB初心者、MATLAB芸人を志す人、最適化芸人を志す人

〇この記事で覚えて頂きたいこと
・最適化の基礎の関数 fminsearch はツールボックス不要なこと
・最適化における出力関数の扱い方
・グローバル変数と永続変数の使い分け(参考までに)

以下元コードを細切れにして解説を行うので、ちゃんとしたコードが欲しい方はGitHubよりダウンロード下さい。

元ネタ

最適化を介して曲線近似
https://jp.mathworks.com/help/matlab/math/example-curve-fitting-via-optimization.html
これを起点として、色々足し算しています。

①グローバル変数の宣言

最適化の途中経過を逐次画面に表示するには、出力関数内でこれら処理を行う必要があります。
なのでメインの関数から出力関数に渡す変数をグローバル宣言します。

clear all;
close all;

%% グローバル変数の宣言
global xdata; %曲線近似を行う対象データのx
global fun; %曲線近似を行う関数
global vidObj; %動画作成用のビデオオブジェクト

②曲線近似対象データの作成、プロット

グローバル宣言したxdata にxを、
ydataにyを代入します。曲線近似を行う関数はサンプルなのでなんでもＯＫですが、qiitaで見て面白かった最適速度模型の関数にしました。
（ウェブカメラ動画からx,yのデータを抽出⇒モデル化したら面白そう）

y=\tanh(ax-b)+\tanh(b)

a=4,b=0.8のもとydataを作成し、rand()でノイズを与えます。

xdata = [0:0.01:1];
fun = @(a,b,x) tanh(a*x-b) + tanh(b);
ydata = fun(4,0.8,xdata) + ( rand(size(xdata)) - 0.5 ) /10;

scatter(xdata,ydata);
axis([0 1 0 2]);grid on;
hold on;

このデータに対し最適化を用いてa,bの導出を行います。

③ビデオ作成の準備

MATLABのビデオ作成機能ですが、標準出力はaviです。
MATLAB芸人としてはtwitterに投稿したいのでmp4にします。
(qiitaに貼りたい場合はgif変換する)

%% ビデオ作成の準備
vidObj = VideoWriter('result.mp4','MPEG-4');
vidObj.FrameRate = 1/0.1;
open(vidObj);

単純に可視化のみでよくて、動画化は不要の場合はvidObjに関わる処理は不要です。

④最適化を実施、同時に可視化＆動画化

%% 最適化を介して曲線近似
% 目的関数の定義 入力は単一のベクトルのみしか許可されない点注意
efun = @(x) sum( ( fun(x(1),x(2),xdata) - ydata ).^2 );

% 最適化オプションの設定 出力関数outfunは後述する
options = optimset('Display','iter','OutputFcn',@outfun);

x0 = [0 0]; %最適化するa,bの初期値
bestx = fminsearch(efun , x0, options); %最適化の実行

%% ビデオ作成の完了
close(vidObj);

メインの関数はここまでとなります。
可視化＆動画化については出力関数にて行うのがミソです。

⑤出力関数にて可視化＆動画化

グローバル変数と永続変数を使います。
永続変数は、関数の呼び出し間で破棄して欲しくない変数に対し宣言します。
今回の例だと、プロット及びテキストに関するハンドルを永続変数にしています。これは可視化における近似曲線を１本だけにする　＆　テキストを逐次書き換えるためで、
各ハンドルをdeleteから呼ぶことで前回の線およびテキストを消しています。
(テキストは内容を書き換える処理でもいいが、面倒なので消して書き直している)

function stop = outfun(x, optimValues, state)

%% グローバル変数の宣言
    global xdata;
    global fun;
    global vidObj;

 %% 永続変数の定義
    persistent p;
    persistent t;
 
 %% 出力関数の中身
    stop = false;
    delete(p);delete(t); %近似曲線、テキストデータを図から削除
    
    %近似曲線の作成、プロット
    y = fun(x(1),x(2),xdata); 
    p = plot(xdata,y,'-r');
    
    %テキストデータの作成、プロット
    str = ['iter=' num2str(optimValues.iteration) newline];
    str = [str 'a=' num2str(x(1)) newline];
    str = [str 'b=' num2str(x(2))];
    t = text(0.7,0.4, str,'FontSize',14); 
    
    writeVideo(vidObj, getframe(gcf));
    drawnow
end

ちなみにですがgcfはget current figureで、今アクティブになっているfigureのハンドルが取れます。gcaならget current axis。
昔、何のことだか分からなくて呪文だと思っていた時期があります。

また出力関数ですが、別ファイルに分けて書くことも出来ます。私はファイルが増えるとメンテが面倒なのでメイン関数と同ファイルに書いちゃいます。（MATLABのこういういい加減なとこ好きです）

結果を再掲。頑張ってる感じがしていいでしょう？
最適化結果もちゃんとa=4,b=0.8に近づいていて文句なしです。

おまけ

iterationの可視化を併記すると更に感動的に！コードはこちら
なおiterationの可視化のみなら、最適化オプション設定のみで可能です。

これに更に仮装大賞のBGMでも付ければ、立派なMATLAB芸になることでしょう。