Home > Programing

2013.06.03

EmacsとgdbでVisual Studioライクなデバッグ環境

C言語を勉強するとき,文法や標準関数をマスターすることに並ぶくらい重要なのが,デバッガを使いこなせるようになることです.C言語を習いたての頃は,書いたコードが予期せぬ動作をしたり,ポインタの使い方を間違えてプロセスを落とすこと(セグメンテーション違反)なんてザラです.そういう失敗をたくさん経験して,自分の書いたコードの間違いに気づき,修正していくことでコーディング力が養われていきます.しかし,自分の書いたプログラムの何処が悪いのか,皆目見当もつかないような状況に陥ってしまうと,プログラミングが嫌いになってしまいます.

デバッガを文字通りに解釈すると「バグを取る」ためのツールで,プログラマにとってはコードの間違いを発見するための検証手段として欠かせません.これに加えて,正しいコードをデバッガ上で動かしてみることも,C言語の理解促進のために重要であると私は考えています.ループや条件分岐でプログラムの流れがどのように変わるのか,ポインタは具体的にどんな値を取っているのか,main関数のargcやargvにはどんな値がセットされているのか,関数呼び出しでどのようにスタックフレームができるのか,auto変数・static変数でメモリ空間がどのように使われているのか等,デバッガを使ってプログラムの動きを追ってみることで,C言語やコンピュータ・アーキテクチャの勉強になります.

かくいう私も20数年前,Microsoft Quick C 2.0という開発環境でC言語を勉強しました.20年以上前の製品ですが,Quick CはC言語の初心者に優しい逸品でした.MS-DOS上で動作するコンソール・アプリケーションながら開発統合環境として動作し,エディタ,コンパイラ,デバッガの機能が一通り揃っていました.特にデバッガが秀逸で,ブレークポイント,ステップ実行,ウォッチ変数などがエディタ上から操作できました.Quick Cのデバッガが使いやすかったお蔭で,挫折せずにC言語の勉強が出来たと言っても過言ではありません.この頃に習得したデバッグの技(操作体系)から離れられず,C/C++のコーディングにはMicrosoft Visual Studioを用いてきました.研究用のC言語のプログラムは,Windows上で開発してからlinux上に移植しています.

今回,東北大学のプログラミング演習Aという講義でデバッガの使い方を教えようと思い立ち,色々調べていくと,emacsとgdbでMicrosoft Visual Studioに近いデバッグ環境が実現できることを知りました.元々emacs+gcc使いの人には常識なのかもしれませんが,”emacs gdb”で検索して上位に出てくるようなページを見るだけでは,いまいちピンと来ていませんでした.あまりにも感動したので(本当は講義のためですが),emacsでC言語のプログラムをデバッグする例を紹介するムービーを作っちゃいました.サクサクとデバッグできる感じを伝えたかったので,ゆっくりと説明するような作りになっておらず,字幕を読むのはキツイかもしれません.そんなときは,ムービーを止めながら観て頂ければと思います.文字が滲んで見づらい場合は画質を上げて下さい(オリジナルは1080pのHD画質).

