Artificial Bee Colony アルゴリズムを使って関数の最小値を求めてみた

#群知能アルゴリズムとは
群知能アルゴリズムとは生物のふるまいをアルゴリズム化したものです。よくこのアルゴリズムを用いて関数の最適化を行ったりします。
このアルゴリズムで有名なのが、鳥の群れのふるまいをアルゴリズム化した粒子群最適化法（PSO）と呼ばれているものだと思います。他にも様々にあり、アリのふるまいをアルゴリズム化した人工アリコニー法(ACO)、ミツバチの行動をアルゴリズム化した人工蜂コロニー法(ABC法)、かっこうのふるまいをアルゴリズム化したもの、魚のふるまいをアルゴリズム化したものなどなどいろいろあります。ここではその中でもミツバチのふるまいをアルゴリズム化した手法であるArtificial Bee Colony(ABC)法を紹介します。この手法、アルゴリズムがそんなに複雑でないわりになかなか良い性能を発揮し、粒子群最適化法よりも比較的良い性能を出します。
以下アルゴリズムです。

#人工蜂コロニーアルゴリズム

##ステップ1　初期化
まず、各探索点の初期位置をランダムに設定します。
設定後最良個体を求めて$x_{best}$とします。次の式となります。

\begin{align}
x_{ij} &= \psi_{ij} \\
i_{b} &= \arg \min_{i \in N} f_{i}(x_{i}) \\
x_{best} &= x_{i_{b}}
\end{align}

$i$は個体番号, $i_{b}$は最良個体番号, $x_{best}$及び$x_{i_{b}}$は最良個体位置を表します。 $\psi_{ij}$は$[-C, C]$の範囲で発生させる一様乱数とします. $C$は定数で、評価関数の探索領域を設定します。

##ステップ2　Employ Bee

以下の式により探索を実施します。

\begin{align}
 v^{g}_{ij} &= x^{g}_{ij} + \phi^{g}_{ij}(x^{g}_{ij} - x^{g}_{mj}) \\
 x^{g}_{ij} &= 
	\begin{cases}
	  v^{g}_{ij} & f_{i}(v^{g}_{i}) \leq f_{i}(x^{g}_{i}) \\
	  x^{g}_{i} & \text{それ以外}
        \end{cases} \\
 C_{i} &= 
	\begin{cases}
	0 & f_{i}(v^{g}_{i}) \leq f_{i}(x^{g}_{i}) \\
	C_{i}+1 & \text{それ以外}
	\end{cases}
\end{align}

$i$は個体番号, $j$はランダムに決定された次元番号, $g$は世代数, $m$はランダムに決定された個体番号を表します. $\phi$は$[-1, 1]$の範囲で発生する一様乱数, $x_{i}$は個体$i$の位置,$v_{ij}$は
$x_{ij}$の更新候補個体, $f_{i}$は個体$i$の適応度, $C_{i}$は連続して個体$i$の位置が更新されなかった回数を表します.

##ステップ3　Onlooker Bee

個々体の相対適応度を算出し, ルーレット戦略により更新候補個体を選択し, 以下の式により位置を更新します.

\begin{align}
 p^{g}_{i} &= \dfrac{f^{g}_{i}}{\sum^{N}_{n=1}f^{g}_{n}} \\
 P^{g}_{i} &= \sum_{j=1}^{i}p^{g}_{j} \\
 c &= \displaystyle{\arg_{i \in N}}( P^{g}_{i} \leq \zeta \leq P^{g}_{i+1}) \\
 v^{g}_{cj} &= x^{g}_{cj} + \phi^{g}_{cj}(x^{g}_{cj} - x^{g}_{mj}) \\
 x^{g}_{cj} &= 
    \begin{cases}
	v^{g}_{cj} & f_{c}(v^{g}_{c}) \leq f_{c}(x^{g}_{c}) \\
	x^{g}_{c} & \text{それ以外}
    \end{cases} \\
 C_{c} &= 
	\begin{cases}
	    0 & f_{c}(v^{g}_{c}) \leq f_{c}(x^{g}_{c}) \\
	    C_{c}+1 & \text{それ以外}
	\end{cases}
