symbolic演算と数値化（対数尤度のHessian計算を題材にして）

パラメトリックな対数尤度関数を定義し、これの微分やヘシアンが欲しいが、手計算でやる気がしないとき、symbolic math toolbox を使って、 symbolic な微分（数値微分ではない解析的な微分）をすることができる。

準備

loglikelihood function が、パラメタ q,r,s,t,u,v,a の関数で書けているものとする。

loglikelihood.m

function L = loglikelihood( theta, Data )
q = theta(1);
r = theta(2);
s = theta(3);
t = theta(4);
u = theta(5);
v = theta(6);
a = theta(7);

L=0;
for i=4:5
  for j=1:10
     L = L + Data(i,j) * log( q * (r-1) + .......

以下のようにして上記を一つの symbolic な数式にまとめあげることができる。

script_hoge.m

syms q,r,s,t,u,v,a;
theta = [q,r,s,t,u,v,a];
Data = sym('Data', [4,5]);
L = loglikelihood( theta, Data );

gradient や Hessian を得る

7次元パラメタベクトル theta=[q,r,s,t,u,v,a] に関する gradient や Hessian 行列は以下の様に得られる。

script_hoge.mの続き

grad = gradient( loglikelihood, theta );
H = hessian( loglikelihood, theta );

gradient や Hessian の数値を得る

データ Data やパラメタ theta　の数値が、それぞれ DataVaue や thetaValue の数値ベクトルの形で与えられているとき、これを代入するには以下のようにする。

script_hoge.mの続き

grad_val = subs( grad, Data, DataValue );
grad_val = subs( grad_val, theta, thetaValue );
H_val = subs( H, Data, DataValue );
H_val = subs( H_val, theta, thetaValue );

二度に分けて代入する必要があることに注意。しかし theta の中身である q,r,... に対していちいち explicit に代入する必要はない。

なお、Dataやthetaの次元が大きくなってくると計算時間がかかりすぎることもあるかもしれません。筆者は Data として 4元分割表をあらわす4×4行列を入れて theta には7次元のパラメタを入れたところ計算はほぼ一瞬でしたが、 H symbolic な表示には数秒かかりました。

半日かけて計算ノートを数ページをつぶしながらあとでバグを恐れる事態を避けることができる幸せを皆様にも！