MATLABでレプリカ交換法を試してみた

 ベイジアンにとってMCMCは非常に強力な手法です。例えば観測データを元にパラメータの事後分布を推定したいときに使えます。しかし、事後分布の形状によっては十分なサンプルサイズを得るまでに非常に大きな計算量を必要とする欠点があります。具体的には、分布が多峰性を持つ場合と、峰が細長い場合です。PRML(パターン認識機械学習)では後者に対する解決策としてHamiltonian MCという手法が紹介されていますが、モデルが微分可能でなければならないという限界がありますし、多峰性に対する解決策ではありません。ここではPRMLで触れられていない方法のひとつとしてレプリカ交換法を試してみます。

準備

 まず、多峰性を持つ分布としてOld Faithful間欠泉のデータをEMアルゴリズムで混合ガウス分布にしたものを用意しました。1行目で読み込んでいるfaithful.txtは
http://www.math.pku.edu.cn/teachers/jjia/Computational%20Statistics/faithful.txt
ここから持ってきました。スペース区切りのデータを2列の行列として読みます。

M = dlmread('faithful.txt');
G = fitgmdist(M,2);
[X,Y] = meshgrid(linspace(1.5,5.5),linspace(40,100));
Z = pdf(G,[X(:),Y(:)]); % 20151214追記
% Z = pdf(G,[reshape(X,[numel(X),1]),reshape(Y,[numel(Y),1])]);
fig = figure(1);
clf(fig);
fig.Position = [100,100,1200,900];
hold on;
scatter(M(:,1),M(:,2),'o');
contour(X,Y,reshape(Z,size(X)));

PRML下巻153ページと同じfigureが得られました。
f:id:eukaryo:20151213185025p:plain

Metropolis法

 準備した混合ガウス分布からMCMCでサンプリングすることで、多峰性に対するMCMCの振る舞いを観察します。まず普通のMetripolis法を試してみました。
コードは以下の通りです。

function sample = Metropolis_MCMC(GMM_model, data, time)
sample = zeros(time,numel(data(1,:)));
sample(1,:) = data(randi(size(data,1)),:);
for t = 2:time
    diff = randn(1,2).*[2.0,20.0];
    candidate = sample(t-1,:) + diff;
    L_now = log(pdf(GMM_model,sample(t-1,:)));
    L_next = log(pdf(GMM_model,candidate));
    if rand < exp(L_next - L_now)
        sample(t,:) = candidate;
    else
        sample(t,:) = sample(t-1,:);
    end
end
end

 試すためのコードは以下の通りです。

sample = Metropolis_MCMC(G,M,100000);
sample(1:10000,:) = []; % burn-in
scatter(sample(:,1),sample(:,2),'rx');
fig = figure(2);
clf(fig);
autocorr(sample(:,1),50);

 結果として、正しくサンプリングできているように見えますが、別の峰に移るのが困難であるため自己相関が高くなっています。
f:id:eukaryo:20151213185147p:plain
f:id:eukaryo:20151213185157p:plain

レプリカ交換法

 別の峰に移るのが困難であるのは、峰の尤度に対して鞍の尤度が低すぎるからです。そこで、混合ガウス分布対数尤度に一定値を掛ける、言い換えると精度行列に一定値を掛けることで鞍を越えやすいようにした分布を考えます。例えば、最初のコードの4行目を

T_inv = 0.1;
Z = exp(log(pdf(G,[X(:),Y(:)]))*T_inv); % 20151214追記
% Z = exp(log(pdf(G,[reshape(X,[numel(X),1]),reshape(Y,[numel(Y),1])]))*T_inv);

とすると、等高線は
f:id:eukaryo:20151213181816p:plain
となります。掛けている値を逆温度と呼びます。この分布から通常のMetropolis法でサンプル列を得ると、当然ながら自己相関は低くなります。(散布図は省きます)
f:id:eukaryo:20151213182441p:plain
 ここで、異なる逆温度のサンプル列を同時並行的に得つつ、双方のサンプル列の詳細釣り合い条件を保ったままサンプルを交換することができます。具体的には
レプリカ交換法 - Wikipediaなどにありますが、交換する確率を
\min \left( 1, \exp\left\{{(E_i - E_j) \left( \frac{1}{kT_i} - \frac{1}{kT_j}\right)}  \right\}\right)
とします。Eは逆温度を考えないときの対数尤度です。逆温度が低い側の尤度が高いならば必ず交換します。
 逆温度が低いサンプル列で鞍を越え、逆温度が高いサンプル列と交換することによって自己相関を低くできるわけです。