\end{align}

$i$は個体番号, $j$はランダムに決定された次元番号, $g$は世代数, $m$はランダムに決定された個体番号を表す. $\phi$は$[-1, 1]$の範囲で発生させる一様乱数, $x_{i}$は個体$i$の位置,
$v_{ij}$は$x_{ij}$の更新候補個体, $f_{i}$は個体$i$の適応度, $C_{c}$は連続して個体$c$の位置が更新されなかった回数を表す.
$c$はルーレット選択により選択した個体番号を表す. $\zeta$は$[0,1]$の一様乱数とする.
$p_{i}$は世代$g$における相対適応度を表し, $P_{i}$はルーレット戦略における選択確率を表す. $f_{c}$は世代$g$における適応度を表す.

##ステップ4　Scout Bee

更新が一定回数なされなかった個体に対して以下の式に基づいて新たな位置を探索します.

\begin{equation}
 x_{ij} = 
	\begin{cases}
		x_{min}+\psi(x_{maxj}-x_{minj}) & C_{limit} \leq C_{i} \\
		x_{ij} & \text{それ以外}
	\end{cases}
\end{equation}

$x_{maxj}$は探索個体の中で最も最大となる位置を表し, $x_{minj}$はもっとも小さい値を表します. $\psi$は$[0, 1]$までの一様乱数を表す. $C_{i}$は連続して個体$i$の位置が更新されなかった回数を表し, $C_{limit}$はその閾値とします.

##ステップ5　最適値の更新

現世代と現世代までに得られた最良個体を比較して, 現時点での最良個体を更新します.最良個体の適応度が最適解に収束したならば終了し, していないならばステップ2へ戻ります.

\begin{align}
	i_{b} &= \arg \min_{i \in N} f(x_{i}) \\
	x_{best} &= 
	\begin{cases}
		x_{i_{b}} & f_{i_{b}}(x_{i_{b}}) \leq f_{best}(x_{best}) \\
		x_{best} & \text{それ以外}
	\end{cases}
\end{align}

$i$は個体番号, $i_{b}$は最良個体番号, $x_{i_{b}}$は現世代での最良個体, $x_{best}$は最良個体位置を表します. $x_{best}$は前世代までの最良個体を表します.

#実際にやってみた結果

簡単のため, 一番簡単なsphere関数で試してみました.
sphere関数とは以下の関数です.
最適値は$f(x) = 0$で,
最適値のときの値は$(x_{1}, x_{2}, …, x_{N}) = (0, 0, …, 0)$です.

f(x) = \sum^{N}_{i=1}x_{i}^{2}

10次元と, 30次元, 50次元で試してみました。収束状況は以下の通りです.だいたいこのくらいで収束するみたいですね。グラフの縦軸が関数の値で, 横軸が繰り返し計算回数です.

10次元の結果

30次元の結果

50次元の結果

#ソースコード(Java)
計算に使用したソースコードは以下の通りです。乱数に関してはメルセンヌツイスターのライブラリである、
Sfmt.javaを使用しています。詳細はこちら。
http://www001.upp.so-net.ne.jp/isaku/rand2.html

AbcTestSimple.java

public class AbcTestSimple 
{
	public static void main( String[] args)
	{
		int iGenerationNum = 100000;
		int iAbcDataNum = 60;
		int iAbcSearchNum = 30;
		int iAbcVectorDimNum = 300;
		int iAbcLimitCount = 180;
		double lfSolveRange = 5.12;
		
		// 初期化を行います。
		Abc abc =  new Abc( iGenerationNum, iAbcDataNum, iAbcVectorDimNum, iAbcSearchNum, iAbcLimitCount, lfSolveRange);
		
		for( int i = 0;i < iGenerationNum; i++ )
		{
			// ABC法の計算を行います。
			abc.vAbc();
			// 現時点での関数の最小値を出力します。
			abc.vOutputGlobalMinAbcDataConstFuncValue();
		}
	}
}

