映像奮闘記

Numpy.hpp - .npyファイルを手軽に扱うためのヘッダライブラリ

2014-03-02T21:36:00.001+09:00

Introduction

ご存知の通り，C++は画像処理などの大規模計算に適した言語であり，Pythonはこれらの結果を簡単にscipyやmatplotlibを使って可視化できます．そのため，多くの研究者はコアの計算をC++で行う一方で，データ整理や実験結果の可視化にはPythonを用いるアプローチが一般的でした．

その際に重要となってくるのは，Numpy標準のファイル形式である.npyファイルを読み書きするための，C++ライブラリの存在です．これまでにもlibnpyなどに代表されるいくつかのライブラリが存在していましたが，予めコンパイルする必要があるため気軽に使えない，またMSVCなどの環境で使用できないなどのいくつかの欠点がありました．

Numpy.hppはこれらの問題を解決するために作られたライブラリです．Numpy.hppの特徴は次の2つです．

Header-only: Numpy.hppはコンパイルする必要がありません．Numpy.hppを使用するためにあなたがすることは，ただ#include "Numpy.hpp"をあなたのソースファイルに加えるだけです．
外部ライブラリに依存しない: C++03のSTL環境だけで完結しているので，MSVCを含むどの環境下でも容易に動作します．

Numpy.hppは以下のURLよりダウンロードできます．
https://gist.github.com/rezoo/5656056

使用方法

このセクションでは，Numpy.hppの使い方について簡単に説明します．

`.npy`ファイルの読み込み

template<typename Scalar>
void LoadArrayFromNumpy(
    const std::string& filename, std::vector<int>& shape,
    std::vector<Scalar>& data);

.npyファイルを読み込みたい場合は，aoba::LoadArrayFromNumpy()関数を使用してください．

const std::string& filename: 入力先のファイル名
std::vector<int>& shape: テンソルの階数とその長さを格納するベクトルの出力先
std::vector<Scalar>& data: データ本体の出力先

std::vector<T>の型Tは指定された.npyファイルの型と同一でなければなりません．もし異なる型を指定した場合には，Numpy.hppは例外を送出します．

aoba::LoadArrayFromNumpy<T>()関数の簡単なサンプルは次の通りです:

#include <iostream>
#include "Numpy.hpp"

// Now I assume np.load("file.npy").shape == (4, 5)
std::vector<int> s;
std::vector<double> data;
aoba::LoadArrayFromNumpy("/path/to/file.npy", s, data);
std::cout << s[0] << " " << s[1] << std::endl; // 4 5
std::cout << data.size() << std::endl; // 20

Numpy.hppはまた，入力先のテンソルの階数が必ず固定である場合に使えるいくつかのショートカット関数を提供しています．

#include "Numpy.hpp"

int rows, cols;
std::vector<double> data;
aoba::LoadArrayFromNumpy("/path/to/file1.npy", rows, cols, data);
std::cout << rows << " " << cols << std::endl;

aoba::LoadArrayFromNumpy(
    "/path/to/file2.npy", a0, data); // 1st-order tensor
aoba::LoadArrayFromNumpy(
    "/path/to/file3.npy", b0, b1, b2, data); // 3rd-order tensor
aoba::LoadArrayFromNumpy(
    "/path/to/file4.npy", c0, c1, c2, c3, data); // 4th-order tensor

`.npy`ファイルの保存

template<typename Scalar>
void SaveArrayAsNumpy(
    const std::string& filename, bool fortran_order,
    int n_dims, const int shape[], const Scalar* data);

.npyファイルを保存したい場合は，aoba::SaveArrayAsNumpy()関数を使用してください．

const std::string& filename: .npyファイルの出力先
bool fortran_order: このデータの格納順がcolumn-majorである場合は，このフラグをtrueに設定してください．
int n_dims: テンソルの階数
const int shape[]: テンソルの長さを表す配列のポインタ
const Scalar* data: テンソルの実データを指すポインタ

aoba::SaveArrayAsNumpy<T>()関数の簡単なサンプルは次の通りです．

