テキトーなメモ帳

テキトーなメモ帳

BFGSを実装してみた

BFGSを実装してみました。
BFGSとは準ニュートン法におけるHesse行列の逆行列の近似手法の一つです。


って言われてもよくわからないかもしれません。
機械学習のほとんど(CRF,MaximumEntropyなど)で行われているのは目的関数の最適化なのですが、その際の偏微分に用いられているアルゴリズムの一つです。
(といっても広く流布されているような機械学習ツールで使われているのはBFGSではなく、ちょっとこれの空間計算量が良くなったLBFGSだと思います)


プログラム:
github.com



参考にしたのは下記ページです

Spaghetti Source - 無制約非線型最適化 (準ニュートン法)
非線形計画法
BFGS 法 - 数理計画用語集


では、非線形計画法 のページの説明と今回自分が書いたコードの対応を取って説明していきます。

1. 初期点 x(0) を与え,k = 0,H(0) = I (単位行列)とします.

は下記に対応します。
テスト
https://github.com/titsuki/algorithm-quasi-newton/blob/master/t/calc.t

my $init_param = Math::MatrixReal->new_from_cols([ [0,0,0] ]); #初期点 x(0) を与え のところ

とプログラム本体の

    my ($rows, $columns) = $self->{x}->dim();
    my $B = Math::MatrixReal->new_diag([map { 1;} @{ [1..$rows ] }]); # "H(0) = I"のHは$Bと対応します。紛らわしいですが、Hesse行列の逆行列の近似 == $B == Hです。

2. g(k) = ∇f(x(k)) を計算し,∥g(k)∥ < ε ならば,計算を終了します.

は下記に対応します。

    my $prev_fx;
    do {
	$prev_fx = $fx;
        ...
	$fx = $self->f->($self->{x});
        ...
    } while(abs($fx - $prev_fx) > $EPS);

3. p(k) = -H(k)g(k), x(k+1) = x(k) + α(k)p(k) とします( α(k) は 1 次元探索で決める).

は下記に対応します。
「1 次元探索」とはここでは黄金分割探索のことを指しています。

sub golden_section_search {
    my ($self,$gradient_direction) = @_;

    my $p = $self->{x}; # $pは"x(k+1) = x(k) + α(k)p(k)"のx(k)
    ...
sub golden_section_search {
    my ($self,$gradient_direction) = @_; # $gradient_directionというのが"p(k) = -H(k)g(k)"のp(k)のことです。

    my $p = $self->{x};
    my $r = 2 / (3 + sqrt(5));
    my $a = 0;
    my $b = 1;
    my $t = $r * ($b - $a);
    my $c = $a + $t;
    my $d = $b - $t;
    my $fc;
    my $fd;
    
    $fc = $self->f->($p + $c * $gradient_direction);
    $fd = $self->f->($p + $d * $gradient_direction);
    
    while($d - $c > $EPS){
	if($fc > $fd){
	    $a = $c;
	    $c = $d;
	    $d = $b - $r * ($b - $a);
	    $fc = $fd;
	    $fd = $self->f->($p + $d * $gradient_direction);
	} else {
	    $b = $d;
	    $d = $c;
	    $c = $a + $r * ($b - $a);
	    $fd = $fc;
	    $fc = $self->f->($p + $c * $gradient_direction);
	}
    }
    return $p + $d * $gradient_direction; #ここの$dがα(k)。黄金分割探索により求められる。
}

4. 行列 H(k+1) を計算します.計算方法として,以下のような方法があります.ただし,各方法において,
   s(k) ≡ x(k+1) - x(k), y(k) ≡ g(k+1) - g(k)
とします.

BFGS公式(出典:BFGS 法 - 数理計画用語集


\begin{eqnarray}
H_{k+1}=H_k+\left(1+\frac{y_k^\top H_k y_k}{y_k^\top s_k}\right)\frac{s_k s_k^\top}{y_k^\top s_k}-\frac{H_k y_k s_k^\top + s_k y_k^\top H_k^\top}{y_k^\top s_k}
\end{eqnarray}

は下記に対応します。

    ...
    do {
    ...
	my $y = $g - $prev_g; # y(k) ≡ g(k+1) - g(k)はここ
	my $s = $self->{x} - $prev_x; # s(k) ≡ x(k+1) - x(k)はここ

        ### BFGS公式はじまり ####
	my $ys = (~$y * $s)->element(1,1);
	my $yBy = (~$y * $B * $y)->element(1,1);
	$B = $B
	    + (1.0 + $yBy / $ys) * (1.0 / $ys) * ($s * ~$s)
	    - (1.0 / $ys) * ($B * $y * ~$s + $s * ~$y * ~$B);
        ### BFGS公式おわり ####
    } while(abs($fx - $prev_fx) > $EPS);
    ...
}

以上です。

遊び方は、テストケース参照
https://github.com/okaoka/algorithm-quasi-newton/blob/master/t/calc.t

このテストケースにおける、dfというサブルーチンの返り値は"勾配ベクトル"です。
勾配ベクトルがよくわからない場合は、下記ページがとても参考になると思います。
勾配ベクトルの意味と例題 | 高校数学の美しい物語