Abc.java

public class Abc 
{
	private int iGenerationNumber;				// 計算回数
	private int iAbcDataNum;					// コロニーの数
	private int iAbcVectorDimNum;				// 蜂の特徴ベクトル
	private int iAbcSearchNum;					// 探索点数
	private int iAbcLimitCount;					// 更新しなかった回数
	private double[][] pplfAbcData;				// ABCデータ配列
	private double lfGlobalMaxAbcData;			// 大域的最適値
	private double lfGlobalMinAbcData;			// 大域的最小値
	private double[] plfGlobalMaxAbcData;		// 大域的最適解を表す粒子のデータ
	private double[] plfGlobalMinAbcData;		// 大域的最適解を表す粒子のデータ
	private double[][] pplfVelocityData;		// ABCデータ速度配列
	private double[] plfFitProb;				// 適合度相対確率
	private double[] plfFit;					// 適合度
	private int[] piNonUpdateCount;				// 更新しない回数
	private Sfmt rnd;
	private double lfSolveRange;				// 探索範囲

	/**
	 * <PRE>
	 * 　人口蜂コロニーの初期化を実行します。
	 * </PRE>
	 * @param iGenCount       計算回数
	 * @param iGenNum         コロニー数
	 * @param iGenVectorDim   探索点の次元数
	 * @param iSearchNum  　  探索点の数
	 * @param iLimitCountData 更新しなかった回数
	 * @param iIntervalMinNum 最低反復回数
	 * @param iAlphaData      探索点上位数
	 * @param lfDrData        収束状況判定
	 * @param lfBoundData     適合度許容限界値
	 * @param lfAccuracyData  適合度評価精度
	 */
	public Abc( int iGenCount, int iGenNum, int iGenVectorDim, int iSearchNum, int iLimitCountData, double lfSolveRangeData )
	{
		int i,j;

		iGenerationNumber	= iGenCount;
		iAbcDataNum			= iGenNum;
		iAbcVectorDimNum	= iGenVectorDim;
		iAbcSearchNum		= iSearchNum;
		iAbcLimitCount		= iLimitCountData;
		lfSolveRange		= lfSolveRangeData;
		lfGlobalMaxAbcData = -Double.MAX_VALUE;			// 大域的最適値
		lfGlobalMinAbcData = Double.MAX_VALUE;			// 大域的最小値

		long seed;
		seed = System.currentTimeMillis();
		rnd = new Sfmt( (int)seed );

		pplfAbcData 		= new double[iAbcDataNum][iAbcVectorDimNum];
		pplfVelocityData	= new double[iAbcDataNum][iAbcVectorDimNum];
		plfFit				= new double[iAbcSearchNum];
		plfFitProb 			= new double[iAbcSearchNum];
		piNonUpdateCount 	= new int[iAbcSearchNum];
		plfGlobalMaxAbcData = new double[iAbcVectorDimNum];
		plfGlobalMinAbcData = new double[iAbcVectorDimNum];

		for( i= 0;i < iAbcDataNum; i++ )
		{
			for(j = 0;j < iAbcVectorDimNum; j++ )
			{
				pplfAbcData[i][j] = 0.0;
				pplfVelocityData[i][j] = 0.0;
			}
		}
		for( i = 0;i < iAbcSearchNum ; i++ )
		{
			plfFit[i] = 0.0;
			plfFitProb[i] = 0.0;
			piNonUpdateCount[i] = 0;
		}
		for( i = 0;i < iAbcVectorDimNum; i++ )
		{
			plfGlobalMaxAbcData[i] = 0.0;
			plfGlobalMinAbcData[i] = 0.0;
		}
		
		for( i= 0;i < iAbcDataNum; i++ )
		{
			for(j = 0;j < iAbcVectorDimNum; j++ )
			{
				pplfAbcData[i][j] = lfSolveRange*(2.0*rnd.NextUnif()-1.0);;
			}
		}
	}