コードと結果は以下の通りです。

function [sample, rate] = Replica_Exchange_MCMC(...
    GMM_model, data, time, replica_num, d_temp)

T_inv = 1.0./(d_temp.^(0:replica_num-1));

sample = zeros(time,numel(data(1,:)),replica_num);
for n = 1:replica_num
    sample(1,:,n) = data(randi(size(data,1)),:);
end

rate = [0,0];
for t = 2:time
    for n = 1:replica_num
        diff = randn(1,2).*[2.0,20.0];
        candidate = sample(t-1,:,n) + diff;
        L_now = log(pdf(GMM_model,sample(t-1,:,n)));
        L_next = log(pdf(GMM_model,candidate));
        if rand < exp((L_next - L_now)*T_inv(n))
            sample(t,:,n) = candidate;
        else
            sample(t,:,n) = sample(t-1,:,n);
        end
    end
    for n = 1+rem(t,2):2:replica_num-1
        L_lowtemp = log(pdf(GMM_model,sample(t,:,n)));
        L_hightemp = log(pdf(GMM_model,sample(t,:,n+1)));
        if rand < exp((L_hightemp - L_lowtemp)*(T_inv(n)-T_inv(n+1)))
            tmp = sample(t,:,n);
            sample(t,:,n) = sample(t,:,n+1);
            sample(t,:,n+1) = tmp;
            rate(1) = rate(1) + 1;
        else
            rate(2) = rate(2) + 1;
        end
    end

end
sample = sample(:,:,1);
rate = rate(1) / sum(rate);
end
sample = Replica_Exchange_MCMC(G,M,100000,6,1.6);
sample(1:10000,:) = [];
scatter(sample(:,1),sample(:,2),'rx');
fig = figure(2);
clf(fig);
autocorr(sample(:,1),50);

f:id:eukaryo:20151213180911p:plain
f:id:eukaryo:20151213180922p:plain

まとめ

 レプリカ交換法を紹介し、分布が多峰性を持つ場合に役に立つ手法であることを示しました。

www.amazon.co.jp
www.amazon.co.jp

ビットプレーン構造を使ってDNAの検索をちょっと高速化する方法

 文字列の全文検索を(Suffix Arrayを使って)実装する際に、「文字列の先頭からn文字の中に文字aがいくつあるか」を数えたいことがあります。文字列が自然言語の場合はウェーブレット行列というデータ構造を使うのが良いようです。ウェーブレット行列なら文字の種類数kに対してlog_{2}k段の操作で済みます。しかし、DNAの塩基配列を文字列とみなして"全文検索"を行いたい場合、4文字しかないのでウェーブレット行列はあまり美味くないと思われます。(実際に試してはいないので推測です)
 ウェーブレット行列と共に用いられるデータ構造に完備辞書というものがあります。今回はDNAの塩基配列のBWTをビットプレーン構造で表現したうえで、完備辞書の考え方を参考にしつつ文字の数を素早く数える方法を紹介します。

方法

 例えば、{3,2,1,0,2,2,1,1}という配列中に2がいくつあるかを数えたいとします。シンプルな方法としてはfor文で走査してif文でカウントしていくものがあります。

	uint32_t a[8] = { 3, 2, 1, 0, 2, 2, 1, 1 };
	uint32_t num = 0;
	for (uint32_t i = 0; i < 8; i++)if (a[i] == 2)num++;

 この方法はbranch misprediction penaltyのために結構遅くなります。文字が4種類しかないので、if (a[i] == 2)が割と頻繁にtrueになるからです。
 そこで、{3,2,1,0,2,2,1,1}を2進数で縦書きして横に読むことを考えます。

元の数列 3 2 1 0 2 2 1 1
2の位 1 1 0 0 1 1 0 0  ←51
1の位 1 0 1 0 0 0 1 1  ←197

 2進数の0b00110011は10進数の51であり、同じく0b11000101は197です。表と逆向きに読んでいる理由は後述します。このように、同じビット位置のビットを集めたデータ構造をビットプレーン構造と言います。この計算を具体的に書くと以下のようになります。

	uint32_t a[8] = { 3, 2, 1, 0, 2, 2, 1, 1 };
	uint32_t b[2] = { 0, 0 };
	for (uint32_t i = 0; i < 8; i++)
	{
		if(a[i] & 1)b[0] += 1 << i;
		if(a[i] & 2)b[1] += 1 << i;
	}
	assert(b[0] == 197 && b[1] == 51);

 元の文字列中に2がいくつ出るかを求めるには、2==0b10であることから

	uint32_t num = __popcnt(~b[0] & b[1]);

