映像奮闘記

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

2011-11-24T20:11:00.001+09:00

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

さて、最近ようやっと一仕事が終わりまして、長ったらしい記事が書けるようになりました。ですので、今回は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倍以上の圧縮効率です。このように、少数のコードブックサイズで効率よくベクトルを表現できることが直積量子化の利点です。

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

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

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

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

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

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

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

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

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

2011-09-24T01:44:00.000+09:00

"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を採用したほうが使い勝手としては向上するよねという、まぁそういった言い訳です。

OpenCVの画素値アクセスについて

2011-09-16T14:45:00.000+09:00

"Fire King" by A Continuous Lean

何度も何度も繰り返されている議題ですが、今回は画像処理に良くある『全画素を変換する』という処理のパフォーマンスについて考察します。
あなたはOpenCVのcv::Mat_<T>を使って、全ての画素値をある関数funcで変換するという処理を行いたいとします。それで、たぶんなんですけど、普通に何も考えず組むとしたらこんな感じのコードになると思います:

for(int y=0; y<src.rows; ++y) {
    for(int x=0; x<src.cols; ++x) {
        dst(y, x) = func(src(y, x));
    }
}

なんの変哲もない普通のループです。このコードは普通に動きます。おわり。ちゃんちゃん。

…そういきたいのですが、パフォーマンスを気にする方々はここでいろいろとした疑問が沸き起こってくるわけです。それじゃあどこが悪いのか。1つずつ疑問点を上げ、その上で上記のコードがある程度の合理性を持っているコードであることを示していきます。

疑問1: 関数呼び出しを行っているからその分遅い？
src(y, x)という処理は結局cv::Mat<T>::operator()という関数を呼び出していることに相当します。なので、関数呼び出しを数十万のピクセルに対して行っているため、関数呼び出しにかかるコストが無視できないレベルで発生していくのではないか？？

回答: 遅くならない
OpenCVのcv::Mat_<T>::operator()はインライン展開されます。なので、関数の呼び出しに対して発生する速度的なコストは発生しません。安心して使いましょう。

疑問2: 掛け算を余計に行なっているために遅い？
結局cv::Mat_<T>::operator()の行っていることはコードに直すとこんな感じです。

src(y, x) == ((T*)(src.data + src.step.p[0]*y))[x]

でも、この項をよく見てみると、xのループに対してyの値は固定されているため、第二項の計算は正直不要でしょう。なのに掛け算の計算でコストが余分にかかってしまうため、速度が遅くなるのでは？？

回答: 遅くなるけどそんなに問題じゃない
確かに遅くなりますが、そこまで問題となる箇所ではありません。時代は3GHzに突入しているのです。整数の掛け算にかかるコストなんてたかが知れています。

そんなところよりも、別の箇所を改善したほうが大抵の場合劇的に向上します。仮に、funcの中にstd::expやstd::logが含まれていたとすると…

dst(y, x) = std::exp(-(src(y, x) - mu)*(src(y, x) - mu)); // ぐぎゃー
dst(y, x) = -src(y, x)*std::log(src(y, x)); // おもいー

…それだけで大分重くなってしまいます。sqrt, sin, cosはまだ我慢できるにしても(それでも掛け算よりは結構重い)、指数、対数関数は普通に使うとかなり重いんです。CV系統などの精度を求めないような計算でしたら、スプライン関数などの近似式に置き換えることを検討しましょう。

at<T>(), operator()は範囲チェックを行ってくれる
昔のOpenCVにはat<T>()と同じような処理をしてくれるCV_IMAGE_ELEMマクロなんて言うものもありましたが、これはもう時代遅れの代物です。なぜなら、cv::Matのat<T>(), cv::Mat_<T>のoperator()は与えられたx, y座標が範囲内のものであるかどうか、きちんとチェックを行ってくれるからです。範囲外である場合には警告を送ってくれます。そういった意味でもCV_IMAGE_ELEMよりatやoperator()を利用していくべきでしょう。

『余計な範囲チェックを行っているために遅くならないのか？』という疑問もあるかもしれませんが、この範囲チェックはReleaseモードで消去されます。なので実際の動作には全く影響されないのです。

cv::Mat::at<T>(), cv::Mat_<T>::operator()は早い
at<T>(), operator()はきちんと考えられているメンバ関数です。ポインタを用いてがりがり書いたコードよりは少し遅いかもしれませんが、この遅さは全体のアルゴリズムに影響するような遅さではありません。ほっと胸を撫で下ろし、安心して使いましょう。

でもやっぱり速度が気になる
とはいっても、きちんと書かれていないという神経質な方もいるかもしれません。at<T>()なんて負けた気がして使いたくない……その場合のコードは大体以下のようになります。

for(int y=0; y<src.rows; ++y) {
    uchar* src_ptr = src.ptr<uchar>(y);
    uchar* dst_ptr = dst.ptr<uchar>(y);
    for(int x=0; x<src.cols; ++x) {
        dst_ptr[x] = func(src_ptr[x]);
    }
}

これだけじゃあ高速化されたといっても微々たるものなので、せっかくですからOpenMPによる並列化を行ってみましょう。さらにassertでサイズが合っているかチェックを行うようにして…そうすると最終的なループはこんな感じになります。

assert(src.size == dst.size);
#pragma omp parallel for
for(int y=0; y<src.rows; ++y) {
    uchar* src_ptr = src.ptr<uchar>(y);
    uchar* dst_ptr = dst.ptr<uchar>(y);
    for(int x=0; x<src.cols; ++x) {
        dst_ptr[x] = func(src_ptr[x]);
    }
}

…大分行数が増えてしまいました。ただ自分がやりたいことは全画素に対してfuncを適用させたい、ただそれだけなのに、ここまで行数を書く必要があるのです。しかもぱっと見で何をしたいのか良く分からない。面倒ですし分かりにくい。

ここで、C++を少しかじったことのある人でしたら、ユニタリー関数を適用させる表式のアルゴリズム関数を適用すればいいんじゃね？といった考えに至るのは自然な発想でしょう。こんなふうに:

template<typename SrcType,
         typename DstType,
         typename UnaryFunction>
void transform_image(const cv::Mat_<SrcType>& src,
                     cv::Mat_<DstType>& dst,
                     UnaryFunction unary_op) {
    assert(src.size == dst.size);
    #pragma omp parallel for
    for(int y=0; y<src.rows; ++y) {
        // equivalent to ptr<T>(y)
        const SrcType* src_ptr = src[y];
        DstType* dst_ptr = dst[y];
        for(int x=0; x<src.cols; ++x) {
            dst_ptr[x] = unary_op(src_ptr[x]);
        }
    }
}

上の形の式をヘッダファイルあたりに定義してやれば、後は以下のようなラムダ式を利用して書くだけで自動的に最適なループ処理を行ってくれます:

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

この表式ですと、どんな変数がなんの変数に影響を及ぼすのか、全画素に対してどのような処理を行うのか一目瞭然です。また、引数に渡すunary_opはコンパイラの最適化処理によって自動的にインライン展開されるため速度的なオーバーヘッドはありませんし、OpenMPによる並列化も行ってくれます。また任意の関数オブジェクトを作用させることができるので、非常に使い回しの効く関数となっています。

画素云々について言及している方はこれまで数多く見かけましたが、このような画像に特化したアルゴリズム関数について言及している方は、何故か誰もいないというのが現状です。なのでまぁ自分用にざっくり書いているコードなのでいろいろと不備はありますけど、上記のようなアルゴリズムを詰め合わせた補助ライブラリを近いうちに上げていきたいと思います。

追記(2011/09/24): 記事のとおり、補助ライブラリであるcvutilをgithubに上げました。詳細は別記事『OpenCVアルゴリズムを高速に記述できる補助ライブラリcvutilを公開します』を参照してください。

補足
『イテレータを使って標準のアルゴリズムを適用すればいいだけなのでは？』と思う方もいるかもしれません。つまりこんな感じですね:

std::transform(src.begin(), src.end(), dst.begin(), [](uchar x) -> uchar {
    return cv::saturate_cast<uchar>(x*2);
});

確かにこの式でも動くんですけど、速度は遅いです。なぜならこのイテレータの指しているデータは、メモリ上に連続していないからです。言い換えると、任意の自然数Nに対して『&(*(iter + N)) == &(*iter) + N』であることが保障されていないため、きちんと対象のデータを指し示すように補正しなければならない。その補正処理に余分な負担がかかって遅くなってしまう。ただ、もしメモリ上に連続であるならば直接ポインタをイテレータとして入れてやることで、(きちんと動作する保障はないですが)高速に動作してくれるでしょう。並列化はされていませんけれど。

また、上記のtransform_imageは並列化は行われておりますが、SIMDによる高速化は行われておりません。これをジェネリックに行うというのは現状としてなかなか厳しいのですが、Boost.SIMDの登場によってある程度の希望の光は見えてくるのではないかと思っています。

追記:

(2011/09/16)コードをコンパイルが通るように若干の修正。ucharを暗黙的に使うのではなく明示的に指定するようにした
(2011/09/17)記事を書いていて初めて知ったんですけど、ptr<T>(y)はcv::Mat_<T>だとoperator[]でオーバーロードされていたんですね。恐らくこれを用いることがOpenCV側の意図しているところだろうと判断し、上の式を書き換えました。

boost.gilを利用して透過画像を作る

2011-06-26T01:20:00.006+09:00

『対象の画像から、特定の色の箇所だけを透明にした透過画像を作る』

こういった課題に一ヶ月近くかかって出来なかったとかいう話を聞いて、流石にそれはかかり過ぎだろうと脊髄反射で言ってしまったわけですけど、それじゃあ自分だったらどうだろうとか後で思ってちょっとやってみることにしました。好き勝手言って出来なかったとか恥ずかしすぎですから、やっぱりそこらへんは実際にやってみる責任みたいなものは感じているわけです。というわけでboost.gilを使って書いてみた結果が以下。

#include <boost/gil/gil_all.hpp>
#include <boost/gil/extension/io/png_io.hpp>

using namespace boost::gil;

template<typename Pixel>
struct transparent_functor {
public:
    explicit transparent_functor(Pixel bg_pixel): bg_color_(bg_pixel) {};
    Pixel operator()(const Pixel& src) {
        if(src == bg_color_) {
            Pixel dst = src;
            get_color(dst, alpha_t()) = 0;
            return dst;
        } else {
            return src;
        }
    }
private:
    Pixel bg_color_;
};

int main() {
    const char* filename = "src.png";
    const rgba8_pixel_t background_color(0, 255, 0, 255);

    rgba8_image_t src_image;
    png_read_image(filename, src_image);
    rgba8_image_t dst_image(src_image.dimensions());

    transform_pixels(
        const_view(src_image),
        view(dst_image),
        transparent_functor<rgba8_pixel_t>(background_color));
    
    png_write_view("dst.png", view(dst_image));

    return 0;
}

結果としてはこんな感じです。上が元の画像で、下が結果の画像。正常に抜き出されていることが分かります:

余談

プログラムより文章を紡ぎ出すのに時間がかかって、改めて自分の文章力の無さにうんざりしています
というか最近暑い！自分は寝たいのに梅雨ですごい蒸し暑いから寝付けず気になってこんな不毛なことやっちゃってるわけです。今度こそ寝ます！！！

Thrustから見る、reduceを用いたアルゴリズム実装

2011-06-02T00:40:00.002+09:00

"The Great Day of His Wrath" by John Martin

はじめに

世はまさに大並列時代。火を手にした人類が飛躍的な進化を遂げたように、並列化のパラダイムは、今まで到底不可能と思われていた速さで計算を行うことができるようになりました。

ところが、どんな手法にも欠点は存在するもので、実際に実装しようとすると非常に難しい。何故かというと、並列には並列化特有の問題が存在しており、愚直に実装してしまうとCPUより早くなってしまうどころか遅くなってしまうことだってあり得るのです。これを回避するにはGPUの内部構造についてきちんと理解をした上で、実装したいアルゴリズムそれぞれの場合に特化したコーディングを行う必要がある。

しかしよくよく考えるとおかしな話です。私たちが実装したいのはあくまで手法であり、ハードウェアではありません。なぜこのような詳細について把握する必要があるのでしょう。これはきちんとした抽象化がなされていない結果生み出された、言ってみれば妖です。この妖によって余計な手間が生み出され、コーディング量は増大し、余計なバグが生成され、絶望した開発者はうつになり、自殺者は増え、世界は暗雲に包まれるー

