遥かへのスピードランナー

シリコンバレーでAndroidアプリの開発してます。コンピュータービジョン・3D・アルゴリズム界隈にもたまに出現します。

DMC(Dynamic Markov Coding)のによるデータ圧縮プログラムを書いてみた

MG

最近Managing Gigabytes勉強会に参加しているのでせっかくなので、この本に載っているアルゴリズムを使ってプログラムを組んでみました。

今回実装したのは、「2.5 SYMBOLWISE MODELS」の後半で説明されている「Dynamic Markov Coding(DMC)」です。書籍の他に、元論文「G. Cormack and R. Horspool, "Data compression using dynamic Markov modelling,"」を参考にしました。

実装はC++で行い、ソースはgithubに置きました。（CentOS 5.2+gcc 4.1.2で動作確認済）
http://github.com/thorikawa/MG/tree/master/dmc

以下、アルゴリズムの概要と実装上の工夫などをまとめてみます。
意見・指摘などは絶賛大歓迎です。

DMC(Dynamic Markov Coding)とは

DMCは適応型の算術符号化の一種で、データの流れ全体をマルコフ連鎖で表します。符号化の過程で、マルコフ連鎖の構造と確率を変化させていき、各状態における次状態への遷移確率を元に符号化します。

上の図は元論文から抜粋したものですが、状態Aから出ている矢印とその添字は、状態Aで読み込む次のビットが0である確率が50%、1である確率が50%であることを示しています。状態Aで0のビットを読み込むと、50%という確率に基づいて符号化を行った後に、状態Aから遷移する確率を更新して、状態Bに遷移します。

圧縮率を高めるためには、各状態における遷移確率が偏るようにマルコフ連鎖の構造を変化させていく必要があります。そのため、DMCではクローニングという仕組みを導入しています。クローニングを示した図が以下です。（元論文より抜粋）

(a)がクローニングを行う前で、AとBという２つの状態から状態Cに遷移していますが、ある一定の条件を満たしたときに、状態CはCとC’に分離します。ある条件というのは以下のソースを参照してほしいですが、元論文では２つのパラメタ閾値１と閾値２を設定して、AからCへのの遷移カウントが閾値１を超え、かつAからCへの遷移カウントとA以外からCへの遷移カウントとの差が閾値２を超えた場合にクローニングが行われます。

  // クローニングを行う。
  void DoCloning (uchar bit) {
    uint trans_count = current_state_->GetTransCount(bit);
    State *next_state = current_state_->GetNextState(bit);
    uint next_count = next_state->GetTotalCount();
    
    if (trans_count > cloning_threshold1_ && (next_count - trans_count) > cloning_threshold2_) {
      State *new_state = new State(last_state_++);
      current_state_->SetNextState(bit, new_state);
      double ratio = (double) (trans_count + 1) / (next_count + 2);
      // ratioの正規化(CloningはCloning後の次状態のtrans_count <= 現状態のtrans_countを保証していないので、ratioが1を超える場合がある。
      if (ratio > 1) {
        ratio = 1;
      }
      for (int i=0; i<=1; i++) {
        uint next_count_bit = next_state->GetTransCount(i);
        uint new_next_count_bit = (uint)(ratio * next_count_bit);
        new_state->SetNextState(i, next_state->GetNextState(i));
        new_state->SetTransCount(i, new_next_count_bit);
        next_state->SetTransCount(i, next_count_bit - new_next_count_bit);
      }
      cloning_count_++;
    }
  }

クローニングが行われないと、A→C→Dと、B→C→Eの２ケースのビット遷移が極端に多い場合でも、Cに到達した時点でAからきたかBからきたかの情報が失われてしまい、C→D・C→Eの遷移確率が均一化されてしまうので、効率のよい圧縮はできません。クローニングを行うことにより、出現頻度の多いパターンがたどる経路を分離させることができ、各状態における遷移確率を偏らせ、効率のよい圧縮が出来る、と理解しています。

Guazzoによる算術符号化アルゴリズム

次に、各状態における遷移確率から符号化を行う部分です。
Managing GigabytesにはDMCに特化したこの部分のロジック詳細は記載されていませんが、今回の実装では元論文でも紹介されているGuazzoによる算術符号化アルゴリズムを利用します。

Guazzoのアルゴリズムでは、整数値を取る数直線上をデータの出現確率に応じて分割します。今回はデータの入出力をビット単位で扱うこととしまうので、0と1に対応する二区間に数直線を分割していけばよいです。
数直線の幅はNを任意の整数値として[0x0,(1<<(N+1))-1]で初期化し、データを1ビット読み込むたびに、新たな上限と下限を計算していきます。そして計算し直した上限と下限を比較して両者の上位ビットが一致する場合は、一致するビットを確定ビットとしてアウトプットします。このビットシフトによって、上限と下限の最上位ビットは常に異なっている状態になり、数直線の幅を確保することができます。