とすることで計算できます。__popcntは立っているビットの数を求めて返すx86の組み込み関数です。
 配列aの先頭からx文字(x<-[1..8])に限って数えたい場合、(2^x)-1との論理積を取れば求められます。表と逆向きに読んでいたのはこのためです。

	uint32_t num = __popcnt(~b[0] & b[1] & ((1 << x) - 1));

 文字列の長さnがプロセッサのワード長wを超える場合、返すべき値をw文字間隔で予め計算しテーブルとして保持しておけば、上述の計算を1度行うことで任意のnに対応できます。

ソースコード

 まず、BWTや補助データを作る関数です。8文字ごとでなく64文字ごとになっています。また、終端文字を含めて5文字から成ります。

void make_BWT(
const uint8_t* src, const uint32_t target_len,
uint32_t* sa, uint64_t* bwt,
uint32_t* C, uint32_t* Occ)
{
	//srcが元文字列とする。srcのBWTを行って結果をbwtに返す。
	//文字列は2文字以上であり、最後の文字はゼロ(=$)であり、
	//かつ他の文字はATGC(=1,2,3,4)のいずれかであることを仮定する。

	assert(2 <= target_len);
	assert(src[target_len - 1] == 0);
	for (uint32_t i = 0; i <= target_len - 2; i++)
	{
		assert(1 <= src[i] && src[i] <= 4);
	}

	//srcの接尾辞配列を構築してsaに格納する。
	make_suffix_array(src, target_len, sa);

	//LUTを初期化する。
	for (uint32_t i = 0; i < 5; i++)C[i] = 0;
	const uint32_t occ_len = (target_len / 64) * 4;
	for (uint32_t i = 0; i < occ_len; i++)Occ[i] = 0;
	const uint32_t bwt_len = ((target_len + 63) / 64) * 3;
	for (uint32_t i = 0; i < bwt_len; i++)bwt[i] = 0;

	//BWTを行いつつCを構築する。C[i]にはiより小さい文字の総出現回数を格納する。
	for (uint32_t i = 0; i < target_len; i++)
	{
		const uint32_t n = sa[i] ? src[sa[i] - 1] : 0;
		C[n]++;

		const uint64_t mask = 1ULL << (i % 64);

		const uint32_t index = (i / 64) * 3;
		if (n & 1)bwt[index + 0] |= mask;
		if (n & 2)bwt[index + 1] |= mask;
		if (n & 4)bwt[index + 2] |= mask;
	}
	for (uint32_t i = 1; i <= 4; i++)C[i] += C[i - 1];
	for (uint32_t i = 4; 1 <= i; i--)C[i] = C[i - 1];
	C[0] = 0;

	//Occを構築する。Occ[i*4+x] = Occ(x, i)とし、
	//Occ(x, i)にはbwt[0]~[i*64+63]における文字xの出現回数を格納する。
	for (uint32_t i = 0; i < occ_len; i += 4)
	{
		const uint32_t i_prev = i == 0 ? i : i - 4;
		const uint32_t j = (i / 4) * 3;

		Occ[i + 0] = Occ[i_prev + 0] + __popcnt64(bwt[j + 0] & ~bwt[j + 1]);
		Occ[i + 1] = Occ[i_prev + 1] + __popcnt64(~bwt[j + 0] & bwt[j + 1]);
		Occ[i + 2] = Occ[i_prev + 2] + __popcnt64(bwt[j + 0] & bwt[j + 1]);
		Occ[i + 3] = Occ[i_prev + 3] + __popcnt64(bwt[j + 2]);
	}
}

 次に、ビットプレーン構造からOccを補間して返す関数です。

uint32_t interpolate_Occ(
const uint64_t* bwt, const uint32_t target_len,
const uint32_t* Occ,
const uint8_t c, const uint32_t num)
{
	//Occ(c, num)を計算して返す。
	//cは0,1,2,3のいずれかであると仮定する。ATGC=0123とする。

	if (target_len < num)return 0;//Occ(x, -1) = 0

	//0<=x<=num/64なる任意のxについて、Occ[x*4+c]==Occ(c, x*64+63)が成り立つ。
	//すなわち、numを64で割った余りが63なら、表の値をそのまま返して良い。
	if (num % 64 == 63)return Occ[(num / 64) * 4 + c];

	const uint32_t prev = num < 64 ? 0 : Occ[((num - 64) / 64) * 4 + c];
	const uint64_t t = 1ULL << (num % 64);
	const uint64_t mask = t | (t - 1);
	const uint32_t index = (num / 64) * 3;

	switch (c)
	{
	case 0: return prev + __popcnt64(bwt[index + 0] & ~bwt[index + 1] & mask);
	case 1: return prev + __popcnt64(~bwt[index + 0] & bwt[index + 1] & mask);
	case 2: return prev + __popcnt64(bwt[index + 0] & bwt[index + 1] & mask);
	case 3: return prev + __popcnt64(bwt[index + 2] & mask);
	default: assert(0);
	}
	return -1;
}

 最後に、Occを使って実際に文字列検索を行う関数です。

