粒子群最適化

最適化手法について調べていたところ，wikipediaに粒子群最適化という方法があったので，実装してみた．

wikipediaの記述どおりにざっくりと実装．パラメタの次元数をテンプレート引数にしてある．

//[ParticleSwarmOpt.h]
//	粒子群最適化

#pragma once

#include <vector>
#include <array>

//"粒子群最適化"による目的関数の最小化
//	N_PARAM : 最適化対象パラメタの個数(次元数)
//
//[使い方]
//	コーディング：
//		これを継承して，仮想関数を実装する．
//	最適化計算：
//		1) Initialize()で初期化
//		2) Update()を任意回コールして最適化
//		3) GetCurrBestParam()で結果パラメタ値を取得
template< size_t N_PARAM >
class PSO
{
public:
	//パラメタベクトルXおよび速度Vの型
	typedef std::array<double,N_PARAM> Vec_t;
protected:
	//粒子情報(CreateInitialState()引数用)
	struct ParticleSetting
	{
		Vec_t X0;	//初期位置
		Vec_t V0;	//初期速度
		std::vector< unsigned int > GroupIndexesBelong;	//粒子が所属グループ群のindex
	};
private:
	//粒子
	struct Particle
	{
		Vec_t X;	//現在位置
		Vec_t V;	//速度
		Vec_t BestX;	//この粒子の，過去最良目的関数値になった位置
		double BestOFV;	//この粒子，過去最良目的関数値
		std::vector< struct Group* > Groups;	//この粒子が所属するグループ

		Particle(){	BestOFV = std::numeric_limits<double>::max();	}
	};

	//グループ
	struct Group
	{
		Vec_t BestX;	//このグループの，過去最良目的関数値になった位置
		double BestOFV;	//このグループの，過去最良目的関数値
		
		Group(){	BestOFV = std::numeric_limits<double>::max();	}
	};

	//-----------------------------------------------------
public:
	PSO()
	{
		SetParam( 1.0, 1.0, 1.0 );
		m_GlobalBestOFV = std::numeric_limits<double>::max();
	}

public:
	//処理パラメタ値のSet,Get
	void SetParam( double W, double C1, double C2 )
	{	m_W=W;	m_C1=C1;	m_C2=C2;	}

	double W() const {	return m_W;	}
	double C1() const {	return m_C1;	}
	double C2() const {	return m_C2;	}

	//初期状態の生成.
	//最適化処理の最初の状態を生成する．
	//Update()の繰り返しを行う前に，コールする必要がある
	bool Initialize();

	//最適化処理の1stepを実行．
	//Initialize()に成功したら，気が済むまでこれを繰り返しコールすることで最適化処理を進行させる．
	//[Ret]
	//	最良パラメタ（GetCurrBestParam()が返す値）が更新された場合はtrue．
	bool Update();

	//結果取得
	Vec_t GetCurrBestParam() const {	return m_GlobalBestX;	}
	double GetObjectiveFuncVal_at_CurrBestParam() const {	return m_GlobalBestScore;	}


	//粒子位置情報取得(状態見る用)
	void GetCurrParamsOfPartilces( std::vector<Vec_t> &rDst )
	{
		rDst.resize( m_Particles.size() );
		for( size_t i=0; i<m_Particles.size(); ++i )
		{	rDst[i] = m_Particles[i].X;	}
	}
	//粒子速度情報取得(状態見る用)
	void GetCurrVelocityOfPartilces( std::vector<Vec_t> &rDst )
	{
		rDst.resize( m_Particles.size() );
		for( size_t i=0; i<m_Particles.size(); ++i )
		{	rDst[i] = m_Particles[i].V;	}
	}

	//-----------------------------------------------------
protected:
	//乱数で 0<=ret<=1 な値を返す
	virtual double Rand0to1() = 0;

