Coding Memorandum

プログラミングに関する備忘録

スポンサーサイト

上記の広告は1ヶ月以上更新のないブログに表示されています。
新しい記事を書く事で広告が消せます。

SVMコード

SVM(SMO法)のコードを載せておきます。理論についてはこちらをどうぞ。


2011.9.18追記:
ここでも書いたのですが,引数 target は -1 or 1 と設定頂く方が正しいはずです。(理論的にも正しいはず)
初出時は簡単なテストデータで動作確認を行っていたためか 0 or 1 でも動作確認が取れていたよう記憶しています。本格的に学習させる場合には 0 or 1 では学習が収束しませんでした。
大まかな流れ

SMOではラグランジュ未定乗数が教師データの数量分作成されます。これらを下記の処理フローで更新していきます。

ラグランジュ乗数の更新評価は1つずつ順番に行い,全ての評価を終えたところで実際に更新が行われたかをチェックします。1件でも更新したものがあれば,更新の作業を続けていきます。

下記のコードでは,個々のラグランジュ未定乗数の更新チェックはexaminUpdate()で行っています。ここでは,KKT条件を満たすかどうかのチェックを行います。
KKT条件を満たさないとき,ラグランジュ乗数を更新させます。SMOでは2点をペアとして更新しますので,update()でペアとなるラグランジュ乗数を探し,stepSMO()で2点を更新しています。

実装上の考慮点

SMOを実装する上では,次の点の考慮が必要です。

▼ Loose KKT Conditions
Platt の本に書かれているのですが,KKT条件のチェックは緩く行います。下記のコードでは,eps,tolerance変数がその役割を担っています。

▼ 誤差値のキャッシュの更新
ラグランジュ乗数の更新時には,キャッシュしている誤差値の更新を行います。
Platt の本を含め,多くの資料で更新対象a[i]に対応する誤差は 0 にすると書かれているのですが,a[i] がclippingされるときは,必ずしも 0 にはなりません。下記コードでは,clippingされたときには実際に誤差値を求めることとしています。(Line.274)

a[i]とペアを組むa[j]に関する誤差は,基本的には誤差更新の式(Line.264)で更新可能なようですが,何かの条件で実際の誤差と合わなくなることがあるようです(条件は分かっていませんが)。下記コードでは,a[j]の誤差は更新式を使わず誤差値を実測することとしました。(Line.279)

コード

下記にSVMのコードを示します。本コードは,クラス分けは2値(0 or 1),学習/テストデータの各属性値は[0,1]に正規化した値で動作させました。その他の条件で動作するかは未評価です。

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <vector>
#include <algorithm>
using namespace std;

class SVM
{
public:
	// 学習
	int learning( vector< vector<double> >& train_set,	// 教師データ群
				  vector<double>& target				// 正解出力 (-1 or 1)
				);
	// 判別
	double discriminate( vector<double> test			// テストデータ
					   );
private:
	int examinUpdate(int i);		// a[i]の更新評価(KKT条件のチェック)
	int update(int i);				// a[i]との更新ペアa[j]を探す
	int stepSMO( int i, int j );	// a[i],a[j]を更新する
	double f(int i);
	double kernel( const vector<double> p, const vector<double> q );
private:
	vector<double> a;				// ラグランジュ乗数
	vector<double> w;				// 重み  a x target
	vector<int> sv_index;			// サポートベクトル(0でないa[])のインデックス
	vector<double> err_cache;		// エラー値のキャッシュ
	vector< vector<double> > m_train_set;	// 教師データ
	vector<double> m_y;			// 正解データ
	int train_set_size;				// 教師データの数量
	int data_size;					// 教師データの要素数
	double b;						// 閾値
	double C;						// 
	double eps;						// ラグランジュ乗数評価時の余裕値
	double tolerance;				// KKT条件評価時の余裕値
	double Ei, Ej;					// エラー値
};

