はじめての多重解像度解析

はじめに

離散ウェーブレット変換による多重解像度解析について興味があったのだが、教科書や解説を読んでも説明が一般的、抽象的過ぎてよくわからない。個人的に躓いたのは

1. スケーリング関数とウェーブレット関数の二種類が出て来るのはなぜだ？
2. 結局、基底を張ってるのはどっちだ？
3. 出て来るのはほとんどウェーブレット関数なのに、最後に一個だけスケーリング関数が残るのはなぜだ？

といった点。このあたり、わかってからテキストを見直すとすごく簡単なのだが、わかるまでにちょっとギャップがある。なので本稿では、多重解像度解析の一番簡単な例である一次元のハールウェーブレットによる変換を例に、そのあたりを大づかみに理解することを目的とする¹。

簡単な例

以下、簡単のため1次元の離散空間を考える。$N=2^m$個の格子点があり、それぞれに数値が与えられているとしよう。例えば、$m=3$、すなわち8点のデータとして、こんな数列が与えられたとする。

data = [5,3,7,5,5,3,1,5]

棒グラフっぽく図示するとこんな感じ。

これを、半分の解像度、すなわち$m=2$、$N=4$点のデータで近似することを考える。最も安直には、隣り合うデータの平均値を用いることだろう。こんな感じ。

これで、もともと8点のデータ[5,3,7,5,5,3,1,3]が、4点のデータ[4,6,4,2]で近似された。

さて、8点を4点で近似したのであるから当然情報は失われているのであるが、ここでもともとの解像度で隣り合っていたデータの差分値を覚えておくことにする。

例えば、データの一つ目と2つめの「5,3」というデータに対して、平均値「4」と、右から見て左の値の差分の半分である「1」を覚えておく。

すると、平均値と差分の情報があるから、「4,1」という情報から「5,3」という元のデータを復元できることになる。そこで、隣り合う数字のうち、左側に平均値、右側に「右からみた左」の差分の半分の値を記録する。そうするとこうなる。

水色の部分が平均値、赤色の部分が差分値となる。右の図から元のグラフ(左)を再現できるため、情報は全く失われていない。

さて、元々8点のデータを、4点の近似データと差分に置き換えた。新たにできた4点の近似データを、さらに2点で近似することにする。

例えば、隣あう水色のデータ「4,6」を、平均値「5」と、差分値「-1」で置き換える。そして、左側の水色のデータに平均を、右側の水色のデータに差分値を保存する。すると、水色のデータがさらに半分に減る。

最後に残った2つの水色の棒も、その平均と差分に置き換えてしまおう。

最終的に、水色のデータ(全データの平均値)が一つと、差分データ7つの情報に変換された。

これは、最も簡単なウェーブレットであるハールウェーブレット(Haar Wavelet)による多重解像度解析に他ならない。各プロセスで情報は全く失われていないから、各工程は逆にたどることができ、最後から最初まで戻ることができる。すなわち、この変換は可逆変換である。

スケーリング関数とウェーブレット

$N=2^m$の一次元格子点を考える。格子点の数だけ自由度があれば、当然のことながらその格子点で表現できる全ての情報を保持することができる。例えば$m=3$の時、空間は8次元空間となるが$(1,0,0,0,0,0,0,0)$、$(0,1,0,0,0,0,0,0,0,0)$・・・$(0,0,0,0,0,0,0,1)$の8つの基底があればこの空間の情報全てを表現できる。いま、$2^m$個の格子をレベル$m$の格子と呼ぼう。この空間を張る基底をレベル$m$のスケーリング関数と呼ぼう。レベル$m$の、$n$番目の基底を$\phi_{m,n}$と表現する。例えばレベル$3$のスケーリング関数は8個ある。

さて、いまレベル$m$の基底で張られた空間でデータが表現されていたとする。このデータをレベル$m-1$のスケーリング関数で表現することを考える。

これは、例えば

(1,0,0,0,0,0,0,0)\\
(0,1,0,0,0,0,0,0)\\
(0,0,1,0,0,0,0,0)\\
(0,0,0,1,0,0,0,0)\\
(0,0,0,0,1,0,0,0)\\
(0,0,0,0,0,1,0,0)\\
(0,0,0,0,0,0,1,0)\\
(0,0,0,0,0,0,0,1)

の8本のベクトルで張られていた空間を、