#include "Numpy.hpp"

// Now I assume s.size() == 2
aoba::SaveArrayAsNumpy("/path/to/file.npy", false, 2, &s[0], &data[0]);

aoba::LoadArrayFromNumpy<T>()と同様に，Numpy.hppはまたいくつかのショートカット関数を提供しています．

// the following functions assume that fortran_order == false.
aoba::SaveArrayAsNumpy("/path/to/file.npy", rows, cols, &data[0]);

aoba::SaveArrayAsNumpy(
    "/path/to/file2.npy", a0, &data[0]); // 1st-order tensor
aoba::SaveArrayAsNumpy(
    "/path/to/file3.npy", b0, b1, b2 &data[0]); // 3nd-order tensor
aoba::SaveArrayAsNumpy(
    "/path/to/file4.npy", c0, c1, c2, c3, &data[0]); // 4th-order tensor

もし何か質問などありましたら，お気軽に私(rezoolab at gmail.com)に連絡ください :) Enjoy!

Numpy.hpp - Simple Header-only Library for handling .npy files

2014-03-02T16:36:00.001+09:00

Introduction

As you know, C++ is a suitable language for computing the large-scale problems, and python is a great language in terms of plotting and visualizing these results by using several attractive libraries such as scipy and matplotlib. For that reason, it is quite a common approach for researchers to use both languages – C++ and python – for the large-scale computation and the visualization.

The important thing in this case is an existence of a C++ library for reading/writing .npy files (a default file extension of numpy). Although several libraries have been proposed to handle such .npy files, they require precompiling and cannot work well in several environments (e.g. MSVC).

To handle these problems, I reimplemented a new C++ library called Numpy.hpp. Comparing with the existing libraries, Numpy.hpp has two advantages:

Header-only: you do not need to precompile numpy.hpp. All you have to do is just add #include "Numpy.hpp" into your source file.
No external dependencies: numpy.hpp uses only STL libraries in C++03. This means that numpy.hpp can work well in any environments including MSVC.

You can download Numpy.hpp from the following URL:
https://gist.github.com/rezoo/5656056

Usage

In this section I will briefly introduce the usage of Numpy.hpp.

Load `.npy` files

template<typename Scalar>
void LoadArrayFromNumpy(
    const std::string& filename, std::vector<int>& shape,
    std::vector<Scalar>& data);

If you wish to load .npy files in C++, please use aoba::LoadArrayFromNumpy() functions.

const std::string& filename: the filename of the input.
std::vector<int>& shape: the output of the tensor’s rank and its lengths.
std::vector<Scalar>& data: the destination of the data.

Please note that the type of std::vector<T> must be equivalent to the type of the specified .npy file. Numpy.hpp will send an exception if you specify the different type.

A simple example of aoba::LoadArrayFromNumpy<T> function is as follows:

#include <iostream>
#include "Numpy.hpp"

// Now I assume np.load("file.npy").shape == (4, 5)
std::vector<int> s;
std::vector<double> data;
aoba::LoadArrayFromNumpy("/path/to/file.npy", s, data);
std::cout << s[0] << " " << s[1] << std::endl; // 4 5
std::cout << data.size() << std::endl; // 20

Numpy.hpp also provides several shortcuts that can be used where the order of tensor is fixed.

#include "Numpy.hpp"

int rows, cols;
std::vector<double> data;
aoba::LoadArrayFromNumpy("/path/to/file1.npy", rows, cols, data);
std::cout << rows << " " << cols << std::endl;

aoba::LoadArrayFromNumpy(
    "/path/to/file2.npy", a0, data); // 1st-order tensor
aoba::LoadArrayFromNumpy(
    "/path/to/file3.npy", b0, b1, b2, data); // 3rd-order tensor
aoba::LoadArrayFromNumpy(
    "/path/to/file4.npy", c0, c1, c2, c3, data); // 4th-order tensor

Save `.npy` files

template<typename Scalar>
void SaveArrayAsNumpy(
    const std::string& filename, bool fortran_order,
    int n_dims, const int shape[], const Scalar* data);