int SVM::learning( vector< vector<double> >& train_set, vector<double>& target ){
	// 初期化
	a.resize(train_set.size());
	w.resize(train_set.size());
	err_cache.resize(train_set.size());
	sv_index.clear();
	m_train_set = train_set;
	m_y = target;
	train_set_size = train_set.size();
	if( train_set_size < 1 )
		return 0;
	data_size = train_set[0].size();
	vector< vector<double> >::iterator ts_it;
	for( ts_it = train_set.begin(); ts_it != train_set.end(); ts_it++ )
		if( ts_it->size() != data_size )
			return 0;
	if( target.size() != train_set_size )
		return 0;

	// パラメータ設定
	C = 1000.0;
	eps = 0.001;
	tolerance = 0.001;
	// 初期値設定
	fill(a.begin(), a.end(), 0.0);
	fill(err_cache.begin(), err_cache.end(), 0.0);
	b = 0.0;

	bool alldata = true;
	int changed;
	int loop = 0;
	while(1){
		changed = 0;
		if( loop > 500 )
			break;
		for( int i = 0; i < train_set_size; i++ ){
			if( alldata || (a[i] > eps && a[i] < (C-eps)) ){
				changed += examinUpdate(i);
			}
		}
		if( alldata ){
			alldata = false;
			if( changed == 0 ){
				break;
			}
		}else{
			if(changed == 0){
				alldata = true;
			}
		}
		loop++;
		printf("loop %d : changed %d\n", loop, changed);
	}

	for( int i = 0; i < train_set_size; i++ ){
		if( a[i] != 0.0 ){
			sv_index.push_back(i);
		}
	}
	// 重みw 計算
	vector<int>::iterator it,it0;
	for( it = sv_index.begin(); it != sv_index.end(); it++ ){
		w[*it] = m_y[*it] * a[*it];
	}
	return 1;
}

int SVM::examinUpdate(int i){
	double yFi;
	if( a[i] > eps && a[i] < (C-eps) ){
		Ei = err_cache[i];	// f(x)-y
	}else{
		Ei = f(i) - m_y[i];
	}
	yFi = Ei * m_y[i];		// yf(x)-1

	// KKT条件のチェック
	if( (a[i] < (C-eps) && yFi < -tolerance) || (a[i] > eps && yFi > tolerance) ){
		return update(i);
	}
	return 0;
}

int SVM::update(int i){ // a[i]
	// a[j] の決定
	double max_Ej = 0.0, Ej;
	int max_j = -1;
	// 1
	int offset = (int)(((double)rand()/(double)RAND_MAX) * (double)(train_set_size-1));
	for( int j = 0; j < train_set_size; j++ ){
		int pos = (j+ offset) % train_set_size;
		if( a[pos] > eps && a[pos] < (C-eps) ){
			Ej = err_cache[pos];
			if( fabs(Ej-Ei) > max_Ej ){
				max_Ej = fabs(Ej-Ei);
				max_j = pos;
			}
		}
	}
	if( max_j >= 0 ){
		if( stepSMO(i,max_j) == 1 ){
			return 1;
		}
	}
	// 2
	offset = (int)(((double)rand()/(double)RAND_MAX) * (double)(train_set_size-1));
	for( int j = 0; j < train_set_size; j++ ){
		int pos = (j+ offset) % train_set_size;
		if( a[pos] > eps && a[pos] < (C-eps) ){
			if( stepSMO( i, pos ) == 1 ){
				return 1;
			}
		}
	}
	// 3
	offset = (int)(((double)rand()/(double)RAND_MAX) * (double)(train_set_size-1));
	for( int j = 0; j < train_set_size; j++ ){
		int pos = (j+ offset) % train_set_size;
		if( !(a[pos] > eps && a[pos] < (C-eps)) ){
			if( stepSMO( i, pos ) == 1 ){
				return 1;
			}
		}
	}
	return 0;
}