(1,1,0,0,0,0,0,0)\\
(0,0,1,1,0,0,0,0)\\
(0,0,0,0,1,1,0,0)\\
(0,0,0,0,0,0,1,1)

の4本のベクトルで表現することに対応する²。この4本のベクトルが、レベル2のスケーリング関数$\phi_{2,0}, \phi_{2,1}, \phi_{2,2}, \phi_{2,3}$に対応する。

当然のことながら、8次元の空間を4次元で表現しようとしているのだから、情報が落ちてしまう。その落ちた情報を拾うのが以下のウェーブレット関数である。

(1,-1,0,0,0,0,0,0)\\
(0,0,1,-1,0,0,0,0)\\
(0,0,0,0,1,-1,0,0)\\
(0,0,0,0,0,0,1,-1)

これがレベル$2$のウェーブレット関数$\psi_{2,0}, \psi_{2,1}, \psi_{2,2}, \psi_{2,3}$である。

要するに、Haar基底によるウェーブレット変換とは、もともと

(1,0,0,0,0,0,0,0)\\
(0,1,0,0,0,0,0,0)\\
(0,0,1,0,0,0,0,0)\\
(0,0,0,1,0,0,0,0)\\
(0,0,0,0,1,0,0,0)\\
(0,0,0,0,0,1,0,0)\\
(0,0,0,0,0,0,1,0)\\
(0,0,0,0,0,0,0,1)

という基底で張られていた情報を、

(1,1,0,0,0,0,0,0)\\
(0,0,1,1,0,0,0,0)\\
(0,0,0,0,1,1,0,0)\\
(0,0,0,0,0,0,1,1)\\
(1,-1,0,0,0,0,0,0)\\
(0,0,1,-1,0,0,0,0)\\
(0,0,0,0,1,-1,0,0)\\
(0,0,0,0,0,0,1,-1)

と、基底の取り直しをしたに過ぎない。つまり、

1. レベル$m$のスケーリング関数$\{\phi_{m,n}\}$は、レベル$m$の格子上の情報を完全に表現できる(基底を張る)、
2. レベル$m-1$のスケーリング関数でレベル$m$の情報を表現しようとすると、情報が失われる。
3. その失われた情報を補完するのがレベル$m-1$のウェーブレット関数$\{\psi_{m-1,n}\}$である。

という関係にある。

さて、レベル$m$のスケーリング関数$\phi_{m,n}$は$2^m$個あるが、それを$2^{m-1}$個のスケーリング関数$\phi_{m-1,n}$とウェーブレット関数$\psi_{m-1,n}$に基底を取り直す。残った$2^{m-1}$個のスケーリング関数を、さらに$m-2$レベルのスケーリング関数とウェーブレット関数で表現し・・・ということを再帰的に繰り返すと、最終的にレベル0のスケーリング関数$\phi_{0,0}$と、様々なレベルのウェーブレット関数$\psi_{m,n}$に分解される。

これが、多重解像度解析で「ほとんどウェーブレット関数で展開されるが、最後に一つだけスケーリング関数が残る」理由である。

プログラム例

以下は、与えられたレベルmの一次元配列dataをハールウェーブレットで変換するRubyスクリプトである。

def transform(data,m)
  m.times do |l|
    (2**(m-l-1)).times do |i|
      i1 = i*(2**(l+1))
      i2 = i1 + 2**l
      s = (data[i1] + data[i2])*0.5
      d = (data[i1] - data[i2])*0.5
      data[i1] = s
      data[i2] = d
    end
  end
end

単純に、隣り合うデータの平均値を左に、差分を右に保存する処理を再帰的に行っている³。

元データとして、レベル8(つまり256点)の、こんな$\tanh$を食わせて見る。

M = 8
N = 2 ** M
data = Array.new(N) do |i|
  Math::tanh((i.to_f-N.to_f/2.0)/(N.to_f*0.1))
end

これをウェーブレット変換したデータはこうなる。

これのデータを、逆変換するのは簡単。隣り合うデータに対して、差分を足したものを左に、引いたものを右に入れれば良い。

def inv_transform(data,m)
  m.times do |l2|
    l = m - l2 -1
    (2**(m-l-1)).times do |i|
      i1 = i*(2**(l+1))
      i2 = i1 + 2**l
      s = (data[i1] + data[i2])
      d = (data[i1] - data[i2])
      data[i1] = s
      data[i2] = d
    end
  end
end

先程のデータを逆変換すると元に戻る。

データ圧縮