そのような混沌の中、ある一筋の光、Thrustが現れました。ThrustはCUDAを抽象化し、C++のパラダイムを用いることで本質にのみ注力できるよう開発されたライブラリです。ユーザはただIteratorに対してアルゴリズムを適用するかの如くdevice_vector<T>と各種アルゴリズムを駆使することで、流れるようなプログラムを組むことが可能となります。このような抽象を用いることで、ユーザは詳細な内部構造について頭を悩ませる必要が無くなりました。誰が書いているのか分からないから使いたくない？ThrustはNVIDIAに所属している開発者Jared HoberockとNathan Bell氏によって書かれています。2人とも並列化のプロフェッショナルなので、普通の人が愚直に実装するよりも大抵の場合効率的でしょう。

さて、Thrustには並列化を行うための多種多様なアルゴリズムが存在しています。そして、Thrustはオープンソースです。これはつまり、このソースコードを読めば並列化のプロがどのようにして各種アルゴリズムの並列化を実装しているのかが分かるということです。

Reductions

ということで、最近はThrustのソースコードをちょくちょく読んでいるところです。まだ読み始めた最中なのですが、特に感動したのはReductionsの項でした。

数々あるSTLのアルゴリズムを普通に実装しようとすると、一つ一つのアルゴリズムすべてに対して、コアレッシングなどのGPGPU特有の問題を回避したコードを記述しなければならず、非常に面倒です。Thrustはこの問題をreduce関数によって解決しました。この関数は並列化が可能なので、これを用いて書かれた関数は自動的に並列化されます。つまり、Thrustはいわば『全から一を返す』アルゴリズムすべてをreduceを用いて記述することによって、冗長な記述を回避しているのです。

このreduceを用いた並列化の実装について簡単に説明していくことが、今回の記事の目的です。reduceを初めて聞いた方は前回の記事『並列環境におけるreduceとscanアルゴリズム』を参照してください。

ここで、本来のThrustはC++によって記述しているのですが、そのままC++で書くと非常に冗長で解り辛くなってしまうので、今回は本質のみを抽出するためHaskellを用いて説明します。あまりHaskellに慣れていないので不適切な書き方等あるかもしれませんが、そのへんはご指摘頂けると幸いです。また、今回のreduceにはfoldr関数を代わりに用いています。分かると思いますが一応。

実装

max_element, min_element, minmax_element

その名の通りリストの中の最大値、最小値、最大と最小値のタプルを返します。

*Main> max_element [1,3,2,4,3,1]
4
*Main> min_element [1,3,2,4,3,1]
1
*Main> minmax_element [1,3,2,4,3,1]
(1,4)

ここらへんはまぁ普通に分かります:

max_element :: Ord a => [a] -> a
max_element list = foldr max (head list) list

min_element :: Ord a => [a] -> a
min_element list = foldr min (head list) list

minmax_element :: Ord a => [a] -> (a, a)
minmax_element list = foldr minmax_func (first, first) $ zip list list where
    first = head list
    minmax_func (min1, max1) (min2, max2) = (min min1 min2, max max1 max2)

all_of, any_of, none_of

all_of, any_of, none_of関数は対象のリストと真偽値を判定する関数を引数に取り、

すべて当てはまっている
一つ以上当てはまっている
すべて当てはまっていない

場合にはTrueを返す関数です。pythonではall(), any(), C++ではまんまのやつがありますね。

*Main> all_of (\x -> x > 3) [1,3,2,4,3,1]
False
*Main> any_of (\x -> x > 3) [1,3,2,4,3,1]
True
*Main> none_of (\x -> x > 3) [1,3,2,4,3,1]
False

これも実装的にはらくちん:

all_of :: (a -> Bool) -> [a] -> Bool
all_of binary_pred list = foldr (&&) True $ map binary_pred list

any_of :: (a -> Bool) -> [a] -> Bool
any_of binary_pred list = foldr (||) False $ map binary_pred list

none_of :: (a -> Bool) -> [a] -> Bool
none_of binary_pred list = all_of (not . binary_pred) list

関数合成を用いてエレガントに解決。

inner_product, transform_reduce

inner_productは2つのリストに二項演算子を作用させてから、その結果からなるリストをreduceする関数です。transform_reduceはリストをいったんユニタリー関数を用いて変形させてから、reduceを実行する関数です。

*Main> inner_product (*) (+) 0 [0,1,2] [3,4,5]
14 -- = 0*3 + 1*4 + 2*5

*Main> transform_reduce (\x -> x*2) (+) 0 [1,2,3]
12 -- = 1*2 + 2*2 + 3*2

C++だといろいろ冗長に書かないといけないですけど、Haskellだったら余裕です:

inner_product :: (a -> b -> c) -> (c -> c -> c) -> c -> [a] -> [b] -> c
inner_product binary_op1 binary_op2 init list1 list2 =
    foldr binary_op2 init $ map operate $ zip list1 list2 where
        operate (x, y) = binary_op1 x y

transform_reduce :: (a -> b) -> (b -> b -> b) -> b -> [a] -> b
transform_reduce unary_op binary_op init list =
    foldr binary_op init $ map unary_op list

count_if, count

count_ifは条件に合致している要素数を返す関数で、countは指定した値と同じ要素の数を返す関数です。

*Main> count_if (\x -> x < 3) [1,2,3,4,5]
2
*Main> count 3 [1,3,2,2,3]
2

count_if :: (a -> Bool) -> [a] -> Int
count_if binary_pred list = foldr (+) 0 count_list where
    count_list = map transform_func list where
        transform_func x = if binary_pred x then 1 else 0

count :: Eq a => a -> [a] -> Int
count value list = count_if (== value) list

ThrustではTrueだったら1, Falseだったら0となるリストを作成して、それを足していく手法をとっていました。
countはcount_ifが思いつけば楽勝でしょう。

find_if, find_if_not

find_ifはリストの中で条件に合致している『一番最初の要素』を示すイテレータを返します。find_if_notは逆で、条件に合致していない『一番最初の要素』を示すイテレータを返します。今回は単純に要素のラベル番号を返すようにしました。なかったら最後のラベル+1が帰ってきます。

*Main> find_if (\x -> x > 3) [1,3,2,5,4]
3
*Main> find_if (\x -> x > 7) [1,3,2,5,4]
5
*Main> find_if_not (\x -> x > 3) [1,3,2,5,4]
0

これは微妙に入り組んでいます:

find_if :: (a -> Bool) -> [a] -> Int
find_if binary_pred list = snd result where
    result = foldr find (False, length list) tuple_list where
        find (b1, i1) (b2, i2)
            | b1 && b2  = (True, min i1 i2)
            | b1        = (b1, i1)
            | b2        = (b2, i2)
            | otherwise = (False, max i1 i2)
        tuple_list = zip (map binary_pred list) [0..]

find_if_not :: (a -> Bool) -> [a] -> Int
find_if_not binary_pred list = find_if (not . binary_pred) list

ただ、find_ifはちょっと本家と違います。インデックスが小さい項が常に左の引数となるため省略してるんだとは思いますけど、対称性を満たしていないのがちょっと気持ち悪いので、ここでは交換則を満たすように関数を変形しています。

実際の実装と課題
ここではhaskellを用いて実装しましたが、もちろん実際の実装は異なっており、zipはzip_iterator, mapはtransform_iterator, タプルはtupleを用いて実装しています。これらのイテレータの評価は遅延的に行われるので、複数のメモリアクセスが発生せず、アルゴリズムの高速化に貢献します。

これらのfancy iteratorを使用した戦略は非常に賢いとは思いますが、同時に欠点もあります。複数のイテレータを大量に組み合わせるために、通常のSTL algorithmのようにbeginとendを使用してしまうと、両方に対してfancy iteratorを組み合わせなければならないため、2倍の手間がかかってしまうのです。

これらの解決案としてはboost::rangeのようなRange構造を実装することです。これをパイプ演算子で繋げていくような形をとっていけば、明快で分かりやすく、速度も十分保証されたプログラムを記述できることでしょう。

Wiresharkを用いたHTMLの初歩的なスクレイピング

2011-05-18T01:55:00.004+09:00

前回はぽんとスクリプトを上げて「これで画像がダウンロードできますよ」とやったわけだけど、もちろん実際にGoogle Art Projectの資料が上げられているわけでもないので、きちんとHTMLやJavascriptのソースを見ながら読解していくことになる。

簡単な動作機構だったら単にHTMLのソース見るだけで良いのだけど、Google Art Projectのような複雑な機構だとそうもいかない。行数が多いと面倒だし大抵の場合難読化されてたりで一筋縄ではいかないからだ。そこで、今回はパケット解析というもうひとつの武器を利用することで内部構造を別の側面から理解することにした。思考の流れを書き綴っていくので、HTMLのスクレイピングをやってみたい方にとって何か参考になれば幸いです。

目標: 高画質な絵画をGoogle Art Projectを利用して手に入れる

まずGoogle Art Projectのサイトへアクセス。絵画ごとに固有のURLが割り当てられており、このURLを解析すれば絵画を入手できると判断。見るからにソースコードからの解析は面倒くさそうだったので、パケット解析ソフトウェアWiresharkを用いて解析を行った。

本質のみに注力するためhttpプロトコルのみキャプチャ。適当なURLにアクセスし、ズームを1回行い、パケットキャプチャを中断。解析を開始。

画像本体らしきURLにアクセスしているパケットを発見。複数のアクセスがあるってことはたぶん画像は複数に分割されていてそれを読み込んでいるんだろう。調べると/unique-id=xi-yj-zkの形でアクセスしていたので、対応するxy座標とズーム値zによって対応する画素値にアクセスすることができると推測。試しにURLにアクセスしてみると、問題なくダウンロードができた。

また、x,yの値を変更してのアクセスを行った結果も問題なくアクセスできた。ただ、どうやら範囲外のxy値にアクセスすると404が返されるらしい。当たり前といえば当たり前。アクセス禁止などのペナルティはないっぽいので、(ちょっと礼儀悪いけど)404が返されたらそれが範囲外と解釈することにする。

一旦まとめてみる。画像は分割されており、これらのアクセスはなんかよく分からないunique-idという、おそらく絵画について表しているパラメータとx,y,z座標の2つだけでダウンロードができる。ついでに、リファラー, cookie, 認証などの面倒な作業をせずとも直接のダウンロードが許されていることも分かる。そこらへんはGoogleさん寛容なんだなぁと思い、まぁAPIを提供してるくらいだからなぁ、とかいろいろ考える。

他の画像に対してのホストを確認してみると、どうやら[lh3〜lh6]までの複数のホストへとアクセスしているらしい。このホストアクセスがランダムか、それとも固有のホストを呼び出さなければならないのか確認してみたかったので、試しにホストだけを変更して、その他のパラメータは固定したままでのアクセスを試みる。結果としては問題なくダウンロードできた。たぶん負荷を分散するという単純な理由から、ランダムにJavascript側で割り振っているんだろう。ここで初めてHTMLのソースを閲覧。"cmn/js/content.js"のソースを開き(難読化はされてないようだ)"ggpht"を検索。

リストに入れているところからすると、画像へのアクセスはたぶんこのリストの中だけに限られているんじゃないかなと思う。試しにlh7へアクセス。404が正常に返された。的中。この指針で間違いないので次へ進む。

次に気になるのはunique-idだ。一体何処から仕入れているのか知らないと全く手のつけようがない。考えられるのは2つで、

アクセス毎に固有のunique-idを発行し、一定期間のみアクセスできるようにしている
HTMLかどこかに参照しやすいように書かれており、それを抜き出してアクセスする

という発想だ。自分は前者を疑っていたので、それっぽいパケットがあるかどうか監視。結果としてはまったく存在していなかった。というと後者かなと思い、対象のHTMLからunique-id(今回の場合だと"vtzd432xqzOpi_3X8JAJ_buly36QKI0n9zGeeWNgBcwsYJTsAKK5mbM9Pgg")を検索。

ビンゴ！これで駒はすべて揃った。まとめに入る。対象の絵画をダウンロードするためには、

取得したい絵画のアドレスへアクセス
<div id="microscope" data-microId="unique-id">となっている項を取得、unique-idを得る
http://lh3.ggpht.com/ - http://lh6.ggpht.com/のいづれかにアクセス。
http://lh3.ggpht.com/unique-id=xi-yj-zk のような形で画像を取得
取得した複数の画像を繋ぎ合わせる

後は非常に簡単でただこの単純なアルゴリズムに従って書き下していくだけ。その産物がこいつ。こんなに早く仕上がるってことは、たぶんGoogleさんは本格的にダウンロードを禁止したいわけじゃないんだろう。

スクレイピングとしては楽な部類で、普通だともうちょっと複雑な認証をかましてくるので自動化しにくい。Captcha入ると流石に無理だけど、おそらく現在の文字認識技術をきちんと応用すれば解けるんじゃないかと思っている。

余談

こういう解析は楽しい。どんな感じの楽しさかっていうと、与えられたパズルを解いていく感覚にちょっと近い
自分はだいたい1時間くらいで解いたけど、早い人はたぶんもっと早いんだろうなって思う。まだまだ精進が足りない