	//パラメタ値Xに対する目的関数値f(X)の計算．
	//[Args]
	//	X : 位置(パラメタ値)
	//[Ret]
	//	位置Xにおける目的関数値を返す．
	virtual double ObjectiveFunc( const Vec_t &X ) = 0;

	//初期化．
	//粒子群の初期状態と，グループの個数を引数に設定して返す．
	//[Args]
	//	rDstParticles : 粒子群初期状態受取
	//	rDst_nGroup : グループ個数受取．（正常な値は1以上）
	//[Ret]
	//	成功時はtrue，失敗時はfalse
	virtual bool CreateInitialState( std::vector<ParticleSetting> &rDstParticles, unsigned int &rDst_nGroup ) = 0;

	//-----------------------------------------------------
private:
	//粒子の移動
	void MoveParticle( Particle &P )
	{
		//粒子が属するグループ群の中で，最も良い結果を探索
		Vec_t GroupX;
		{
			double OFV = std::numeric_limits<double>::max();
			for( auto pG : P.Groups )
			{
				if( pG->BestOFV < OFV )
				{
					OFV = pG->BestOFV;
					GroupX = pG->BestX;
				}
			}
		}
		//速度の更新と移動
		for( size_t i=0; i<N_PARAM; ++i )
		{
			P.V[i] *= m_W;
			P.V[i] += m_C1 * Rand0to1() * ( P.BestX[i] - P.X[i] );
			P.V[i] += m_C2 * Rand0to1() * ( GroupX[i] - P.X[i] );

			P.X[i] += P.V[i];
		}
	}

	//目的関数値の更新
	//グローバルな最良値が更新されたらtrueを返す
	bool UpdateOFV( Particle &P )
	{
		double OFV = ObjectiveFunc( P.X );
		if( OFV >= P.BestOFV )
		{	return false;	}

		P.BestOFV = OFV;
		P.BestX = P.X;
		bool ret = false;

		for( auto pG : P.Groups )
		{
			if( OFV < pG->BestOFV )
			{
				pG->BestOFV = OFV;
				pG->BestX = P.BestX;

				if( OFV < m_GlobalBestOFV )
				{
					m_GlobalBestOFV = OFV;
					m_GlobalBestX = P.BestX;
					ret = true;
				}
			}
		}
		return ret;
	}

private:
	std::vector< Particle > m_Particles;	//粒子群
	std::vector< Group > m_Groups;	//グループ群
	Vec_t m_GlobalBestX;	//粒子群全体の，過去最良目的関数値になった位置
	double m_GlobalBestOFV;	//粒子群全体の，過去最良目的関数値

	//速度Vの更新式のパラメタ
	double m_W;	//慣性定数
	double m_C1, m_C2;	//粒子が自身のベスト位置方向にいく度合，グループのベスト位置に行く度合
};

//Initialize()
template< size_t N_PARAM >
bool PSO<N_PARAM>::Initialize()
{
	//初期情報の取得
	std::vector<ParticleSetting> PSs;
	unsigned int nGroups = 0;
	if( !CreateInitialState( PSs, nGroups ) )
	{	return false;	}

	//返ってきたデータのチェック
	if( nGroups==0 || PSs.empty() )
	{	return false;	}

	for( const auto &PS : PSs )
	{
		if( PS.GroupIndexesBelong.empty() )
		{	return false;	}

		for( unsigned int g : PS.GroupIndexesBelong )
		{
			if( g >= nGroups )
			{	return false;	}
		}
	}

	//既存データ破棄
	m_Groups.clear();
	m_Particles.clear();
	m_GlobalBestOFV = std::numeric_limits<double>::max();

	//データ生成
	m_Groups.resize( nGroups );
	m_Particles.resize( PSs.size() );

	auto iP = m_Particles.begin();
	for( const auto &PS : PSs )
	{
		iP->X = iP->BestX = PS.X0;
		iP->V = PS.V0;
		iP->Groups.reserve( PS.GroupIndexesBelong.size() );
		for( unsigned int g : PS.GroupIndexesBelong )
		{
			iP->Groups.push_back( &(m_Groups[g]) );
		}

		UpdateOFV( *iP );
		iP++;
	}

	return true;
}