int SVM::stepSMO( int i, int j ){
	if( i == j )
		return 0;

	double ai_old = a[i], ai_new;
	double aj_old = a[j], aj_new;
	double U,V;
	if( m_y[i] != m_y[j] ){
		U = max(0.0, ai_old - aj_old);
		V = min(C, C+ai_old - aj_old);
	}else{
		U = max(0.0, ai_old + aj_old - C);
		V = min(C,   ai_old + aj_old);
	}
	if( U == V )
		return 0;

	double kii = kernel(m_train_set[i], m_train_set[i]); 
	double kjj = kernel(m_train_set[j], m_train_set[j]);
	double kij = kernel(m_train_set[i], m_train_set[j]);
	double k =  kii + kjj - 2.0*kij;
	if( a[j] > eps && a[j] < (C-eps) ){
		Ej = err_cache[j];
	}else{
		Ej = f(j) - m_y[j];
	}

	bool bClip = false;
	if( k <= 0.0 ){
		// ai = U のときの目的関数の値
		ai_new = U;
		aj_new = aj_old + m_y[i] * m_y[j] * (ai_old - ai_new);
		a[i] = ai_new; // 仮置き
		a[j] = aj_new;
		double v1 = f(j) + b - m_y[j] * aj_old * kjj - m_y[i] * ai_old * kij;
		double v2 = f(i) + b - m_y[j] * aj_old * kij - m_y[i] * ai_old * kii;
		double Lobj = aj_new + ai_new - kjj * aj_new * aj_new / 2.0 - kii * ai_new * ai_new / 2.0 
					 -m_y[j] * m_y[i] * kij * aj_new * ai_new
					 -m_y[j] * aj_new * v1 - m_y[i] * ai_new * v2;
		// ai = V のときの目的関数の値
		ai_new = V;
		aj_new = aj_old + m_y[i] * m_y[j] * (ai_old - ai_new);
		a[i] = ai_new; // 仮置き
		a[j] = aj_new;
		v1 = f(j) + b - m_y[j] * aj_old * kjj - m_y[i] * ai_old * kij;
		v2 = f(i) + b - m_y[j] * aj_old * kij - m_y[i] * ai_old * kii;
		double Hobj = aj_new + ai_new - kjj * aj_new * aj_new / 2.0 - kii * ai_new * ai_new / 2.0 
					 -m_y[j] * m_y[i] * kij * aj_new * ai_new
					 -m_y[j] * aj_new * v1 - m_y[i] * ai_new * v2;

		if( Lobj > Hobj + eps ){
			bClip = true;
			ai_new = U;
		}else if( Lobj < Hobj - eps ){
			bClip = true;
			ai_new = V;
		}else{
			bClip = true;
			ai_new = ai_old;
		}
		a[i] = ai_old; // 元に戻す
		a[j] = aj_old;
	}else{
		ai_new = ai_old + (m_y[i] * (Ej-Ei) / k);
		if( ai_new > V ){
			bClip = true;
			ai_new = V;
		}else if( ai_new < U ){
			bClip = true;
			ai_new = U;
		}
	}
	if( fabs(ai_new - ai_old) < eps * (ai_new+ai_old+eps) ){
		return 0;
	}

	// a[j]更新
	aj_new = aj_old + m_y[i] * m_y[j] * (ai_old - ai_new);
	// b更新
	double old_b = b;
	if( a[i] > eps && a[i] < (C-eps) ){
		b += Ei + (ai_new - ai_old) * m_y[i] * kii +
				  (aj_new - aj_old) * m_y[j] * kij;
	}else if( a[j] > eps && a[j] < (C-eps) ){
		b += Ej + (ai_new - ai_old) * m_y[i] * kij +
				  (aj_new - aj_old) * m_y[j] * kjj;
	}else{
		b += (Ei + (ai_new - ai_old) * m_y[i] * kii +
			       (aj_new - aj_old) * m_y[j] * kij +
		      Ej + (ai_new - ai_old) * m_y[i] * kij +
			       (aj_new - aj_old) * m_y[j] * kjj ) / 2.0;
	}
	// err更新
	for( int m = 0; m < train_set_size; m++ ){
		if( m == i || m == j ){
			continue;
		}else if( a[m] > eps && a[m] < (C-eps) ){
			err_cache[m] = err_cache[m] + m_y[j] * (aj_new - aj_old) * kernel( m_train_set[j], m_train_set[m] )
										+ m_y[i] * (ai_new - ai_old) * kernel( m_train_set[i], m_train_set[m] )
										+ old_b - b;
		}
	}

	a[i] = ai_new;
	a[j] = aj_new;
	if( bClip  ){
		if( ai_new > eps && ai_new < (C-eps) ){
			err_cache[i] = f(i) - m_y[i];
		}
	}else{
		err_cache[i] = 0.0;
	}
	err_cache[j] =  f(j) - m_y[j];

	return 1;
}