Google Art Projectを利用して高画質の絵画を手に入れる

2011-05-16T01:20:00.002+09:00

"The Starry Night" by Vincent van Gogh

2011年2月2日、GoogleさんはGoogle Art Projectというサイトをリリースしました。どういうサイトかっていうとMoMAニューヨーク近代美術館やTate Britainなど、世界中の美術館にある有名な絵画を渡り歩くことができるという趣旨のサイトで、誰もが一度は見たことがある絵画がー家に居ながらー非常にリアリスティックな解像度で見ることができる。恐ろしい時代になったもんです。

ところがこのサイト、高解像度で見られるのは非常にいいことなのですが、肝心のダウンロードができない。せっかく高解像度なんですから保存できてもいいはずなのに、だけどできない。Google Mapなんかは常に最新の情報に更新しなければならないので、ダウンロードに制限をかけているのは合理性があるんですけど、絵画とか全く更新する必要がないでしょう。著作権も切れてるのに。よく分からないです。

ということで、なければ自分でなんとかしようという精神で、1時間くらいでちゃっちゃとPythonでスクリプト組んでダウンロードするようにしてみました。即興で組んだのでいろいろ粗い面はありますけど、そこら辺は勘弁してください。

#!/usr/bin/env python
#-*- coding: utf-8 -*-

import sys
import re
import urllib2
import random
from PIL import Image

def usage():
    print "%s [URL]" % sys.argv[0]
    return 1

def get_id(url):
    response = urllib2.urlopen(url)
    if response.code != 200:
        raise TypeError
    data = response.read()
    return re.search(
        r'<div id="microscope" data-microId="(.+)">', data).group(1)

def get_filename(x, y):
    return "%i_%i.jpg" % (x, y)

def download_image(target_id, x, y, z):
    urllist = [
        "http://lh3.ggpht.com/",
        "http://lh4.ggpht.com/",
        "http://lh5.ggpht.com/",
        "http://lh6.ggpht.com/"
    ]
    base_url = random.choice(urllist)
    url = "%s%s=x%i-y%i-z%i" % (base_url, target_id, x, y, z)

    try:
        response = urllib2.urlopen(url)
    except urllib2.HTTPError:
        return False
    with file(get_filename(x, y), "wb") as fp:
        fp.write(response.read())
    return True

def download_and_get_max_x(target_id, x, z):
    result = download_image(target_id, x, 0, z)
    if result:
        return download_and_get_max_x(target_id, x+1, z)
    else:
        return x

def download_and_get_max_y(target_id, y, z):
    result = download_image(target_id, 0, y, z)
    if result:
        return download_and_get_max_y(target_id, y+1, z)
    else:
        return y

def combine_images(name, max_x, max_y):
    x_size = 512
    y_size = 512

    result = Image.new("RGB", (max_x*x_size, max_y*y_size))
    for y in range(max_y):
        for x in range(max_x):
            image = Image.open(get_filename(x, y))
            result.paste(image, (x*x_size, y*y_size))
    result.save(name)

def main():
    if len(sys.argv) == 1:
        return usage()
    depth = 3
    target_addr = sys.argv[1]
    target_id = get_id(target_addr)
    max_x = download_and_get_max_x(target_id, 0, depth)
    max_y = download_and_get_max_y(target_id, 1, depth)
    for y in range(1, max_y):
        for x in range(1, max_x):
            download_image(target_id, x, y, depth)
    combine_images("result.png", max_x, max_y)

if __name__ == "__main__":
    sys.exit(main())

使い方は簡単で、

$ python artproject.py http://www.googleartproject.com/museums/tate/mariana-248

のようにURLを指定してあげると勝手に取得して結合してくれます。こんな感じで:

"Mariana" by John Everett Millais

あぁそれにしてもお美しい……巷ではオフィーリアちゃんが有名だけど、やっぱりマリアナさんのほうが綺麗だよ。絶対。

余談

もともとの動機としては高解像度な絵画の壁紙が欲しくて、いろいろと探し回ったわけですが、どうもGoogle Art Projectしか高解像度の絵画が載っていないので仕方なくといった次第です。
利用はどうやらGoogle's Terms of Serviceに準じているようです。ざっと見した感じだとダウンロードを禁止する条項はないですし著作権も切れているので全く問題ないように思えるんですが、もし警告とかきたらチキンですから遠慮無く取り下げます。ご了承ください。
このネタを題材に時間があればHTMLのスクレイピング、パケット解析についての入門記事を書いてみたいです。
なんだかんだ言ってちゃんと実物見たほうが綺麗。ミレイの作品は以前東京に行ったときにBunkamuraで見ました。Bunkamuraはそういう催し物よくやっていて、そこらへんが東京羨ましかったり

並列環境におけるreduceとscanアルゴリズム

2011-04-19T01:28:00.001+09:00

"fireking" by cmyk1219

1. 概要
reduceとscanはGPGPUなどの並列環境下において重要な役割を果たす関数だけど、あまりこれらの関数、特にscanについて言及している記事はあまり見かけない。なので今回はreduceとscanについて調べた結果をまとめていきたいと思う。

2. reduce, scan関数の概要
2.1. reduce
scanについてはreduceが分かればより理解しやすくなるので、ひとまずはreduceから始めることにする。

reduceとは、端的に言うと、対象のリストを与えられた二項演算子で纏め上げる関数だ。とは言っても、いきなりそんなことを言われてもよく分からない。まずは例を見ていこう:

>>> import operator
>>> reduce(operator.add, [1, 2, 3, 4, 5])
15
>>> reduce(operator.sub, [1, 2, 3, 4, 5])
-13
>>> reduce(operator.mul, [1, 2, 3, 4, 5])
120

reduceにおいて重要なのは二項演算子であり、上記の場合、reduceは以下のような演算を行う。

(((1 + 2) + 3) + 4) + 5 = 15
(((1 - 2) - 3) - 4) - 5 = -13
(((1 * 2) * 3) * 4) * 5 = 120

これにより、最初の演算はsum関数、3番目の関数はprod関数と等価であることが分かる。
リストの型は何も数値に限られている訳ではない:

>>> reduce(operator.and_, [True, False, True])
False
>>> reduce(operator.or_, [True, False, True])
True

これより、上記2つの演算はall, any関数と機能的に等価であることが分かるのだけれど、all, any関数のほうが効率的なので、実際はちゃんとall, any関数を使うべきだ。

他にもreduceを用いて様々な関数、アルゴリズムを記述できるが、あまりにも多いので今回は後回しにして、scanを重点的に解説していく。

2.2. scan
scan関数はreduceと比べて文献が少ない。言及している文章もあまり見かけないが、並列処理においてscanは重要な関数だ。残念ながらpythonで実装されていないけれど、仮に実装されていたとしたら以下のような働きをする:

>>> scan(operator.add, [1, 2, 3, 4, 5])
[1, 3, 6, 10, 15]
>>> scan(operator.sub, [1, 2, 3, 4, 5])
[1, -1, -4, -8, -13]
>>> scan(operator.mul, [1, 2, 3, 4, 5])
[1, 2, 6, 24, 120]

きちんと書き下せばどのような働きをしているのか分かりやすい:

[1, 1+2, (1+2)+3, ((1+2)+3)+4, (((1+2)+3)+4)+5]

結局のところ、scanは先頭からreduceしていった結果をリストにまとめ、それを返す関数となる。scanは別名prefix sumとも呼ばれている。
ちなみに、haskellではscanl1(scanl)という名前で実装されている。

Prelude> :type scanl1
scanl1 :: (a -> a -> a) -> [a] -> [a]
Prelude> scanl1 (+) [1,2,3,4,5]
[1,3,6,10,15]

また、CUDAライブラリであるthrustにはinclusive_scan(exclusive_scan)という名前で実装されている。

#include <thrust/scan.h>
int data[6] = {1, 0, 2, 2, 1, 3};
thrust::inclusive_scan(data, data + 6, data); // in-place scan
// data is now {1, 1, 3, 5, 6, 9}

3. 並列による制限
何故並列下においてreduce, scanが重要になってくるのか、それは、両方ともメニーコアな環境において効率的に動作するアルゴリズムが複数提唱されているからだ。reduce, scanは並列と相性がいい。つまり、既存のアルゴリズムをどうにかしてreduce, scanに落とし込むことができれば、並列化の恩恵を受けれられると。

ただし、並列環境においては、演算にいくつかの制限が生じてくる。具体的には、集合Gと演算子"・"の組において、

演算子"・"は交換則"a・b = b・a"を満たしていなければならない
演算子"・"は結合則"(a・b)・c = a・(b・c)"を満たしていなければならない

後者に関してはあたり前の話で、これが成り立っていなければ並列に計算を行うことができない。前者については多少曖昧で、実際には満たさなくてもよいアルゴリズムも作ることはできるけど、コアレッシングなどの問題を考えると満たしていることが強く求められる。実際、thrustのreduceは交換則を満たしていることを前提としている。

reduceの場合はこの2つの制限で十分だけれど、scanの場合は上記2つの制限に加えてさらに2つの制限を満たしておくと都合がいい。具体的には、集合Gと演算子"・"の組において、

ただ一つの単位元"e"が存在する
演算子"・"は任意の元"a"に対して"a・a^{-1 =}a^-1・a = e"を満たす、ただ一つの逆元"a^-1"が存在する

これら4つの制限を満たしている集合Gと演算子"・"の組を『アーベル群』という。要するに、scan関数は、対象の型から成る集合と二項演算子の組がアーベル群を成していれば都合がいいということになる。

具体例を出そう。整数と和の組(Integer, (+))は明らかに上記4つの制限を満たす。よって、この組は並列にreduce, scanオペレーションを適用でき、かつ都合が良い。

1 + 2 == 2 + 1
(1 + 2) + 3 == 1 + (2 + 3)
1 + (-1) == (-1) + 1 == 0

真偽値と論理積の組(Bool, (&&))は上記2つの制限を満たすため並列にreduce, scanオペレーションを実行可能だが、後半2つの制限は満たしていないため都合は良くない。

3次元の回転行列と積の組(Matrix, (*))は交換則を満たしていないため―少なくともthrustでは―並列にreduce, scanオペレーションを実行することはできない。

3.1. 都合が良いとは？
とは言っても、一体アーベル群であればどこがどう都合がいいのか？それは、リスト中において局所的にreduceを行いたい場合に役に立つ。

リスト[a₀, ... , a_n]にscanオペレーションを実行して

を手に入れたとする。この場合、i < jにおいてb_j+b_i^-1を計算すると、

となる。つまり、i+1からjまでの要素についてのreduceを僅かO(1)で求めることができる。

これには大きな利点がある。例えば、大規模な時系列データを、その周りのデータから平均化(平滑化)したいとする。各要素ごとに計算を行うのは明らかに非効率であるので、まずscanオペレーションを(+)に対して適用し、その後で差分を計算し、正規化を行ってやればよい。結局、scan (+)のやっていることはf(x)の不定積分F(x)を求めることと本質的に等価となる。この考え方はCVでは、積分画像として様々な箇所で用いられる。

そしてさらに、この操作はアーベル群であれば何でも良い……とは言っても、アーベル群であることは結構厳しい制限となる。

例を出すと、リスト中の任意の範囲ーiからjの成分がすべてある条件式に適合しているかどうか確かめたいとする。
関数型の考え方だと、まずはmap関数を用いて(C++でならtransform_iteratorを用いて)リストの成分をBoolに変え、その後でscanを利用することを思いつくかもしれない。しかし、残念ながら(Bool, (&&))はアーベル群でないので、そのまま愚直に適用することができない。

Prelude> scanl1 (&&) $ map (<3) [1,2,3,2,1]
[True,True,False,False,False]
-- 最後の2要素については条件を満たしているのに、Falseで埋もれてしまって分からない！

この場合は、単純にアーベル群である(Integer, (+))をscanに適用してやればよい。

Prelude> scanl1 (+) $ map (\x->if x<3 then 1 else 0) [1,2,3,2,1]
[1,2,2,3,4]
-- きちんと第3要素のみ条件から外れていることが分かる

Pythonっぽく書くとたぶんこんな感じ:

>>> import operator
>>> f = lambda x: 1 if x < 3 else 0
>>> scan(operator.add, [f(x) for x in [1,2,3,2,1]])

4. より詳しく知りたい方は
scanアルゴリズムのGPGPU実装に関する詳細については、現在のところ『Parallel Prefix Sum (Scan) with CUDA(Mark Harris, NVIDIA)』が一番分かりやすい。より詳しく知りたい方はそちらをどうぞ。

というより、基本NVIDIAの出しているテキストは分かりやすい。すごいことです。

Moving

