[プログラミング] ビット並列アルゴリズムを使った編集距離

ふと、ビット並列アルゴリズムを使った編集距離を計算するアルゴリズムを書きたくなったので書いてみた。
まず、通常の編集距離であるLevenshtein Distanceを求めるアルゴリズムは以下のように書ける

	int levenshteinDistance(String A, String B) {
		int m = A.length();
		int n = B.length();
		int dp[] = new int[n + 1];
		int next[] = new int[n + 1];
		for (int i = 0; i <= n; i++)
			dp[i] = i;
		for (int i = 1; i <= m; i++) {
			next[0] = i;
			for (int j = 1; j <= n; j++) {
				if (A.charAt(i - 1) == B.charAt(j - 1)) {
					next[j] = dp[j - 1];
				} else {
					int r = dp[j - 1];
					r = min(r, dp[j]);
					r = min(r, next[j - 1]);
					next[j] = r + 1;
				}
			}
			int tmp[] = dp;
			dp = next;
			next = tmp;
		}
		return dp[n];
	}

このアルゴリズムの時間計算量はO(nm)、空間計算量はO(n)となっている。(ここでn,mはそれぞれ文字列の長さ)

この問題に対して、1999年にMyersがビット並列手法を用いてO(floor(m/w) n)の方法を提案した[1]。また2003年にHyyröがO(floor(m/w) d)(dは2つの文字列の編集距離)の方法を提案している[2]。

ビット並列手法は通常のCPUが32ビットや64ビットの数値を1命令で扱えることを利用して複数の計算を一度に実行する手法である。編集距離の問題においてはDPの各値が列および行が1つしかずれていないときに高々1しかずれていないということを利用してやることができる。これに関しては文献[3]などに詳しい。また[4]のHyyröの論文リストにはビット並列手法を利用した文字列照合アルゴリズムに関して多くの論文がある。

以下に文献[2]のFigure 3を実装したコードを示す。下のコードではm <= 64を仮定している。これに関しては例えばスペルチェックなどの応用においては64文字以下で十分なことが多く、場合によっては以下のコードでも十分である。

	void createPM(String query, long PM[]) {
		for (int i = 0; i < query.length(); i++) {
			char c = query.charAt(i);
			PM[c] |= 1l << i;
		}
	}
	
	int levenshteinDistance(int m, long PMs[], String B) {
		int D = m;
		long D0, HP, HN, VP, VN;
		long top = 1l << (m - 1);
		VP = ~0; VN = 0;
		for (char c : B.toCharArray()) {
			long PM = PMs[c];
			D0 = (((PM & VP) + VP) ^ VP) | PM | VN;
			HP = VN | ~(D0 | VP);
			HN = D0 & VP;
			if ((HP & top) != 0) {
				++D;
			} else if ((HN & top) != 0) {
				--D;
			}
			VP = (HN << 1) | ~(D0 | ((HP << 1) | 1));
			VN = D0 & ((HP << 1) | 1);
		}
		return D;
	}

実験結果

今回実験としてランダムな長さ50の文字列1000個に対して、ランダムな固定長のクエリを100個生成し、編集距離の近いものを見つけるのにかかる時間を測定した。クエリの長さは4から64まで4刻みで変化させた。

結果は以下の図のようになり、クエリ長がある程度長いときにはビット並列化の効果が高くなることが分かる。

その他

また、longを2個利用して128ビットまで対応したものを以下に示す。また、多倍長演算を実装したクラスを作成し、任意の長さの文字列を扱えるプログラムも作成したがプリミティブ型を利用したものに比べ10倍程度遅く、現在速度を改善できないかと考えているところである。

	void createPM(String query, long PMl[] , long PMh[]) {
		for (int i = 0; i < query.length(); i++) {
			char c = query.charAt(i);
			if(i < 64){
				PMl[c] |= 1l << i;				
			}else{
				PMh[c] |= 1l << (i - 64);								
			}
		}
	}

	int levenshteinDistance128(int m, long PMl[] , long PMh[], String B) {
		if(m > 128)return -1;
		int D = m;
		long D0, HP, HN, VP, VN;
		long D0_h, HP_h, HN_h, VP_h, VN_h;
		long top = 1l << (m - 1);
		if(m > 64){
			top = 1l << (m - 65);
		}
		VP = ~0; VN = 0;
		VP_h = ~0; VN_h = 0;
		for (char c : B.toCharArray()) {
			long PM = PMl[c];
			long PM_h = PMh[c];
			long x = (PM & VP); 
			//check x + VP is overflow
			//(cf. Hacker's Delight 2-12)
			long carry = ((x & VP) | ((x | VP) & ~(x + VP))) >>> 63;
			D0 = ((x + VP) ^ VP) | PM | VN;
			D0_h = (((PM_h & VP_h) + VP_h + carry) ^ VP_h) | PM_h | VN_h;
			HP = VN | ~(D0 | VP);  HP_h = VN_h | ~(D0_h | VP_h);
			HN = D0 & VP; HN_h = D0_h & VP_h;
			if(m <= 64){
				if ((HP & top) != 0)++D;
				else if ((HN & top) != 0)--D;				
			}else{
				if ((HP_h & top) != 0)++D;
				else if ((HN_h & top) != 0)--D;	
			}
			long H_h = (HP_h << 1) | (HP >>> 63);
			long H2_h = (HN_h << 1) | (HN >>> 63);
			VP = (HN << 1) | ~(D0 | ((HP << 1) | 1)); VP_h = H2_h | ~(D0_h | H_h);
			VN = D0 & ((HP << 1) | 1); VN_h = D0_h & H_h;
		}
		return D;
	}

参考文献

[1] G. Myers: A fast bit-vector algorithm for approximate string matching based on dynamic progamming. Journal of the ACM, 46(3): 395-415, 1999.
[2] Heikki Hyyrö: A Bit-Vector Algorithm for Computing Levenshtein and Damerau Edit Distances, In the proceedings of the Prague Stringology Conference 2002 (PSC 2002)
[3] Heikki Hyyrö: Extending and Explaining the Bit-Parallel Approximate String Matching Algorithm of Myers, Technical report A2001-10 of the Department of Computer and Information Sciences, University of Tampere.
[4] http://www.cs.uta.fi/~helmu/pubs/pubs.html