2011/04/19

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

"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を行いたい場合に役に立つ。

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

を手に入れたとする。この場合、i < jにおいてbj+bi-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の出しているテキストは分かりやすい。すごいことです。

2011/03/31

Moving

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

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

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

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

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

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

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

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

2011/03/11

自分は無事です

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

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

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

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

2011/02/16

OMakeのContributionsに載ったった


やったね!

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

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

2011/02/15

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

"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++と画像処理の両方に秀でていなければ記述できないので、使っている人が殆どいない。
  • 英語のサイトですら文献が殆どない
正直なところ、短所の部分がかなり厳しい。かなりC++をやりこんでいないと使えないということは、それだけでユーザ層を大分狭めます。仕方がないことですけど。

Viewの概要
GILでは、画像に対する操作はすべて『View』を起点にして行われます。Viewは以下のような構造を持ち、ジェネリックな操作を行えるよう工夫してあります。


View
|-- Locator
|    `-- Pixel
|         |-- Channel
|         `-- Layout
|              |-- Color Space
|              `-- Channel Mapping
`-- dimensions

(上図は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.