2011-03-31T19:17:00.004+09:00

本当は4/1に書く予定でしたが、そういえば4/1はエイプリルフールでしたので今日書きます。

えーと、本日、東北大学理学部物理学科を卒業し、明日から東北大学大学院情報科学研究科に入学します。所属する研究室はコンピュータビジョンなど、主に視覚的な情報について研究している研究室ですので、たぶんこれから書く記事はそれ系統の内容が多くなると思います。このブログのタイトルに少し近づきました。

元々このブログは映像を作り、その技術的な情報についてログを残せるよう立ち上げたブログでした。確かに現在は映像の「え」の字もないブログに仕上がっていますが、もちろん今でも映像は好きですし、好きなプロダクションさんの映像はいつも観に行くようにしています。そして全く役に立たなかったかと言われればそういうわけではなくて、個人的に仕事を頼まれることもありましたし、Webサイトや印刷物のデザインをする際にも役立ちましたし、ここで得た知見のお陰で様々な方とお話することもできました。

物理学科に入ったのも元々は相対論や量子力学の深遠な世界について知りたいために入りました。正直に言って完全に理解できたかというと怪しいところでありますし、理学から工学という他分野に移るわけで、一見すると今までの知識が役立たなくなるように見えます。が、実はそうでもなく、統計力学(*1)や相対論で使っている手法をコンピュータビジョンに輸入している論文があったり、案外様々な場面で役立ったりしているわけで(*2)。

全然違う領域に見えても絡んでくる領域というものは少なからずあるものです。ということで、映像と物理がお互いに絡み合っている、コンピュータビジョンの道に進むことを決めました。

これからもどうぞよろしくお願いします。

p.s. 家族は無事でした、よかったよかった

*1 ここで宣伝ですが、田崎先生の『統計力学I, II』はお薦めです。この本はぜひ理学系だけでなく工学系の方々にも読んでもらいたい本です。値段も2冊併せて7,000円くらいと安い。新品のゲームソフト1本買うくらいならこれを買いましょう
*2 もちろん線形代数や基本的な偏微分も数式や概念を学んだりする上で役立ちました。一番学んで良かったのはブラケット記法とヒルベルト空間、あれは素晴らしい

自分は無事です

2011-03-11T20:51:00.001+09:00

自分は出張で京都にいましたので、無事です。もし身内の方が見ていましたらご安心ください。
ただ暫くは京都に滞在することになると思います。早く交通機関が復旧してほしいです。

ただ家族がどうなっているのか全く分からない。怖いです、不安です。

[文句]
iPhoneではSoftbankの災害用伝言ダイヤルに『登録できない』！！！

明らかにsoftbank側の怠慢、糞だ。糞すぎる

OMakeのContributionsに載ったった

2011-02-16T02:22:00.001+09:00

The OMake build system: User Contributions

やったね！

余談
翻訳をすることは勉強にもなって非常にいいですし皆さんにもぜひやって欲しいのですが、
その分、なんというか、日本の中へと篭ってしまうことにも繋がっている気がします。

英語に臆することなく積極的にアウトプットできるよう、自分自身の心構えをもう少しきちんとしなければならんですね。

boost.GILで組む、ジェネリックな画像処理のアルゴリズム

2011-02-15T03:54:00.012+09:00

"Image dans le néant" by gelinh

boost.GILは凄い！開発者の頭の良さがビシビシと伝わってくる！
ということで、今回はGILに関しての紹介記事を書こうと思います。

概要
あなたは画像処理のエキスパート。顧客の依頼で、8bitのRGB画像を処理するアルゴリズムを記述していたとします。ところが対象となるデバイスの仕様を調べていた際に、実はRGBA画像にも対応させなければいけないことが分かりました。面倒だと思いながら書いていた矢先、さらにBGRやABGR,さらには16や24bitにも対応したアルゴリズムを記述しなければならないことが判明しました。なんということでしょう…これらの画像すべてに対してアルゴリズムを書くなんて、とてもじゃないですがやってられない。やめてくれ！って感じです。

boostに付いてくるGILを用いることで、画像に対する操作をよりジェネリックに行うことが出来ます。具体的に言うと、あなたはGILの作法に従ってアルゴリズムを書き下すだけで、メモリに展開されている画像データの色空間、色深度、レイアウト、配置順に『関わらず』、任意のデータ配置に対して適用できるアルゴリズムを作ることができるのです。最初に提示した問題が一気に解決しちゃいました。

長所と短所
長所

画像処理に対してジェネリックなアルゴリズムを記述でき、少ない記述量で全てをカバーできる。処理速度も普通に書き下した場合と同等(静的なポリモーフィズムしか使っていないので、コンパイル時に必要な計算はすべて終わっている)。
様々なアダプタが存在している(e.g. subimage_view, rotated180_view)ので、ばんばん遅延評価して一気に処理を行うことが可能。
普通の人が普通に書くよりも、GILのアルゴリズムに従って書いたほうが大抵の場合早い
ソースコードが比較的読みやすい。コード量も短くなる。
誰が作ったかわかんないもの利用したくねーという方も安心のAdobe製。画像処理の専門家が作ってるのでとても合理的。

短所

テンプレートを多用しており『非常に難解』。
全てを抽象化しているために『非常に難解』。
C++と画像処理の両方に秀でていなければ記述できないので、使っている人が殆どいない。
英語のサイトですら文献が殆どない

(上図はMemory basedで、画素のメモリ配置がInterleavedな場合のView構造)

以下のリストで、それぞれの名称がどのような役割を持っているのか説明します。

View: (x, y)座標における画素値の他、特定のX、Y方向におけるイテレータや全ピクセルにわたって捜査できるxyイテレータなど、「様々な手法によって画素にアクセスできる手段」を提供する、一種の抽象的なクラスを表しています。画像の幅、高さを表すdimensionsと、後述するLocatorを持ちます。(注: 画像のデータを実際に保持している『わけではありません』。Viewはあくまで、対象の画素値へとアクセスするための『手段』を提供しているだけに過ぎないのです)
Locator: 画像の2次元的な『位置』を表します。ランダムアクセスイテレータの2次元版のようなものと考えれば理解が早いと思います。具体的に言いますと、 *loc あるいは *(loc + point2(dx, dy)) などで現在の画素値や(dx, dy)だけ離れた位置の画素値を取得することは出来ますが、Locatorの現在地(xy座標)や画像の大きさ(width, height)を取得することはできません。
Pixel: Locatorによって返される画素値を表します。後述するChannelとLayoutを持ちます。
Channel: 画素の色深度を表します(8bit, 16bit, etc…)。各型の取りうる最大値、最小値も取得できます。
Layout: 色空間のレイアウトを表します。ちょっと分かりにくいので、下の2つを読んだほうが早いでしょう。

Color Space: 色空間を表します(e.g. RGB, CMYK, HSV, Gray, etc…)。
Channel Mapping: 実メモリ上での配置順を表します(e.g. RGB, BGR, etc…)。ただし、この情報はInterleave画像のみ用います(Planar画像だと配置順は関係ないため)。

アルゴリズムを用いて書いてみる
それでは実際にGILに付属しているアルゴリズムを用いて実際に書き下してみます。以下の2つのアルゴリズムを用いてみましょう。

boost::gil::transform_pixel(const View& src, const View& dst, F functor)
対象のピクセルを受け取って別のピクセルを返す関数を、全ピクセルに対して適用します。この関数は、画素値と出力値が1対1対応である場合に用います(e.g. 明度、コントラスト, レベル補正, トーンカーブ, etc...)。

boost::gil::transform_pixel_positions(const View& src, const View& dst, F functor)
周囲のピクセルの値を元に結果となるピクセルを返す関数を、全ピクセルに対して適用します。transform_pixelとは、受け取る引数がロケータであるという点が異なっており、周囲の画素値を元に実際の出力値を求める場合に有効です(e.g. Blur, Sovel, Laplacian, etc...)。

実際にやってみれば、どんな感じか分かるでしょう:

#include <boost/gil/gil_all.hpp>
#include <boost/gil/extension/io/jpeg_io.hpp>
#include <opencv2/imgproc/imgproc.hpp>

using namespace boost::gil;

struct change_brightness {
public:
    change_brightness(): weight_(1.0) {};
    explicit change_brightness(float weight): weight_(weight) {};

    template<typename Pixel>
    Pixel operator()(Pixel src) {
        typedef typename channel_type<Pixel>::type channel_type;
        Pixel dst;
        for(int c=0; c<num_channels<Pixel>::value; ++c) {
             dst[c] = cv::saturate_cast<channel_type>(
                     static_cast<float>(src[c] * weight_));
        }
        return dst;
    }
private:
    float weight_;
};

struct change_brightness_locator {
public:
    change_brightness_locator(): weight_(1.0) {};
    explicit change_brightness_locator(float weight): weight_(weight) {};

    template<typename Locator>
    typename Locator::value_type operator()(Locator loc) {
        typedef typename Locator::value_type pixel_type;
        typedef typename channel_type<Locator>::type channel_type;
        pixel_type src = loc(0, 0);
        pixel_type dst;
        for(int c=0; c<num_channels<Locator>::value; ++c) {
            dst[c] = cv::saturate_cast<channel_type>(
                    static_cast<float>(src[c]) * weight_);
        }
        return dst;
    }
private:
    float weight_;
};

template<typename Locator>
struct x_gradient_functor {
public:
    explicit x_gradient_functor(Locator loc)
    : left_(loc.cache_location(-1, 0)),
      right_(loc.cache_location(1, 0)) {};

    typename Locator::value_type operator()(Locator loc) {
        typedef typename Locator::value_type pixel_type;
        typedef typename channel_type<Locator>::type channel_type;
        pixel_type src_left = loc[left_];
        pixel_type src_right = loc[right_];
        pixel_type dst;
        for(int c=0; c<num_channels<Locator>::value; ++c)
            dst[c] = (src_left[c] - src_right[c])/2;
        return dst;
    }

private:
    typename Locator::cached_location_t left_;
    typename Locator::cached_location_t right_;
};

template <typename SrcView, typename DstView>
void x_gradient(const SrcView& src, const DstView& dst) {
    assert(src.dimensions() == dst.dimensions());
    transform_pixel_positions(
            subimage_view(src, 1, 0, src.width()-2, src.height()),
            subimage_view(dst, 1, 0, src.width()-2, src.height()),
            x_gradient_functor<typename SrcView::locator>(src.xy_at(0, 0)));
}

int main() {
    const char* filename = "lena.jpg";
    rgb8_image_t src_image;

    jpeg_read_image(filename, src_image);
    rgb8_image_t dst_image(src_image.dimensions());

    transform_pixels(
            const_view(src_image),
            view(dst_image),
            change_brightness(1.5));
    jpeg_write_view("result_normal.jpg", view(dst_image));

    transform_pixel_positions(
            const_view(src_image),
            view(dst_image),
            change_brightness_locator(2.0));
    jpeg_write_view("result_locator.jpg", view(dst_image));

    x_gradient(const_view(src_image), view(dst_image));
    jpeg_write_view("result_gradient.jpg", view(dst_image));

    return 0;
}

また、画像の出力例を下記に示します:

左から、元画像, change_brightness, change_brightness_locator, x_gradientの結果

change_brightness_locator()については、for_each_pixel用のアルゴリズムをlocatorを用いて適用させた場合にどうなるか確認するためだけに実装しましたので、実務的ではありません。また、値を型の持つ最大、最小値に丸めるためにOpenCVのcv::saturate_castを使用しています。

main()についてはどんなことをやってるのか大体感覚的に理解できるはずですので、とりあえず3つのファンクタに関して説明を加えていきます。

typedef typename channel_type<Pixel>::type channel_type;

channel_type<T>::typeで実際に使われているチャネルの型を取得できます。Tにはピクセル、ロケータ、ビューのどれを指定しても構いません。今回の場合ですと8bitのチャネルを指定しているので、型にはunsigned charが入ります。

num_channels<Pixel>::value

ピクセルを構成しているチャネルの数を取得できます。この場合ですと色空間がRGBなので値には3が入ります。この値はコンパイル時に決定されるので、最適化を用いることでループ展開が行われ、速度的なオーバーヘッドはありません。

typedef typename Locator::value_type pixel_type;

ロケータが返すピクセルの型を取得できます。この場合ですとrgb8_pixel_tが入ります。

pixel_type src = loc(0, 0);

()演算子によって現在のピクセル値を取得しています。他にも*演算子でイテレータっぽくアクセスすることもできます。

explicit x_gradient_functor(Locator loc)
    : left_(loc.cache_location(-1, 0)),
      right_(loc.cache_location(1, 0)) {};