//Update()
template< size_t N_PARAM >
bool PSO<N_PARAM>::Update()
{
	//粒子群の移動
	for( auto &P : m_Particles )
	{	MoveParticle( P );	}

	//目的関数計算と，｛グループ，グローバル｝情報更新
	bool bGlobalUpdated = false;
	for( auto &P : m_Particles )
	{
		bGlobalUpdated |= UpdateOFV( P );
	}
	return bGlobalUpdated;
}

				
で，２次元パラメタ空間で動作テストしてみる．つまらないテストコードは以下．
				
#include "ParticleSwarmOpt.h"

#include "opencv2/core/core.hpp"
#include "opencv2/highgui/highgui.hpp"
#include "opencv2/imgproc/imgproc.hpp"

#include <random>
#include <iostream>

class TestOpt : public PSO<2>
{
public:
	TestOpt()
		: m_RandomEngine( std::random_device()() )
	{
	}

public:
	void Show()
	{
		cv::cvtColor( m_FuncMap, m_ShowImg, CV_GRAY2BGR );

		std::vector< Vec_t > Xs, Vs;
		GetCurrParamsOfPartilces( Xs );
		GetCurrVelocityOfPartilces( Vs );
		
		size_t n = Xs.size();
		for( size_t i=0; i<n; ++i )
		{
			cv::Point P0( cvRound(Xs[i][0]), cvRound(Xs[i][1]) );
			cv::circle( m_ShowImg, P0, 2, cv::Scalar(0,255,32), -1 );

			cv::Point P1( cvRound(Xs[i][0]+Vs[i][0]), cvRound(Xs[i][1]+Vs[i][1]) );
			cv::line( m_ShowImg, P0,P1, cv::Scalar(255,32,32) );
		}

		Vec_t BestX = GetCurrBestParam();
		cv::circle( m_ShowImg, cv::Point( cvRound(BestX[0]),cvRound(BestX[1]) ), 3, cv::Scalar(0,0,255), -1 );

		cv::imshow( "PSO<2>", m_ShowImg );
	}

protected:
	//
	virtual double ObjectiveFunc( const Vec_t &X ) override
	{
		int x = cvRound( X[0] );
		int y = cvRound( X[1] );
		if( x<0 || x>=m_FuncMap.cols || y<0 || y>=m_FuncMap.rows )
		{	return 1000;	}	//※適当にでかい値

		return m_FuncMap.at<unsigned char>( y,x );
	}

	//
	virtual double Rand0to1() override
	{	return std::uniform_real_distribution<>(0.0,1.0)( m_RandomEngine );	}

	//
	virtual bool CreateInitialState( std::vector<ParticleSetting> &rDstParticles, unsigned int &rDst_nGroup ) override
	{
		const int W = 320;
		const int H = 240;
		{//目的関数の具合をてきとーにつくる
			m_FuncMap.create( H,W, CV_8UC1 );
			for( int y=0; y<H; ++y )
			{
				unsigned char *p = m_FuncMap.ptr<unsigned char>( y );
				for( int x=0; x<W; ++x, ++p )
				{
					*p = cvRound( Rand0to1()*128 );
				}
			}

			int x = cvRound( W*0.1 + W*0.8*Rand0to1() );
			int y = cvRound( H*0.1 + H*0.8*Rand0to1() );
			cv::circle( m_FuncMap, cv::Point(x,y), 7, cv::Scalar(0), -1 );
			//m_FuncMap.at<unsigned char>(y,x) = 0;

			cv::medianBlur( m_FuncMap, m_FuncMap, 15 );

			std::cout << "正解の箇所 = (" << x << ", " << y << ")" << std::endl;
		}

		//
		const unsigned int N = 100;
		rDst_nGroup = N;
		rDstParticles.resize( N );

		const double PI = acos(-1.0);
		for( unsigned int i=0; i<N; ++i )
		{
			ParticleSetting &P = rDstParticles[i];

			P.X0[0] = Rand0to1() * (W-1);
			P.X0[1] = Rand0to1() * (H-1);
			double V = Rand0to1()*5;
			double theta = 2*PI*Rand0to1();
			P.V0[0] = V * cos( theta );
			P.V0[1] = V * sin( theta );

			P.GroupIndexesBelong.reserve( 3 );
			P.GroupIndexesBelong.push_back( i>1  ?  i-1  :  N-1 );
			P.GroupIndexesBelong.push_back( i );
			P.GroupIndexesBelong.push_back( i+1<N  ?  i+1  :  0 );
		}

		//
		m_ShowImg.create( H,W, CV_8UC3 );
		return true;
	}

private:
	std::mt19937 m_RandomEngine;
	cv::Mat m_FuncMap;
	cv::Mat m_ShowImg;
};