	/**
	 * <PRE>
	 * 　人工蜂コロニー最適化法を実行します。
	 * </PRE>
	 */
	public void vAbc()
	{
		// employee bee の動作
		vEmployBeeOrigin();

		// onlookers beeの動作
		vOnlookerBeeOrigin();

		// scout bee の実行
		vScoutBeeOrigin();

		// 更新します。
		vGetGlobalMaxMin();
	}

	/**
	 * <PRE>
	 *	 　Employ Beeを実行します。
	 * </PRE>
	 */
	private void vEmployBeeOrigin()
	{
		int m,h;
		int i,j;
		double lfRand = 0.0;
		double lfFunc1 = 0.0;
		double lfFunc2 = 0.0;
	// employee bee の動作
	
		// 更新点候補を算出します。
		m = rnd.NextInt(iAbcSearchNum);
		h = rnd.NextInt(iAbcVectorDimNum);
	
		for( i = 0;i < iAbcSearchNum; i++ )
		{
			lfRand = 2*rnd.NextUnif()-1;
			for( j = 0; j < iAbcVectorDimNum; j++ )
				pplfVelocityData[i][j] = pplfAbcData[i][j];
			pplfVelocityData[i][h] = pplfAbcData[i][h] + lfRand*( pplfAbcData[i][h] - pplfAbcData[m][h] );
	
			// 各探索点と更新しなかった回数を格納する変数を更新します。
			lfFunc1 = lfSphere( pplfVelocityData[i] );
			lfFunc2 = lfSphere( pplfAbcData[i] );
	
			if( lfFunc1 < lfFunc2 )
			{
				pplfAbcData[i][h] = pplfVelocityData[i][h];
				piNonUpdateCount[i] = 0;
			}
			else	piNonUpdateCount[i] = piNonUpdateCount[i] + 1;
		}
	}

	/**
	 * <PRE>
	 *   Onlooker Beeを実行します。
	 * </PRE>
	 */
	private void vOnlookerBeeOrigin()
	{
		int i,j;
		int c,m,h;
		double lfRes = 0.0;
		double lfRand = 0.0;
		double lfFitProb = 0.0;
		double lfProb = 0.0;
		double lfPrevProb = 0.0;
		double lfFunc1 = 0.0;
		double lfFunc2 = 0.0;
		double lfObjFunc = 0.0;

		lfRes = 0.0;
		for(j = 0;j < iAbcSearchNum; j++ )
		{
			// 適応度の算出
			lfObjFunc = lfSphere( pplfAbcData[j] );
			if( lfObjFunc >= 0.0 )	lfFitProb = 1.0/( 1.0+lfObjFunc );
			else			lfFitProb = 1.0+Math.abs( lfObjFunc );
			lfRes += lfFitProb;
			plfFit[j] = lfFitProb;
		}
		// 適応度の正規化
		for( j = 0;j < iAbcSearchNum; j++ )	plfFitProb[j] = plfFit[j]/lfRes;
		// ルーレット戦略を実行
		lfProb = lfPrevProb = 0.0;
		lfRand = rnd.NextUnif();
		c = 0;
		for( j = 0;j < iAbcSearchNum; j++ )
		{
			lfProb += plfFitProb[j];
			if( lfPrevProb <= lfRand && lfRand <= lfProb )	c = j;
			lfPrevProb = lfProb;
		}
		// ルーレット選択した探索点に対して更新候補を求めます。

		// 更新点候補を算出します。
		// 更新点候補を乱数により決定します。
		m = rnd.NextInt(iAbcSearchNum);
		h = rnd.NextInt(iAbcVectorDimNum);

		lfRand = 2*rnd.NextUnif()-1;
		for( j = 0; j < iAbcVectorDimNum; j++ )
			pplfVelocityData[c][j] = pplfAbcData[c][j];
		pplfVelocityData[c][h] = pplfAbcData[c][h] + lfRand*( pplfAbcData[c][h] - pplfAbcData[m][h] );

		// 更新点候補を次のように更新します。
		lfFunc1 = lfSphere( pplfVelocityData[c] );
		lfFunc2 = lfSphere( pplfAbcData[c] );
		if( lfFunc1 < lfFunc2 )
		{
			for( j = 0;j < iAbcVectorDimNum; j++ )
				pplfAbcData[c][j] = pplfVelocityData[c][j];
			piNonUpdateCount[c] = 0;
		}
		else	piNonUpdateCount[c] = piNonUpdateCount[c] + 1;
	}