ロケータの変動値をキャッシュしておくことで、高速化を図ります。Locator::cached_location_tによってロケータのキャッシュ型を取得します。キャッシュを使ったアクセスには[]演算子を用います。
仰々しく言っていますが、要はstd::ptrdiff_tみたいなものです。移動する値は同じなんだから記憶しておきましょうということです。

subimage_view(src, 1, 0, src.width()-2, src.height())

subimage_view()によって画像の領域を切り出しています。これはアダプタで、実際の評価は遅延的に行われます。

まとめ
大体こんな感じでアルゴリズムを組んでいくのがGIL流です。遅延評価とアルゴリズムを武器に、ジェネリックなアルゴリズムを作っていきましょう。

余談
C++に関しての記事を投稿するのは正直に言って『少々怖い』です。大筋では合っていると思いますが、完全に合っているという自信があまりないからです。ですが、GILに関しては本当に、本当に参考となる文献が少ない。サンプルコードも多くない。こういった中でこの記事が少しでも参考程度になっていただければ幸いです。

それと、ざっくり組んでしまったのでこうしたほうがいいよとか、間違い等あれば遠慮なく言ってもらえると助かります。

(2/16) 追記: ソースコード若干の修正、説明や画像例の追加。
(3/31) typo: e.q. -> e.g.

画像をぼかす: いろんなボケ、リアルなボケ

2010-12-18T18:25:00.010+09:00

『画像をぼかす』と一口に言っても、その方法にはいろんなものがあります。一番有名なのはPhotoshopの『ガウスぼかし』に見られるようなGaussian Blurですが、残念ながらこれはカメラに見られるようなリアルなぼかしと比べるとどうしても見劣りします。原因はいくつか考えられて、一つはGaussian Blurに使われているガウス関数の勾配が穏やかなため、現実のレンズのシチュエーションを考慮していないという点があります。

それじゃあ勾配を急にした関数ー例えばフェルミ分布関数のようなーをカーネルに採用すれば良いのかというとそういうわけではありません。有名なLena画像で試してみましょう。

左が元画像、中央がカーネルにガウス関数を用いた画像、右がフェルミ分布関数を用いた画像。
確かに中央よりは綺麗だけど、それでもまだリアルとは言えない

それじゃあ他に何が違うのでしょう？端的に言うと、黒も白も等価に扱って平均してしまうのがマズい。

画像をぼかすという行為は要は畳込み積分で、周囲のピクセル値を特定のウェイトで平均してやってることになります。で、普通にブラーしてしまうとピクセルの輝度に関わらず常に一定の割合で平均してしまいますので、なんか濁ったような画像が生成されてしまうと。

通常のブラー。カーネルにはよりリアルに見せるため正六角形を用いている

それではどのようなアプローチを用いればいいのでしょうか？レンズブラーを観察してみると、明るい箇所が強調されてぼかされていることが分かります。言い換えると、輝度の高いピクセルが優先的に混合されている、つまり予めRGBの全ピクセル値を適当な単調増加関数で写像してから畳込み積分を行って、その後対応する逆関数で戻してやればいい。要はhirax.netとやっていることは同じで、expしてからボカしてlogすると、こういうわけです。

結果の画像。比較的リアルにボカされていることが分かる

並べてみれば、違いがはっきり(クリックで拡大)

Lenaさんだけでは少し分かりづらいので、適当な背景写真に適用してみましょう。

森林画像

まず、通常の手法でぼかしてしまうとどうなるでしょう？

森林画像(通常のブラー)

お世辞にも綺麗なボケとは言えません。それじゃあ今回の手法でぼかしてみます。

森林画像(リアルブラー)

なんということでしょう…リアルなボケ画像ができました。

比較画像(クリックで拡大)

今回組んだOpenCVのソースコードを下記に示します。2.2でコンパイルしましたが、多分2.0以上だったら普通に通ると思います。

#include <algorithm>
#include <cmath>
#include <string>
#include <opencv2/imgproc/imgproc.hpp>
#include <opencv2/highgui/highgui.hpp>
#include <opencv2/core/operations.hpp>

struct exp_each_element
{
public:
    exp_each_element(float factor): factor_(factor/256.0){};
    cv::Vec3f operator()(cv::Vec3b &color) {
        return cv::Vec3f(
                std::exp(static_cast<float>(color[0])*factor_),
                std::exp(static_cast<float>(color[1])*factor_),
                std::exp(static_cast<float>(color[2])*factor_));
    }
private:
    const float factor_;
};

struct log_each_element
{
public:
    log_each_element(float factor): factor_(factor/256.0){};
    cv::Vec3b operator()(cv::Vec3f &color) {
        return cv::Vec3b(
                cv::saturate_cast<uchar>(std::log(color[0])/factor_),
                cv::saturate_cast<uchar>(std::log(color[1])/factor_),
                cv::saturate_cast<uchar>(std::log(color[2])/factor_));
    }
private:
    const float factor_;
};

cv::Mat make_fermi_kernel(float radius, float range) {
    const int size = static_cast<int>(2.0*range + radius)*2 + 5;
    const float mu = static_cast<float>(size/2);
    const float range_inversed = 1.0 / range;
    cv::Mat dst(size, size, CV_32F);
    for(int i=0; i<size; ++i) {
        for(int j=0; j<size; ++j) {
            const float u = static_cast<float>(i - mu);
            const float v = static_cast<float>(j - mu);
            const float r = std::sqrt(u*u + v*v);
            dst.at<float>(i, j) =
                    1.0/(std::exp(range_inversed*(r - radius)) + 1.0);
        }
    }
    cv::normalize(dst, dst, 1.0, 0, CV_L1);
    return dst;
}

cv::Mat make_regular_polygon_kernel(int number, float radius) {
    const int size = static_cast<int>(2.0*radius) + 1;
    const float center = static_cast<float>(size)/2.0;
    std::vector<cv::Point> dst_vector;
    for(int i=0; i<number; ++i) {
        const float theta = 2.0*M_PI*i/number;
        const cv::Point point = cv::Point(
                center + radius*std::cos(theta),
                center + radius*std::sin(theta));
        dst_vector.push_back(point);
    }
    cv::Mat dst(size, size, CV_32F);
    const cv::Scalar color(1.0, 1.0, 1.0, 1.0);
    cv::fillConvexPoly(dst, &dst_vector[0], number, color);
    cv::normalize(dst, dst, 1.0, 0, CV_L1);
    return dst;
}

int main(int argc, char** argv) {
    const std::string source = argc >= 2 ? argv[1] : "lena_std.tif";
    const std::string destination = argc >= 3 ? argv[2] : "out.jpg";
    const float factor = 6.0;

    cv::Mat image = cv::imread(source);
    CV_Assert(image.type() == CV_8UC3);
    cv::Mat kernel = make_regular_polygon_kernel(6, 8.0);

    cv::Mat exp_image(image.size(), CV_32FC3);
    std::transform(image.begin<cv::Vec3b>(), image.end<cv::Vec3b>(),
            exp_image.begin<cv::Vec3f>(), exp_each_element(factor));

    cv::filter2D(exp_image, exp_image, -1, kernel);

    std::transform(exp_image.begin<cv::Vec3f>(), exp_image.end<cv::Vec3f>(),
            image.begin<cv::Vec3b>(), log_each_element(factor));

    cv::imwrite(destination, image);
    return 0;
}

余談

実質80行くらいで完成したので嬉しいです。可読性も多分結構読みやすい部類だと思う
C++のコードが何やってるか分からない方は、まず関数オブジェクトとSTLについて調べてみましょう。あれはいいもんです

追記(11/03/23)

どんなデバイスであるか、たとえば、それが入力デバイスであるのか、出力でバイスであるのか、はたまた、計算処理デバイスであるのか次第で「輝度値の単位系」をどのように保持（処理）するべきかは違うことでしょう。
大切なことは、取り扱う「値」がどんな単位であるのかを見失わないことなのか、と思います。

hirax.net::「コンピュータ画像処理」と「輝度値の単位系」

自分がこのような『任意の単調増加関数で良い』という具合で筆を進めたのには2つ理由があります。

まず1つ目が、『全てのデジタルカメラにおいて、ピクセルは輝度の対数の次元を取っているのだろうか？』という単純な疑問です。

きちんと物理的/工学的な意味合いを持つ画像処理にするためには、格納されているピクセルデータの次元がきちんと揃えられている必要があります。自分自身も普通に考えて輝度の対数が格納されているんだろうということで、任意ではなく指数関数に制限された関数系を使いましょうと運びたかったのですが、そうするための証拠が不足していた。

もしかしたら自分の予想もつかない方法で処理されている可能性があるわけで、そういった証拠が十分にないことから、単調増加関数と少し曖昧とした意味合いを含ませたのです。

もう一つの理由は『今回の画像処理は物理的に正確である必要がない』という理由です。これが大きかった。

そもそも今回の画像処理はどういう目的があるのでしょうか？ぼけて普通よりもリアルに見えて面白い、それで終わりです。物理的に計測をするわけではなく、ざっくり『それっぽく』見えればそれでいいという類のものです。そのような目的の下では、下手に物理的正確さを求めて処理時間を長くするよりは、計算速度の早いフェイクのほうがより適しています。

だから何も指数関数に制限する理由はあまりないわけで、より早くてよりそれっぽく見えれば、そちらのほうが優れていると、そう言えるわけです(ただ、今回の計算量からすると最も大きなオーダーとなっているのは畳み込み積分の箇所で、単調増加関数自体のコストはそれに比べれば微々たるものでしょう)。

OMakeマニュアル日本語訳 PDFバージョンをリリースしました

2010-10-27T02:03:00.002+09:00

"Aint they cute?" by Nicki Varkevisser

題名の通りで、以前から取り組んでおりました『OMakeマニュアル日本語訳』の全ページを一つに纏めたPDFバージョンをリリースしました。

http://omake-japanese.sourceforge.jp/OMake.pdf

よりダウンロードできます。電子書籍でまとめてご覧になりたい方や、職場などの環境で使いたいけれどネットが使えないのでごりごり印刷したい…といった場合に便利です。

元々は渋川さんの『LaTeX経由でのPDF作成』という記事を拝見しまして、試験的にPDF版を作成してみたところ、思いのほか綺麗に製本されていることからきちんとリリースしてみようと思い立った次第です。殆ど素の状態から手を加えていないので暫定的ですが、後はきちんとスタイルを組んでいって、もうちょっと読みやすいレイアウトにしていけばいいかなって感じです。

思えば個人的な暇潰しから始まったこのプロジェクトも現在はPDFにして約200ページにもなり、きちんと一つの事柄に集中して取り組めたことは、学部時代のいい経験になったと考えています。まあ、これからも暇があればちょくちょくレイアウト、文章含めて修正していきますので、どうぞ今後ともよろしくお願いします。

確率微分方程式を解くー幾何ブラウン運動、バシチェックモデルの場合

2010-10-04T01:58:00.004+09:00

"New York Fashion Week Fall 2007: Doo Ri" by Art Comments

確率微分方程式という分野がある。どういうものかっていうと、要は確率過程を微分方程式によって表そうという試みであり、微分方程式の中に微小の確率変数が含まれている。具体的にはこんな感じ。

このdWはブラウン運動を表していて、時間がたつにつれてうねうねと確率的に動く。期待値と標準偏差はそれぞれ0と√tで、時間がたつにつれて動く場所もどんどん広がっていく。この運動はランダムウオークとか酔歩とも呼ばれていて、酔っ払ったおっさんがうねうねとあちこちに歩いているように見えることからこういう名前がついたそうな。上のような式で表すことができる確率過程を一般には『伊藤過程』と言って、他にもストラトノビッチ確率解析とかいうのもあるらしいけど、詳しくやっていないので分かりません。

これの何が便利かっていうと確率過程を通常の微分方程式のように扱えるところで、確率変数をまるで微分方程式を解くかのように求めることができる。とは言っても通常の微分方程式にはない規則がいくつかあって、初めはすごく奇妙に思えるかもしれない。

具体的には、以下のように確率変数x(t)とy(t)を定めた場合、

このy(t)もまた伊藤過程に従い、確率微分方程式dyは以下のルールに従う。

普通だったら消えるはずの2乗の項が残ってしまう。それじゃあ(dx)^2はどのようにして計算するのかというと、以下の公式に従う。

なんと、(dW)^2がdtに変化してしまう！これを『伊藤の公式』という。なんでそうなるのか？それを知りたい方は専門書を漁ればいくらでも載っているので、そこらへんを参照すればいいと思う。ここが通常の微分方程式と違うところで、逆を言えばそこさえ守ればきちんと積分によって求めることができる。専門書だとこの箇所がすごく分かり辛い表現で説明されているけど、詳細をはしょればこんな感じなはず。

それじゃあ、まず手始めに幾何ブラウン運動と呼ばれるモデル