以下が符号化部分のコードです。

  void Encode (uchar bit) {
    ulong mp = CalculateMP();
    Encode(bit, mp);
  }
  void Encode (uchar bit, ulong mp) {
    if (bit == 1) {
      lower_bound_ = mp;
    } else {
      upper_bound_ = mp - 1;
    }
    
    // 上限と下限を比較して確定しているbitがあればencoded bitとしてバッファに格納する。
    // 確定しているbit分だけ、上限・下限をシフトする。
    while ((lower_bound_ & kMSBIT) == (upper_bound_ & kMSBIT)) {
      unsigned char outbit = lower_bound_ >> (kN-1);
      AddBuffer(outbit);
      lower_bound_ = (lower_bound_ << 1) & kMSMASK;
      upper_bound_ = ((upper_bound_ << 1) & kMSMASK) | 1;
    }
    
    // クローニングを行う。
    DoCloning(bit);
 
    // マルコフ連鎖を更新する。
    current_state_->SetTransCount(bit, current_state_->GetTransCount(bit)+1);
    current_state_ = current_state_->GetNextState(bit);
  }

  // 分離点(mp)を計算する。
  ulong CalculateMP () {
    double p0, p1;
    ulong mp;
    uint totalCount = current_state_->GetTotalCount();
    p0 = (double) (current_state_->GetTransCount(0) + 1) / (totalCount + 2);
    p1 = (double) (current_state_->GetTransCount(1) + 1) / (totalCount + 2);
    mp = (ulong) ((p1 * lower_bound_ + p0 * upper_bound_) / (p0 + p1));
    if (mp <= lower_bound_) {
      mp = lower_bound_ + 1;
    }
    // mpの最右bitを1に変更
    mp = mp | 1;
    // 最右bitを1に変更した結果、upper_bound_を超える可能性がある。その場合、mp=upper_bound_とする。
    if (mp > upper_bound_) {
      mp = upper_bound_;
    }
    return mp;
  }
Guazzoアルゴリズムの実践における課題

元論文ではGuazzoのアルゴリズムを実装するにあたり２つの問題があると述べています。

  • １つ目の問題は、整数値の数直線を上限と下限を分数で分割していくため、過程が進むにつれてより大きな精度が必要になる、という問題。
  • ２つ目の問題は、符号化するデータは有限であるという制約のもと、過不足のないデータを復元しなくてはいけないという問題。

まず一つ目の問題について、Guazzoは、正確に境界値を求める必要性を弱めることで解決しているそうです。上記のソースコードではCalculateMPの中で、確率からmp（数直線上の分岐点）を求める際にdouble型から整数値に丸めたり、mpの末尾bitを1に固定化したり（理由は後述）しています。このようにmpが正確な値ではなくても、圧縮時・復元時ともに同一の基準（論文ではfidelity criterion（忠実度基準）と呼んでいます）に従って計算が行われていれば、圧縮・復元に際しては大した問題にはなりません。

二つ目の問題については、かなり手ごわい問題で、たとえば1ビットが圧縮後に2ビットになった場合、圧縮後の2ビットがインプットされても元の2ビットを復元できる保障がない点に注意を払う必要があります。例えば、区間が[0x0000...0000, 0x1111...1111]として、入力ビット"0"を圧縮するとします。はじめの0で分岐点が"0x0010....0001"となり、区間が[0x0000...0000, 0x0010....0000]に縮小された場合、符号化された出力ビットは"00"になります。同様の状態でデコーダは分岐点の"0x0010....0001"以上かどうかで復元後のビットが0なのか1なのかを判断しますが、入力ビット"00"を復元しようとした場合、今後のビットの出現次第で、"0x0000....0000"や"0x0011....1111"のように分岐点以上にも以下にもなる可能性があります。そのため、この時点では復元することはできないのです。

すなわちaビットをbビットに符号化したとして、bビットあれば元のaビットを過不足なく復元できるとは限らないわけです。
この問題を解消するために、符号化時に二つの工夫をしています。

それは「mpの最右ビットを必ず1にしていること」と「符号化の最後にmpの最右ビットを除くN-1ビットを出力していること」です。これだけで過不足なく元データを復元することが出来るようになります。
理由を説明します。まずmpの最右ビットが1であるという条件から、数直線上の区間は必ず（ビットシフトの演算をのぞき）縮小されていきます。そのため、最後の入力ビットを符号化したときのmpのN-1ビット目までと最後に出力されるmpのN-1ビット目までは必ず異なるものになります。最後に出力されるmpのN-1ビットをインプットとすれば、最後の入力ビットを符号化したときのmp以上か未満かを必ず決定できる、つまり最後の入力ビットを決定できるのです。

