C++で2種類のコンテナを同時にソートする

C++でデータを扱っていると，2種類のコンテナを同時にソートしたい場合が時々出てくる．どういうことかというと，

int main() {
    int keys[]   = {5, 3, 1, 2, 4};
    int values[] = {1, 2, 3, 4, 5};
    aoba::sort_by_key(keys, keys+5, values);
    for(int i=0; i<5; ++i) {
        std::cout << keys[i] << " ";   // 1 2 3 4 5
    }
    std::cout << std::endl;
    for(int i=0; i<5; ++i) {
        std::cout << values[i] << " "; // 3 4 2 5 1
    }
    std::cout << std::endl;
    return 0;
}

と，片方のソート結果を元にもう片方をソートするような感じ．これはCUDAライブラリであるthrustに標準搭載されており，sort_by_keyという名前がついている．

さて，こういうときに一からソートアルゴリズムを組み直すなんてことは馬鹿らしい．じゃあどうするか．基本的には既存のstd::sortに渡すイテレータを変形することで，なんとか解決しようとする．そうなると，boostにzip_iteratorが存在することに気付けば，イテレータをboost::tupleの組で表現することによって，まとめてソートさせればいいような気がしてくる．こんな風に:

std::sort(
    boost::make_zip_iteator(
        boost::make_tuple(
            keys_first, values_first)),
    boost::make_zip_iterator(
        boost::make_tuple(
            keys_last, values_last)),
    compare_functor());

しかし，そうは問屋が降ろさない．何故なら，boost::zip_iteratorはWritableではない．これはソートアルゴリズムの要件を満たしていないため，適用させようとしてもコンパイルエラーが発生してしまう！

それじゃあどのようにして解決するかというと，一番楽な方法としてはWritableなコンセプトを満たしたzipイテレータを新しく作ってやって，それを用いてstd::sortを行うのがいい．とはいっても，C++でそれをやるのは若干手間のかかる作業になる．具体的にはこんな感じで実装した:
実装したといっても，実際には『Sorting two arrays simultaneously』の内容を若干手直しするだけなので基本的にはそんなに大変じゃないし，ただ利用するぶんには単純にsort_by_key()を呼び出すだけでよい．

しかし，なんでこんなにC++は行数が多くなってしまうのか．haskellなんか一行で済むのに．

超一様分布列をBoost.Randomを用いて実装する

ご存知の通り、C++のライブラリには高性能な擬似乱数生成エンジンBoost.Randomが存在している。Boost.Randomは擬似乱数を生成するエンジン部分と、一様分布、ガウス分布、ベルーヌーイ分布など好みの分布を生成する出力部分の2つで構成されている。

しかし、ことモンテカルロ法においては主に収束速度を早める目的として、擬似乱数ではなく低食い違い量列(Low-Discrepancy Sequence, LDS)と呼ばれる数列が用いられる場合がある。これを用いると低次元での積分計算の際に収束速度を大きく改善することができ、具体的には従来のモンテカルロ法で(1/√N)だった収束速度が約(1/N)くらいにまで落とすことができる。

さて、このような低食い違い量列をBoost.Randomで用いるためにはエンジン部分を製作すればいいだけで、実際には以下のような形で実装すれば良い。今回は簡単なHalton列を用いて実装を行った:

具体的には、対象のファンクタがUniform Random Number Generatorのコンセプトを満たすためには以下の記述が必要となる:

• X::result_typeにoperator()が返す型を指定
• operator()にエンジン部分を記述
• min(), max()にエンジンが生成する分布の最小値、最大値を指定

さらに、対象のファンクタが上述の内容に加えPseudo-Random Number Generatorのコンセプトを満たすためには:

• X()にデフォルトの状態で生成されるようなコンストラクタを指定
• X(...)にユーザ指定の状態で生成されるコンストラクタを指定
• seed(...)にエンジンの状態を変更するためのメンバ関数を指定

するだけで良い。後はBoost側が良きに計らってくれる。非常に便利だと思う。

余談
C++0xには同等のstd::randomが付属しているため、これを利用すればいいのではという声もあるとは思うが、std::uniform_real_distribution<>にこのLDSを食わせたところ、どうみてもLDSではない値が出力される結果となった。恐らくresult_typeがdoubleという一般の擬似乱数ではない型を指定しているために起こる現象だと思うが、ソースを見ていないので一先ずは従来のBoost.Randomの手法を用いることにする。

Restricted Boltzmann Machineの学習手法についての簡単なまとめ

