Coding Memorandum

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

スポンサーサイト

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

超覚えやすいグラフの探査アルゴリズム

グラフの経路探査で有名なダイクストラ法,A*探索を深さ優先探索(DFS)・幅優先探索(BFS)とのアナロジーで理解できることに気づきました。これを覚えておけば,スクラッチから何も参照せずにダイクストラ法,A*探索のコードを書けるようになります。
(以下の説明は,ある程度アルゴリズムを理解していることを前提としています。「超分かりやすい」でない点に注意ください)

グラフ探査の抽象化コード
グラフ探査を抽象化したコード(C++風擬似コード)は次のようになります。
Search()
{
  Container que;                      // 抽象コンテナ ・・・(1)

  // スタート地点をキューに入れる
  Vertex s(スタート地点);
  que.push(s);

  while( !que.empty() ){
    Vertex v = que.pop();
    if( isGoal(v) )                   // ゴールに到達したら終了
      return;
    while( Vertex w = adjacent(v) ){  // vに隣接する頂点wを順番に得る
      if( Condition(w) ){             // 探索候補選択の条件 ・・・(2)
        que.push(w);
      }
    }
  }
}
グラフ探索のエッセンスだけを抜き取るとDFS,BFS,ダイクストラ法もA*法もこのようなコードで表すことができます。 (実際には,経路・訪問済み記録等のコードが足されます)

「(1)コンテナ」と「(2)探索候補選択の条件」が各アルゴリズムで異なってきます。

DFS,BFS

DFSの場合は,コンテナに「スタック(LIFO)」を用います。BFSの場合は,「キュー(FIFO)」を用います。

探索候補の条件は,DFS・BFSともに「wが未訪問である」という条件になります。

A*

A*探索の場合は,コンテナに「優先度付きキュー(priority queue)」を用います。優先順位は次の式で決まります。

f*(v) = g(v) + h*(v)

g(v) はスタート地点から頂点vまでの距離(ステップ数,辺の重みの総和等),h*(v) は頂点vからゴール地点までの見積距離を示します。"*"表記はヒューリスティック(発見的)関数であることを示しています。

h*(v) は 実際のゴールまでの最短距離 h(v) に対して,次の式を満たすと,

0 <= h*(v) <= h(v)

A*探索は最短の経路を見つけることを保証します。

探索候補の条件は,

「wが未訪問 or  (prior f*(w) > new f*(w))」

となります。頂点wが訪問済みの場合,それまでのwの見積もり値(prior f*(w))に対して,新しい経路(v→w)の見積もり値(new f*(w))が小さい場合には再度探索候補として採択します。

ダイクストラ法

ダイクストラ法の場合は,コンテナに「上書き可能な優先度付きキュー」を用います。優先順位は「スタートから頂点vまでの距離」となります。
コンテナへの格納(push)は,初めて頂点vへ(経路候補として)到達したときにそのときの距離で格納します。その後の探索で,頂点vへのより短い距離の経路候補が見つかった場合には,その距離値で上書きします。このため,「上書き可能な優先度付きキュー」が必要となります。

探索候補の条件は,

「wが未訪問 or (prior D(w) > new D(w))」

となります。prior D(w)はこれまで見つかっている頂点wへの最短距離,new D(w)は新しい経路(v→w)による距離を示します。
各頂点vのD(v)の初期値をINT_MAX等にすることによって,「wが未訪問」という条件チェックを省くことが可能です。

前述のようにコンテナへのpush()は,すでに頂点wの要素が格納済みの場合は,新しい距離値で上書きとなります。C++ STLでは,このようなコンテナはないので,set で erase(),insert()を組み合わせる,または自前でコンテナを実装する等の処理が必要です。

経路の記録

DFS,BFS,ダイクストラ法では,queからpop()された段階で頂点vへの経路が確定します。

また,DFS,BFSは未訪問点のみpush()しますので,push()する段階で頂点vへの経路が確定します。ダイクストラ法の場合は,上書きでpush()しますので,(pop()される前の)最後のpush()で頂点vへの経路が確定します。

DFS,BFS,ダイクストラ法での各頂点への経路情報は,頂点に対応する一次元配列を用意し,push()時にどの頂点から遷移したかを記録しておけば十分です。一つ前の頂点を順に辿ることでスタート地点までの経路を得ることができます。