bool search_BWT(
const uint8_t* string, const uint32_t string_len,
const uint64_t* bwt, const uint32_t target_len,
const uint32_t* C, const uint32_t* Occ,
uint32_t* lb, uint32_t* ub)
{
	//元の文字列の中にstringがあるかを探し、suffix_arrayの下限と上限とをそれぞれlbとubとに格納して返す。
	//返り値は stringがある⇔true とする。
	//stringの各文字は1,2,3,4のいずれかであると仮定する。ATGC=1234とする。

	for (uint32_t i = 0; i < string_len; i++)
	{
		assert(1 <= string[i] && string[i] <= 4);
	}

	(*lb) = 0;
	(*ub) = target_len - 1;
	for (uint32_t i = string_len - 1; i < string_len; i--)
	{
		(*lb) = C[string[i]] + interpolate_Occ(bwt, target_len, Occ, string[i] - 1, (*lb) - 1);
		(*ub) = C[string[i]] + interpolate_Occ(bwt, target_len, Occ, string[i] - 1, (*ub)) - 1;
		if ((*ub) < (*lb))return false;
	}
	return true;
}

結論

 ビット演算は楽しい!


waifu2xを画像可逆圧縮に使ってみた(…が、うまくいかなかった)

 waifu2xは不可逆圧縮されたり縮小されたりした画像を元に戻すためのツールです。deep learningの一種であるCNN(Convolutional Neural Network)が使われており、良い感じに復元されると評判です。原著論文はarXivにあるので自由に読むことができます。
[1501.00092] Image Super-Resolution Using Deep Convolutional Networks
 元画像を縮小してjpgで圧縮した画像をwaifu2xで復元したら一般には誤差が残るわけですが、その誤差とjpg画像を合わせれば可逆圧縮画像とみなせるのではないかと思い、実際に試してみました。

方法

 大まかな手順としては、
(1)テスト画像をMSペイントで開き、縦横50%に縮小してjpgで保存する。
(2)そのjpg画像をwaifu2xで元の大きさに復元する。
(3)テスト画像と復元後の画像との誤差を求め、ビットマップ画像として書き出す。
(4)そのビットマップ画像を7zでLZMA圧縮する。
 となります。前回の記事ではjpgで公開されている画像をテストに使いましたが、今回はpngで公開されている画像としてこれを使ってみました。ただし、縦を切り詰めて1150ピクセルにしました。縦横のピクセル数がともに偶数だと、復元後の画像と元画像との大きさが同じになって実装しやすくなるためです。

www.pixiv.net
 また、waifu2xには派生ソフトがいろいろあるのですが、今回は簡単に操作できるこれを使いました。inatsuka.com
テスト画像と復元後の画像の誤差を求めるコードは以下のようになります。ここ以外は前回の記事のコードとほとんど一緒です。

	//imgが元画像で、img2がimgをjpg圧縮してからwaifu2xで復元した画像とする。
	//imgとimg2の差分を求めてimgを更新する。
	for (int y = 0; y < img->height; y++)for (int x = 0; x < img->width * 3; x++)
	{
		const int pos = y * img->widthStep + x;
		const unsigned char d = (unsigned char)(img2->imageData[pos] - img->imageData[pos]);
		img->imageData[pos] = (char)(d < 128 ? 255 - 2 * d : 2 * d - 256);
	}

結果

 以下のようになりました。waifu2xにはノイズ除去設定が3通りあるのですが、念の為全て試してみました。2列目は誤差画像をLZMA圧縮したファイルの容量であり、jpg画像の容量は含んでいません。jpg画像は175 KBでした。

ノイズ除去 容量(MB)
なし 2.12
1.92
1.93

 元画像をPNGで圧縮すると1.62 MB、JPEG 2000で圧縮すると1.46 MBになりますから、PNGと比較してもかなり悪い結果です。