double SVM::discriminate( vector<double> test ){
	if( test.size() != data_size ){
		return 0.0;
	}
	vector<int>::iterator it;
	double eval = 0.0;
	for( it = sv_index.begin(); it != sv_index.end(); it++ ){
		eval += w[*it] * kernel(m_train_set[*it], test);
	}
	eval -= b;

	return eval;
}

double SVM::f(int i){
	double F = 0.0;
	for( int j = 0; j < train_set_size; j++ ){
		if( a[j] == 0.0 )
			continue;
		F += a[j] * m_y[j] * kernel(m_train_set[j], m_train_set[i]);
	}
	F -= b;
	return F;
}

#define GAUSSIAN

double SVM::kernel( const vector<double> p, const vector<double> q )
{
	double r = 0.0;
#ifndef GAUSSIAN
	// 多項式カーネル
	double p = 4.0;		// Tuning Parameter
	r = 1;				// Tuning Parameter
	for( int i = 0; i < data_size; i++ ){
		r += p[i] * q[i];
	}
	r = pow( r, p );
#else
	// ガウシアンカーネル
	double delta = 1.0;		// Tuning Parameter
	for( int i = 0; i < data_size; i++ ){
		r += (p[i] - q[i]) * (p[i] - q[i]);
	}
	r = -r / (2*delta*delta);
	r = exp(r);
#endif
	return r;
}
スポンサーサイト

サポートベクターマシン (SVM)

サポートベクターマシンについて,情報源に関するメモを纏めておきます。

▼ 参考書
Marahon Match用に次の参考書籍を購入しました。
        サポートベクターマシン入門

理論から実践(インプリ)まで一通り書かれていて良い本だと思います。機械学習の基礎となるところから書かれていますので,機械学習の入門書としても良いように思います。
また,各章の最後に「さらなる文献と話題」の節が設けられ,参考文献の紹介が充実しています。

▼ Web公開の資料
サポートベクターマシンは旬な技術なようで,数多くの資料が公開されています。(特に大学の研究室での公開資料が多いですね)

今回は,以下の資料を参考にさせて頂きました。

SVMを調べていると「2次の目的関数が得られ,これを最適化すれば良い」というところまではすぐに到達できると思います。しかしながら,この最適化問題をどう解くのかという部分を見つけるのには,わりと時間がかかりました。(結局,この部分の詳細を知ることができず参考書籍を購入することになったのですが)

いくつかの手法がある中で,逐次最小最適化(Sequential Minimal Optimization,SMO)という手法が優れているようです。

▼ SMO実装に関する参考資料
SMOの実装については前述の参考書籍にも触れられているのですが,実際にコードを書き起こすにはいくつかの情報が足りていません。コードを書くためには本家の資料を参照する必要があります。
        J. Platt. Fast training of support vector machines using sequential minimal optimization.
        In B. Schölkopf, C. J. C. Burges, and A. J. Smola, editors, Advances in Kernel Methods
        -- Support Vector Learning
, pages 185--208. MIT  Press, 1999.

 ここから入手できます。