	/**
	 * <PRE>
	 * Scout Beeを実行します。
	 * </PRE>
	 */
	private void vScoutBeeOrigin()
	{
		int i,j,k;
		double lfRand = 0.0;
		// 新たな探索点を求めて探索を実行します。
		for( i = 0;i < iAbcSearchNum; i++ )
		{
			if( piNonUpdateCount[i] > iAbcLimitCount )
			{
				for( k = 0;k < iAbcVectorDimNum; k++ )
				{
					lfRand = rnd.NextUnif();
					pplfAbcData[i][k] = lfSolveRange*(2.0*lfRand-1.0);
				}
			}
		}
	}
	
	/**
	 * <PRE>
	 * 現時点での目的関数の全体を通しての最大、最小値を求めます。
	 * </PRE>
	 */
	private void vGetGlobalMaxMin()
	{
		int i,j;
		int iMinLoc = 0;
		double lfObjFunc = 0.0;

		for( i = 0;i < iAbcSearchNum; i++ )
		{
			lfObjFunc = lfSphere( pplfAbcData[i] );
			if( lfGlobalMinAbcData >= lfObjFunc )
			{
				iMinLoc = i;
				lfGlobalMinAbcData = lfObjFunc;
				for( j = 0; j < iAbcVectorDimNum; j++ )
					plfGlobalMinAbcData[j] = pplfAbcData[i][j];
			}
			if( lfGlobalMaxAbcData <= lfObjFunc )
			{
				iMinLoc = i;
				lfGlobalMaxAbcData = lfObjFunc;
				for( j = 0; j < iAbcVectorDimNum; j++ )
					plfGlobalMaxAbcData[j] = pplfAbcData[i][j];
			}
		}
	}

	/**
	 * <PRE>
	 * 　現時点でのもっともよい粒子の位置を出力します。
	 * </PRE>
	 */
	public void vOutputGlobalMinAbcData()
	{
		int i;
		// 現時点での各粒子の目的関数の値を出力します。
		for( i = 0; i < iAbcVectorDimNum; i++ )
		{
			System.out.print( plfGlobalMinAbcData[i] + "," );
		}
		System.out.println("");
	}

	/**
	 * <PRE>
	 * 　現時点でのもっともよい粒子の目的関数値を出力します。
	 * </PRE>
	 * @return 関数値
	 */
	public void vOutputGlobalMinAbcDataConstFuncValue()
	{
		// TODO 自動生成されたメソッド・スタブ
		System.out.println( lfSphere( plfGlobalMinAbcData ) + "," );
	}
	/**
	 * <PRE>
	 *   Sphere関数を実行します。
	 * </PRE>
	 * @param plfData ベクトルデータ
	 * @return Sphere関数の値
	 */
	private double lfSphere( double[] plfData )
	{
		double lfRes = 0.0;
		for(int i = 0;i < plfData.length; i++ )
			lfRes += plfData[i]*plfData[i];
		return lfRes;
	}
}	

#参考
もっとも初期に提案されたものに近い手法を紹介しましたが, 研究が進んでいて
より高速にかつ高次元に最適値に到達できる改良されたABCの手法が提案されています.

#参考文献
宇谷昭秀, 長島淳也, 牛膓 隆太, 山本尚生, "Artificial Bee Colony (ABC) アルゴリズムの高次元問題に対する解探索性能の強化", 電子情報通信学会論文誌 D, Vol. J94-D No.2, 2011, pp.425-438.