考察

 waifu2xは確かに滑らかに拡大してくれていますし、jpg特有のノイズは除かれているのですが、jpg圧縮で失われたディテールを完全には復元できていません。実際に復元後の画像と元画像を比較してみると細部のディテールが結構失われていることがわかります。(興味のある人は試してみると良いです) そのため、誤差画像が予想していたより黒っぽくなり、圧縮率が上がりませんでした。あと、分かってはいたのですが、waifu2xはコーデックに使うにはあまりにも計算量が大きすぎます。最初はいくつかの元画像で試そうかと思っていたのですが、計算が重すぎて1枚でギブアップしました。

LOCO-I予測+LZMAで画像を可逆圧縮してみた

 画像を可逆圧縮する方法にはいくつかありますが、多くの場合は(1)ある画素の輝度値を周りの画素から予測する→(2)予測誤差を符号化する という2ステップから成り立っています。今回は、(1)にLOCO-I予測を使い、(2)にLZMAを使って実際に可逆圧縮を行ってみました。

原理

 まず、LOCO-I予測に関してはwikipediaの以下の記事を参考にしました。
Lossless JPEG - Wikipedia
 LOCO-I予測は上・左・左上の3つの画素の値を元に当画素の値を予測するものです。具体的な式はwikipediaの記事にあります。なぜ上・左・左上の三箇所なのかというと、復号するときに左上から順番に復号していくので、当画素を復号しようとする時点では右や下の画素値は未知だからです。
 LZMAアルゴリズムに関しては英語版wikipediaに細かい解説があります。
Lempel–Ziv–Markov chain algorithm - Wikipedia, the free encyclopedia
 LZMAはzipと同じ辞書式圧縮ですが、エントロピー符号化にハフマン符号化でなくRange Coderを使っているなどの理由から、zipより高い圧縮率を期待できます。ちなみにPNGはzipと同じハフマン符号化を使っています。(というか、zipそのものを使っています)

方法

 下記のプログラムを使ってテスト画像をLOCO-I予測し、得たビットマップ画像を7zというフリーソフトLZMA圧縮しました。テスト画像には ご注文はうさぎですか? ホワイトデー特製PC用壁紙(1920x1200)の5人分を試してみました。
スペシャル -TVアニメ「ご注文はうさぎですか?」公式サイト-

結果

名前 PNG (MB) JP2 (MB) 提案手法 (MB)
ココア 3.82 2.29 2.50
チノ 3.92 2.38 2.57
リゼ 4.00 2.38 2.61
シャロ 3.93 2.38 2.58
千夜 3.78 2.25 2.42

 2列目の"PNG"は、元画像をMSペイントで開きPNGで保存した場合の大きさです。PNGはソフトによって圧縮率が変わります。3列目の"JP2"は、GIMPを使って元画像をJPEG2000 losslessに変換したときの大きさです。
 変換後の画像はこんな感じになります。(縮小してあります)
f:id:eukaryo:20150702003648p:plain

考察

 全ての画像の圧縮率でJPEG2000のほうが優れています。LOCO-I予測はJPEG-LSというJPEG2000より古い規格で使われていた方法ですし、エントロピー符号化に関してもJPEG2000はより洗練された方法を使っています。

ソースコード

 C++OpenCVで書きました。

#include <opencv2/core/core.hpp>
#include <opencv2/opencv.hpp>
#include <opencv2/highgui/highgui.hpp>
#include <stdlib.h>
#include <iostream>

using namespace cv;