次のサイトには,SMOの擬似コードが掲載されています。
       SMO Pseudo-Code

 コードについてはこちらを参照ください。

TopCoder Marathon Match AgentMatching

Marathon Match AgentMatching に参戦しました。

Event Calenderにも予定がなく,急に始まったこのMarathon Matchでしたが,折りしも夏休みに入ったところであり,また今回の問題は機械学習ということで学生のころの研究分野でしたので,何とかなるかもと淡い期待を持ちつつ参加することにしました。(賞金付きというのも魅力的でしたが)

しかしながら,コトはそんなに上手く運ばず,結構な時間をかけながら紆余曲折の末,何とかそれなりのスコアに到達できたというところで終わりました・・・

本問題の解法として,まずは機械学習ということでニューラルネットかID3(C45) かというのが頭に浮かんでいました。教師データが離散値であるのでID3で行けるのではと思い,試してみました。
ID3についてはWikipediaを参照してみてください。(そんなにメジャーなものではないと思っていたので,項目があることに驚きでした)
    http://ja.wikipedia.org/wiki/ID3

この方法での 結果はほとんどスコアが伸びず,ある程度手を入れてみたものの,わりと早めに見切りを付けることとなりました。教師データにノイズが乗っていることからも難しいだろうなとの予測もありましたので。
 

続いてはサポートベクターマシンにトライしてみました。サポートベクターマシンは,ここ最近よく目にするようになって,一度はさわってみたいなと思っていたので丁度良い機会でした。非線形の分類ができて,ノイズにも(ある程度)対応するとのことでしたので,結構の良い結果を得られるのではと期待を持って取り組んでみました。
サポートベクターマシンについては実装に関する資料がほとんどなく(理論の解説は多く見受けられますが),インプリにはかなり苦労させられました。この知見については,別記事でご紹介できればと思っています。

しかしながら,こちらの結果も芳しくなく,しかもID3よりも悪いという想定外の結果となってしまいました。誤差の許容パラメータが厳しすぎると学習が収束せず,緩くすると全て「真」(あるいは「偽」)となってしまい,収拾がつかない状況に。
このような結果を受け,今回の学習データは矛盾するデータが多く,分類するという考え方がそぐわないのだろうとの考えに至りました。
 

残り時間も少なくなり,残り2日しかない中で確率ベースのアルゴリズムを考えてみました。
学習データの重要そうな項目について,その各項目の値ごとに「真」となる確率を求めておきます。テストデータに対しては,各値ごとに「真」となる確率を引きだし,その確率値の重み付き和を求め,出力としてみました。(重みはシミュレーションで学習させる)

これが(やっつけアルゴリズムでありながら)意外なことに上手く行き,何とか満足できる位置までスコアアップしました。
もっと早くこれに到達できていれば,もう少し上位を狙えたような気がします。
 

TopCoder SRM 446

久々のSRM。
このところ仕事が忙しく,なかなか時間が合いませんでした。

久しぶりに参加した今回のSRMはケアレスミスでスコアを大きく落とす展開となり,散々な結果となってしまいました。

Div2 500のCubeWalkingの問題(Div1 250と同じ)では,進行方向をx,y座標の増分dx,dyで管理してみたのですが,進行方向の90℃回転で次のようなコードを書いてしまいました・・・

if( *it == 'L' ){
    dx = dy;
    dy = -dx;
}else{  // *it=='R'
    dx = -dy;
    dy = dx;
}
dyの代入の前にdxを壊しているという凡ミス。これでExample Testが通ってしまうところも不幸でした。

250問題の方は終了間際でミスに気付き,再submitで最低点。

1000点問題は結構なコード量になりそうでしたので,心が折れました。

動的計画法 その1

イントロダクション

動的計画法(Dynamic Programing:以下「DP」と記述)は,ある種の最適化問題を効率的に(少ない計算時間で)解くための方法である。