一方,A*探索の場合は,ある頂点vについて(別の経路で)何度も評価される(何度もpop()される)可能性があります。すなわち,pop()された段階でも頂点vへの最短経路は確定しません。このため,経路情報はコンテナに入れる要素(Vertex)に記録して持ちまわる必要があります。

まとめ

まとめると下表のようになります。

 

コンテナ

探索候補の条件

優先順位

DFSスタック未訪問
BFSキュー未訪問
A*優先度付きキュー未訪問 or
  旧f*(w) > 新f*(w)
f*(v)
  =g(v)+h*(v)
ダイクストラ上書き
優先度付きキュー
未訪問 or
  旧D(w) > 新D(w)
D(v)
スポンサーサイト

プログラミングコンテストの数学

Topcoder SRM等の問題を見ると,ときおり数学的知識を要する問題が出てきます。
学業から離れて数年も経つと,仕事に関係しない知識はどんどんこぼれるもので,不意に数学の知識を問う問題に出会うと戸惑います。

この辺りのところを,今後のコンテストのためにメモにまとめておこうと思います。

素数

まずは素数判定から。素数判定には様々な手法が提案されていますが,小さな数の素数判定であれば次の単純な判定処理で十分ではないかと思います。

bool isPrime(unsigned int n)
{
    if( n == 2 ) return true;
    if( (n&1) == 0 || n == 1 ) return false;// 1と,2以外の偶数は素数でない
    int upper = (int)sqrt((double)n);
    for( int i = 3; i <= upper; i += 2 ){// 3以上の奇数を評価
        if( n % i == 0 ) return false;
    }
    return true;
}

判定対象nより小さい数で割り切れるかを見ます。割り切れた場合は素数でありません。割り切れるかの判定は&radicnまでで十分です。
ループ回数が入力値の平方根/2ですので,UINT_MAX(=4G)であってもたかだか30000回ほどのループ回数で済みます。


素数のリストが欲しい場合は,「エラストテネスのふるい」を用いるのが一般的です。

void sieve( int n, vector<int>& primes )
{
    if( n < 2 ) return;
    primes.push_back(2);
    if( n == 2 ) return;

    vector<bool> p(n/2,true);// p[x] => 2*x+3が素数であるかを示す
    int upper = (int)sqrt((double)n);
    upper = (upper-3)/2;

    for( int i = 0; i <= upper; i++ ){
        int t = i*2+3;
        if( p[i] ){
            primes.push_back(t);
            p[i] = false;
            for( int j = (t+1)*i; j < n/2; j += t ){
                p[j] = false;
            }
        }
    }
    for( int i = upper; i < (n-1)/2; i++ )
        if( p[i] )
            primes.push_back(i*2+3);
}

このコード例では,nまでの素数をvector primes に格納します。

素数判定用のワーク配列 pn/2 のサイズで用意し,素数であることを示すtrueで初期化します。コード例では,配列 p を3以上の奇数を示すインデックス値の配列として考えています。
「エラストテネスのふるい」では次の処理を繰り返します。

  1. 配列 p の中で,trueである最も小さい数を primes へ入れる
  2. 項1で選択した数の倍数をfalseにする

ふるいをかける上限は,素数の判定と同様に &radicn までで十分です。また,素数 pi の倍数をfalseにする部分では pi×pi 以降の数について処理すれば十分となります。これ以前の数は pi より小さい素数の倍数処理で既に false にされています。

最大公約数

規約分数を求めるときに最大公約数が必要となります。 最大公約数は,ユークリッドの互助法で求められます。

unsigned int gcd( unsigned int a, unsigned int b )
{
    unsigned int t;
    if( a < b ){
        t = a; a = b; b = t;
    }
    while( b != 0 ){
        t = a % b;
        a = b;
        b = t;
    }
    return a;
}

ユークリッド互助法は,下記の式が( b≠0 のとき)成り立つことを使っています。

この関係式を繰り返し適用することで最小公倍数を取得します。

3個以上の数の最大公約数では次の関係式が成り立ちます。


参考情報として「最小公倍数」は次の式で求められます。
unsigned int lcm( unsigned int a, unsigned int b )
{
    return a*b/gcd(a,b);
}

基数変換

10進数と n進数の相互変換を考えます。10進数 dn進数値の次の関係式を用います。



▼n進数への変換

void toNAdic( unsigned int decimal, vector<int>& coef, int base )
{
    do{
        unsigned int c = decimal % base;
        coef.push_back(c);
        decimal /= base;
    }while( decimal != 0 );
}

n進数の各桁の数値がvector coef に格納されます。

▼10進数への変換