int encode()
{
	//画像を読み込み、このプログラムに対応している形式の画像かどうか確認する。
	IplImage* img = cvLoadImage("test.png", CV_LOAD_IMAGE_ANYCOLOR | CV_LOAD_IMAGE_ANYDEPTH);
	if (img == NULL)return 0;
	if (img->depth != IPL_DEPTH_8S && img->nChannels != 3)
	{
		cvReleaseImage(&img);
		return 0;
	}

	//LOCO-I予測を行い、残差画像をdestに格納する。
	unsigned char* dest = (unsigned char*)calloc(img->widthStep * img->height * 3, sizeof(unsigned char));
	for (int y = 0; y < img->height; y++)for (int x = 0; x < img->width * 3; x++)//各ピクセルの各色について
	{
		//左、上、左上の画素について、それがあるならその輝度値、ないならゼロを得る。
		const int left = (x < 3) ? 0 : (unsigned char)img->imageData[y * img->widthStep + (x - 3)];
		const int upper = (y == 0) ? 0 : (unsigned char)img->imageData[(y - 1) * img->widthStep + x];
		const int upperleft = (x < 3 || y == 0) ? 0 : (unsigned char)img->imageData[(y - 1) * img->widthStep + (x - 3)];

		//LOCO-I予測を行い、predictに格納する。
		int predict;
		if (max(left, upper) <= upperleft)predict = min(left, upper);
		else if (upperleft <= min(left, upper))predict = max(left, upper);
		else predict = left + upper - upperleft;

		//実際の画素値との差を求めてdestに格納する。
		predict -= (unsigned char)img->imageData[y * img->widthStep + x];
		dest[y * img->widthStep + x] = (unsigned char)((predict + 256) % 256);
	}

	//destの中身をimgに全部コピーする。見やすいように、全単射の写像をかませる。
	for (int y = 0; y < img->height; y++)for (int x = 0; x < img->width * 3; x++)
	{
		const unsigned char a = dest[y * img->widthStep + x];
		img->imageData[y * img->widthStep + x] = a < 128 ? 255 - 2 * a : 2 * a - 256;
	}

	//ウィンドウを生成して画像を表示する。
	cvNamedWindow("test_loco_i", CV_WINDOW_AUTOSIZE);
	cvShowImage("test_loco_i", img);

	//何かのキーが押されたら保存して終了。
	cvWaitKey(0);
	int param[] = { CV_IMWRITE_PNG_COMPRESSION, 9 };
	cvSaveImage("test_loco_i.png", img, param);
	cvSaveImage("test_loco_i.bmp", img);
	cvDestroyWindow("test_loco_i");
	cvReleaseImage(&img);
	free(dest);

	return 0;
}
int decode()
{
	//画像を読み込み、このプログラムに対応している形式の画像かどうか確認する。
	IplImage* img = cvLoadImage("test_loco_i.png", CV_LOAD_IMAGE_ANYCOLOR | CV_LOAD_IMAGE_ANYDEPTH);
	if (img == NULL)return 0;
	if (img->depth != IPL_DEPTH_8S && img->nChannels != 3)
	{
		cvReleaseImage(&img);
		return 0;
	}

	//LOCO-I予測を行い、元画像を復元する。
	for (int y = 0; y < img->height; y++)for (int x = 0; x < img->width * 3; x++)//各ピクセルの各色について
	{
		//左、上、左上の画素について、それがあるならその輝度値、ないならゼロを得る。
		const int left = (x < 3) ? 0 : (unsigned char)img->imageData[y * img->widthStep + (x - 3)];
		const int upper = (y == 0) ? 0 : (unsigned char)img->imageData[(y - 1) * img->widthStep + x];
		const int upperleft = (x < 3 || y == 0) ? 0 : (unsigned char)img->imageData[(y - 1) * img->widthStep + (x - 3)];

		//LOCO-I予測を行い、predictに格納する。
		int predict;
		if (max(left, upper) <= upperleft)predict = min(left, upper);
		else if (upperleft <= min(left, upper))predict = max(left, upper);
		else predict = left + upper - upperleft;

		//かませてあった写像の逆変換を行い、予測誤差を求める。
		int e = (unsigned char)img->imageData[y * img->widthStep + x];
		e = e % 2 ? (255 - e) / 2 : (256 + e) / 2;

		//実際の画素値を求めてimgを更新する。
		predict -= e;
		img->imageData[y * img->widthStep + x] = (unsigned char)((predict + 256) % 256);
	}

	//ウィンドウを生成して画像を表示する。
	cvNamedWindow("test_decode", CV_WINDOW_AUTOSIZE);
	cvShowImage("test_decode", img);

	//何かのキーが押されたら終了。
	cvWaitKey(0);
	cvDestroyWindow("test_decode");
	cvReleaseImage(&img);

	return 0;
}
int main()
{
	////予測誤差を濃淡に変換する写像が全単射であることを確認する。
	//int f[256] = { 0 };
	//for (int a = 0; a < 256; a++)
	//{
	//	int e = a < 128 ? 255 - 2 * a : 2 * a - 256;
	//	int r = e % 2 ? (255 - e) / 2 : (256 + e) / 2;
	//	assert(a == r);
	//	assert(f[e]++ == 0);
	//}

	encode();
	decode();

	return 0;
}

SOR法の加速パラメータωを最適化してみた

 SOR法とは連立一次方程式Ax=bを反復的に解く方法の一つです。SOR法はハイパーパラメータとして加速パラメータωというものを必要とします。教科書的にはωは予め定数として与えておくものらしいですが、今回は誤差が最小になるようなωを各ステップごとに求めながら計算していく方法を実装してみました。