アルゴリズム関連の参考書籍では,次のように説明される。

 k 個の要素だけをとったときの最適解が表の形で与えられているとする。これに要素を1個付け加えたときに,最適解がどのように変わるかをこの表をもとに計算し,表を書き換える。
 この操作を続けて,n 個の要素すべてを使った最適解が求まればよい。表の添え字としては,問題の性格を規定するパラメータのうち要素の個数以外のものを選ぶ。
アルゴリズムとデータ構造 (岩波講座 ソフトウェア科学) より
問題を細かく小さい問題に分けて,それぞれにベストの解答を見つけよ。それらの回答を張り合わせて,全体の問題の回答を導け。
プログラマのための論理パズル 難題を突破する論理思考トレーニング より

DPをインターネットで検索して得られる情報も概ねこのような説明であると思う。これらの説明と,いくつかの例題を見ることでDPのやりたいこと(方法)を(ある程度)理解することはできるであろう

しかしながら,実際に問題の解法としてDPを適用しようとした場合,「何を“小さい問題”として定義するのか?」,「何を“(k 個の)要素”とするのか?」を考え付くことが難しく(そもそもDPで解くべき問題であるかも見抜けないことも多い),解を得るまでに至らないのである。

ここでは,DPを数学面から見直し*1,各種例題を数学的に捉えなおすことでDPの理解を深める。「DPで考える」ための考え方を習得することを目指す。

DPの数学表現

動的計画法は「離散値を取る変数で表される関数の最大値または最小値(最適値)」を効率的に求める方法である。

変数 xi が 離散値{ a1,a2,・・・,an }をとるとき,

式1
の最大値,最小値を求める問題を考える。

一般的には,この関数の最適値を求めるには xi の全ての組み合わせについて関数の評価を行わなければならない。しかしながら,関数 f が次のような形で表されるとき,全件評価を行わずに最適値を求める方法が存在する。

式2
ここでは,この関数の最大値を求めることを考える。

x1 に着目すると, x1f1(x1),h1(x1,x2) にのみに関係する。変数 x2 に対して,f1(x1)+h1(x1,x2)の最大値を値とする関数f2(x2)を考える。

式3
式(1)が最大値をとるとき,次式
式4
も最大値をとり,式(1)の最大値と同じ値となる。

続いて式(2)について同様に,変数x3に対して,f2(x2)+h2(x2,x3)の最大値を値とする関数f3(x3)を考える。

式5
式(2)が最大値をとるとき,次式
式6
も最大値をとり,同じ値となる。

以下,同様の式の置きかけを続けることで,最終的には関数fn(xn)を次のように導くことができる。

式7
fn(xn)の最大値は式(1),(2),(3)の最大値と同じである。 このようにして部分ごとの最大値の計算を積み重ねることで,最終的に式(1)の最大値を求めることができる。

上述の説明は“最大値/max”を“最小値/min”に言い換えても成立する。(関数値の最小値を求めるときにも使える)

DPの計算手順は次の通りとなる。

関数f1(x1),hi(xi,xi+1) (i=1,n-1)が与えられたとき
1. f2 = max(f1(x1)+h1(x1,x2))を求める。
    x1, x2は離散値なので,f2(x2)も離散値となる。==> 配列へ保存

f2(x2)の値を用いて,
2. f3 = max(f2(x2)+h2(x2,x3))を求める。

以下,f4f5と求め,最終的には fnを求める。

fn の最大値が式(1)の最大値となる。

DPの適用

変数xi, xi+1を“問題”で定義される何らかの状態を表す値として捉えると,hi(xi,xi+1) は現在の状態 xi から次の状態 xi+1 へ遷移すると得られる値(得点)として考えられる。

DPは,この状態を進めることで得られる値(得点)の和について,その最適値を求める問題に適用することができる。

最適値に対応する各変数の値

式(1)が最適値をとるときの変数 xi の値が必要な場合には,次のように求める。

変数 xi+1 の関数,