If you with to save .npy files, please use aoba::SaveArrayAsNumpy() functions.

const std::string& filename: the destination of the .npy file.
bool fortran_order: if you store the data as column-major order, please set this flag as true.
int n_dims: the rank of the tensor.
const int shape[]: the array that contains the tensor’s lengths.
const Scalar* data: the raw data of the tensor.

A simple example of aoba::SaveArrayAsNumpy<T> function is as follows:

#include "Numpy.hpp"

// Now I assume s.size() == 2
aoba::SaveArrayAsNumpy("/path/to/file.npy", false, 2, &s[0], &data[0]);

As with aoba::LoadArrayFromNumpy<T>(), Numpy.hpp also provides several shortcut functions.

// the following functions assume that fortran_order == false.
aoba::SaveArrayAsNumpy("/path/to/file.npy", rows, cols, &data[0]);

aoba::SaveArrayAsNumpy(
    "/path/to/file2.npy", a0, &data[0]); // 1st-order tensor
aoba::SaveArrayAsNumpy(
    "/path/to/file3.npy", b0, b1, b2 &data[0]); // 3nd-order tensor
aoba::SaveArrayAsNumpy(
    "/path/to/file4.npy", c0, c1, c2, c3, &data[0]); // 4th-order tensor

If you have any questions, do not hesitate to ask me (rezoolab at gmail.com) :) Enjoy!

ノンパラメトリック平均場近似と確率伝播法の提案

2013-12-31T23:06:00.000+09:00

2013年も残り少なくなってきました．そこで，今回は研究の締めくくりとして，今年のCVPRで発表した論文について紹介したいと思います．

僕の研究は，大雑把に言うとグラフィカルモデルの最適化手法に関するものです．本研究の貢献を簡単にまとめると，次のようになります．

グラフィカルモデル内の確率変数が連続的である場合，計算上の理由から，変分ベイズ法ならびに確率伝播法で扱える分布のモデルには制限があります．提案手法を使うことで，非指数型の分布族に関しても変分ベイズ法を効率的に適用できるようになりました．
さらに，同様の議論を確率伝播法にも適用することで，確率伝播法をノンパラメトリックに扱うこともできるようになりました．

実際の論文では対象のグラフィカルモデルを一階のマルコフ確率場に限定していますが，変分ベイズ法などの別の領域にも容易に適用できます．また，議論を簡単にするため実際の論文では確率変数の離散化を焦点に当てた形で書いていますが，その流れに沿って書いてもあまり面白く無いので，ここではどのような過程でこのような研究に取り組んだのか簡単に書いていきたいと思います．

グラフィカルモデルは機械学習の多くの分野で使われてきた確率モデルです．この確率モデルを用いて対象の問題を高精度に解くためには，確率変数の周辺化が重要となってきます．周辺化の計算を直接行うことは困難であるため，大きく分けて2つの方法論で周辺分布を近似的に求めることが一般的です．一つはギブスサンプリングに代表される，大量のサンプルをばらまいてモンテカルロ的に推定解を求める方法，もう1つは平均場近似(変分ベイズ法)や確率伝播法に代表される，変分原理に基づいた近似的推論法を使う方法です．後者は得られる推定解の精度が比較的悪いものの，特定の反復方程式に従って周辺分布を更新するだけで良いため，高速に周辺分布の推定解が求められる利点を持ちます．本研究は後者の推定手法に関するものです．

変分原理に基づいた周辺分布の反復計算は，グラフィカルモデル内の確率変数が連続的であるか離散的であるかで異なります．変数が連続的である場合，周辺分布の密度関数は少数のパラメータからなるパラメトリックな連続分布として表現されます．そして，平均場近似，確率伝播法の両方とも，実際の最適化には変分原理に従って導出された特定の反復方程式に従って，対象のパラメータを更新していきます．しかしながら，このような更新は更新後の分布が更新前の分布と同じパラメータで表現される必要があるため，扱える分布モデルの範囲は非常に限られていました．以上の背景から，僕は平均場近似ならびに確率伝播法をノンパラメトリックに扱えるような方法論を新しく提案すれば，良いインパクトになるのではないかと考えました．