方法

 SOR法のアルゴリズムは簡単なので説明は控えます。
SOR法 - Wikipedia
今回の実装では(式が簡単になるので)ガウス・ザイデル法ではなくヤコビ法を元にして実装しました。よって漸化式は
x_i^{(k+1)} = x_i^{(k)} + \omega\frac{1}{a_{ii}} \left( b_i-\sum_{j=1}^n a_{ij}x_j^{(k)} \right)
となります。ここで、ωを変数としたままx^{(k+1)}を求めると各要素がωの一次式になるので、誤差の二乗ノルム\left\| Ax^{(k+1)}-b\right\|はωの二次式で表されます。その二次式の最小値を\omega_{opt}としてx^{(k+1)}の各要素に代入すればよいというわけです。

結果

 結果として、ωを決め打ちするSOR法と比べて反復回数は少なくなりました(当たり前ですが)。ただ、\left\| Ax^{(k+1)}-b\right\|を求める際に行列ベクトル積を2回計算する必要があるのですが、そのデメリットを補えるほど加速されるかは場合によります。私が試したところ、非対角成分を小さな負の値にしたときに\omega_{opt}が2より遥かに大きくなるステップが見つかり、その場合はこの方法のほうが結果として速くなりました。

#include <stdio.h>
#include <conio.h>
#include <assert.h>
#include <random>