このムービーで使っているバグ入りコード(sum.c)とemacsの設定ファイル(gdb.el)です.この設定ファイルでは,gdb-many-window, gdb-use-separate-io-buffer, gud-tooltip-modeと,以下のようなショートカットキーを採用し,快適なデバッグ環境を実現しています.

  • [F1]: カーソル位置(もしくはリージョン)の変数の内容を表示(gud-print
  • [Shift]+[F1]: カーソル位置の変数をウォッチに追加(gud-watch
  • [F5]: プログラムの実行を再開(gud-cont
  • [Shift]+[F5]: デバッガを終了
  • [F6]: カーソル行まで実行(gud-until
  • [F9]: カーソル行のブレークポイントを設定/解除(gud-breakもしくはgud-remove
  • [Shift]+[F9]: main関数の先頭行にブレークポイントを設定(gdbに”break main”と入力するのと等価)
  • [F10]: ステップオーバー(現在行を実行,関数の中には入らない)(gud-next
  • [F11]: ステップイン(現在行を実行,関数の中にも入る)(gud-step
  • [Shift]+[F11]: ステップアウト(現在の関数を抜けるまで実行)(gud-finish

(GUDのデフォルトでは,ステップオーバーやステップインのために3ストロークが必要のようですが,それって疲れないのかなぁ)

2011.08.12

CRFツールのベンチマーク

CRFsuite 0.12のリリースに合わせて,チャンキングタスクによるベンチマークを更新しました.見所は,0.11→0.12でどのくらい速くなったのか,高速だと謳っているWapitiがどのくらい速いかです.その他,各学習アルゴリズムによる性能(学習速度,精度)差も,興味深い点かと思います.比較したツールは,

実験では,L2正則化(C=1),L-BFGSの終了条件は直近10回の反復で目的関数の改善率が1e-5を下回ったとき,平均化パーセプトロンの反復回数は50に固定,という条件にしました.

まず,CRFsuite 0.12で実装された平均化パーセプトロン(AP)は,非常に高速でありながら,L-BFGSやSGDに迫るタグ付け性能を示しています.精度はL-BFGSやSGDと比較して0.1~0.3ポイントくらい低下するようですが,学習データを1回(iteration)処理するのに要する時間は,1.7-1.9倍高速です.1回の反復に要する時間が短くなるのは,各事例に対してViterbiアルゴリズムを適用するだけで学習を進めることができ,手間のかかる対数尤度や勾配の計算を省くことができるからです.さらに,今回の実験では8回の反復でピークの性能を示したので,反復の回数を開発データセットなどでうまく調整できれば,非常に有用な学習アルゴリズムと言えます.

比較したソフトウェアの中で後発のWapitiはなかなか速いですが,一回の反復に要する時間を比較するとCRFsuiteより遅いようです.Wapitiのデフォルトのパラメータでは,L-BFGSの収束判定が非常に甘いため,30回強の反復で最適化が終了してしまいます(このため収束判定条件を統一しました).L-BFGSのiteration回数と精度のグラフを見ていると,30回付近では精度の改善が続いているため,この段階で最適化を止めてしまうのは早すぎます.WapitiはML compで長い間最速を維持していると謳っていますが,収束判定が甘いのが主要因だとすると,利用には注意が必要かもしれません.

CRFsuiteもML compにエントリーしたいのですが,系列ラベリングタスクの素性生成は各ソフトウェアに任されており,しかも,各ソフトウェアはどのタスクのデータが来たのか分からないため,有益な比較がしづらい状況にあります.つまり,与えられたデータの各列の意味がソフトウェア側からは(トリッキーな方法を使わない限り)分からないので,品詞タグ付けやチャンキングで同じ素性セットを使うことになってしまいます.さらに,素性数や学習アルゴリズム,パラメータなどで学習時間は大きく変わるので,純粋な実装の「良さ」をML compの実験環境で見積もるのは難しいと思います.

CRFsuite 0.11と0.12の比較(L-BFGS学習時)では,sparse素性で1.38倍,dense素性で1.46倍の速度向上が見られました.大した速度向上ではないかもしれませんが,0.12ではソースコードの可読性・再利用性を劇的に改善したので,この結果には満足しています.

2011.08.11

CRFsuite 0.12 released

CRFsuiteのバージョン0.12をリリースしました.このバージョンは変更点がてんこ盛りです.張り切りすぎたせいで,ドキュメントの更新に手間取り,コードが出来てからリリースまで1年くらい経過してしまいました.変更点はChangeLogの通りですが,このブログではいくつか補足しながら紹介します.

ソースコードの大がかりな再構成を行いました.特に,グラフィカルモデルの処理に関する部分(対数尤度や勾配の計算など)と,学習アルゴリズムに関する部分を分離しました.今回のソースコードの再構成により,新しい学習アルゴリズムの追加や,異なる形状のグラフィカルモデル(例えば属性とラベルバイグラムに対する素性や二次マルコフ素性)を追加しやすくしました.ソースコードの再構成をやるモチベーションは前々からあったのですが,@yuutatさんの「Fast Newton-CG Method for Batch Learning of Conditional Random Fields」の研究に混ぜて頂いたのをきっかけに,コードの大がかりな修正を行いました.新しいグラフィカルモデルを追加したり,ヘッシアン・ベクトル積のDPによる計算を追加することには,まだ未着手ですが,そのための基礎を固めました.

ソースコードを綺麗にしたついでに,平均化パーセプトロンPassive-AggressiveAdaptive Regularization of Weights (AROW) の学習アルゴリズムを追加しました.学習アルゴリズムは,-aオプションで指定することになりました.学習アルゴリズム毎のパラメータの一覧は,-Hオプションで見ることができます.

ソースコードの見直しを行っているときに,「あー昔は理解が足りず無駄な計算をしているなぁ」と思うところがいくつかあったので,その部分を修正しました.日記で書いていたSSEのexp計算ルーチンも取り込まれています.これらの修正により,だいたい1.4~1.5倍速くなりました.新しいベンチマーク結果を載せておきました.

系列の開始(BOS)と終了(EOS)の素性を自動的に生成するのを廃止しました.以前のCRFsuiteでは,BOS,EOSという特殊なステートを定義し,そのステートから各ラベルへの遷移素性(バイグラム素性)を作っていました.考えてみたら,系列の先頭と末尾にEOS/BOSを表す属性を追加すれば十分だったので,BOS/EOSに対する特殊なコードを削除しました.BOS/EOS素性を作りたいときは,各インスタンスの先頭と末尾のアイテムに,”__BOS__”や”__EOS__”みたいな属性を追加してください.

いくつかの学習パラメータの名前が変更になりました.重要なところでは,”regularization.sigma”が”c1″及び”c2″になりました(c1 = 1.0 / regularization.sigma; c2 = 0.5 / regularization.sigma^2).これにより,L1正則化,L2正則化が同時に使えるようになりました.L1正則化,L2正則化をOFFにするときは,それぞれ”c1″,”c2″を0にセットしてください(sigmaによる設定だと無限大にする必要があったので,cによる設定に変更しました).

交差検定(cross-validation)やテストセットによる評価(holdout-evaluation)を,データセットにグループ番号を割り振ることで実装しました.学習と同時にテストセットで評価を行うには,複数のデータセットをCRFsuiteに入力し,どのデータセットで評価を行うのか,グループ番号を指定します.交差検定を行うときは,N個のデータセットを明示的に与えるか,-gオプションを使い入力されたデータセットを自動的にN分割すればOKです.

データセットのフォーマットに関するドキュメントを書きました.データフォーマットはCRFsuiteの特徴の一つでもあるのですが,ユーザからの質問が多いので,きちんと書くことにしました.

これまでは,#でコメントを書くことを許していたのですが,自然言語のテキストから属性を作るときにうっかりと”#”を使ってしまうトラブルが起きそうなので,廃止しました.

予測されたラベル系列の条件付き確率(-pオプション)や,各ラベルの周辺確率(-iオプション)を出力できるようにしました.これらの機能には多くの要望が寄せられました.

CRFsuiteのC言語APIを大幅に改訂しました.ライブラリの名前はlibcrfからlibcrfsuiteに変更になりました.ライブラリが公開している関数名は,”crfsuite_”という文字列から始まるようになりました.C言語のAPIドキュメントを書きました.

CRFsuiteのライブラリを簡単に使うために,C++のラッパー(crfsuite_api.hppcrfsuite.hpp)を書きました(C++/SWIGのAPIドキュメント).さらに,C++のAPIをSWIGを経由して他の言語から呼び出せるようにしました.手始めにPythonのモジュールを構築し,学習とタグ付けがPythonから行えるようになりました.モジュールのビルド方法はREADMEを参照してください.

CRFsuiteに与えるデータセットを生成するためのサンプルスクリプトを改良しました.チャンキングのサンプル(chunking.py)は,わかりやすい属性名を出力するようになりました.また,CRFsuiteのPythonモジュールを利用して直接タグ付けを行う機能も実装されました.CRF++互換のテンプレートを処理するためのスクリプト(template.py)を追加しました.その他,固有表現抽出(ner.py),品詞タグ付け(pos.py)のサンプルも追加しました.これらのサンプルを改変するだけで,CRFsuiteを新しいタスクに適用できると思います.

メモリリークを何点か修正し,メモリ使用量を若干削減しました.

モデルファイルが存在しないときに落ちる問題を修正しました.

2010.09.02

exp(x)の高速計算 ~実験編~

最後に,e^xの計算がどのくらい高速化できたのか,実験してみました.アルゴリズムとして,以下のものを比較しました.

  • libcのexp関数
  • パデ近似を用いたもの(cephesの実装に近いです)
  • 5次~13次でテーラー展開を0を中心として行ったもの
  • 5次~13次の最良多項式近似を[-0.5, +0.5]の範囲で行ったもの
  • 5次~13次の最良多項式近似を[0, \log 2]の範囲で行ったもの
  • 5次~11次の最良多項式近似を[0, \log 2]の範囲で行ったものをSSE2で実装したもの

SSE2の実装で11次で打ち止めしているのは,11次と13次で演算精度の差が殆どないことが分かり,単に実装するのが面倒くさかったからだけです.0を中心としたテーラー展開や,最良多項式を[-0.5, +0.5]の範囲で構成する手法は,nの値を切り下げではなく,四捨五入(つまり0.5を足してから切り下げる)で求めます(Wapitisse_mathfun.hのexp計算ルーチンは四捨五入を使っています).

実験では,0を平均としたガウス分布で乱数を10,000,000個発生させ,その値のexpを全て計算するのに要する時間と,libcの演算結果を正しいものと仮定したときの誤差を測定しました.実験環境は,Intel(R) Xeon(R) CPU 5140 (2.33GHz) 上でDebian Linuxを動作させたアプリケーション・サーバーです.gccのバージョンは4.1.2で,コンパイルオプションは「-O3 -fomit-frame-pointer -msse2 -mfpmath=sse -ffast-math -lm」です.実験に用いたコードは,こちらからダウンロードもしくは閲覧できます.

横軸に計算時間,縦軸に誤差のRMS(二乗平均平方根)をプロットしたものを示します.Performance of fast exp computation

縦に一本線が入っているのが,libcの計算時間です.これを基準に誤差を測定しているので,エラーは0です.この線よりも左側に入っていないと,高速化の意味が無いことになります.各アルゴリズムは,左上が5次の多項式近似で,右下に向かって7次,9次,11次,13次と,精密かつ遅くなっていきます.だいたい,10^{-16}あたりで計算誤差の改善が止まっています.これは倍精度浮動小数点の有効桁数である16桁と大体一致しています.

計算精度に着目すると,テイラー展開 < 最良近似式 [-0.5, +0.5] < 最良近似式 [0, log 2] という順で,どの近似式も次数を上げていくと,10^{-16}付近まで誤差が改善します.ただ,アルゴリズムをCで実装してしまうと,libcと比較しても大した速度向上が得られません.5次まで近似精度を落としたとしても,2倍程度しか速くなりません.libcのexpの実装も相当に高速化されているので,一筋縄ではいかないようです.

これに対し,アルゴリズムをSSE2で実装したものは,4倍~8倍程度高速化されています.しかも,計算精度はSSE2化していないものとぴったり一致しているので,計算精度を犠牲にすることなく,高速化が達成できたことになります.e^xを大量に計算する状況では,11次の最良近似式を用いたSSE2ルーチンを使えば,演算精度を落とすことなく高速化できますし,近似式の精度が低くて良ければ,次数を下げてさらに高速化できます.次期バージョンのCRFsuiteでは,その他の改良も含めて,SSE2ルーチンを搭載する予定です.

exp(x)の高速計算 ~SSE2実装編~

exp(x)の高速計算 ~理論編~ の内容を元に,SSE2で実装してみます.ここからは,SSE2に関する基礎知識が必要になるのですが,すべて書くのは面倒なので,簡単にまとめると,

  • SSE2では,倍精度浮動小数点の2個の値に対して同じ演算をすると速い.
  • 現在のIntel CPUでは,SSE専用のレジスタ(128bit)が8個ある.
  • 加減乗除の算術演算や,ビットシフトなど,殆どの命令がSSE専用に準備されている.
  • メモリとSSE専用のレジスタ間でデータ転送を行うときは,メモリアドレスが16バイトでアライメントされている必要がある(厳密には「すべき」という言い方が正しいが).
  • C言語のソースコードの中でSSE用のコードを手っ取り早く書く方法は,コンパイラのintrinsicを使う方法.Microsoft Visual C++とgccでほぼ完全な互換性があるのが嬉しい.
  • コンパイラintrinsicは,アセンブラではないので,自分が意図したとおりの最適化コードが生成されるとは限らない.一応,コンパイラが最適化して速いコードを生成することになっているが,実際はそう上手くいかないことの方が多いので,コンパイル後のアセンブラコードを確認する必要がある.
  • SSEは,かなり昔にへるみさんのページで勉強しました.大体感触を掴んだら,MSDNのドキュメントなどを見ながら,コンパイラintrinsicを書けば,そんなに難しくないと思います.
  • SSEの命令を並べるときは,命令のロード・演算・ストアのパイプライン処理を頭の中で考えるべき.

だいたいこんな所ですかね.先ほどのexp(x)の計算ルーチン(5次最良多項式近似)をSSE2で書くと,こういうコードになります.

#include < emmintrin.h >

#ifdef _MSC_VER
#define MIE_ALIGN(x) __declspec(align(x))
#else
#define MIE_ALIGN(x) __attribute__((aligned(x)))
#endif

#define CONST_128D(var, val) \
    MIE_ALIGN(16) static const double var[2] = {(val), (val)}

void remez5_0_log2_sse(double *values, int num)
{
    int i;
    CONST_128D(one, 1.);
    CONST_128D(log2e, 1.4426950408889634073599);
    CONST_128D(c1, 6.93145751953125E-1);
    CONST_128D(c2, 1.42860682030941723212E-6);
    CONST_128D(w5, 1.185268231308989403584147407056378360798378534739e-2);
    CONST_128D(w4, 3.87412011356070379615759057344100690905653320886699e-2);
    CONST_128D(w3, 0.16775408658617866431779970932853611481292418818223);
    CONST_128D(w2, 0.49981934577169208735732248650232562589934399402426);
    CONST_128D(w1, 1.00001092396453942157124178508842412412025643386873);
    CONST_128D(w0, 0.99999989311082729779536722205742989232069120354073);
    const __m128i offset = _mm_setr_epi32(1023, 1023, 0, 0);

    for (i = 0;i < num;i += 4) {
        __m128i k1, k2;
        __m128d p1, p2;
        __m128d a1, a2;
        __m128d xmm0, xmm1;
        __m128d x1, x2;

        /* Load four double values. */
        x1 = _mm_load_pd(values+i);
        x2 = _mm_load_pd(values+i+2);

        /* a = x / log2; */
        xmm0 = _mm_load_pd(log2e);
        xmm1 = _mm_setzero_pd();
        a1 = _mm_mul_pd(x1, xmm0);
        a2 = _mm_mul_pd(x2, xmm0);

        /* k = (int)floor(a); p = (float)k; */
        p1 = _mm_cmplt_pd(a1, xmm1);
        p2 = _mm_cmplt_pd(a2, xmm1);
        xmm0 = _mm_load_pd(one);
        p1 = _mm_and_pd(p1, xmm0);
        p2 = _mm_and_pd(p2, xmm0);
        a1 = _mm_sub_pd(a1, p1);
        a2 = _mm_sub_pd(a2, p2);
        k1 = _mm_cvttpd_epi32(a1);
        k2 = _mm_cvttpd_epi32(a2);
        p1 = _mm_cvtepi32_pd(k1);
        p2 = _mm_cvtepi32_pd(k2);

        /* x -= p * log2; */
        xmm0 = _mm_load_pd(c1);
        xmm1 = _mm_load_pd(c2);
        a1 = _mm_mul_pd(p1, xmm0);
        a2 = _mm_mul_pd(p2, xmm0);
        x1 = _mm_sub_pd(x1, a1);
        x2 = _mm_sub_pd(x2, a2);
        a1 = _mm_mul_pd(p1, xmm1);
        a2 = _mm_mul_pd(p2, xmm1);
        x1 = _mm_sub_pd(x1, a1);
        x2 = _mm_sub_pd(x2, a2);

        /* Compute e^x using a polynomial approximation. */
        xmm0 = _mm_load_pd(w5);
        xmm1 = _mm_load_pd(w4);
        a1 = _mm_mul_pd(x1, xmm0);
        a2 = _mm_mul_pd(x2, xmm0);
        a1 = _mm_add_pd(a1, xmm1);
        a2 = _mm_add_pd(a2, xmm1);

        xmm0 = _mm_load_pd(w3);
        xmm1 = _mm_load_pd(w2);
        a1 = _mm_mul_pd(a1, x1);
        a2 = _mm_mul_pd(a2, x2);
        a1 = _mm_add_pd(a1, xmm0);
        a2 = _mm_add_pd(a2, xmm0);
        a1 = _mm_mul_pd(a1, x1);
        a2 = _mm_mul_pd(a2, x2);
        a1 = _mm_add_pd(a1, xmm1);
        a2 = _mm_add_pd(a2, xmm1);

        xmm0 = _mm_load_pd(w1);
        xmm1 = _mm_load_pd(w0);
        a1 = _mm_mul_pd(a1, x1);
        a2 = _mm_mul_pd(a2, x2);
        a1 = _mm_add_pd(a1, xmm0);
        a2 = _mm_add_pd(a2, xmm0);
        a1 = _mm_mul_pd(a1, x1);
        a2 = _mm_mul_pd(a2, x2);
        a1 = _mm_add_pd(a1, xmm1);
        a2 = _mm_add_pd(a2, xmm1);

        /* p = 2^k; */
        k1 = _mm_add_epi32(k1, offset);
        k2 = _mm_add_epi32(k2, offset);
        k1 = _mm_slli_epi32(k1, 20);
        k2 = _mm_slli_epi32(k2, 20);
        k1 = _mm_shuffle_epi32(k1, _MM_SHUFFLE(1,3,0,2));
        k2 = _mm_shuffle_epi32(k2, _MM_SHUFFLE(1,3,0,2));
        p1 = _mm_castsi128_pd(k1);
        p2 = _mm_castsi128_pd(k2);

        /* a *= 2^k. */
        a1 = _mm_mul_pd(a1, p1);
        a2 = _mm_mul_pd(a2, p2);

        /* Store the results. */
        _mm_store_pd(values+i, a1);
        _mm_store_pd(values+i+2, a2);
    }
}

SSE2のintrinsicに書き下しただけという感じですが,いくつか注釈を付けるとすると,

  • CONST_128Dマクロは,倍精度浮動小数点の定数を2個並べた定数を定義し,そのまま128bitのSSEレジスタにロードできるようにしている.
  • MIE_ALIGNマクロは,difference between gcc and vcのページを参照.他にも,Visual C++とgccの差異に関してまとめられていて,このページの情報はすごく有用です.
  • 128bit(倍精度浮動小数点×2)の演算を2つ並べて,パイプライン処理になりやすいように配慮している.そのため,引数nは4の倍数でなければならない.

exp(x)の高速計算 ~理論編~

ロジスティック回帰やCRFなどの対数線形モデルの学習でよく出てくるのが,expの計算です.これをSSE2を使って高速化するのが,今回のテーマです.まずは,背景の理論を説明します.

まず,指数関数e^xを2の指数関数2^aに変換することを考えます(なぜ2の指数関数かはいずれ分かります).

e^x = 2^a

両辺の自然対数をとり,aについて解くと,a = x/\log 2 .IEEE754など,2を基数とした指数部を採用している浮動小数点形式では,整数nに対して2^nを容易に構築できるので,上式の実数解aの代わりに整数,n = \left\lfloor x/\log 2 \right\rfloor を用い,e^xの大まかな値を計算することを考えます.ただし,\lfloor a \rflooraを超えない最大の整数を表します.

さて,anで近似したときの誤差(a - n)の範囲は,0 \leq a - n < 1 .誤差(a - n)\log 2を乗じたものを,bと定義すると,

b = (a - n) \log 2 = x - n \log2

上式の両辺の指数をとると,

e^x = 2^n e^b

ここで,naを超えない最大の整数なので,bの値域は,0 \leq b < \log 2

これらのことから,e^xは以下のステップで計算出来ることが分かります.

  1. n = \left\lfloor x/\log 2 \right\rfloor を計算する.
  2. b = x - n \log2 を計算する.
  3. e^bを定義域(0 \leq b < \log 2)において,テーラー展開もしくは最良近似多項式などで近似計算する.
  4. 浮動小数点形式を利用して2^nを構築し,e^bとの積を計算する.

ステップ1は簡単そうに見えますが,負の数を整数にキャストすると,0に近づいてしまうことに注意が必要です.つまり,doubleからintへのキャストは,正の値は切り下げ,負の値は切り上げになってしまいます.このことに注意すると,このステップは次のように書けます(a < 0で,ifによる分岐を作らないのが賢いやり方).

a = LOG2E * x;
a -= (a < 0);
n = (int)a;

ステップ2はそのまま書けばOKです.

ステップ3では,e^bを多項式に近似します.多項式近似というとテーラー展開が一般的ですが,「・・・点の周りについてテーラー展開すると」という表現の通り,ある点の近傍の近似は良いのですが,そこから離れた点では近似精度があんまりよくないという欠点があります.Cephesというライブラリでは,パデ近似(Padé approximant)を採用していますが,除算を1回だけ行う必要があります.最近のCPUは速くなったとはいえ,SSE2で倍精度のパック除算(DIVPD)をやるのに必要なレイテンシ・スループットは,70・70で,倍精度のパック乗算(MULPD)の7・6と比べると,除算が圧倒的に遅いことになります.そこで,今回は最良近似式を用います.なにやら難しそうに見えますが,Sollya というソフトウェアを使えば,簡単に最良多項式が求まります.以下の例では,[0, \log 2]の範囲で,5次の最良多項式を求めています.

> remez(exp(x), 5, [0, log(2)]);
0.99999989311082729779536722205742989232069120354073 + x *
(1.00001092396453942157124178508842412412025643386873 + x *
(0.49981934577169208735732248650232562589934399402426 + x *
(0.16775408658617866431779970932853611481292418818223 + x *
(3.87412011356070379615759057344100690905653320886699e-2 + x *
1.185268231308989403584147407056378360798378534739e-2))))

多項式がホーナー形式(horner scheme)で出てくるので,至れり尽くせりです.これをそのままC言語に書き直せば,ステップ3の実装は完了です.

y = 1.185268231308989403584147407056378360798378534739e-2;
y *= b;
y += 3.87412011356070379615759057344100690905653320886699e-2;
y *= b;
y += 0.16775408658617866431779970932853611481292418818223;
y *= b;
y += 0.49981934577169208735732248650232562589934399402426;
y *= b;
y += 1.00001092396453942157124178508842412412025643386873;
y *= b;
y += 0.99999989311082729779536722205742989232069120354073;

最後のステップ4ですが,もちろんpow(2, n)なんてやってしまったら,せっかくの高速化が台無しです.浮動小数点形式IEEE754では,2を基数として指数部が構成されているので,2^nはちょっとした工夫で一瞬で作れます.あんまり移植性がないコードですが,2^nは次のように作ります.

typedef union {
    double d;
    unsigned short s[4];
} ieee754;

ieee754 u;
u.d = 0;
u.s[3] = (unsigned short)(((n + 1023) << 4) & 0x7FF0);

こうすると,u.dで2^nの値にアクセスできます.なぜそうなるかは,IEEE754の仕様と睨めっこすればすぐに分かると思います.あとは,u.dとyの積を計算すれば,e^xの計算が出来たことになります.

2010.03.07

SimString (類似文字列検索ライブラリ) 1.0 released

SimStringという類似文字列検索ライブラリをBSDライセンスでリリースしました.類似文字列検索とは,文字列集合(データベース)の中から,クエリ文字列と似ているものを見つけ出す処理です.コンピュータは,正確に一致する文字列を探すのは得意ですが,表記揺れに出くわすと,途端に対応できなくなります.例えば,「スパゲティ」に対して,レストラン情報などを返すサービスにおいて,「スパゲッティ」や「スパゲティー」などの表記揺れが検索クエリに与えられると,通常のデータベースでは情報を提示することが出来ません.類似文字列検索を用いると,表記揺れが検索クエリに与えられても,「スパゲティ」という既知語を代替クエリとして提案したり,「スパゲティ」の情報をダイレクトに引き出すことができるようになります.

似てる語を探す技術って,文字列処理の基本中の基本で,自然言語処理では当たり前のように使われていてもおかしくないと思うのですが,研究ではほとんど見かけません.2つの文字列が与えられたときに,似てるか似てないかの識別する研究は,色々あるのですが,ある文字列が与えられたときに,どの文字列が近いかを高速に収集する方法は,あまり見かけません.検索エンジンを提供している会社では,検索クエリログに基づくクエリ訂正などを行っていますが,検索クエリログは誰でも入手できるわけではありません.類似文字列検索は,データベースやデータマイニングの分野で盛んに研究がされていますが,単語や用語レベルの文字列が簡単に扱えるツールが見つからなかったので,自分で作ってみました.

SimStringは,「文字列集合の中で,検索文字列との類似度がある閾値以上のものをすべて返す」という処理を高速に行えるように設計されています.類似度関数としては,コサイン係数,ジャッカード係数,ダイス係数,オーバーラップ係数をサポートしています.類似度を計算するときの文字列の特徴として,文字Nグラムを採用しています.

例えば,3グラム(N=3)を用いたとき,「スパゲティ」は「$$ス」「$スパ」「スパゲ」「パゲテ」「ゲティ」「ティ$」「ィ$$」と表現されます.このとき,文字列の先頭と末尾を表す特殊な記号「$」が挿入されていますが,これを用いるか用いないかはカスタマイズ可能です(デフォルトでは用いないことになっている).また,Nグラムの単位(すなわちNの値)も,自由に設定できます.また,日本語などのマルチバイト文字でも正確にNグラムが計算できるように,ユニコードをサポートしています.

SimStringは,類似文字列検索を正確に,かつ高速に行えるように設計されています.Google Web 1Tコーパス の英単語(13,588,391文字列)から,コサイン類似度が0.7以上の文字列を検索するのに必要な時間は,1クエリあたり平均 1.10 [ms] 程度です.

実装はC++で,ヘッダファイルをインクルードするだけで他のアプリケーションに組み込めるようになっています.また,PythonとRubyのSimStringモジュールも,ビルドできるようになっています.他言語のモジュールはSWIGを用いてラッパを自動生成しているので,Perlなどの言語にも対応できるはずです(ビルド報告,サンプルプログラムなど歓迎致します).Pythonで類似文字列が一瞬で検索出来るなんて,楽しい世界が広がると思いませんか?

$ python
Python 2.5.1 (r251:54863, Oct 15 2007, 23:38:19)
[GCC 3.4.6 20060404 (Red Hat 3.4.6-8)] on linux2
Type "help", "copyright", "credits" or "license" for more information.
>>> import simstring
>>> db = simstring.reader('web1tuni/web1tuni.db')
>>> db.measure = simstring.cosine
>>> db.threshold = 0.9
>>> db.retrieve('approximate')
('approximat', 'pproximate', 'approximate', 'approximated', 'approximatel',
'approximates', 'approximatey', 'anapproximate', 'approximatelt', 'approxima
tely', 'reapproximate', 'toapproximate')
>>> db = simstring.reader('web1tja/unigrams.db')
>>> db.measure = simstring.cosine
>>> db.threshold = 0.7
>>> print(' '.join(db.retrieve('スパゲッティー')))
スパゲッティ スパゲッテー スパゲティー スパッティー スパゲッティー スパゲッ
ティーニ スパゲッティー・ スパッゲッティー スパゲッティーニ・ ・・スパゲッテ
ィー スープスパゲッティー スパゲッティーモンスター
>>>

いろんなクエリを入れながら遊んでいるだけで,研究ネタが浮かんできそうです.この例では,単にコサイン類似度で類似文字列を検索しているだけなので,表記揺れ,語尾変化,派生語,その他別の概念を指すものが混ざって獲得されます.ここからルールや機械学習を用いて,所望の語を選別するのは,別タスクになります.

アルゴリズムの詳細については,今週東大で開催される情報処理学会の全国大会で発表する予定です.学会初日の朝一番の発表(1C-1)ですが,興味のある方は,是非遊びに来て下さい.

2009.10.07

SQLite3のINSERT文を高速化

こちらのFAQを参照.小ネタですが,以前に1回調べて,また忘れていたので,メモ.

2009.09.27

Classias 1.0 released

Classiasという分類のための機械学習アルゴリズムの実装を公開しました.今のところ,L1/L2正則化ロジスティック回帰(最大エントロピー法),L1/L2正則化L1損失線形カーネルサポートベクトルマシン(SVM),平均化パーセプトロンをサポートしています.学習アルゴリズムとしては,平均化パーセプトロン,L-BFGS法,OWL-QN法,Pegasos,Truncated Gradient(L1-FOLOS)を実装してあります.カーネルは使えませんが,線形識別モデルを高速に学習できるようになっています.二値分類,多クラス分類,候補選択(明示的に与えられた候補の中からスコア最大のものを選ぶタスク)をサポートしています(SVMは今のところ二値分類のみ).

このツールはもともと,最大エントロピー法を自分で使うために実装したもので,作り始めてからもう2年くらい経過しています.去年のColingやEMNLPの論文の機械学習の箇所は,このツールで実験しています.文字列による素性,自動分割交叉検定など,自然言語処理の実験に便利な機能に重点を置いて作ったつもりです.詳しくは,使い方を参照してください.

ソースコードも,インスタンスのデータ構造,損失関数,学習アルゴリズムなどのコンポーネントを,C++テンプレートクラスとして実装するという設計になっています.L-BFGS法を実装しているlibLBFGSを除けば,ヘッダファイルをインクルードするだけでClassiasを利用したプログラムが書けるようになっています.クラスリファレンスやサンプルコードはドキュメントを参照してください.新しい学習アルゴリズムも,簡単に追加できるようになっているのですが,ドキュメントをもう少し充実させていかなければなりません・・・.

LIBLINEARやLIBSVMと一緒にrcv1.binaryでパフォーマンスを測ってみました.まず,分類精度に関してはロジスティック回帰とL1損失SVMに大差がなく,平均化パーセプトロンはちょっと悪いという予想通りの結果で,Classias,LIBLINEAR,LIBSVMのどれも同程度の分類精度を出しています.

学習速度に関しては,LIBLINEAR速すぎです.L1正則化のロジスティック回帰ではそれほどスピードの差はありませんが,L2正則化のロジスティック回帰では,LIBLINEARの方が4倍くらい速くなっています.アルゴリズム的には,L-BFGS法と信頼空間ニュートン法の戦いで,後者はヘッシアンを使っているので,速くなるのも頷けます.ただ,収束の判定条件が違うと思うので,後できちんと検証してみようと考えています.

SVMに関しては,LIBLINEARの圧勝(LIBLINEARの方が10~20倍高速)でした.ClassiasのPegasosも健闘していますが,学習率や収束判定の調整で10~20倍の速度差を埋められるかは疑問です.Dual coordinate descentはやっぱり強いですね.こちらも後で実装して遊んでみたいです.

実装するといろいろ分かることがあるので,これからも実験しながら育てていきたいと考えています.

2009.09.24

CRFsuite 0.9 released

CRFsuite 0.9をリリースしました.libLBFGS 1.8と組み合わせてビルドすることが出来なくなっていたので,その問題の修正のみです.また,試験的にLinuxバイナリの配布を始めてみました.libLBFGSと静的にリンクしてあるので,このバイナリを落として実行するだけで動くはずです.

いろいろやっていたら,もう2ヶ月も経っていた・・・.

Next »