を求めてみる。これは株価の変動を表す最も簡単なモデルなんだけど、何故これが株価の変動を表しているのかっていうのは、とりあえず解いてみれば判明するはず。

さて、普通の考えだったら『変数分離法を用いてxを分離させて、そのまま積分すればよくね？』って発想になるんだけど、何度も言っているようにこれは確率微分方程式なので、まずは伊藤の補題を用いて変数変換を行い、xを分離させてみることにする。y = ln(x)とおくと、伊藤の補題より

となる。(-1/2 σ^2)の項が出てくるところがミソで、普通に変数変換してしまうと出現しない。これで定数部分だけに分離できたので、そのまま愚直に積分を行うと、

となる。これは対数正規分布の形となっているため、株価の変動を表す簡単なモデルとなっていることは容易に理解できると思う。

それでは、次にVasicekモデルと呼ばれるモデル

について求めてみる。これは主に国債金利の利子率の時間的変動を表すモデルで、いったいこいつの何が重要なんじゃいっていうと、金融商品のリスクを推定する際に使われる。金融商品の価格は利子率によって変動するので、利子率の変動を推定することは金融商品のリスクを推定することに繋がる(*1)。つまり、今の自分のポートフォリオがどれくらいのリスクを持っているのか判断できる。そもそもこの記事を書くきっかけになったのはこいつが原因で、wikipediaさんにはただ『積分するとこれになりますよ』って書いてあるだけで、どうやって積分するのかについては書いていなかった。しょうがないので自分で手を動かして、ようやっと理解できたので公開してみようと。

さて、まずはu=b-rとおいて変数変換すると、

となる。ただ普通にChain Ruleやってるように見えるけど、たまたま結果が同じだっただけで、どんな変数変換でもきちんと伊藤の補題を使わなきゃいけない。そんで、物理やってる人はピンとくるかもしれないけど、これはランジュバン方程式の形となっているので、積分には同様の手法を利用できる。まぁ自分はピンとこなかったんですけど。それはさておき、v=e^(at)uとおいてさらに変数変換してやると、

となり、いい感じに余計な項が消えてブラウン運動だけになった。よってこれを積分してやると

あとはごちゃごちゃ整理してやると

これでwikipediaさんにあるのとおんなじ形になった。めでたしめでたし。

なんかブログに書くと導出するのは楽ちんなように思えるけど、結構悩んだし大変だった。どんなもんでも答えを見るとすぐ分かるし思いつきそうなんだけど、実際になにもない地点から思いつくかどうかっていうのはまた別な話で、それを考えると新しい学問を生み出すっていうのは凄まじいほどの労力が必要だ。実際、未だにどういう発想で量子力学を思いついたのか本気で理解できない。あんなのヒントがいくらあっても無理だ。そういうことを考えた一日でありましたとさ。

*1: もちろんこの説明は適当でないが、詳しく説明してもあんまり意味がないのでこう書くことにする

追記(2010/10/04)
1/2の項が抜けていたので修正。

マルコフ連鎖モンテカルロ法を実装してみよう

2010-09-28T18:31:00.002+09:00

"Chain Bridge, Budapest" by szeke

約2ヶ月ぶりということで、ご無沙汰しております。書きたいネタというのは結構あって書こう書こうとは思っているんですが、なにより書くとなるとあれもこれも伝えなきゃいけないみたいな思いになって、結局分量の多さから諦めてしまうというのが結構続いています。もう少し気を張らずに更新していきたいものです。

さて、最近の自分はマルコフ連鎖モンテカルロ法(Markov Chain Monte Carlo, MCMC)についておりました。何にも知らない方ににとってはよく分からないかもしれませんが、この手法はマルコフ連鎖が持つ簡明さとモンテカルロ法が持つ実用性が合わさった凄まじい手法なんです。そしてなによりとてもエレガント。

とりあえず知らない人のためにてきとーに解説しますと、与えられた適当な関数から確率変数をサンプリングするための公式です。べつにこれじゃなくても確率変数のサンプリングにはいろんなアルゴリズムがあるわけですが、これを使うと関数の計算コストがやけに高い場合だとか、ピーク部分を集中的にサンプリングしてくれたりとかで結構便利で、金融や工学、物理の分野でも非常に役立っております。計算物理をやる人にとっては必須に近い手法といっても過言ではありません。

そんなMCMCですが、実際に提案されている手法には様々なものがあります。不変分布から遷移確率を逆算するため、自由度が増えるから遷移確率は一意じゃなくなってしまうんです。今回は最も有名なMetropolis-Hastingsアルゴリズム(*1)を用いて、任意の関数に従って乱数のサンプリングを行うコードを練習がてらpythonで書いてみました。

#!/usr/bin/env python
#-*- coding:utf-8 -*-

import sys
import math
import random

def main():
   f = lambda x: math.exp(-x*x / 10.0)
   for x, fx in metropolis(f):
       print x, fx

def metropolis(func, start=0.0, delta=2.0, burn_in=100, seed=1):
   metro_iter = _metropolis(func, start, delta, seed)
   for i, value in enumerate(metro_iter):
       if i > burn_in:
           yield value
           break
   for value in metro_iter:
       yield value

def _metropolis(func, start, delta, seed):
   random.seed(seed)
   initial_x = start
   initial_fx = func(initial_x)
   proposed_x, proposed_fx = 0.0, 0.0
   while True:
       proposed_x = initial_x + random.uniform(-delta, delta)
       proposed_fx = func(proposed_x)
       ratio = proposed_fx / initial_fx
       if ratio > 1.0 or ratio > random.random():
           initial_x, initial_fx = proposed_x, proposed_fx
           yield proposed_x, proposed_fx
       else:
           yield initial_x, initial_fx

if __name__ == '__main__':
   sys.exit(main())

うーん・・・きれいだ(*2)・・・

こいつのすごいところは、アルゴリズム自体は非常にシンプルで実装しやすいというところにあります。行数も上を見ると分かるのですがとっても少ないし、コア部分はたった25行くらいで構成されています。

なのに使っている学問は結構複雑で、定常分布から遷移行列を導出するために詳細釣り合いの条件からMetropolis的制約を用いて導出、サンプリングを行っております。このアルゴリズム自体は結構知ってる人が多いかもしれませんが、実際に理論方面で理解している人はどうなんでしょう。

んで、正規分布のガウス関数からサンプリングを行った結果がこれ。

赤線が元の確率分布。きちんとサンプリングが行われているようです。今回は正規分布で計算を行いましたが、非負の関数であればなんでも構いません。

あとはこれを研究テーマへと落とし込むことができれば、とりあえず一区切りつけるかなって感じです。

*1: 提案確率関数が対称であるという制約を持つのがMetropolis。それを非対称の場合にも拡張したのがHastings。個人的にはMetropolis-HastingよりMetropolisのほうが綺麗で好きです。というか提案確率密度関数が非対称じゃないといけないようなシチュエーションがよく分からない
*2: yieldとジェネレータって便利ですよね。計算速度も効率的だし、慣れてくるとすべて遅延評価で行わせたいくらい

ドルコスト平均法は賢い投資方法と言えるのか？補足

2010-07-23T01:49:00.009+09:00

"Pyramid of Capitalism" by Warren Noronha

以前ブログの記事を読んでいただいた友人から感想を伺うと、「なんかよく分からない」といった感想が帰ってきました。基本的に説明が不足しているなぁと思った箇所は実際いくつもあって、金融の確率論を少しでもやると大体感覚的に理解できるのですが、普通の人はそんなのやらないのが当たり前ですので、基本的な概念について少し説明したいと思います。この分野に興味ない人でもちょこちょこと面白そうな話を挟んでいくので、暇つぶしに見てやってください。

1. 確率変数、期待値、分散

確率変数というのは試行によって出てくる値が確率的に異なってくる変数のことです。これは抽象的なので、もうすこし具体的にいきましょう。例えば、以下のコイントスゲームを考えます。最初あなたは1万円を所持していたとして、表であれば所持金が2倍となり、一方で裏であれば所持金が0.5倍、つまり半分になります。このような試行を3回行ったとして、最終的な資金額をXと表すと、このゲームは以下の式で表せます。

ここでのX₁はそれぞれの試行の確率変数に相当し、50%の確率で0.5, 50%の確率で2.0の値をとります。
で、これだけではただ定式化しただけで、何にも面白いことは分かりません。ですがこいつの期待値と分散を求めてやることで、このゲームは得なのか損なのか、またどれくらいのリスクがあるゲームなのか調べてやることができるのです。

例えば、Xの期待値をE[X]としてやりましょう。ここで、期待値の性質として以下を確認しておきましょう。まず、XとYが確率変数とすると、その和の期待値は単純にそれぞれの期待値の和で表せます。

また、XとYが独立(それぞれの試行がもう一方の試行に影響を及ぼさない)であるとすると、その積も同じような性質を持ちます。

この性質を用いると、Xの期待値は以下のように表せます。

つまり、このゲームを行うことによってあなたは9500円くらいの利益を得ることが期待できます。もうこれはやるっきゃないと。

このゲームが得であることは理解できましたが、ただこのゲームにどれくらいリスクがあるのかは不透明です。期待値が10000円プラスになるからといって、40%の確率で全財産を失ってしまうみたいなゲームですと、普通は行いたくないですよね。このリスクについて求めるためには分散を計算してやればいいのですが、これはそれなりに面倒です(*1)。前回ごっちゃごっちゃやってたのはこの面倒さによるものが大半ですので、今回は説明しません。とりあえず分散がリスクに対応しているということが大事です。

世の中の投資行動には常にリスクとリターンが付きまといます。一般的にリスクとリターンは紙一重と言われており、リスクの高い行動はよりリターンも大きいというのが感覚的な理解です。ただ、このリスクとリターンを確率変数を用いて計算することによって、どの行動が一番リスクとリターンのバランスが取れている行動であるのか、判断することができます。

で、どの戦略を用いれば合理的に資産を増やすことができるのかという問題ですが、実はこの戦略は既に知られており、ある戦略を用いればリスクをリターンに変化させることで指数関数的に資産を増やすことができます。その戦略について興味がある方は適当に調べると出てくると思います。

2. なぜ乗法モデルを仮定しているのか？

まず、株価の変動は確率的であるから確率変数として表せるという話は納得できると思います。では、株価などの資産モデルには様々なものがあるのに、なぜそれぞれの確率変数が独立と仮定した乗法モデル

をわざわざ採用するのでしょうか？別に、例えば

のような加法モデルでも構わないわけです。何故なんでしょう？

まず独立と仮定した理由については、単に計算が楽だからです。正直独立と仮定しないと、ここまで綺麗な形がでるのか、そもそも解析的に解けるのかどうかすら分かりません。ただこれを仮定したことによって、ごちゃごちゃ計算できて結果的にそれなりに綺麗な数式へと展開することができたという次第です。

乗法モデルを仮定した理由については、金融業界の指数関数的な性格によるものです。例えば、ローンなんかは年利15%ーフリーローンですとこんくらいでしょうかーとありますが、これは1年間につきローンが15%増加していくので、ふるまいとしては指数的になります。同様に、国債についても年利で計算するので指数的な計算となります。金利が低いといってほいほい借金をすると危険だという理由もこれによるものです。逆に、金持ちがより金持ちになる理由の一つでもあります。

これは個人的な意見ですが、金融の世界で乗法モデルが成り立っている理由は、人間が金に対して貧欲であるからだと思っています。はっきりと言いますが、世の中の結構な割合の人間は指数的なふるまいを理解できません。ですからマイホームを数十年ローンでぽんと買ってしまったり、リボ払いで決済をしたりしてしまいます。

ちなみに、これらの金額はすべて年金公式と呼ばれる公式で算出できます。1期間後に支払いを始め、n期間にわたって金額Aを支払う場合、金利をrとおくと借入額Pは

もしくは

と表せます。これはすごく役立つ公式なので、店員の言葉に騙されないでおくこともできるかもしれません。ついでに、ローンを借りて全額返した場合に損をする差額はnA-Pで計算できて、

となります。こういう情報ってあんまり表沙汰にならないんですけれど。

追記(10/08/23)
上の式は正直汚いからどんな振る舞いをするのか感覚的に理解し辛いですけど、利率が少ない、つまりrが微小のときは (1+r)^n ～ 1 + nr で近似できますから損失額は

となります。ですから利率が低いときは、分割回数をちょっと多くしても線形的な負担にかならない。確かに損失額が増えることは確かですけど、一括払いで生活に支障がでるよりは、3, 4回くらいローンで払うっていう選択肢もありなのかもしれません(*2)。

*1: 1次元のising模型の解析解を求めたり、極座標のラプラシアン導出するのよりは面倒じゃない
*2: ただそこまでして本当に欲しいものかそれっていう疑問はありますし、30年ローンとかで家を購入するような場合には全く通用しません。そもそもこのご時勢に長期のローンを組むなんて狂気の沙汰です。狂っている。