近年の機械学習ではDeep Learningと呼ばれる分野が一世を風靡しています．コンピュータビジョンや自然言語処理，音声認識などの分野では何らかの問題を解こうとした際に，まず対象の入力データからSIFTやケプストラムといった何らかのアルゴリズムを用いて特徴ベクトルを抽出し，ごりごりと判別していくといった流れが一般的です．しかし，その特徴ベクトルを生成するという生のデータから本質となる部分を抽出するアルゴリズム自体は研究者が一生懸命考えながら作るのが普通でした．

Deep Learningの分野で最も有名な手法の一つであるDeep Belief Nets(DBN) [Hinton06]は，研究者がアルゴリズムを作るのではなく，それ自体も機械学習にやらせましょうという動機で生まれたアルゴリズムです．DBNではまるで一昔前にやたら流行ったニューラルネットワークのように各ノードを層状に配置し，それらのパラメータをRestricted Boltzmann Machine(RBM)と呼ばれるモデルを用いて学習することによって自動的に特徴ベクトルを生成します．これが既存手法を大きく超える精度を出してしまったためにDeep Learningは一躍有名となり，ニューラルネットの再来か，機械学習は人工知能の夢を見るか，みたいに言われているような状況です．最近ではGoogleがネコのニューロンを自動的に学習したとかでニュースになっていました．

さて，そのような状況にもかかわらず，そのDeep Belief Netsやその改良版であるDeep Boltzmann Machineを学習するためには，まずその基礎となっているRBMについて学ぶ必要があります．しかしそのRBMですらお世辞にもわかり易いものとは言えず，しかもその最終的な更新式やその導出過程についてきちんと書いている文献は今まで確認したところ存在していません(*1)．「チュートリアルや論文だけでは良く分からないので，実際に手を動かして確認したい．だがその計算過程は公開されていないため，結局良く分からないままで終わる」ことが多いように思われました．

ということで，RBMの計算過程について書いた個人用のpdfファイルを公開してみることにします．外部に見せるということで多少文章を詳しくしましたが，これ単体で学ぶのではなく他の資料と比較しながら見るとベターかと思います．多分数式は間違ってないと思いますが，何かありましたらご一報頂けると幸いです．

http://dl.dropbox.com/u/2048288/RestrictedBoltzmannMachine.pdf

*1 Deeplearning.netの式は比較的詳しいですが，実際に計算したところ微妙に間違えている箇所を数ヶ所発見しています

直積量子化(Product Quantization)を用いた近似最近傍探索についての簡単な解説

"aka motsu-nabe" by chatani

概要
冬の寒さも一段と厳しくなってまいりました。おでんや鍋が恋しくなる季節です。

さて、最近ようやっと一仕事が終わりまして、長ったらしい記事が書けるようになりました。ですので、今回は2011年にTPAMIで発表された、近似最近傍探索についての論文『Product quantization for nearest neighbor search』について簡単に紹介したいと思います。

この論文は2011年に発表された、最近傍探索アルゴリズムの決定打です。シンプルな理論でありながら既存手法を打ち破るほどの強力な性能を有し、速度も非常に高速、かつ省メモリなのでスマートフォンに載せ、リアルタイムで動作させることも可能です。

以前この手法はCV勉強会@関東で紹介されたらしいのですが、具体的に紹介しているページは(最近すぎるので当たり前ですが)現在のところあまり見かけていません。ただ個人的にこの手法はあと5年くらいは本線で戦えるものと思っておりますので、きちんと解説した記事を出すことはある程度有意義なものかなぁと、そういう風に考えています。

最近傍探索の課題
最近傍探索は最も古典的な手法でありながら、現在においても頻繁に利用されているアルゴリズムの一つです。近年のプロセッサの性能向上によって大規模なデータ処理が家庭のパソコンでもできるようになり、コンピュータビジョンにおいても数億件規模(10^8~9)の高次元ベクトルデータを最近傍探索を用いて処理することで、良い精度の結果を得ることが当たり前の時代になってきました。

ただそうは言っても、最近傍探索をそのまま用いることは現実的ではなく、いくつかの課題が浮上してきます。具体的には以下のとおりです:

• いかに圧縮した状態でデータを格納するか: ベクトルデータをそのままの状態で保持することはメモリ容量の圧迫へ繋がります。例えば128要素のベクトルを考えた場合、floatを使用した場合でも32×128=4096bitを1エントリで消費してしまいます。IOアクセスが生じてしまうと速度は一気に低下してしまうため、ベクトルデータはすべてメモリ上に展開しておくことが望ましいー
  しかし、今回我々が扱いたいエントリ数というのは数億とかそういったオーダーです。そのオーダーのベクトルをそのまま格納するのは現在のメモリ容量だとなかなか厳しい。ですので、情報量を保ったままベクトルをいかに圧縮できるかが、実用的な面において重要な課題となります。