また、最後のmpのN-1ビット目で止めている、というのは余分な1ビットが符号化されてしまうのを防止する意味もあります。最後に出力されるmpのNビット目を出力してしまうと、mp以上が確定してしまい、余分な"1"が符号化されてしまうことになります。しかし、N-1ビット目までの情報だと、デコーダは最後に出力されるmp以上か未満かを判別することはできません。つまり余分な1ビットは符号化されないのです。

Extra 7bits

さらに考慮すべき点があります。それはこのロジックでは符号化がビット単位になされていますが、実際にはデータはバイト単位だということです。符号化した後のビット数が8に達するたびにバイト単位で書き出しを行っていますが、最後に8に達しなかったあまりのビットは出力されることはありません。この状態のままでは完全な復元が出来なくなってしまいます。

そこで元論文で提案しているのは、符号化の最後に余分な7ビット(Extra 7 bits)を符号化することです。

  // encodedされた最終bitが、バイトとして出力されるようにダミーの7bitをencodeする。
  // また、encode対象の最終bitが正常にencodeされるように、mpの末尾1bitを除いた全てのbit列をencoded bitとして出力する。
  void EncodeFinish () {
    ulong mp;
    for (int i=0; i<7; i++) {
      mp = CalculateMP();
      // ダミーのbitとして、encoded bitが最低1bit確定されるbitを選択し、encodeする。
      if ((lower_bound_ & kMSBIT) == (mp & kMSBIT)) {
        // lower_bound_の最上位bitとmpの最上位bitが一致する場合は、ダミーのbitとして0をencodeする。
        Encode(0, mp);
      } else {
        // upper_bound_の最上位bitとmpの最上位bitが一致する場合は、ダミーのbitとして1をencodeする。
        Encode(1, mp);
      }
    }
    mp = CalculateMP();
    // mpの最後の1bitを除き、出力する。
    while (mp != kMSBIT) {
      uchar outbit = (mp >> (kN-1));
      AddBuffer(outbit);
      mp = (mp << 1) & kMSMASK;
    }
  }
};

ここで余分な7ビットを符号化しています。7ビットはなんでもよいかといえばそうではなく、それまでに符号化したビット全体がバイト列として出力されるためには、本体の符号化後ビット数が8N+1ビットだった場合を考慮して、最低でも7ビット分余分に出力を行う必要があります。そのため、0と1のうち、符号化した後に最低でも1ビット符号化される方のビットを選択して出力していきます。

これでバイト単位に丸めることによって、ビットが切り捨てられてしまったとしても、データ本体を符号化したビットは生き残るので、完全な形で復元することができます。

ベンチマーク

最後に、この圧縮手法で、圧縮率のベンチマークをとってみます。
テスト対象としては、定番所でThe Canterbury Corpusを利用しています。
まずはファイル種類別の圧縮率の相関図です。

それぞれのデータ特性についてあまり詳しく分析できていないのですが、テキストファイルではおおよそ35〜50%くらいの圧縮性能が実現できています。id:naoya氏が最近日記に書かれているPerl で Range Coder (再挑戦)のRange Coderでも、alice29.txtの圧縮率は50%をクリアできていなかったので、マルコフ連鎖を利用するだけで、大幅な性能改善を図ることが出来た、と理解しています。

次にCloningの閾値（上記ソースのcloning_threshold1=cloning_threshold2の値）を変化させていったときの圧縮率の変化です。他の圧縮アルゴリズムの比較として右端にZIPでの圧縮結果をプロットしています。

元論文では閾値が2のときがおおむね最も圧縮率が高い、と記述されていましたが、実際には上記グラフの通り、閾値が8〜32くらいのときに圧縮率がピークになるようです。直観的にも、閾値が低すぎると、各状態における状態遷移確率が偏る前に、すぐにクローニングしてしまうので、、圧縮率も低下するものと理解できます。

また、ZIPと比較すると、ptt5（FAXファイル）やEXCELファイル以外ではまだまだ圧縮率で負けていることがわかります。

まとめ

DMCを使って圧縮を行うと、エントロピーを最小にするように動的にマルコフ連鎖を最適化しながら圧縮を行うため、全体の出現頻度表をベースに圧縮するよりも効率のよい圧縮化を行うことができます。他のアルゴリズムとの比較など、考察が足りない部分は多々ありますが、それはおいおいということで。あらかじめ次数を決めてかかるPPMなどよりも、圧縮率はよいものと思われます。

参考文献

[1] Ian H. Witten, Alistair Moffat, Timothy C. Bell, Managing Gigabytes: Compressing and Indexing Documents and Images, Morgan Kaufmann Publisher, 1999
[2] G. Cormack and R. Horspool, "Data compression using dynamic Markov modelling," Computer J., 30 No. 6, 541-550, 1987.
[3] 高速な算術圧縮を実現する「Range Coder」(記事), 2006, CodeZine, 翔泳社