void init_matrix(double* src, const int dim, const bool diag_domi, const double lowerbound, const double upperbound)
{
	//srcは[dim*dim]の配列であると仮定する。これをdim次正方行列とみなして、
	//対角成分が全て正であるような対称行列になるように適当な値を代入して返す。
	//diag_domiがtrueなら狭義対角優位にする。

	assert(lowerbound < upperbound);

	std::mt19937 gen(1729);
	std::uniform_real_distribution<> ra(lowerbound, upperbound);

	for (int i = 0; i < dim; i++)
	{
		for (int j = 0; j < i; j++)src[i * dim + j] = src[j * dim + i];//対称行列なので
		for (int j = i; j < dim; j++)src[i * dim + j] = ra(gen);
	}

	//対角成分の処理。対角優位にしない場合でも、対角成分で割る処理があるので非ゼロの値にしておく。
	for (int i = 0; i < dim; i++)
	{
		double sum = 0.0;
		for (int j = 0; j < dim; j++)sum += abs(src[i * dim + j]);
		src[i * dim + i] += diag_domi ? sum + ((double)dim) * 0.01 : abs(src[i * dim + i]) + 1.5;
	}
}
void init_vector(double* src, const int dim)
{
	std::mt19937 gen(1729);
	std::uniform_real_distribution<> ra(0.0, 10.0);

	for (int i = 0; i < dim; i++)src[i] = ra(gen);
}
void mul(const double* mat, const double* vec, const int dim, double* dest)
{
	for (int i = 0; i < dim; i++)
	{
		dest[i] = 0.0;
		for (int j = 0; j < dim; j++)dest[i] += mat[i * dim + j] * vec[j];
	}
}
void sub(const double* vec1, const double* vec2, const int dim, double* dest)
{
	for (int i = 0; i < dim; i++)dest[i] = vec1[i] - vec2[i];
}
void mad(const double* vec1, const double p, const double* vec2, const int dim, double* dest)
{
	for (int i = 0; i < dim; i++)dest[i] = vec1[i] * p + vec2[i];
}
void Jacobi(const int dim, const bool diag_domi, const double epsilon, const double lowerbound, const double upperbound)
{
	assert(0 < dim);

	double* A = (double*)malloc(sizeof(double) * dim * dim);
	double* x = (double*)malloc(sizeof(double) * dim);
	double* b = (double*)malloc(sizeof(double) * dim);

	init_matrix(A, dim, diag_domi, lowerbound, upperbound);
	init_vector(b, dim);
	for (int i = 0; i < dim; i++)x[i] = b[i];

	double* x1 = (double*)malloc(sizeof(double) * dim);
	double* v1 = (double*)malloc(sizeof(double) * dim);
	double* v2 = (double*)malloc(sizeof(double) * dim);

	//連立一次方程式Ax=bをJacobi法で解く。

	for (int loop = 0; loop < 1000; loop++)
	{
		//この時点で、近似解はxに格納されている。

		//ヤコビ法を適用する。
		for (int i = 0; i < dim; i++)
		{
			x1[i] = b[i];
			for (int j = 0; j < dim; j++)if (i != j)x1[i] -= A[i * dim + j] * x[j];
			x1[i] /= A[i * dim + i];
		}
		for (int i = 0; i < dim; i++)x[i] = x1[i];

		//誤差ノルム||Ax-b||を求めてerrorに格納する。
		mul(A, x, dim, v1);
		sub(v1, b, dim, v2);
		double error = 0.0;
		for (int i = 0; i < dim; i++)error += v2[i] * v2[i];

		printf("反復%d回目  Jacobi  error=%f\n", loop, error);
		_getch();

		//新しい近似解が目標の精度に達しているなら終了。
		if (error < epsilon)break;
	}

	free(A);
	free(x);
	free(b);
	free(x1);
	free(v1);
	free(v2);
}
void SOR(const int dim, const bool diag_domi, const double epsilon, const double lowerbound, const double upperbound, const double omega = -1.0)
{
	assert(0 < dim);

	double* A = (double*)malloc(sizeof(double) * dim * dim);
	double* x = (double*)malloc(sizeof(double) * dim);
	double* b = (double*)malloc(sizeof(double) * dim);

	init_matrix(A, dim, diag_domi, lowerbound, upperbound);
	init_vector(b, dim);
	for (int i = 0; i < dim; i++)x[i] = b[i];

	double* x1 = (double*)malloc(sizeof(double) * dim);
	double* v1 = (double*)malloc(sizeof(double) * dim);
	double* v2 = (double*)malloc(sizeof(double) * dim);
	double* v3 = (double*)malloc(sizeof(double) * dim);

	//連立一次方程式Ax=bをSOR法で解く。

	for (int loop = 0; loop < 1000; loop++)
	{
		//この時点で、近似解はxに格納されている。

		//加速パラメータωとしてSOR法を適用後の新しい近似解ベクトルx'の値を求める。x'=x1*ω+xとする。
		for (int i = 0; i < dim; i++)
		{
			x1[i] = b[i];
			for (int j = 0; j < dim; j++)x1[i] -= A[i * dim + j] * x[j];
			x1[i] /= A[i * dim + i];
		}

		double error = 0.0;
		double omega_using = omega;

		if (0.0 < omega && omega < 2.0)//引数omegaに有効な値が指定されているならそれを使う。
		{
			for (int i = 0; i < dim; i++)x[i] += omega * x1[i];

			//誤差ノルム||Ax-b||を求めてerrorに格納する。
			mul(A, x, dim, v1);
			sub(v1, b, dim, v2);
			for (int i = 0; i < dim; i++)error += v2[i] * v2[i];
		}
		else//最適なomegaを求める。
		{
			//ベクトルx'の各成分はωの1次式であるが、ここで||Ax-b||を求める。
			//||Ax-b||はスカラーでありωの2次式で表される。||Ax-b||=p2*ω^2+p1*ω+p0とする。
			mul(A, x1, dim, v1);
			mul(A, x, dim, v2);
			sub(v2, b, dim, v3);
			double p2 = 0.0, p1 = 0.0, p0 = 0.0;
			for (int i = 0; i < dim; i++)
			{
				p2 += v1[i] * v1[i];
				p1 += v1[i] * v3[i] * 2.0;
				p0 += v3[i] * v3[i];
			}

			//||Ax-b||が最小となるようなωの値を求める。
			omega_using = -p1 / (2.0 * p2);

			//そのωのときの||Ax-b||を求める。
			error = p0 - ((p1 * p1) / (4.0 * p2));

			//そのωを使って新しい近似解を求め、xに格納する。
			mad(x1, 1, x, dim, v1);
			for (int i = 0; i < dim; i++)x[i] = v1[i];

		}

		printf("反復%d回目  ω=%f  error=%f\n", loop, omega_using, error);
		_getch();

		//新しい近似解が目標の精度に達しているなら終了。
		if (error < epsilon)break;
	}

	free(A);
	free(x);
	free(b);
	free(x1);
	free(v1);
	free(v2);
	free(v3);
}
int main()
{
	printf("JacobiとSOR(ω=1.0)が同じ結果になることを確認する。\n");
	Jacobi(100, true, 0.00001, -1.0, 1.0);
	SOR(100, true, 0.00001, -1.0, 1.0, 1.0);

	printf("ωを最適化させれば速く収束することを確認する。\n");
	SOR(100, true, 0.00001, -1.0, 1.0);

	printf("ω最適化がうまく働く例。\n");
	SOR(100, true, 0.00001, -1.0, 0.0);

	printf("おわり");
	_getch();
	return 0;
}