このような元々の平均場近似や確率伝播法では扱えない分布をどうにかして扱えるよう拡張できないだろうかという要望は多く，これまでにいくつかの試みが提案されています．これらの試みはすべて次の混合分布の仮定を下にしています．すなわち，平均場近似や確率伝播法で扱う近似分布やメッセージを，次の混合ガウス分布で仮定することです．

もし上の仮定の下で混合ガウス分布のパラメータに関する反復方程式を導出できたのであれば，対象のグラフィカルモデルがどのような分布であったとしても平均場近似や確率伝播法を適用できる．つまり，平均場近似や確率伝播法の可能性を一気に広げることができます．

しかし，このアイデアを直接適用してもあまり上手くいかないことが知られていました．なぜなら，変分原理は計算の段階で混合ガウス分布のエントロピーを計算する必要があるためです．混合ガウス分布のエントロピーは計算困難であるため，このままでは反復方程式を解析的に導出できません．

このような問題に対応するための方法として，これまで大きく分けて2つのアプローチが提案されてきました．1つは解析解を出すために，混合分布のエントロピーにJensenの不等式を用いてさらなる近似を行い，解析解を求める方法です[3]．しかしながら，この方法論は平均場近似や確率伝播法で導入した近似分布の他に，さらなる近似をおくことになります．もう1つはモンテカルロ法を用いて，混合ガウス分布のパラメータの反復解を力技で求める方法です[1][2]．しかしながら，この方法論ではサンプル数に依存せず高速に推定解を求められる，確率伝播法の利点を失うことになります．

これ以外にエントロピーの計算をどうにかして解決できないだろうかと思い，思考を巡らせました．前述の通り，そもそもの問題は混合ガウス分布のエントロピーが計算困難であることでした．ガウス分布の裾が別のガウス分布に干渉してしまうため，それぞれ独立した分布として扱うことができないのです．そうであれば，最初の出発点に混合ガウス分布を選ぶのではなく，互いに干渉しない別の分布を選べば良いのではないでしょうか？例えばヒストグラムのような，互いに干渉せず，エントロピーが容易に計算できる混合矩形分布を採用するのは？

ここで，hは変数iにおける，s番目の矩形分布を表します．本研究で新しく提案したアイデアは基本的にはこれだけです．この矩形分布の位置とサイズは，重なっていなければ変数空間のどの箇所に配置しても構いません．すなわち，この表現を用いて変数空間の重要な箇所を密に離散化し，重要でない箇所を疎に離散化することもできます．さらに，変数毎に配置する矩形分布の個数は異なっていても構いません．すなわち，変数の重要性に従って，変数の離散化の度合いを変えることもできます．

次に，上の表現を用いて本来の連続的な最適化問題をパラメータαの最適化問題へ変換し，パラメータに関する反復方程式を導出しました．詳しくは言及しませんが，最終的に導出した反復方程式の形は対象の変数が離散的な場合の平均場近似，確率伝播法の形と非常に似通っています．唯一の違いはエネルギー関数に加わる補正項の存在です．この補正項は変数空間の離散化の「非一様性」を補正する役割をもちます．すなわち，提案手法によって，グラフィカルモデルの連続的な変数空間を非一様に離散化し，ノンパラメトリックに扱えるようになったのです．

以上，急ぎですが本研究の解説記事を書かせていただきました．興味を持っていただけた方は，次のURLからダウンロードしていただければ幸いです．

http://www.vision.is.tohoku.ac.jp/index.php/download_file/view/35/177/167/

参考文献
[1] A. Ihler et. al., Particle Belief Propagation, AISTATS, 2009
[2] E. B. Sudderth et. al., Nonparametric Belief Propagation, CVPR, 2003
[3] S. J. Gershman et. al., Nonparametric Variational Inference, ICML, 2012

条件付き確率場の推論と学習

2013-04-26T00:03:00.002+09:00