ウェーブレット変換は、$N$個のデータを$N$個の異なるデータに変換するもので、この変換では情報は落ちていないから可逆変換である。しかし、せっかくウェーブレット変換したので、データを圧縮することを考えよう。

まず、先程の変換では平均と差分を保存していた変換に$\sqrt{2}$をかけることにする。それに対応して、逆変換は$\sqrt{2}$で割らなければならない。

def transform(data,m)
  m.times do |l|
    (2**(m-l-1)).times do |i|
      i1 = i*(2**(l+1))
      i2 = i1 + 2**l
      s = (data[i1] + data[i2])/Math.sqrt(2.0)
      d = (data[i1] - data[i2])/Math.sqrt(2.0)
      data[i1] = s
      data[i2] = d
    end
  end
end

def inv_transform(data,m)
  m.times do |l2|
    l = m - l2 -1
    (2**(m-l-1)).times do |i|
      i1 = i*(2**(l+1))
      i2 = i1 + 2**l
      s = (data[i1] + data[i2])/Math.sqrt(2.0)
      d = (data[i1] - data[i2])/Math.sqrt(2.0)
      data[i1] = s
      data[i2] = d
    end
  end
end

この状態で、ウェーブレットの自乗重みについて「上位30%まで」残し、残りは0としてしまおう⁴。

transform(data,M)

data2 = data.map{|x| x**2}.sort.reverse
th = data2[N*0.3] # 自乗重みの上位30%をスレッショルドに設定

data.map!{|x| x**2 < th ? 0 : x } # スレッショルド以下のウェーブレットについては重みをゼロに

inv_transform(data,M)

すると、差分が小さいところは無視される。

ややデータがガタガタしたのがわかるかと思う。上位10%だけ残すことにすると、よりガタガタする。

これは、基底の数を固定した状態で、もともとのデータとの平均二乗誤差をもっとも小さくするように残す基底を選んだことに対応する⁵。ハールウェーブレットでは、平なところがあるデータほど、効率よく圧縮できる。

まとめ

もっとも簡単なウェーブレットであるハールウェーブレットによる変換、逆変換、圧縮を実装してみた。人によって躓きポイントは異なるだろうが、僕の場合は、

1. 同じ解像度の情報を別の基底で表現すること(レベル$m$の情報は、レベル$m$のスケーリング関数か、レベル$m-1$のスケーリング関数とウェーブレットで二通りの表現ができる)
2. スケーリング関数は再帰的に「より粗い」スケーリング関数とウェーブレットに展開され、最後に「一番粗い」スケーリング関数が残る

というところを理解するのに時間がかかった。

ウェーブレット変換(というか多重解像度解析)の例として、絵が4分割された図(なんか左上に縮小された図があって、残りの3パネルに黒っぽいデータがあるやつ)がよく出てくるが、これも1次元のデータのうち、平均値を左に寄せ、差分を右に寄せた図に対応すると思うと理解しやすいと思う。

ウェーブレットについて本格的に勉強したければ、例えば教科書としてドブシーの「ウェーブレット10講」なんかを読むと良いんじゃないでしょうか。

本稿で利用したコードを

に置いておく。

1. 「はじめての多重解像度解析」と書いたが、これは「まったくはじめての人にもわかるように書いた」という意味ではなく、単に僕が「はじめて多重解像度解析してみたよ」という意味。 ↩

2. 面倒なのでここでは規格化定数は無視した。 ↩

3. ここで、展開した基底を正規直交基底だと思うためには係数に$\sqrt{2}$とか出て来るが、フーリエ変換の$2\pi$と同様に、逆変換と辻褄があっていればどうでも良い。 ↩

4. ちなみに「どの場所のどのレベルの重みであったか」も合わせて記録する必要があるから、30%しか重みを残さなくても、データが30%に圧縮されるわけではない。といっても、場所によってレベルも自動的に決まるから、「(場所, 重み)」のペアを保存すればよく、前者を32bit整数、後者を倍精度実数で保存することにすれば、上位30%だけ残したデータは、もとのデータに比べて45%の圧縮になるはず(たぶん)。 ↩

5. 気持ちとしては、2倍粗いウェーブレットが、自乗した時に2倍の重みを持つように$\sqrt{2}$をかけている。これが自乗重みが小さい奴から順番に切った時に、元の関数からの平均自乗誤差を最も小さくする形になっていると思うが、あんまり自信がない。間違ってたらごめんなさい。 ↩