• いかに高速に解を見つけられるか: 次元数の多いベクトルの距離計算は基本的に低速です。さらにそのまま実装してしまうとすべてのエントリに対して距離の計算をする必要があるため速度はO(N)となり、実用的ではありません。1エントリあたりの計算コストを削減したり、そもそも明らかに異なる(クエリより離れすぎている)エントリを計算しないことで、計算を高速に行うことも同じく重要となります。

これらの課題を克服するため、近年では近似最近傍探索と呼ばれる手法が盛んに研究されております。その中でも有名なアルゴリズムはLSH(Locality-Sensitive Hashing)であり、ハッシュ関数によって抽出される近傍エントリのみを計算することで、高速に最近傍探索を行います。ただこの手法は計算用のハッシュを別に格納する構造であるため、全体のサイズは減少するのではなく増大してしまいます。

また、Preferred Researchで紹介されているMinHashは、ベクトルの要素数が定まっておらず、成分が{0,1}の2値しかとらないような特殊なベクトルを計算する際においては有効な手法でありますが、そうでない場合においてはいったんベクトルをバイナリ表現に落とし込まなければならないので一般的に精度は低下してしまいます。今回紹介する直積量子化は、対象のベクトルが高次元であり、かつユークリッド距離による最近傍探索を行いたい場合において威力を発揮する手法です。

スカラ量子化
直積量子化について説明を行う前に、まず量子化とは何かについて説明していきます。物理をかじったことのある方ですと『量子化 = 物理量を演算子へ変換する』ものと思っちゃいますけど、そっちの量子化ではなくある数を特定のビット列で表現することの量子化です。

分割された領域に、特定のビットを割り当てる

例えば、0〜1の間をとる数を2ビットで表現しろという問題を考えた場合、普通は上図のように0〜1を4分割し、それぞれの数値を割り当てることで表現を行います。このように、少ないビット数で対象の数値を表現することを一般に『スカラ量子化(Scalar Quantization)』といいます。スカラ量子化によってfloat型(32bit)の数値は、精度を犠牲にすることによって僅か2bitで表現することができるようになりました。

ベクトル量子化
ところで、このスカラ量子化は対象のデータが数値であった場合には効率的ですが、ベクトルの場合だと効率的とは言えません。

スカラ量子化では、ベクトルを効率良く表現することができない

例えば上の図について考えてみましょう。この点群を1エントリ辺り1ビットで表現しようとした場合、スカラ量子化によって各成分を量子化してしまっては最低でも2ビットが必要となりますが、物事を空間的に考えることで発想を転換します。つまり、それぞれのクラスタを代表する点を予めコードブックとして保持しておき、各点は一番近い代表点のインデックスを保持しておくのです。これによってそれぞれの点は僅か1ビットで表現することができるようになりました。このような手法のことを一般に『ベクトル量子化(Vector Quantization)』といいます(*1)。

ベクトル量子化によって、ベクトルを2つの代表点によって表せるようになった

*1: 本手法のクラスタリングにはk-means algorithmを使用しています。k-meansに関しては『クラスタリングの定番アルゴリズム「K-means法」をビジュアライズしてみた』が視覚的に分り易いのでそちらを参照してください。

直積量子化
ベクトル量子化はスカラ量子化よりも少ないビット数で効率良くベクトルを表現できるのですが、同時に欠点もあります。次元が増大するにつれて使用するコードブックが指数的に増えてしまうのです。

これは次元の呪いに起因するものです。数次元〜数十次元くらいだったらそんなに問題とはならないんですけど、数百とか数千次元になってくるととんでもない量のコードブックが必要となり現実的でない。『直積量子化(Product Quantization)』はこの次元の呪いを解決するために用いられる、スカラ量子化とベクトル量子化の中間を行く手法です。

高次元のベクトルを直積表現によって分解する

具体的な手法は下記の通りです。まず、N次元のベクトルをM分割し、N/M次元のサブベクトルを生成します。そして、そのサブベクトルをそれぞれベクトル量子化することによって、インデックスを(i1, i2, ..., iM)のタプルで表現します。そうすれば次元の呪いは発生せず、データ表現も比較的低ビットで抑えられます。

例として128次元のベクトルを考えてみましょう。まず、このベクトルを16分割することで8次元のサブベクトルへ落とし込みます。そして、そのサブベクトルをベクトル量子化によって8bit表現に持ってきた場合、対象のベクトルは8×16=128bitで表現でき、かつ使用するコードブックサイズも16×256=4096と少数で済みます。float型のサイズで考えると30倍以上の圧縮効率です。このように、少数のコードブックサイズで効率よくベクトルを表現できることが直積量子化の利点です。