条件付き確率場の推論と学習 from Masaki Saito

　前回の研究室ゼミで話した，条件付き確率場についての簡単な紹介スライドを公開します．内容は基本的なPairwise CRFについての話と，それが実際にどのような問題に対して応用がなされているか，またもっとも使われる最適化手法の一つである，平均場近似と確率伝搬法についての簡単な解説と，CRFの初歩的な学習法といった流れです．

補足:
　CRFの学習は自然言語処理など他の分野でもよく使われている学習法なのでおそらく知っている方は多いと思うのですが，普通はダイレクトに対数尤度を最大化することでパラメータを求めています．しかし，今回はKL距離という確率分布間の近さを図る指標を用いて対数尤度を拡張し，パラメータを求めています．なぜわざわざそうするのかというと，だいたいの最適化手法がKL距離の最小化という議論に帰着できて，すっきりするからです．KL距離を使えば平均場近似，確率伝搬法，最尤推定，あと時間の都合上話せなかったTRWやGeneralized Belief Propagationなどの最適化手法がすべて統一的に理解できます．これによって覚えることが少なくなって楽できるだけでなく，このKL距離を使った新しい最適化手法を提案することもできます．

近況報告

2013-04-25T02:48:00.005+09:00

お久しぶりです．細々とした近況はTwitterに書いているのですが，あくまでTwitterはつぶやきであって，記事でない．定期的に自分の中でまとめた内容をブログに吐き出す必要がありますね．

修士論文を提出し，晴れて今年度から博士後期課程の1年生です．不安が全くないと言えば嘘になりますが，それでも毎日好きなことを学び，研究できる環境は非常に楽しい．みなさんもぜひ博士課程に進学しましょう．唯一不満に思うところは頻繁に書類を書く業務が入るあたりですが，仕方のないことです．
後できちんとした内容を書きますが，2013年のCVPR(コンピュータビジョンのトップカンファレンス)に採択されました．去年を含めると今回で2回目の採択ですが，前回と異なるのが，今回はポスターではなくオーラルセッションでの採択というところです．オーラルセッションは非常に狭き門(採択率 60/1870 = 3.2%)でありますので，そのような研究成果を得られたことは非常に嬉しいですし，来年もぜひ通しておきたいところです．
博士後期課程になりましたので，これからは本名メインで活動していきます．プロフィール画像も更新しました．前に使ってた画像は高校時代からのものでしたから，大体7年位使っていたんでしょうか．というかこのブログ，今年で7年目になるのですね．
自身の宣伝もかねて，論文やコードなどはこれから作る自己紹介のページにて積極的に公開していきます．論文についてはできるだけ読んでもらえるよう，簡単な解説も記事にしていきたいと思っております．

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

2012-10-20T21:19:00.000+09:00

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を用いて実装する

2012-08-21T23:58:00.000+09:00

ご存知の通り、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の学習手法についての簡単なまとめ

2012-08-18T22:58:00.001+09:00

近年の機械学習では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)を用いた近似最近傍探索についての簡単な解説

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

"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倍以上の圧縮効率です。このように、少数のコードブックサイズで効率よくベクトルを表現できることが直積量子化の利点です。

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

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

クエリベクトルを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年ローンとかで家を購入するような場合には全く通用しません。そもそもこのご時勢に長期のローンを組むなんて狂気の沙汰です。狂っている。

映像奮闘記

Numpy.hpp - .npyファイルを手軽に扱うためのヘッダライブラリ

Introduction

使用方法

.npyファイルの読み込み

.npyファイルの保存

Numpy.hpp - Simple Header-only Library for handling .npy files

Introduction

Usage

Load .npy files

Save .npy files

ノンパラメトリック平均場近似と確率伝播法の提案

条件付き確率場の推論と学習

近況報告

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

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

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

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

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

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

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

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

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

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

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

Moving

自分は無事です

OMakeのContributionsに載ったった

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

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

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

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

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

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

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

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

`.npy`ファイルの読み込み

`.npy`ファイルの保存

Load `.npy` files

Save `.npy` files

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

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