BLESS by Overture

2010-06-05T14:32:00.009+09:00

BLESS from Overture on Vimeo.

KiraKiraというアイスランドのミュージシャンの音楽PV。作ったのはOvertureというアメリカのユニット。

これの何が凄いのかっていうと、たった2人の夫婦によって作られたっていうところです。やった人は分かると思うんですが、アニメーション制作というのは本当に大変です。それを夫婦の共同作業でやっちゃってるところが仲良さそうで素敵。しかもその作品もとても個性溢れていて、気持ちがいいなぁと思う作品でした。

ブログも良かったらごらんになってください。奥さんが日本人なので片方は日本語で書かれています。

ドルコスト平均法ははたして賢い投資方法と言えるのか？ part2

2010-05-31T10:14:00.010+09:00

"Soviet poster of the transport sector during the Five Year Plan. 1930" by National Library of Scotland

3. 計算

前回の記事で株価が乗法モデル

に従った場合の一般投資戦略とドルコスト戦略の資産はそれぞれ

で与えられることが分かりました。

それではこれよりリスクの定量的評価を行います。私たちが知りたいのは、それぞれの資産のリスクとリターンの値です。この場合リターンは確率変数の期待値、リスクは確率変数の分散(より正確に言うと標準偏差)にそれぞれ対応していることが分かりますので、 X_N の期待値と分散をそれぞれ求めてやればよさそうです。

3.1 期待値の計算

まず簡単な期待値から調べていきましょう。 R_i と R_j は i≠j の場合独立であり、さらに R_i の期待値を μ 、分散を σ² と仮定しているため、 X_N⁽¹⁾ の期待値は

X_N⁽²⁾ の期待値は

となります。

3.2 分散の計算

次に一番知りたい分散についてですが、正直これは期待値に比べて少々面倒くさい。これは分散が期待値のような線形的な関係を持たない上に、ドルコスト平均法の R_i の項が折りたたまっているために生じているのですが、そんなことを憂いていても仕方がないのでさっさと計算していきます。

まず比較的簡単な X_N⁽¹⁾ の分散ですが、

ここで右式の V にのみ注目すると、分散の定義より

関係式

に注意すると、上式は

となりました。よって X_N⁽¹⁾ は

と表現できます。これは比較的何も悩まずに展開できました。

次にドルコスト戦略 X_N⁽²⁾ について同じように求めていきます。

2つの項が現れました。それぞれの項は結構面倒ですので、分けて考えることにします。まず第二項についてですが、

次に第一項について計算します。まず

ここで

についてそれぞれの R_k で乗算し、独立性から期待値を分割して求めるわけですが、きちんと注目してやると

R_k となっている変数は全部で |i-j| 個
R_k² となっている変数は全部で min(N-i+1, N-j+1) 個

となっているため、上式は

と展開できます。

ここで対称性について考えます。当たり前ですが i と j を反転させても上式の値は不変であるため、反転対称性から上式は以下のように展開できます。

以上から、 X_N⁽²⁾ の分散の値は

となります。ここで、数式を見やすくするために

と定義しておくことにし、便宜的にこの変数を『修正バリアンス』とでも名づけておきます。この命名が正しいのかどうか、そもそも名前がつけられているのかについてはよく分かりません。

この分散、和の項が3つ並んできたり少々面倒くさい形になっていますが、よくよく注目してみると σ に関して単調増加の形をしており、 σ=0 を代入すると結果に0が帰ってきたりとちゃんとしたリスク関数になっているのが確認できます。ただこの関数が一般投資戦略に比べてリスクを分散する形になっているかどうかはちょっと判別しにくいですね。

まぁそんなこんなで X_N⁽¹⁾ と X_N⁽²⁾ それぞれにおける期待値と分散を求めることができました。次回はこの関数を利用して定量的な調査を行っていきます。

Note. σに関して単調増加であることと0であることの説明

まず、式の形よりV_N⁽²⁾はσに関して単調増加関数であることは自明。
次にσ=0を代入すると修正バリアンスはμ²になるので、次数を調整した後で合計のラベリングを変更してやると、

となる。

ドルコスト平均法は賢い投資方法と言えるのか？ part1

2010-05-29T03:19:00.016+09:00

"Soviet poster from the Five year plan showing economic investment. 1930" by National Library of Scotland

1. 概要

ドルコスト平均法という投資法があります。一応簡単に説明しておくと、ある区切られた期間ごとに同じ金額だけ対象の商品を購入しておくことで、安定した収益が得られるという戦略です。信託投資やFXをやっている方にとっては馴染み深い戦略かもしれませんし、ファイナンシャルプランナーからお勧めされた方も多いかもしれません。

一定期間同じ金額だけ同じ商品を買うだけというこの戦略、よく初心者にお勧めの投資法として雑誌などでは紹介されていますが、はたして本当にお勧めできる投資法なのでしょうか？またリスクを軽減できるとありますが、一体どれくらいのリスク軽減が望めるのでしょうか？それに対して、リターンはどれくらい減少するのでしょうか？以上を明らかにしていくのが今回の記事の目的です。

定量的に評価するために少々数式が登場してきますが、一応大学の経済学を専攻していたり、統計学について初歩的な知識(期待値や分散など)を持っている方ならばすいすいと読み進めることができるくらいのレベルを目指しています。
金融工学に関しての知識はまったく必要ありません。必要な箇所はすべて簡単に解説していきます。
今回の記事は式の導出までを行います。結果だけ知りたい方は読み飛ばしてもらっても構いません。

2. 2つの戦略

ここではドルコスト戦略についてより詳しく説明した後に、数式を用いて定式化していきたいと思います。今回は比較のため、ドルコスト戦略の他にもう一つ『一般投資』を取り上げていきます。

まず一般投資戦略とは特になんにも考慮することなく、自分の手持ち資産のうち一部分を対象の金融商品に投資して、その後は満期まで放置しておくという戦略です。

一方で、ドルコスト戦略はまず自分の手持ち資産のすべてをN分割し、一定期間毎に対象の金融商品を購入するという行為をN回行って、満期まで保持しておくという戦略です。

この2つの戦略を定式化する前に、まず以下の仮定を行っていきましょう。

時間を離散化してN期間として考え、現在いる時点を0時点、満期をN時点とします。
株価をS_nで表し、S_nは以下の『乗法モデル』に従うとします。

ここで、R_iは確率変数とします。つまり、S_iも確率変数となります。ただし、S₀だけは現時点で株価が判明しているのでただの実数となります。
R_i自体の確率分布は全く仮定していないことに注意してください。つまり、R_iは正規分布に従うとは限りません。
計算を簡単にするため、異なる確率変数R_iとR_jは完全に独立であり、また期待値と分散は同一の値μとσ²を付与するものと仮定します。
もちろん、これは理想化された仮定で、現実的ではありません。実際、R_iが過去の値に影響することは数々の測定から明らかになっています。が、そんなこといちいち考慮に入れるととっても大変ですし、そもそも今回の目的とは趣向が異なります。そこまで正確にやるんでしたら連続時間と仮定して確率微分方程式を解くべきですしね。こういうことするのは専門家だけで結構ですので今回はやりません。

以上の仮定を用いて、実際に2つの戦略を定式化します。まず一般投資戦略の場合、X_nをn時点での手持ち資産と定義すると、以下のように表せます。

ここで、λは手持ち資産のうち投資にかけるウェイトとして定義され、0<λ<1の間をとります。

次に、ドルコスト戦略の場合は、少々複雑ですが以下のように表せます。

ただし、後の計算を楽にするため、

とおきました。

ここで注目しておきたいのが、上式の形です。R_Nの項が最も多く出現し、満期から現時点まで下がっていくにつれて出現回数が減少していくのが分かるでしょう。このことから何が言えるのかというと、『ドルコスト平均戦略は満期近くの株価の変動に最も大きな影響を受ける』ということです。つまり、初期時点で株価が上昇しても、満期近くに株価が大きく下落した場合、ドルコスト戦略を用いると資産がマイナスになってしまうのです。

もちろん、λ=1(全力投資)の場合の投資戦略よりはリスクは分散されていることは式の違いにより明らかです。それは確かなので、ドルコスト戦略はリスクを軽減する戦略であることは間違いありません。ですが、これはλの値を調整することで同じ効果を与えそうです。

ここから分かる結論としては、『ドルコスト戦略は満期近くに株価が大きく上昇すると予想した場合に効力を発揮する戦略である』ということです。ですが通常の場合ですと、満期近くの株価の変動なんて正直予想できるはずがありません(明日の株価ですら予想は難しいというのに)。ドルコスト戦略が初心者にお勧めできるという論調は少々疑問であります。

そうなると結局のところ、フィナンシャルプランナーがドルコスト平均法を進める主な理由は『商品を定期購入してくれる客を作り出す』ためである気がしてなりません。だって長期投資という名目の元で、給料の一定金額を毎月自分たちに差し出してくれて、しかもその回数分だけ手数料をいただける。これほど美味しい客はありません…

…いやいや、そう考えるのはまだ早い。もしかしたらきちんとリスクを定量的に評価した上で勧めているのかもしれません。となると次に気になってくるのが『一体どれくらいリスクが減少するのか、その上でリターンはどれくらい減少するのか』という点です。「どうしてもドルコスト平均法じゃないといけないんですか？最初っから投資額を抑えることもできるんじゃないですか？」という声がどこからか聞こえてきそうです。

疲れたので今日はここまで。次回から実際のリスクとリターンの評価を行っていきます。

2006年に発表されたイグノーベル賞の文学賞が面白い

2010-05-11T20:28:00.009+09:00

"860 Hispanic literature; 869 Portuguese" by Helder da Rocha

2006年に発表されたイグノーベル賞の文学賞が面白いです。どういう論文かっていうと『無駄に難しい単語を多用したとき、そいつは知的に見えるかどうか』っていう論文で、まずタイトルからして "Consequences of Erudite Vernacular Utilized Irrespective of Necessity: Problems with Using Long Words Needlessly." 。完全にふざけてます。Summaryも秀逸で、これ言いたいだけにこの論文書いただろって思うくらい。面白かったので思わずノリでSummary部分を翻訳してしまいました:

必要性と無関係に使用される博学な専門的用語が導く結末：不必要に長い単語を使うことにおける問題
DANIEL M. OPPENHEIMER

執筆者にとって最も奨励されている書き方は、過度に複雑な語彙の使用を避けることだ。しかしながら、大多数の学生は知的な印象を与えようとして、わざと複雑な語彙を使おうとする傾向にある。この論文ではその戦略が効果的であるかどうか、その程度について調査する。実験1-3では語彙の複雑性について扱い、語彙の複雑性と知的に見える効果の間柄には負の関係性があることを調べる。この関係性は元の小論文の質や、小論文の内容がどうであるかや、参加者が最初に期待していた度合いについては関わらない。複雑性の負の影響については、文章の読解速度(fluency)を処理することによって調べた。実験4では直接読解速度について扱い、読解が難しくなる文章はより著者が知的でないと判断されることが判明した。実験5では読解速度を遅くする要因について調査した。以上より、読解が難しくなる文章は明らかに正解であると即座に判断できなくなるため、人々は読解する手がかりの少なさから、文章をあまり信頼しなくなる。事実、文章の本質に関わらない語彙というのは過度に文章を修飾してしまい、結果としてかえって反対の方向へ影響してしまう。これらの影響や、実際の応用についても議論した。Copyright 2005 John Wiley & Sons, Ltd.

適当な翻訳で勘弁！この教授は学生の無駄に難しい言い回しを多用したレポートに飽き飽きして書いたんです。間違いない。

でも、実際こういうのってよくあるよなぁとか思います。『齟齬』『蓋然性』『慚愧』とか使われても、あんま頭良さそうに見えません。最近の例でいいますと対談の際に用いられた『写像』でしょうか。これも別にここで使わなくても、十分に説明できましたし…

やっぱりイグノーベル賞は全体的に面白いです。こういうセンス日本にはないよなぁ。羨ましい限り。

Geant4をCernのリポジトリ経由でインストールする

2010-05-06T20:32:00.014+09:00

"ALICE - A Large Ion Collider Experiment at Cern" by dominikf

Geant4のインストールはソースコードをダウンロードし、ビルドを行うという手法が旧来の日本のページではよく書かれていますが、Debian系列(Ubuntu)のディストリビューションを使っているのであれば、いちいちそんな手間のかかることをしなくてもリポジトリ経由でインストールすることができます。こんな大容量のライブラリですと、ライブラリの依存関係を解決するだけでも大分かかっちゃいますしね。

で、aptを利用してGeant4をインストールするためには、現時点では以下の手順を踏んで行うのが一番楽ちんです。