unsigned int toDecimal( const vector<int>& coef, int base )
{
    unsigned int result = 0, u = 1;
    vector<int>::const_iterator it;
    for( it = coef.begin(); it != coef.end(); it++ ){
        result += *it * u;
        u *= base;
    }
    return result;
}

n進数の各桁の数値vector coef を与えると,10進の数値が返ります。

組み合わせ・順列

憶えていそうであやふやになりがちな「組み合わせ」と「順列」を整理しておきたいと思います。

組み合わせ 
順列 

▼組み合わせ

unsigned int combination( unsigned int n, unsigned int r )
{
    unsigned int numerator = 1, denominator = 1;
    if( n-r < r ) r = n-r;
    if( r == 0 )  return 1;
    for( unsigned int i = n; i > (n-r); i-- )
        numerator *= i;
    for( unsigned int i = 2; i <= r; i++ )
        denominator *= i;
    return numerator/denominator;
}
[long long版]
unsigned long long combination( unsigned long long n, unsigned long long r )
{
    unsigned long long p = 1, a= 1;
    if( n-r < r ) r = n-r;
    if( r == 0 )  return 1;
    for( unsigned long long i = n; i > (n-r); i-- )
        p *= i;
    for( unsigned long long i = 2; i <= r; i++ )
        a *= i;
    return p/a;
}

▼順列

unsigned int permutation( unsigned int n, unsigned int r )
{
    unsigned int result = 1;
    for( unsigned int i = n; i > (n-r); i-- )
        result *= i;
    return result;
}
[long long版]
unsigned long long permutation( unsigned long long n, unsigned long long r )
{
    unsigned long long result = 1;
    for( unsigned long long i = n; i > (n-r); i-- )
        result *= i;
    return result;
}

どちらのコード例もオーバーフローを考慮していませんので,適切な範囲で値を渡してください。

参考までに,階乗計算でオーバフローせずに計算可能な値は次の値までとなります。

long :12! = 479,001,600
long long : 20! = 2,432,902,008,176,640,000

三角関数

C/C++における三角関数の仕様を纏めておきます。

三角関数引数戻り値
sin(r)ラジアンsin値: -1 <= ret <= 1
cos(r)ラジアンcos値: -1 <= ret <= 1
tan(r)ラジアンtan値: π/2で∞であるが,実際には近似値を使うため値が返る。
asin(t)-1 <= t <= 1ラジアン: -π/2 <= ret <= π/2
acos(t)-1 <= t <= 1ラジアン: 0 <= ret <= π
atan(t)任意ラジアン: -π/2 <= ret <= π/2
atan2(y, x)任意ラジアン: -π <= ret <= π  [x=0かつy=0のときは,0が返る]

モジュロ演算

モジュロ演算(剰余演算)では,次の式が成立します。



足し合わせた結果の剰余演算は,それぞれの要素の剰余演算の結果を足し合わせたものと同じ値となります。

グレイ符号・ビットカウント

最後はあまり数学的なトピックではありませんが,グレイ符号とビットカウントについて記しておきます。過去のSRMで調べておかなければと思った項目です。

▼グレイ符号
グレイ符号はビット符号化方式であり,隣り合う値がハミング距離で1しか離れない符合となります。詳細はWikipediaに譲ります。

unsigned int GrayCode(unsigned int val)
{
    return val^(val>>1);
}

▼ビットカウント
数値の2進表現において,1のビットがいくつあるかを数え上げます。cell(spu)だと,cntbですね。

検索してみると次のコードが見つかります。

unsigned int cntb(unsigned int val)
{
    val =  (val & 0x55555555) + (val >> 1 & 0x55555555);
    val =  (val & 0x33333333) + (val >> 2 & 0x33333333);
    val =  (val & 0x0f0f0f0f) + (val >> 4 & 0x0f0f0f0f);
    val =  (val & 0x00ff00ff) + (val >> 8 & 0x00ff00ff);
    return (val & 0x0000ffff) + (val >>16 & 0x0000ffff);
}
最初は1ビット同士の足し算,その結果が2ビットで表されるので続いて2ビット同士の足し算を行っています。これを繰り返し,最後は16ビットの足し算を行って1のビットの合計が求まります。

参考文献

C言語による最新アルゴリズム事典
Mathematics for TopCoders : Topcoder Algorithm Tutorial

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

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

動的計画法 その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ヶ月以上更新のないブログに表示されています。新しい記事を書くことで広告を消せます。