2014/03/02

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

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()関数を使用してください.

  1. const std::string& filename: 入力先のファイル名
  2. std::vector<int>& shape: テンソルの階数とその長さを格納するベクトルの出力先
  3. 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()関数を使用してください.

  1. const std::string& filename: .npyファイルの出力先
  2. bool fortran_order: このデータの格納順がcolumn-majorである場合は,このフラグをtrueに設定してください.
  3. int n_dims: テンソルの階数
  4. const int shape[]: テンソルの長さを表す配列のポインタ
  5. 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

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.

  1. const std::string& filename: the filename of the input.
  2. std::vector<int>& shape: the output of the tensor’s rank and its lengths.
  3. 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.

  1. const std::string& filename: the destination of the .npy file.
  2. bool fortran_order: if you store the data as column-major order, please set this flag as true.
  3. int n_dims: the rank of the tensor.
  4. const int shape[]: the array that contains the tensor’s lengths.
  5. 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/31

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

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/26

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

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

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

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

2013/04/25

近況報告

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