式8
を求めるときに変数 xi+1 の各値に対して,最大値を与える xi を記録しておく。これを mi(xi+1) とする。

fn(xn) の最大値(pnとする)が求まると,mn-1(pn) から xn-1 の値が求まる。これをpn-1とする。 以下,同様にmn-2(pn-1) から pn-2 を求め,pn-3pn-4,・・・を順番に求める。pixi の値である。

例題1

TopCoder Tutorialから
http://www.topcoder.com/tc?module=Static&d1=tutorials&d2=dynProg

このTutorialのIntermediateの問題では,次の例題が挙がっている。

N×Mの表の各マスに,それぞれいくつかのリンゴが置かれている。左上のマスからスタートして,右もしくは下へ進み,通ったマスのリンゴを集めながら右下のマスまで行く。リンゴは最大で何個集められるか?

マスを進む各ステップ(1ステップで1マス移動)の状態(状態=どのマスの上にいるか)を x0, x1, ・・・, xN+M-2 とする。i ステップ目にいるマス xi は,左上のマスを(0,0),右下のマスを(N-1, M-1)とする座標系を考えたとき,(p, q) :p+q=i のマスとなる。 x0は左上のマスであり,xN-1+M-1 は右下のマスである。

fig1

この最適化(最大値)問題は,xi における獲得リンゴ数を f(xi) とすれば,

式1-1
を最大化する問題となる。

ここで,xixi+1 の間には進む方向に関するルールからの依存関係があることに注意が必要である。
xi 上のマス(p, q)は,xi-1 においてマス (p-1, q),もしくは(p, q-1) のときのみ移動可能である。

fig2
xi の位置に関わらず xi+1 上の任意のマスに移動できるのであれば,各ステップでf(xi) の最大値のマスを選択していけば良いこととなる(貪欲法となる)

ここで次の関数h(xi,xi+1)を考える。

式1-3
ここでN/Aは評価されない値である。プログラミングする場合は,この関数値は評価対象から外すこと(行けない経路は評価しない)となる。
進む方向の依存関係を考慮した獲得リンゴ数の式は次のように書くことができる。
式1-4
獲得リンゴ数を最大化するには,この式を最大化すればよい。

この式は前述の式(1)と同等の形式であり,DPを用いて最大値を計算することができる。f1(x1)を次のように定義すると,

式1-5
前述の式の最大値は,次式の最大値と一致する。
式1-6
これを繰り返し fN+M-2(xN+M-2) を求める。この値が求めるべき最大値となる。

例題2

有向グラフの経路最適化問題を考える。

Fig.2-1
SからスタートしてEまで進む経路を考える。進む方向は左から右へのみ進めることとし,i ステップ目にいる頂点を xi で表す。各枝には得点が与えられて,SからEの経路で枝を通過すると得点を獲得する。このとき,得点を最大化する経路を求める。

各枝は2頂点 xixi+1 を結ぶものであるので,枝の得点は次のように表すことができる。

式2-2
この値は,iステップ目の頂点 xi から i+1ステップ目の頂点 xi+1 を結ぶ枝の得点を示す。
SからEまでの経路で得られる得点は次のように書き表せる。
式2-3
この式の最大値は,これまでの説明通りにDPで解くことができる。式(1)に対応する初期値f(x0) は0であると考えればよい。

ここで,グラフの枝ではなく頂点に得点が与えられているケースを考える。この場合は,各ステップ i において,頂点 xi のうち最大値をとる頂点を訪問すれば良いことは自明である。このケースは各ステップで最大となるものを選択していく“貪欲法”となる。

動的計画法は,あるステップの値が1つ前のステップの状態に依存するよう定義されるときに,それらの値の組み合わせ(これまでの例では加算)を最適化するときに用いることができると捉えることができる。


*1参考文献:これなら分かる最適化数学―基礎原理から計算手法まで

FC2Ad

上記広告は1ヶ月以上更新のないブログに表示されています。新しい記事を書くことで広告を消せます。