//
int _tmain(int argc, _TCHAR* argv[])
{
	TestOpt Test;

	Test.SetParam( 0.8, 1, 1 );
	if( !Test.Initialize() )return 0;

	std::cout << "Initial" << std::endl;
	Test.Show();
	cv::waitKey();

	unsigned int iter = 0;
	do
	{
		++iter;
		std::cout << "step " << iter;
		if( Test.Update() ){	std::cout << " *";	}
		std::cout << std::endl;
		Test.Show();
	}
	while( cv::waitKey(0)!='q' );

	return 0;
}

				
適当に乱数でグレースケール画像を作り，２次元パラメタ空間上の目的関数だということにする．黒いところほど関数値が小さい．

				・緑の●が粒子．

				・粒子から出ている青い線が粒子の速度

				・赤●はこれまでに見つけた最良のパラメタの位置．
				


				初期状態の粒子をランダムに配置した様子．

				画像左下付近のブラックホールっぽい個所が目的関数が最小になる場所であるが，初期配置では粒子がそのあたりに存在しないため，ローカルミニマムな箇所に赤●が表示されている．

				問題は粒子のグループ設定をどうするかであるが，全粒子が１つのグループというのだとなんとなくつまらないので，テストコードを見てわかるように，１粒子が３グループに属する形にしてみた．これだと，ある粒子がこれまでより良い場所を見つけても，別グループの粒子にその情報が伝達するのにある程度時間がかかるハズ．

				この状態から最適化処理開始．とりあえず1000step走らせてみた．

				↓

				

				2step目で早くも粒子の一つが偶然最小値付近に到達．

				移動速度がとても大きい粒子が多数存在しているが，元ネタの式に忠実だと多分こうなるのが正解であろう．（用途次第では最大速度制限などしてもよいか？）

				↓

				

				step6．さらに解が更新された．

				情報伝達に時間がかかるため（全粒子にこの赤●の情報がいきわたるには50stepくらい必要になる），まだ赤●の箇所に粒子が集まってくる気配はない．

				↓

				

				step108．だいぶ粒子が赤●付近に集まっているが，他のいくつかのローカルミニマムな箇所にも粒子が集まっている．

				↓

				

				step550．多くの粒子が赤●付近に吸い寄せられた．

				↓

				

				step1000．まだパラメタ空間上を飛びまわって探索している粒子が残っているのが頼もしい．
				
粒子位置の目的関数の勾配等を見ることなく粒子の移動速度が決定されるので，乱数と粒子数頼みではあるが，パラメタ空間の広範囲を探索してくれることが期待できそう．

				目的関数のJacobianやHessianが無くても動かせるので，滑降シンプレックス法と同様，「とりあえずの解」が手早く欲しいときに使えるかも．
				
ただし，終了条件には工夫が要りそうに思う．

				step数や，解の更新度合といった条件だけだと，不十分な状態で止まってしまいそうな気がするので，粒子の分布の広がりのようなものを見るなどすべきかもしれない．
				
by nakiusagi3