距離計算
さて、直積量子化によって各ベクトルのエントリを圧縮した状態で格納したとしましょう。あるクエリベクトルが与えられたとして、そのユークリッド距離が最小となるエントリを見つけるにはどうすればいいでしょうか？

幸運なことに、ユークリッド距離の計算は各要素毎に独立しています。つまり、ただ単にクエリベクトルを各サブベクトルに分け、それぞれの領域で距離を計算し、その和を計算するだけで良いのです。よって、直積量子化を用いて最近傍探索を近似的に行うためには:

1. クエリベクトルをM分割し、N/M次元のサブベクトルに分解する
2. クエリのサブベクトルと、それぞれのコードブックに記された代表ベクトル間の距離を予め計算しておく
3. それぞれのエントリがどのインデックスに位置しているのか確認し、2で計算した距離を参照して足し算を行う
4. 最小距離をもつエントリを候補として提出する

ことで実現できます。これが直積量子化を用いた、近似最近傍探索の核となるアルゴリズムです。

詳細について
この距離計算は非対称性からある程度のバイアスがかかっているため、実際の距離を正確に見積もるためにはこのバイアスを取り除く必要があったりします。また、具体的な距離計算がコードブックだけで済むので省力化されているとはいえ、計算量は依然としてO(N)であり、そこらへんも大規模なデータ数では足を引っ張る原因となります。

んじゃあこれをどうやって解決するのかについては、現論文

http://lear.inrialpes.fr/pubs/2011/JDS11/jegou_searching_with_quantization.pdf

を参照してください。上記2つの懸念についての解決案が記述されており、上の知識があれば多分すらすらと読めるはずです。
本当はこれも解説しなきゃなんないんですけど、書いている内に疲れたちゃったのでとりあえずコアなとこだけということで勘弁してください。気が向いたら更新します。

OpenCVアルゴリズムを高速に記述できる補助ライブラリcvutilを公開します

"The Dunny Collective - Brisbane River" by KayVee.INC

概要
OpenCVを用いてアルゴリズムを組む際、ある画像からある画像へ変換処理を行ったり、関数fを用いて新規に画像を生成したり、画像の結果を1つの値へと集約するなどの定型的な処理は非常に多く行われておりますが、並列処理を行いかつ速度を考慮したプログラムを組もうとするとある程度のコード量を必要とするためそれなりの手間となります。また画像処理によくあるループ文のネストが大量に発生するため、普通に組んでしまうと一見何をやっているのか非常に分かりづらいコードとなってしまいます。

cvutilはこのような冗長な作業を定型化し、複数のアルゴリズム関数とC++0xに新しく搭載されたラムダ式を使用して高速にアルゴリズムを記述することができます。このようにして作られたアルゴリズムは自動的にOpenMPによって並列化されるため、研究用途などの速度をある程度担保しつつもメンテナンスのし易いコードを記述するための手助けとなります。

インストール
cvutilはgithubにて公開しております。ダウンロード先は下記のURLを参照してください。

https://github.com/rezoo/cvutil

cvutilはヘッダオンリーのライブラリです。バージョンがOpenCV 2.x系統であれば、対象のincludeディレクトリにcvutilをコピーするだけで使用できます。
C++0xのラムダ式を用いることで端的に記述することができますが、関数オブジェクトを使用すればある程度古いコンパイラでも正常に動作するでしょう(そこまでして使う利点があるのかどうかは分かりませんが)。

ちなみに、元々は自分用に作られた適当なライブラリなので、バージョンアップに伴ってある程度破壊的な変更を行うこともありますのでご了承ください。その際はバージョンを0.1上げることによって対応します。

アルゴリズム一覧
ここではcvutilに搭載されているアルゴリズム一覧の一部を紹介いたします。これらのアルゴリズムはSTLのそれと似ているため、C++を触ったことのある方は名前から直感的に機能を類推できるでしょう。

transform_image
transform_image関数は対象の画像を、任意の関数を通して出力先へ変換します。cvutilの中では最も汎用性の高い関数となります。

cvutil::transform_image(src, dst, [](uchar x) {
    return cv::saturate_cast<uchar>(x*2); });

この関数は複数の画像を高度に合成するような場合においても利用できます。アリティ(渡される引数の数)によって、最も適したアルゴリズムが自動的に呼び出されます。

cvutil::transform_image(src1, src2, dst, [](uchar x, uchar y) {
    return (x + y)/2; });