Cernのリポジトリを新規に追加します。

deb http://lcg-heppkg.web.cern.ch/lcg-heppkg/debian stable hep
deb-src http://lcg-heppkg.web.cern.ch/lcg-heppkg/debian stable hep

の2つを、Synapticかコマンドライン経由で追加してあげればいいでしょう。
apt-get updateでデータベースを更新します。もちろんそのときこの2つのリポジトリが信用できるのかどうかわかりませんので、aptさんは

W: GPG error: http://lcg-heppkg.web.cern.ch stable Release: 公開鍵を利用できないため、以下の署名は検証できませんでした: NO_PUBKEY F464E3D84061544D

のようなエラーを吐きます。そこで、これらの警告がでないようにしてやります。
コマンドラインを立ち上げて、
sudo apt-get install cern-archive-keyring と入力。警告: 以下のパッケージは認証されていません! cern-archive-keyring 検証なしにこれらのパッケージをインストールしますか [y/N]? と当たり前ですが聞かれるので、yを入力します。信用されていないリポジトリから鍵をインストールして警告を出さないようにするってあたりが、俺を信頼しろ感が出ててとっても素敵。
apt-get update でエラーが出ないことを確認してから、sudo apt-get install geant4でGeant4をインストールして、終了です。

自分の環境(Ubuntu 10.04)ではこの手法で特に問題なくインストールできました。多分Debianでも普通にインストールできると思います。

さて、次に環境変数の設定を行っていきます。ホームディレクトリの".profile"ファイルの末尾に以下を記述します。

export G4WORKDIR=${HOME}/work/g4

これを行うことで、Geant4のワーキングディレクトリを${HOME}/work/g4に設定します。${HOME}以下は各自適当なやつに変えてみてください。

また、Geant4の動作に必要なシェルの変数を書き出すためには、対象のシェル内で

source /usr/share/geant4/env.sh

と入力。きちんとG4WORKDIRが設定されていることを確認してください。

マルコフ連鎖を使って簡単な問題を解いてみよう

2010-04-06T01:49:00.008+09:00

"Cute Gyoza / Meat Dumpling Key Chain" by ♥ Rainbowcatz ♥

はじめに

人工無能、一時期流行りましたねー。

難そうなわりに以外と簡単で、しかも結果が面白いということでRubyやなんかで実装している記事をよく見かけました。そういった中でマルコフ連鎖が注目されていた時期が一時期ありましたが、正直あれはマルコフ連鎖を利用しているというよりは、別のアルゴリズムにマルコフ連鎖という名前をつけたというか、あんまり関連性がないなーという印象でした。

マルコフ連鎖の本質はそういうことではなく、状態の遷移確率を行列によって求めたり、定常分布を固有値問題で計算するところにある、と自分は考えています。マルコフ連鎖についてのより詳しい説明は他のブログにお任せするとして(自分がやるとすぐにボロが出てしまう)、今回はマルコフ連鎖を使って簡単な問題を解き、実際にマルコフ連鎖がどんな感じなのか確かめてみることにします。

もんだい

二つの箱A,Bが用意されている。2つの箱の内部は完全に遮蔽されており、外部からは確認することができない。また、Aの箱には赤玉がN個、Bの箱には白玉がN個入っている。

ここで、両方の箱に腕を突っ込み、それぞれの箱から玉を一個づつ取り出し、交換し、反対側の箱のなかに入れることにしよう。今回の場合、Aの箱には赤玉がN-1、白玉が1個入ることになる。もちろん、Bの箱には赤玉が1個、白玉がN-1個入っている。

この試行をひたすら続けて、そう、例えば無限回行ってみよう。この場合、Aの箱に入っている赤玉の数はどれくらいであろうか？個数に対応した確率分布を求めよ。

回答

1回、2回ならまだしも、無限回なんて計算できるかクソが！そんなふうに考えていた時期が俺にもありました。

まず、Aの箱における赤玉の個数をmとおきます。すると、これは状態mの遷移問題と考えることができ、マルコフ連鎖が使えます。

また、1期間後にとりうることのできるmの遷移はm-1, m, m+1の3通りしか存在しません。以上を利用して、各mにおける状態の遷移確率を求めることができます。

まず、mがm+1に遷移する確率(m->m+1)を求めてみましょう。この場合、確率はAの箱の白玉が選ばれる確率とBの箱の赤玉が選ばれる確率の積に等しいので、

と表せます。

同様にして、m->m, m->m-1 についても表していきます。

mがmに遷移する確率は(i)Aの箱の赤玉が選ばれる確率とBの箱の赤玉が選ばれる確率の積と(ii)Aの箱の白玉が選ばれる確率とBの箱の白玉が選ばれる確率の積の和((i)+(ii))に等しい。
mがm-1に遷移する確率は確率はAの箱の赤玉が選ばれる確率とBの箱の白玉が選ばれる確率に等しい。

以上の議論をまとめると、i->jへの遷移行列pijは以下のように表せます。

ここで、遷移行列はN×Nの次元を持っているものとします。

以上を利用して遷移行列を求め、pythonを用いて定常分布を算出しました。プログラムのソースコードを以下に示します。

#!/usr/bin/env python
#-*- coding:utf-8 -*-
"""
eig.py [N]

定常確率分布を記述したファイルを新しく書き込むプログラムです。
引数には箱に入っている玉の数を指定します。引数が指定されなかった場合には、N=10が用いられます。
"""

import sys

def main():
    number = int(sys.argv[1]) if len(sys.argv) == 2 else 10
    P = getStateMatrix(number)
    result = getFiniteStateSpace(P)
    fp = file("%i.dat" % number, "w")
    fp.write("\n".join([str(abs(num)) for num in result]))
    fp.close()

def getStateMatrix(n):
    """n個における遷移確率行列を返します"""
    from scipy import matrix
    P = [[getState(i, j, n) for j in xrange(n+1)] for i in xrange(n+1)]
    return matrix(P)

def getState(i, j, N):
    if   j == i-1: return float(i*i)/(N*N)
    elif j == i:   return float(2*i*(N-i))/(N*N)
    elif j == i+1: return float((N-i)**2)/(N*N)
    else: return 0.0

def getFiniteStateSpace(mat):
    """
    i->jに遷移する確率を行列(i,j)で表現した遷移確率行列から、定常分布を求めます。
    この関数は行列を転置する必要がない点に注意してください。
    
    mat : N×Nの遷移確率行列(numpyのmatrixでなければならない)
    戻り値 : N次元の定常分布ベクトル
    """
    from scipy import linalg
    la, v = linalg.eig(mat.T)
    result = zip(la, v.T)
    
    # 定常分布を求めるため、固有値が1に近い順で固有ベクトルをソートする
    result_sorted = sorted(result,
                           cmp=lambda x, y: cmp(abs(1.0-x), abs(1.0-y)),
                           key=lambda x: x[0])
    
    eig_value   = result_sorted[0][0]
    # 固有値が1近傍でなかったら例外を送出
    assert abs(eig_value - 1.0) < 1.0e-6, "the eigen value of FSS is apart from 1.0"
    
    fss = result_sorted[0][1]    
    return fss/sum(fss) # 正規化

if __name__ == "__main__":
 sys.exit(main())

実行は対象のソースコードを"eig.py"として保存してから、カレントディレクトリを適切なディレクトリに移動した後で、コマンドライン上から"python eig.py 10"のように指定してあげます。この場合、N=10(赤玉10個, 白玉10個)の際に遷移するであろう定常分布のデータが"10.dat"として書き出されます。

ここで、試しにN=6, 10, 50, 100における確率分布を載せておきます。

N = 6 のとき。離散的でよく分からないが、N=3付近の確率が一番高くなっている。

N = 10, 50 のとき。だんだん滑らかになっていく。

N = 100 のとき。N=50付近のピークが非常に鋭くなっている。

以上より、Nを増やせば増やすほど、N/2近傍により鋭いピークが表れるようになることが確認できます。

どーいうことが言いたいのか

結局のところ、お手玉の数を増やせば増やすほどmは中央付近で安定するわけです。「いや、最初に赤玉と白玉分けられているじゃん」それはそうなんですけど、今回の場合は無限回試行しているので、最初どのように分布してようが全く関係ないわけで。

最初は偏っていた分布がごちゃごちゃやっている内に中央付近に移動し、そこから多少のゆらぎはあるかもしれないが、結果的に中央へと帰ってしまうような、そういうイメージです。

まとめ

確率を無限回にすっ飛ばしたり、ある時点での確率分布を求めたり、モンテカルロ法にも使われているマルコフ連鎖はとっても賢くてすごいやつという記事なのでした

前のセメスターのレポートを微妙に公開してみるという話

2010-04-06T01:27:00.002+09:00

前のセメスターで、『固有値問題』『フーリエ変換』『モンテカルロ法』を利用した問題を作成し、レポートを提出するよう求められた授業がありました。

物理学科で『固有値問題』といいますと、普通は『時間に依存しないSchrodinger e.q.を解いて粒子の波動関数を求める』といった問題が一般的なんですが、そんな問題を解いてもみんな同じことやってるので正直つまらない。別に物理と絡めて解けと指示されたわけではありませんので、そんなら違う問題でもいーかということで自分はマルコフ連鎖を使用して、簡単な問題を解いてみるといったレポートを提出しました。

この問題は自分で設定したんですが、正直簡単なわりになかなかきれいにできてる(うわー自画自賛)ということで、このまま埋もれたままにするのはもったいないなーと思いまして、微妙に一部分を変更して公開してみることにしました。この記事利用してレポート書いたら多分すぐにバレるし、別に問題ないかなぁと。

Bloggerにはてなスターを設置するためには(まとめ)

2010-04-04T23:31:00.004+09:00

はてなのサービス『はてなスター』は、はてなダイアリーだけしか設置できない…というわけではなく、livedoor, blogger, seeseaなど、外部のブログでも導入することができます。

Bloggerにはてなスターを導入するための具体的な方法は『はてなスターをブログに設置するには』を参照すればすべて書いてありますが、この内容はテンプレートがクラシックなど、標準のものしか対応していません。なのでこのブログのようにいろいろとカスタマイズしているブログに設置するためには、いささか面倒な手順を踏まなければいけません。ここではBloggerでの変更手順を、上記の記事を補足する形で説明します。

まず、はてなスターのJavaScriptは、h3タグに囲まれている内容を記事のタイトル、h3タグ内のaタグのリンク先を記事のpermalinkと推測して処理を行っています。よって、ブログのHTML構造がその方式に従っている場合は特に変更する必要はなく、ただheadタグ内直下に

<script type="text/javascript" src="http://s.hatena.com/js/HatenaStar.js"></script>
<script type="text/javascript">
Hatena.Star.Token = 'your-original-token';
</script>

を挿入するだけで良いです。

しかし、h3タグの中にaタグが複数個入っている場合や、ブログタイトルにh2タグなど別のタグを使っている場合には、正解のURLを教えてあげる必要があるため、h3タグの直下に

<h3 class='post-title'>
<a expr:href='data:post.url' style='display:none;' title='permanent link'></a>

とダミーのaタグを挿入し、Hatena.Star.SiteConfigにカスタマイズした辞書を入れてやります。

また、複数のJSを使っている場合、時々JS同士が競合して正常な動作が行われない場合があります。筆者の場合ははてなスターのJSとGoogle AnalyticsのJSが競合してしまったために、Add Starボタンを押してもログインのダイアログが表示され、正常にスターが追加されないというエラーが生じていました(*1)。こういった場合、JavaScriptのロード順を変えることによって解決する場合があります。具体的には、このブログですとbodyタグ終了の直前に以下のようにコードを挿入することによって解決しました。

<script type='text/javascript'>
  var gaJsHost = (("https:" == document.location.protocol) ? "https://ssl." : "http://www.");
  document.write(unescape("%3Cscript src='" + gaJsHost + "google-analytics.com/ga.js' type='text/javascript'%3E%3C/script%3E"));
  </script>
<script type='text/javascript'>
  try {
  var pageTracker = _gat._getTracker("your-UA-code");
  pageTracker._trackPageview();
  } catch(err) {}</script>
<script src='http://s.hatena.ne.jp/js/HatenaStar.js' type='text/javascript'></script>
<script type='text/javascript'>
Hatena.Star.Token = 'your-original-token';
Hatena.Star.SiteConfig = {
entryNodes: {
'div.post': {
uri: 'h3 a',
title: 'h3',
container: 'h3'
}
}
};
</script>

要はGoogle Analyticsの次にはてなスターのJSが読み込まれるように修正したわけです。以上の手順によって、はてなスターの正常な動作を確認しました。この記事が皆様のお役に立てば幸いです。

(*1) おそらく競合の原因はクリックイベントのフックによるものと考えられます。