heap sortを使った上位ランキング取得プログラム
Managing Gigabytes 4.6章で解説されているソートのプログラムを実装してみた。
検索エンジンなどでN個のデータの中から上位r個を取得したい場合、まずN個のデータからなるmax-heapを構成して、ルート（最大値）から順にr個をヒープから取り除くというアプローチが考えられる。しかしN>>rの場合、r個のデータからなるmin-heapを構成して、残りのN-r個のデータをheapのルート（最小値）と順次比較して、ルートよりも大きい場合はルートと入れ替えて、heapを再構成する、という手順を取った方がより計算量が少なくて済む、という話。
10万件のランダムな数値をスペース区切りで出力するプログラム
#include <iostream> #include <stdlib.h> using namespace std; int main (int argc, char *argv[]) { srand((unsigned int)time(0)); for (int i=0; i<100000; i++) { int tmp = rand()/1000; cout << tmp << " "; } cout << endl; return 1; }
10件のmin-heapを構成して、10万件のデータから上位10件を取得・表示するプログラム
#include <iostream> #include <fstream> #define SWAP(a,b) (a ^= b,b = a ^ b,a ^= b) using namespace std; #define HEAPSIZE 10 #define DATAN 100000 int heapsize = HEAPSIZE; int heap[HEAPSIZE]; int datan = DATAN; void make_heap (int root); int main (int argc, char *argv[]) { srand((unsigned int)time(0)); ifstream input_file("testdata"); int testdata[datan]; for (int i=0; i<datan; i++) { input_file >> testdata[i]; } //copy the first r accumlators into the heap for (int i=0; i<heapsize; i++) { heap[i] = testdata[i]; } //convert array into min-heap for (int i=((heapsize+1)/2); i>0; i--) { make_heap(i-1); } for (int i=heapsize; i<datan; i++) { //testdata[i]; if (heap[0] < testdata[i]) { //discard rheap[0] and insert testdata[i] into min-heap heap[0] = testdata[i]; make_heap(0); } } //display top-r number in the ascending order while (heapsize > 0) { cout << heap[0] << endl; heap[0] = heap[heapsize-1]; make_heap(0); heapsize--; } return 1; } void make_heap (int root) { int left = (root+1)*2-1; int right = left+1; if (left >= heapsize) { return; } if (heap[root] <= heap[left]) { if (right >= heapsize || heap[root] <= heap[right]) { //root is min //do nothing } else { //right is min //printf("swap %d and %d\n", root, right); SWAP(heap[root], heap[right]); make_heap(right); } } else { if (right >= heapsize || heap[left] <= heap[right]) { //left is min //printf("swap %d and %d\n", root, left); SWAP(heap[root], heap[left]); make_heap(left); } else { //right is min //printf("swap %d and %d\n", root, right); SWAP(heap[root], heap[right]); make_heap(right); } } }
考察とかは特になしです。他のアルゴリズムとのベンチ比較とかも時間があればやってみたい。
第9回PRML読書会
土曜日はサイボウズ・ラボで行われた第9回PRML読書会に参加しました。
自分は発表者トップバッターでSVMの基本的なところを説明しました。
参加者の方からもいろいろ指摘をいただきました。
- なぜマージンを最大化するとよいのか？の説明で『まず２値に分類された学習データをガウスカーネルでのParzen推定を適用して入力の分布を推定する。誤分類が最小になる分類平面は、ガウスカーネルの分散を→０の極限において、マージンを最大化する分類平面に一致する』とあるが、なぜ分散を０に近づけるのかがわからない。
- そういうものとして理解するしかない？理論的な説明はまだ分からずです。。
- Randomized Algorythmを適用してSVMの計算を高速化する手法がある。
- ちょっとググってみたところこの辺ですかね。いろいろと制約はるみたいですがO(log n)で二次計画問題の近似解が求まる！
- biasをゼロと仮定して、二次計画問題を高速で解く手法が存在する。
- liblinearなどのライブラリでこの手法が利用されている。
多クラスSVMの発表のときにも参加者の間で議論がありましたが、SVMはやはり2クラス分類器で、多クラスに応用する事例というのは実際にも少ないそうです。（PRML7.2で紹介されているRVMはSVMとは名前は似ているが完全な別モノ！）
個人的には2クラス分類、と言われてもあまり実際的な使い道が思いつかないのですが、世間的にSVMがもてはやされてるのはそれだけ2クラス分類をしたい問題が多いということなのか。今後SVMの事例を見る際は、そういう観点からも考えていきたいと思いました。
あと、8月末の第6回読書会でもニューラルネットワークの説明をしたのでいまさらですが資料をあげておきます。
このときにも冒頭にSVMとの比較がありましたが、今回改めて両方を比較してみると、
- SVMはSparseな解になる
- 訓練データの特性にもよるがマージン境界上の訓練データ集合が少ない場合は予測計算量がすくない。
- SVMは局所解問題がない
- NNでは局所解問題があるので、初期値選択を慎重にする必要がある。
- SVMにはマージン最大化による汎化能力がある。
- NNは過学習の問題がある。
- SVMは基本的に二クラス分類器
- NNはどちらかといえば回帰に向いている。
という感じでしょうかな。多クラス分類とか回帰とかだったらSVM以外の選択肢を使った方がよさそうかな、と改めて思いました。
Rでガウス過程による分類を実装
PRMLの6.4.5〜6.4.6の範囲にあるガウス過程による分類をRで実装してみました。
ソースコード全文はgithubにアップしています。
http://github.com/thorikawa/prml/blob/master/gaussian_process_classify.R
ここでは例として、(1,0),(2,0),(3,0)で1、(0,1),(0,2),(0,3)で0の値を取る訓練集合を用いています。
# Training data x=list(c(1,0),c(2,0), c(3,0), c(0,1), c(0,2), c(0,3)) t=c(1,1,1,0,0,0) training_data_num <- length(x)
この訓練集合とカーネル関数をもとに予測分布を導出しています。
ガウス過程においては、訓練集合から予測分布を決める（ほぼ）唯一の要素はカーネル関数になるということなのは分かるのですが、カーネル関数によってどのように予測分布が変わってくるのかという点は良く理解できていません。
この辺をまずは感覚でつかむために、いくつか代表的なカーネルで予測分布を求めてみました。
ガウスカーネル
で与えられるカーネル。
の値を0.2,0,4と変えてみると、id:n_shuyoさんの考察通り、
は訓練データを中心とするガウス分布の尖り具合を表現していることが分かります。
# Gaussian kernel sigma <- 0.2 kernel <- function (x1, x2) { exp(-sum((x1-x2)*(x1-x2))/(2*sigma^2)); } gp(kernel, "gaussian kernel", "sigma=0.2");
# Gaussian kernel sigma <- 0.4 kernel <- function (x1, x2) { exp(-sum((x1-x2)*(x1-x2))/(2*sigma^2)); } gp(kernel, "gaussian kernel", "sigma=0.4");
線形カーネル
で与えられるカーネル。単純な線形分離になります。
# Linear kernel kernel <- function (x1, x2) { sum(x1*x2); } gp(kernel, "linear kernel");
指数カーネル
で与えられるカーネル。
の値を1.0,0.3と変えてみると、ガウス分布と同様、
は指数分布の尖り具合を表していると言えそうです。
# Exponential kernel theta <- 1.0 kernel <- function (x1, x2) { exp(-theta*sqrt(sum((x1-x2)*(x1-x2)))); } gp(kernel, "exponential kernel", "theta=1.0");
# Exponential kernel theta <- 0.3 kernel <- function (x1, x2) { exp(-theta*sqrt(sum((x1-x2)*(x1-x2)))); } gp(kernel, "exponential kernel", "theta=0.3");
多項式カーネル
で与えられるカーネル。
M=2,M=3と次数を変えてみると、、、これは、、、どういう特徴を持つのかいまいち不明です。
# Polynomial kernel kernel <- function (x1, x2) { (sum(x1*x2)+1)^2; } gp(kernel, "polynomial kernel", "degree=2");
# Polynomial kernel kernel <- function (x1, x2) { (sum(x1*x2)+1)^3; } gp(kernel, "polynomial kernel", "degree=3");
まとめ
というわけで理解は浅いのですが、様々なカーネル関数でガウス過程による分類を実装してみました。
カーネルの役割は、訓練集合における任意の２つのデータに対する、出力の類似度をどう評価するかを決定することだと理解しています。
ガウスカーネルや指数カーネルは、RBFであり、中心からの距離に対して反比例する形で値が減少していくため、入力値が近ければ近いほどその出力の値も近い値を示し、予測分布は訓練集合を中心とした分かりやすい分布になります。
線形や多項式カーネルはRBFではないので直感的な理解が現時点では出来ていません。カーネル多変量解析を熟読すればこのあたりも理解できるようになるかなあ。
ニューラルネットワークで画像認識
ニューラルネットワークの簡単な関数近似プログラムを先日書いたので、今は画像認識プログラムを書いてますが、ものすごく簡単なバージョンが出来上がったので晒しておきます。
C++で画像解析部分を作って、Rで訓練データの学習、テストデータの判定をしています。
MNISTの手書き文字から学習データ準備
まずは、インプットとなる画像のデータですが、定番のMNISTの手書き数値データを使います。
元々のIDXフォーマットというフォーマットは再利用性が無さそうなので、http://www.cs.toronto.edu/~roweis/data.htmlから既にJPEG化されたものを引っ張ってきてそれを解析します。
OpenCVを使って2値化＋CSV出力する簡単なプログラム。
#include "cv.h" #include "highgui.h" #include <stdio.h> #include <stdlib.h> const int kDIGIT_W = 28; const int kDIGIT_H = 28; int main (int argc, char* argv[]) { IplImage *src_img, *dest_img; char* imgfile = argv[1]; int answer = atoi(argv[2]); src_img = cvLoadImage(imgfile, CV_LOAD_IMAGE_GRAYSCALE); dest_img = cvCreateImage(cvGetSize(src_img), IPL_DEPTH_8U, 1); // make src image gray scale cvThreshold(src_img, dest_img, 90, 255, CV_THRESH_BINARY); int width = dest_img->width; int height = dest_img->height; int xCount = width / kDIGIT_W; int yCount = height / kDIGIT_H; for (int x=0; x<xCount; x++) { for (int y=0; y<yCount; y++) { int offsetX = x*kDIGIT_W; int offsetY = y*kDIGIT_H; printf("%d", answer); for (int i=0; i<kDIGIT_W; i++) { for (int j=0; j<kDIGIT_H; j++) { int val = dest_img->imageData[(j+offsetY)*dest_img->widthStep+(i+offsetX)]; printf(",%d", val*-1); } } printf("\n"); } } cvReleaseImage(&src_img); cvReleaseImage(&dest_img); return 0; }
これで、以下のように実行すると画像データが格納されたCSVが出来上がります。
1行が1文字を表わしており、1列目が正解の値(数値)、2列目以降が二値化された画像データです。
$ g++ -g -I/usr/include/opencv make_test_data.cpp -o make_test_data -lcxcore -lcv -lhighgui -lcvaux -lml #訓練データの準備 $ ./make_test_data mnist_train0.jpg 0 > mnist_train_all.txt $ ./make_test_data mnist_train1.jpg 1 >> mnist_train_all.txt $ ./make_test_data mnist_train2.jpg 2 >> mnist_train_all.txt .... $ ./make_test_data mnist_train9.jpg 9 >> mnist_train_all.txt #テストデータの準備 $ ./make_test_data mnist_test0.jpg 0 > mnist_test_all.txt $ ./make_test_data mnist_test0.jpg 1 >> mnist_test_all.txt .... $ ./make_test_data mnist_test0.jpg 9 >> mnist_test_all.txt
学習プログラム
次にこのCSVデータをインプットするニューラルネットワークの学習プログラムを組みます。
ニューラルネットの設定は以下です。
- 基本的な構造:入力が784次元、出力が10次元、隠れユニット層を1つ含む2層ネットワーク(PRMLで言うところの)
- 隠れユニット数:30
- 学習パラメータ:0.05
- 学習方法:オンライン学習。ただし、1つのエポックでは60338件の訓練データ集合の中から1000件のデータをランダム抽出し、それを1000エポック繰り返す。
- 出力ユニットの活性化関数:ソフトマックス関数
- 誤差関数:交差エントロピー誤差関数（ただし、プログラム中は活性に対する微分であるyk-tkのみ使われているため、交差エントロピー関数そのものは出てこない。）
#number of hidden unit
s <- 30
#learning rate parameter
l <- 0.05
#count of loop
MAX_EPOCH <- 1000
#read traning data
digit_data <- read.csv("train/mnist_train_all.txt",header=F)
##initialize weight parameter
#weight parameter between input and hidden units
w1 <- matrix(runif(s*length(digit_data), -1, 1), s, length(digit_data))
#weight parameter between hidden units and output
w2 <- matrix(runif(10*(s+1), -1, 1), 10, s+1)
#definition of function
logsumexp <- function (x) {
m <- max(x)
m + log(sum(exp(x-m)))
}
softmax <- function (a) {
sapply(a, function (x) { exp(x-logsumexp(a)) })
}
neuro_func <- function (input) {
a1 <- w1 %*% c(1, input)
z1 <- tanh(a1)
a2 <- w2 %*% c(1, z1)
z2 <- softmax(a2)
return(z2)
}
##train the neural network
for (k in 1:MAX_EPOCH) {
#sample 1000 training data
sample_index = sample(1:nrow(digit_data), 1000)
for (i in sample_index) {
tmp <- digit_data[i,]
#target data
t <- rep(0, 10)
t[as.integer(tmp[1]+1)] <- 1
input <- t(tmp[2:length(tmp)])
#forward propagation
a1 <- w1 %*% c(1, input)
z1 <- tanh(a1)
a2 <- w2 %*% c(1, z1)
z2 <- softmax(a2)
#back propagation
d2 <- z2 - t
w2 <<- w2 - l * d2 %*% t(c(1, z1))
d1 <- d2 %*% w2[,2:(s+1)] * (1-tanh(t(a1))^2)
w1 <<- w1 - l * t(d1) %*% t(c(1, input))
}
}
save(w1, w2, logsumexp, softmax, neuro_func, file="nnet_image.rdata")
テストプログラム
準備したテストデータを読み込んで、ネットワークの出力のうち、最大の出力を持つユニットをネットワークが判断した数値とみなしています。
#read trained parameter
load("nnet_image.rdata")
#read traning data
digit_data <- read.csv("test/mnist_test_all.txt",header=F)
##main function to train the neural network
correct <- 0
wrong <- 0
for (i in 1:nrow(digit_data)) {
tmp <- digit_data[i,]
# target number of class
answer <- tmp[1]
t <- rep(0, 10)
t[as.integer(answer+1)] <- 1
input <- t(tmp[2:length(tmp)])
# forward propagation
z <- neuro_func(input)
zt <- which.max(z) - 1
if (zt == answer) {
correct <- correct + 1
} else {
wrong <- wrong + 1
}
}
cat(paste("correct case=",correct,", wrong case=", wrong, "\n"))
結果とまとめ
テストデータ集合10153件中9360件正解、ということで正解率は約92.5%でした。
（correct case= 9360 , wrong case= 793 ）
自分の手書き文字でも試してみましたが、まあ意地悪いことをしなければ大体正解しているみたいです。
ただ、他の方の結果と比較すると、id:ultraistさんの結果やhttp://yann.lecun.com/exdb/mnist/の例では、2層ネットワークで97%近い正解率まで達しているので、まだまだ。
id:ultraistさんは隠れユニット数が512で、LeCun氏は300とか1000とかなので、僕の隠れユニット数30はちょっと少なすぎたか、という印象。ただ、入力ベクトルが784次元なので、それに対し隠れユニットが1000はいくらなんでも多い、というか一般化できていないんじゃないかそれ？という感じなのだが、どうなんだろう。
とりあえず始めたばかりなので、TODOを列挙。
- 隠れユニット数をいじってみて比較
- 学習パラメータをいじってみて比較
- deskew処理
- オブジェクト抽出処理
- 画素以外の特徴ベクトルを入力にして学習
- 2層ネットワーク以外のやり方と比較
と見えているだけでもまだまだやることてんこ盛りです・・・。地道にやっていきます。
2009/9/13追記
学習時間について
この例ではのべ100万件のデータを学習していますが、全部で丸2〜3日ほどかかりました。
Rでニューラルネットワーク（2変数の関数近似）
で、1変数の関数近似がうまくいったので、調子にのって2変数の関数近似にもチャレンジしてみました。
2変数のsinc関数を、ニューラルネットワークの誤差伝播法を使って近似する（しようとする）ものです。
library(animation) #number of training data N <- 30 #number of hidden unit s <- 5 #learning rate parameter l <- 0.01 #standard deviation of training data sd <- 0.05 #count of loop LOOP <- 10000 #frame interval of animation movie INTERVAL <- 0.1 #total time of animation movie TOTAL <- 20 #initialize weight parameter w1 <- matrix(seq(0.1,length=s*3, by=0.1),s,3) w2 <- matrix(seq(0.1,length=s+1, by=0.1),1,s+1) #generate traning data xt <- seq(-10,10,length=N) yt <- xt f <- function(x, y) { r <- sqrt(x^2+y^2) temp = function (r) { if ( r == 0 ) { return(10) } else { return(10*sin(r)/r) } } sapply(r, temp) } zt <- outer(xt, yt, f) persp(xt,yt,zt,theta = 30, phi = 30, expand = 0.5, col = rainbow(50), border=NA) neuro_func_ <- function (x, y) { a <- w1 %*% c(1, x, y) z <- tanh(a) y <- w2 %*% c(1, z) return(y) } #vectorize neuro_func_ neuro_func <- function (x, y) { mapply(neuro_func_, x, y) } persp(xt, yt, outer(xt, yt, neuro_func), theta = 30, phi = 30, expand = 0.5, col = rainbow(50), border=NA, axes=TRUE, ticktype="detailed") #main function to train the neural network neuro_training <- function () { trigger_count <- as.integer(LOOP / (TOTAL / INTERVAL)) for (k in 1:LOOP) { if (k %% trigger_count == 1) { persp(xt, yt, outer(xt, yt, neuro_func), theta = 30, phi = 30, expand = 0.5, col = rainbow(50), border=NA, axes=TRUE, ticktype="detailed") } for (i in 1:N) { for (j in 1:N) { # forward propagation a <- w1 %*% c(1, xt[i], yt[j]) b <- tanh (a) z <- w2 %*% c(1, b) # back propagation d2 <- z - zt[i,j] w2 <<- w2 - l * d2 * c(1, b) d1 <- d2 %*% (1-tanh(t(a))^2) * t(w2[2:(s+1)]) w1 <<- w1 - l * t(d1) %*% c(1, xt[i], yt[j]) } } } } saveMovie(neuro_training(), interval=INTERVAL, moviename="neural_network_training_movie", movietype="mpg", outdir=getwd(), width=640, height=480) cat("learning end print w1 and w2\n") print(w1) print(w2) persp(xt, yt, outer(xt, yt, neuro_func), theta = 30, phi = 30, expand = 0.5, col = rainbow(50), border=NA, axes=TRUE, ticktype="detailed")
近似したい元関数が以下。
で、結果が以下のムービー。
サンプルデータを10000回ループで読み込ませているんですが、途中から微動だにしなくなってしまいました・・・。
何かパラメータの設定が悪いのか、もっとループを多くしなければいけないのか、まだ原因は分かっていません。。。どなたか知識のある方がこれを読んで頂けたら指摘してもらえると嬉しいです。
Rでニューラルネットワーク（1変数の関数近似）
機械学習・パターン認識方面の勉強初めてから４ヶ月ほど立ちました。最近はnaoya_tさん主催のPRML読書会に参加させて頂いています。
来週末8/29の第6回読書会ではニューラルネットワークの章の発表を担当することになったので、Rを使ってサンプルプログラムを組んでみました。参考にしたのはPRML5.1〜5.3の範囲で、sin関数からサンプリングした点データをニューラルネットワークを使って誤差逆伝播法で学習し、元のsin関数を近似します。
学習前の初期状態が以下の図。赤字が元の関数（線）およびサンプルデータ（点）で青字がニューラルネットワークの出力です。
で、学習後の状態が以下です。
いい感じに再現できています。
以下ソースコード。
library(animation) #number of training data N <- 50 #number of hidden unit s <- 5 #learning rate parameter l <- 0.05 #standard deviation of training data sd <- 0.05 #count of loop LOOP <- 7000 #frame interval of animation movie INTERVAL <- 0.1 #total time of animation movie TOTAL <- 25 #initialize weight parameter w1 <- matrix(seq(0.1,length=s*2, by=0.1),s,2) w2 <- matrix(seq(0.1,length=s+1, by=0.1),1,s+1) #generate traning data xt <- seq(0,1,length=N) yt <- sin(2*pi*xt)+rnorm(N,sd=sd) neuro_func_ <- function (x) { a <- w1 %*% c(1, x) z <- tanh(a) y <- w2 %*% c(1, z) return(y) } #vectorize neuro_func_ neuro_func <- function (x) { sapply(x, neuro_func_) } png(filename = "graphic1.png", width = 480, height = 480, pointsize = 18, bg = "white") plot(function(x) {return(sin(2*pi*x))}, xlab="input", ylab="output", xlim=c(0, 1), ylim=c(-2.5,2.5), sub=paste("S=", s, ",N=", N, ",l=", l), col=2) plot(neuro_func, add=T, col=4, xlim=c(0,1)) points(xt, yt, col=6) #main function to train the neural network neuro_training <- function () { trigger_count <- as.integer(LOOP / (TOTAL / INTERVAL)) for (k in 1:LOOP) { if (k %% trigger_count == 1) { plot(function(x) {return(sin(2*pi*x))}, xlab="input", ylab="output", xlim=c(0, 1), ylim=c(-2.5,2.5), main=paste("LOOP=", k), sub=paste("S=", s, ",N=", N, ",l=", l), col=2) plot(neuro_func, add=T, col=4, xlim=c(0,1)) points(xt, yt, col=6) } for (i in 1:N) { # forward propagation a <- w1 %*% c(1, xt[i]) z <- tanh (a) y <- w2 %*% c(1, z) # back propagation d2 <- y - yt[i] w2 <<- w2 - l * d2 * c(1, z) d1 <- d2 %*% (1-tanh(t(a))^2) * t(w2[2:(s+1)]) w1 <<- w1 - l * t(d1) %*% c(1, xt[i]) } } } saveMovie(neuro_training(), interval=INTERVAL, moviename="neural_network_training_movie", movietype="mpg", outdir=getwd(), width=640, height=480) cat("learning end.print w1 and w2\n") print(w1) print(w2) png(filename = "graphic2.png", width = 480, height = 480, pointsize = 18, bg = "white") plot(function(x) {return(sin(2*pi*x))}, xlab="input", ylab="output", xlim=c(0, 1), ylim=c(-2.5,2.5), sub=paste("S=", s, ",N=", N, ",l=", l), col=2) plot(neuro_func, add=T, col=4, xlim=c(0,1)) points(xt, yt, col=6)
収束する様子をムービーにしてみました。
中央上部に表示されている数値はサンプル点学習のループ回数です。
考察・雑感
学習パラメータは0.05固定としていますが、いろいろといじってみて、この値を大きくしすぎると収束せずに発散したり、振動したりするので、これくらいの値が妥当と判断しました。学習するデータ集合の数を多くした場合は、一つ一つの学習におけるの更新量は小さくするべきなので、これらの学習パラメータも小さくすることが望ましいはずです。
また、いずれにしても精度のよい近似が得られるまではデータ集合を繰り返し学習する必要があります。上記の動画を見ても、視覚的に精度のよい近似に近づくまでは1000回程度の学習ループを必要としています。学習パラメータを、固定ではなく、学習データに基づくパラメータに依存させたりすることでこの近似のスピードをもう少し早めることが出来るのかも知れません。PRMLの後の方で出てくるのかな？
参考文献
R全般に関して、id:syou6162氏のはてダを多いに参考にさせて頂きました。出し惜しみなく情報を晒して頂いてありがとうございます。
あと、Rでアニメーションを作るのは以下のエントリを参考にしました。
やろうと思ったことがほぼ何でもできちゃうR、ステキです！！。
DMC(Dynamic Markov Coding)のによるデータ圧縮プログラムを書いてみた
最近Managing Gigabytes勉強会に参加しているのでせっかくなので、この本に載っているアルゴリズムを使ってプログラムを組んでみました。
今回実装したのは、「2.5 SYMBOLWISE MODELS」の後半で説明されている「Dynamic Markov Coding(DMC)」です。書籍の他に、元論文「G. Cormack and R. Horspool, "Data compression using dynamic Markov modelling,"」を参考にしました。
実装はC++で行い、ソースはgithubに置きました。（CentOS 5.2+gcc 4.1.2で動作確認済）
http://github.com/thorikawa/MG/tree/master/dmc
以下、アルゴリズムの概要と実装上の工夫などをまとめてみます。
意見・指摘などは絶賛大歓迎です。
DMC(Dynamic Markov Coding)とは
DMCは適応型の算術符号化の一種で、データの流れ全体をマルコフ連鎖で表します。符号化の過程で、マルコフ連鎖の構造と確率を変化させていき、各状態における次状態への遷移確率を元に符号化します。
上の図は元論文から抜粋したものですが、状態Aから出ている矢印とその添字は、状態Aで読み込む次のビットが0である確率が50%、1である確率が50%であることを示しています。状態Aで0のビットを読み込むと、50%という確率に基づいて符号化を行った後に、状態Aから遷移する確率を更新して、状態Bに遷移します。
圧縮率を高めるためには、各状態における遷移確率が偏るようにマルコフ連鎖の構造を変化させていく必要があります。そのため、DMCではクローニングという仕組みを導入しています。クローニングを示した図が以下です。（元論文より抜粋）
(a)がクローニングを行う前で、AとBという２つの状態から状態Cに遷移していますが、ある一定の条件を満たしたときに、状態CはCとC’に分離します。ある条件というのは以下のソースを参照してほしいですが、元論文では２つのパラメタ閾値１と閾値２を設定して、AからCへのの遷移カウントが閾値１を超え、かつAからCへの遷移カウントとA以外からCへの遷移カウントとの差が閾値２を超えた場合にクローニングが行われます。
// クローニングを行う。 void DoCloning (uchar bit) { uint trans_count = current_state_->GetTransCount(bit); State *next_state = current_state_->GetNextState(bit); uint next_count = next_state->GetTotalCount(); if (trans_count > cloning_threshold1_ && (next_count - trans_count) > cloning_threshold2_) { State *new_state = new State(last_state_++); current_state_->SetNextState(bit, new_state); double ratio = (double) (trans_count + 1) / (next_count + 2); // ratioの正規化(CloningはCloning後の次状態のtrans_count <= 現状態のtrans_countを保証していないので、ratioが1を超える場合がある。 if (ratio > 1) { ratio = 1; } for (int i=0; i<=1; i++) { uint next_count_bit = next_state->GetTransCount(i); uint new_next_count_bit = (uint)(ratio * next_count_bit); new_state->SetNextState(i, next_state->GetNextState(i)); new_state->SetTransCount(i, new_next_count_bit); next_state->SetTransCount(i, next_count_bit - new_next_count_bit); } cloning_count_++; } }
クローニングが行われないと、A→C→Dと、B→C→Eの２ケースのビット遷移が極端に多い場合でも、Cに到達した時点でAからきたかBからきたかの情報が失われてしまい、C→D・C→Eの遷移確率が均一化されてしまうので、効率のよい圧縮はできません。クローニングを行うことにより、出現頻度の多いパターンがたどる経路を分離させることができ、各状態における遷移確率を偏らせ、効率のよい圧縮が出来る、と理解しています。
Guazzoによる算術符号化アルゴリズム
次に、各状態における遷移確率から符号化を行う部分です。
Managing GigabytesにはDMCに特化したこの部分のロジック詳細は記載されていませんが、今回の実装では元論文でも紹介されているGuazzoによる算術符号化アルゴリズムを利用します。
Guazzoのアルゴリズムでは、整数値を取る数直線上をデータの出現確率に応じて分割します。今回はデータの入出力をビット単位で扱うこととしまうので、0と1に対応する二区間に数直線を分割していけばよいです。
数直線の幅はNを任意の整数値として[0x0,(1<<(N+1))-1]で初期化し、データを1ビット読み込むたびに、新たな上限と下限を計算していきます。そして計算し直した上限と下限を比較して両者の上位ビットが一致する場合は、一致するビットを確定ビットとしてアウトプットします。このビットシフトによって、上限と下限の最上位ビットは常に異なっている状態になり、数直線の幅を確保することができます。
以下が符号化部分のコードです。
void Encode (uchar bit) { ulong mp = CalculateMP(); Encode(bit, mp); } void Encode (uchar bit, ulong mp) { if (bit == 1) { lower_bound_ = mp; } else { upper_bound_ = mp - 1; } // 上限と下限を比較して確定しているbitがあればencoded bitとしてバッファに格納する。 // 確定しているbit分だけ、上限・下限をシフトする。 while ((lower_bound_ & kMSBIT) == (upper_bound_ & kMSBIT)) { unsigned char outbit = lower_bound_ >> (kN-1); AddBuffer(outbit); lower_bound_ = (lower_bound_ << 1) & kMSMASK; upper_bound_ = ((upper_bound_ << 1) & kMSMASK) | 1; } // クローニングを行う。 DoCloning(bit); // マルコフ連鎖を更新する。 current_state_->SetTransCount(bit, current_state_->GetTransCount(bit)+1); current_state_ = current_state_->GetNextState(bit); }
// 分離点(mp)を計算する。 ulong CalculateMP () { double p0, p1; ulong mp; uint totalCount = current_state_->GetTotalCount(); p0 = (double) (current_state_->GetTransCount(0) + 1) / (totalCount + 2); p1 = (double) (current_state_->GetTransCount(1) + 1) / (totalCount + 2); mp = (ulong) ((p1 * lower_bound_ + p0 * upper_bound_) / (p0 + p1)); if (mp <= lower_bound_) { mp = lower_bound_ + 1; } // mpの最右bitを1に変更 mp = mp | 1; // 最右bitを1に変更した結果、upper_bound_を超える可能性がある。その場合、mp=upper_bound_とする。 if (mp > upper_bound_) { mp = upper_bound_; } return mp; }
Guazzoアルゴリズムの実践における課題
元論文ではGuazzoのアルゴリズムを実装するにあたり２つの問題があると述べています。
- １つ目の問題は、整数値の数直線を上限と下限を分数で分割していくため、過程が進むにつれてより大きな精度が必要になる、という問題。
- ２つ目の問題は、符号化するデータは有限であるという制約のもと、過不足のないデータを復元しなくてはいけないという問題。
まず一つ目の問題について、Guazzoは、正確に境界値を求める必要性を弱めることで解決しているそうです。上記のソースコードではCalculateMPの中で、確率からmp（数直線上の分岐点）を求める際にdouble型から整数値に丸めたり、mpの末尾bitを1に固定化したり（理由は後述）しています。このようにmpが正確な値ではなくても、圧縮時・復元時ともに同一の基準（論文ではfidelity criterion（忠実度基準）と呼んでいます）に従って計算が行われていれば、圧縮・復元に際しては大した問題にはなりません。
二つ目の問題については、かなり手ごわい問題で、たとえば1ビットが圧縮後に2ビットになった場合、圧縮後の2ビットがインプットされても元の2ビットを復元できる保障がない点に注意を払う必要があります。例えば、区間が[0x0000...0000, 0x1111...1111]として、入力ビット"0"を圧縮するとします。はじめの0で分岐点が"0x0010....0001"となり、区間が[0x0000...0000, 0x0010....0000]に縮小された場合、符号化された出力ビットは"00"になります。同様の状態でデコーダは分岐点の"0x0010....0001"以上かどうかで復元後のビットが0なのか1なのかを判断しますが、入力ビット"00"を復元しようとした場合、今後のビットの出現次第で、"0x0000....0000"や"0x0011....1111"のように分岐点以上にも以下にもなる可能性があります。そのため、この時点では復元することはできないのです。
すなわちaビットをbビットに符号化したとして、bビットあれば元のaビットを過不足なく復元できるとは限らないわけです。
この問題を解消するために、符号化時に二つの工夫をしています。
それは「mpの最右ビットを必ず1にしていること」と「符号化の最後にmpの最右ビットを除くN-1ビットを出力していること」です。これだけで過不足なく元データを復元することが出来るようになります。
理由を説明します。まずmpの最右ビットが1であるという条件から、数直線上の区間は必ず（ビットシフトの演算をのぞき）縮小されていきます。そのため、最後の入力ビットを符号化したときのmpのN-1ビット目までと最後に出力されるmpのN-1ビット目までは必ず異なるものになります。最後に出力されるmpのN-1ビットをインプットとすれば、最後の入力ビットを符号化したときのmp以上か未満かを必ず決定できる、つまり最後の入力ビットを決定できるのです。
また、最後のmpのN-1ビット目で止めている、というのは余分な1ビットが符号化されてしまうのを防止する意味もあります。最後に出力されるmpのNビット目を出力してしまうと、mp以上が確定してしまい、余分な"1"が符号化されてしまうことになります。しかし、N-1ビット目までの情報だと、デコーダは最後に出力されるmp以上か未満かを判別することはできません。つまり余分な1ビットは符号化されないのです。
Extra 7bits
さらに考慮すべき点があります。それはこのロジックでは符号化がビット単位になされていますが、実際にはデータはバイト単位だということです。符号化した後のビット数が8に達するたびにバイト単位で書き出しを行っていますが、最後に8に達しなかったあまりのビットは出力されることはありません。この状態のままでは完全な復元が出来なくなってしまいます。
そこで元論文で提案しているのは、符号化の最後に余分な7ビット(Extra 7 bits)を符号化することです。
// encodedされた最終bitが、バイトとして出力されるようにダミーの7bitをencodeする。 // また、encode対象の最終bitが正常にencodeされるように、mpの末尾1bitを除いた全てのbit列をencoded bitとして出力する。 void EncodeFinish () { ulong mp; for (int i=0; i<7; i++) { mp = CalculateMP(); // ダミーのbitとして、encoded bitが最低1bit確定されるbitを選択し、encodeする。 if ((lower_bound_ & kMSBIT) == (mp & kMSBIT)) { // lower_bound_の最上位bitとmpの最上位bitが一致する場合は、ダミーのbitとして0をencodeする。 Encode(0, mp); } else { // upper_bound_の最上位bitとmpの最上位bitが一致する場合は、ダミーのbitとして1をencodeする。 Encode(1, mp); } } mp = CalculateMP(); // mpの最後の1bitを除き、出力する。 while (mp != kMSBIT) { uchar outbit = (mp >> (kN-1)); AddBuffer(outbit); mp = (mp << 1) & kMSMASK; } } };
ここで余分な7ビットを符号化しています。7ビットはなんでもよいかといえばそうではなく、それまでに符号化したビット全体がバイト列として出力されるためには、本体の符号化後ビット数が8N+1ビットだった場合を考慮して、最低でも7ビット分余分に出力を行う必要があります。そのため、0と1のうち、符号化した後に最低でも1ビット符号化される方のビットを選択して出力していきます。
これでバイト単位に丸めることによって、ビットが切り捨てられてしまったとしても、データ本体を符号化したビットは生き残るので、完全な形で復元することができます。
ベンチマーク
最後に、この圧縮手法で、圧縮率のベンチマークをとってみます。
テスト対象としては、定番所でThe Canterbury Corpusを利用しています。
まずはファイル種類別の圧縮率の相関図です。
それぞれのデータ特性についてあまり詳しく分析できていないのですが、テキストファイルではおおよそ35〜50%くらいの圧縮性能が実現できています。id:naoya氏が最近日記に書かれているPerl で Range Coder (再挑戦)のRange Coderでも、alice29.txtの圧縮率は50%をクリアできていなかったので、マルコフ連鎖を利用するだけで、大幅な性能改善を図ることが出来た、と理解しています。
次にCloningの閾値（上記ソースのcloning_threshold1=cloning_threshold2の値）を変化させていったときの圧縮率の変化です。他の圧縮アルゴリズムの比較として右端にZIPでの圧縮結果をプロットしています。
元論文では閾値が2のときがおおむね最も圧縮率が高い、と記述されていましたが、実際には上記グラフの通り、閾値が8〜32くらいのときに圧縮率がピークになるようです。直観的にも、閾値が低すぎると、各状態における状態遷移確率が偏る前に、すぐにクローニングしてしまうので、、圧縮率も低下するものと理解できます。
また、ZIPと比較すると、ptt5（FAXファイル）やEXCELファイル以外ではまだまだ圧縮率で負けていることがわかります。
まとめ
DMCを使って圧縮を行うと、エントロピーを最小にするように動的にマルコフ連鎖を最適化しながら圧縮を行うため、全体の出現頻度表をベースに圧縮するよりも効率のよい圧縮化を行うことができます。他のアルゴリズムとの比較など、考察が足りない部分は多々ありますが、それはおいおいということで。あらかじめ次数を決めてかかるPPMなどよりも、圧縮率はよいものと思われます。
参考文献
[1] Ian H. Witten, Alistair Moffat, Timothy C. Bell, Managing Gigabytes: Compressing and Indexing Documents and Images, Morgan Kaufmann Publisher, 1999
[2] G. Cormack and R. Horspool, "Data compression using dynamic Markov modelling," Computer J., 30 No. 6, 541-550, 1987.
[3] 高速な算術圧縮を実現する「Range Coder」(記事), 2006, CodeZine, 翔泳社