generate_image
generate_image関数は任意の関数を用いて、出力先に画像を新しく生成します。generate_imageをそのまま用いることは通常ありませんが、x, y座標を元に画像を生成する関数generate_image_xyや、画像中央に原点がある右手系の座標系u, v座標を元とした関数generate_image_uvは比較的使い勝手が良いと思います。

int width = 512;
int height = 512;
cv::Mat_<cv::Vec3b> dst2(height, width);
cvutil::generate_image_xy(dst2, [=](int x, int y) {
    return cv::Vec3b(0, x*256/width, y*256/height); }); // BGR

reduce_image
reduce_image関数は対象の画像を、任意の関数を通して1つの結果へと集約し、その値を返します。

cv::Mat_<uchar> src3(3, 3);
src3(0, 0) = 0; src3(0, 1) = 1; src3(0, 2) = 2;
src3(1, 0) = 3; src3(1, 1) = 4; src3(1, 2) = 5;
src3(2, 0) = 6; src3(2, 1) = 7; src3(2, 2) = 8;
int result = cvutil::reduce_image(src3, (int)0, std::plus<int>()); // 36

上のサンプルは全要素に対しての足し算を行っています。transfromの並列化は兎も角、reduceの並列化は一般的に面倒なのですが、reduce_imageはそこらへんの処理もきちんと並列化してくれます。

transform_reduce_image
transform_reduce関数は1枚または複数の画像を対象の関数で1つの値へ変換した後に、reduce関数を適用します。2種類の結果を合併した集約結果を返す場合などに役立つと思います。

minmax_element_image
minmax_element_image関数は画像の最小値と最大値を返します。

std::pair<uchar, uchar> minmax_result =
    cvutil::minmax_element_image(src); // (min, max)

integrate_image
integrate_image関数は画像を2次元に積分します。具体的には、inclusive_scan関数を2次元的に適用します。バイナリー関数にstd::plusを適用した場合、この関数は積分画像の生成と等価です。また、画像の型が対象の二項演算子に関してアーベル群を成していた場合、任意の矩形領域に関してのreduceはintegrate_imageによって生成された画像を用いて、定数時間による計算が可能となります(参考: 並列環境におけるreduceとscanアルゴリズム)。

OpenCVにも積分画像を生成する関数は存在しているのですが、この関数は並列化を行いかつ任意の二項演算子を適用できる点が異なります。

cv::Mat_<uchar> src3(3, 3), dst3(3, 3);
src3(0, 0) = 0; src3(0, 1) = 1; src3(0, 2) = 2;
src3(1, 0) = 3; src3(1, 1) = 4; src3(1, 2) = 5;
src3(2, 0) = 6; src3(2, 1) = 7; src3(2, 2) = 8;
cvutil::integrate_image(src3, dst3, std::plus<uchar>());
// [0, 1, 2;     [0,  1,  3;
//  3, 4, 5;  =>  3,  8,  15;
//  6, 7, 8]      9,  21, 36]

cvutilはその他にも色々と関数がありますが、多分名前からどんなことやってるのか大体は理解できると思いますので、とりあえず代表的な関数をいくつか紹介いたしました。

License
自分用に作っているだけなのですが、一応ライセンスはMIT Licenseを適用します。誰でも自由に無制限に使ってかまいません。

とりあえずはそんな感じです。もし何か使っていく上で質問や疑問、あるいはご指摘等ございましたらお気軽にコメントまたはメールをお願いします。

追記
お気づきの通り、本来の文脈ならば、アルゴリズムに渡す引数はcv::Mat_<T>ではなく何かしらの抽象構造Viewを渡すべきなのです。C++のSTL Algorithmはコンテナでなくイテレータという抽象に依存することによって、任意のデータ構造による適用が可能となり、さらにtransformed_iteratorなど、変形を遅延させるような構造を取ることができる。つまり、cvutilにおいても本来はそのような『抽象に依存する』データ構造を起点にアルゴリズムを組んでいくべきなのです。

// こんな感じで書けたら最高なんだけど…
cvutil::transform_image(
    cvutil::make_resized_view(src, dst.size()),
    dst,
    any_functor);

ただこのアナロジーを画像に適用しましょうとなるとなかなか難しく、現状としてはboost.GILしか存在していない。そういった意味で確かにboost.GILは魅力的なライブラリなのですが、GILに適用するalgorithmが殆ど存在していませんし、なによりコンパイル時間がboooooooooooostしてしまう。ぱぱっと書いてぱぱっと書き捨てられるような構造がむしろ研究にとっては重要なわけで、そんな面倒なことするよりも多少汚いOpenCVを採用したほうが使い勝手としては向上するよねという、まぁそういった言い訳です。