明朝仕事で早いので、トンガ王国のタイムゾーンにて日付変更直後に投稿させて頂きます。Machine Learning Advent Calendar 2012１７日目の内容です。
軽い気持ちで手を挙げてみたものの、蓋を開けたら（自分には）高度なテーマが次から次へと。一方で、今回初歩的な内容取り上げてしまったので少々場違い感があります。恐縮気味の@thso83です。どうぞよろしくお願い致します。
さて、12/17の担当回では、反復学習におけるヘッセ行列の役割について、ロジスティック回帰を例にとりながら解説を試みようと思います。12月17日といえばライト兄弟が人類初の動力飛行を成し遂げた日として有名ですが、みなさんも本記事を参考に初めての機械学習を成し遂げてはいかがでしょうか（こじつけ苦しいな・・・）。
さて、ロジスティック回帰は、回帰という名前が紛らわしいですが、2クラス分類問題における一般化線形モデルとしてよく知られています。入力の重み付き線形和に、ロジスティック関数（例えばsigmoid関数）をかませた確率モデルです。ネット上にもたくさん情報が見つかりますし、また、数式に関しては黄色い本が詳しいです。
重みパラメータを適切に設定することにより、あるデータが入力されたときに、そのデータが2クラスのどちらに属するかの確率を予測してくれるのですが、通常その重みパラメータは、反復最適化アルゴリズム（学習）を用いて求められます。
よく出てくる基本的な最適化アルゴリズムの一つとして、最急降下法があります。これはある目的関数（判別問題の場合は、教師データtと予測データyの誤差関数E）を、求めたい変数（重みw)で偏微分（勾配）を求め、関数曲面の最も急な方向にパラメータを近づけていくという手法になります。
例えば、下図のような単純な関数を考えると、勾配を頼りに極値に辿りつけることが直感的に理解できるかもしれません。下側のグラフにおいて、赤い程左手前方向への勾配が急で、逆に青いほど右奥への勾配が急なことを示しています。値が極小になるようにしたければ、勾配が下向きに大きい方向へパラメータを修正すれば、極小点へ効率的にたどり着くことができそうです。
最急降下法は、シンプルで使いやすいアルゴリズムなのですが、学習の効率化（反復回数の削減）において改善の余地が残されています。そこで、今回紹介するヘッセ行列が鍵となります。ヘッセ行列は、目的関数の２階微分を評価するのに用いられ、その固有値が全て正（正定値）なら極小点、固有値が全て負なら極大点、正と負の両方なら鞍点を示すという性質があります。上図で実際に確かめて頂けると思います。
そして、勾配（一次微分）だけではなく、極点での性質（二次微分）も利用した最適化アルゴリズムが、最急降下法よりも効率的なニュートン法となります。実際に２クラスの判別問題を、最急降下法を用いた学習（上）とニュートン法を用いた学習（下）で、それぞれ解いた様子を図に示します。
左図で三角及びバツ印でプロットされているのがそれぞれのクラスに属する学習データです。直線が重みパラメータによってひかれる判別境界線であり、学習回数に応じて濃淡を濃くして、その学習の進捗を表しています。また、右図では誤差関数値が学習回数に応じて減少していく様をプロットしています。
最急降下法（上）が徐々にしか収束しないのに対し、ニュートン法（下）では、最初の数回で効率良く最適解に収束していることが見て取れると思います。以下に最急降下法とニュートン法それぞれの学習を試したコードを掲載しておきます。octaveをインストールして、色々と実験し、収束の仕方の違いなどを確認して頂けると幸いです。
clear all; close all; clc
N = 30;
% data in class 1
x = randn(N, 2)*1.0 + 0.5;
t = ones(N, 1);
% date in class 2
x = [x; randn(N, 2)*0.5 - 0.5];
t = [t; zeros(N, 1)];
% bias term
x = [x, ones(size(x,1), 1)];
[num, dim] = size(x);
% sigmoid function
sigmoid = inline('1.0 ./ (1.0 + exp(-x))');
% number of training iterations
ITER = 60;
% methods for finding approximations
method = {'descent', 'newton'};
weight_default = rand(dim, 1);
%weight_default = zeros(dim, 1);
for m = 1:size(method,2)
% parameters
weight = weight_default;
% history for plot
err_history = 0;
weight_history = weight;
% training iterations
for i = 1:ITER
% Calculate the hypothesis function
y = sigmoid(x * weight);
% cross-entropy error function
err = (1/num) * sum(-t.*log(y) - (1-t).*log(1-y));
% calculate gradient
grad = (1/num) .* x' * (y-t);
% update parameters
switch method{m}
case 'descent'
weight = weight - grad;
case 'newton'
% calculate hessian.
H = (1/num).*x' * diag(y) * diag(1-y) * x;
% H\grad means H^(-1)*grad
weight = weight - H\grad;
end
err_history(i) = err;
weight_history(:, i+1) = weight;
end
% show result
subplot(2, size(method,2), 2*(m-1) + 1)
show(x, t, weight_history)
ylabel(method{m})
subplot(2, size(method,2), 2*(m-1) + 2)
plot(0:size(err_history,2)-1, err_history, '.--')
xlabel('Iteration'); ylabel('E')
drawnow
end
描画用関数
function show(x, t, weight_history)
class1 = find(t);
class2 = find(t == 0);
plot(x(class1,1), x(class1,2), 'ro')
hold on
plot(x(class2,1), x(class2,2), 'bx')
% decision boundary line
for i=1:size(weight_history, 2)
grad = i/size(weight_history, 2);
plot_x = [min(x(:,1)), max(x(:,1))];
if weight_history(2, i) == 0
plot_t = repmat(weight_history(3, i), 1, 2);
else
plot_t = (-1./weight_history(2, i)) .* (weight_history(1, i) .* plot_x + weight_history(3, i));
end
if grad == 1
width = 2;
else
width = 0.5;
end
plot(plot_x, plot_t, 'color', [0, grad, 0], 'linewidth', width)
end
legend('class 1', 'class 2', 'decision boundary')
axis([min(x(:,1)') max(x(:,1)') min(x(:,2)') max(x(:,2)')])
hold off
endfunction
octaveを起動し、main.mを実行すると、サンプルデータを生成し、その学習を行った結果をグラフ表示します。求められたパラメータを利用して、任意データの予測を試すこともできます。以下では、[1 1]と[-1 -1]（３つめの成分は常に１のバイアス）のそれぞれの入力に対して、クラス１に属する確率を求めています。図と見比べて、正しく予測できていることがわかると思います。
> main
> prob = sigmoid([1 1 1]*weight)
prob = 0.99934
> prob = sigmoid([-1 -1 1]*weight)
prob = 0.0011237
なお実際に、より高次元の実問題にヘッセ行列を活用する場合、その計算量（パラメータ数×パラメータ数のサイズの行列計算）が問題となってきます。そこで様々な近似手法が研究されているのですが、それに関しては、PRML復々習レーン#7での発表用に用意した資料がありますので、参考にしていただけると幸いです。
以上となります。
疑問点や間違っている箇所など、コメント頂けると幸いです。
さいごに、このような機会をご用意くださった@naoya_tさんへの感謝で締めくくりたいと思います。
どうもありがとうございました。