はじめに

大学の授業で「実際のデータを使ってスパースモデリングを行いましょう！」という課題が出て、スパースモデリングを使った画像圧縮を考案し()実装したので、こちらにまとめておきます。

概要と結果

要は、JPEGの理論にL1正則化を組み込んで画像圧縮を行ったわけです。考えたときには天才か!?と思いましたが、調べたら研究しつくされていました。ですがもう実装も終えて発表資料もつくっていたので知らないふりして発表しました(「L1正則化画像圧縮」「圧縮センシング」などで検索してください)。

最初に結果だけ示すと、こんな感じです。

圧縮率: 原画像 \lt L1(\lambda =0.01) \lt L1(\lambda =0.1) \lt \lt L1(\lambda =1) \simeq JPEG \lt L1(\lambda =5)

これだけで今回の記事の結果が理解る方もいらっしゃるかもしれませんが、わからなくてもこの図は後で載せる予定ですので根気強く下の記事を見てください。

とにかく、原画像と、それをJPEGで圧縮した画像、JPEGにL1正則化を組み込んだ方法で圧縮した画像です。L1正則化ではパラメータとして $\lambda$ があり、今回の場合これが大きいほど強いスパース性を要求しています。スパース性についてまだ何も触れていないので、 $\lambda$ が大きいと圧縮率が高くなるとお考えください。

さて、どうですか？
L1正則化を用いた圧縮方法は優れているでしょうか？

画像圧縮では２つの観点が重要です。
それは、圧縮率の高さと劣化の少なさです。
上の結果を見てみると、

$\lambda$ が0.01、0.1のとき圧縮率はJPEGより小さいですが劣化も少なく、原画像と同じように見えます。
$\lambda$ が5のとき、圧縮率は最も高いですが、劣化が大きく実用的であるとは言えません。
$\lambda$ が1のとき、圧縮率はJPEGと大体同じくらいになりました。
JPEGの場合、画像にノイズがのっているように見えますが、反対にL1正則化では全体がぼやけたようになっています。

この違いは、各々のアルゴリズムに起因するものですが、どちらがいいと一概に言えるものではないでしょう。この女性の画像に関しても、JPEGのほうが好まれる場合もあるでしょうし、L1正則化のほうが好まれる場合もあるでしょう。

ちなみに、各画像の圧縮率は以下の通りです。

現画像	JPEG	L1( $\lambda$ =0.01)	L1( $\lambda$ =0.1)	L1( $\lambda$ =1)	L1( $\lambda$ =5)
0%	89%	21%	43%	87%	96%
圧縮率の定義を知らないので、今回は

圧縮率 = 100 \times ( 1 - \frac{復元に使用しているデータの個数}{全てのデータの個数} )

という計算式です。すべてのデータの個数は、例えば $256 \times 256$ の画像であれば $256 \times 256 = 65536$ 個です。

$\lambda=5$ のときの圧縮率は96%であり、現画像では65536個のデータで表現していた画像をわずか2396個のデータで表現しています。そう思うとよく表現できているほうなのですかね。。。
※ $\lambda$ の値と劣化の関係は決して一般的なものでなく、画像によって異なります。

スパースモデリングとは

※L1正則化などの知識が既にある人は飛ばしてもらって構いません。

まずはじめに、スパースモデリングとはなにか。

最も大事なことはスパース性です。スパースとは日本語では疎であること、などといいますがデータ分析などの意味合いでは0が多いデータといってもいいかもしれません。

ここで大事なのは0である、ということはそのデータが持っている情報は意味がないもしくは本質的でないということです。

スパースモデリングとは、データはスパースな性質をもっていると仮定し、その上で本質的な情報のみでデータを表現する、ということです。
これはいわゆるオッカムの剃刀や、モーガンの公準という考え方が根底にあります。

機械学習と比較してみると、機械学習は多くのデータからどのようにして結果を表現するかに注目するのに対し、スパースモデリングでは得られた結果からどのデータが要因だったのかを表現しようとするところが大きな違いです。

※オッカムの剃刀の話と過剰適合・過学習を同一視してる意見がありますが、過学習は学習の結果汎化性能が下がってしまうということであり、オッカムの剃刀とは若干意味が異なります。結果として同一のことを主張することもありますが、注意してください。

例：L1正則化

ここで、スパースモデリングの手法の例として線形回帰のL1正則化を挙げます。別名LASSOとも呼ばれています。
詳しく知りたければ他のサイトを閲覧することを強くおすすめしますが、ここでも簡単に触れます。

実数値入力 $x \in R$ から実数値出力 $y \in R$ を結びつける関数

\begin{aligned} y &= \bm{w}^T \bm{\phi} ( x ) + \epsilon \\ &= \sum_{m=1}^M w_m \phi_m (x) + \epsilon \end{aligned}

を考えます。 $\phi (x)$ は特徴量関数と呼ばれ、例えば多項式回帰では

\bm{\phi}(x) = ( \phi_{M-1}(x),\phi_{M-2}(x),...,\phi_{1}(x),\phi_{0}(x))^T = ( x^{M-1},x^{M-2},...,x,1)^T

のようになります。

※多項式回帰では特徴量に $x^n$ が含まれるため非線形のように思えますが、 $\phi_i(x)$ の線形結合により予測を行っているため線形回帰といいます。単回帰や重回帰のみを線形結合としているサイトをたまに見ますが、単回帰は線形回帰の一種と私は認識しています。

ここで、通常の回帰ではある誤差関数 $E(\bm{w})$ に対してその誤差が最も小さくなるようなパラメータ ${\bm{w}_{LS}}$ を最適解とします。代表的な誤差関数として以下のようなものがあります。

E(\bm{w}) = \frac{1}{2}\sum_{n=1}^N\{y-\bm{w}^T\bm{\phi}(x_n)\}^2 \\ \bm{w}_{LS} = \argmin_{\bm{w}}E(\bm{w})

ですが、L1正則化ではこの誤差関数に次のようなペナルティ項を導入し、新たなコスト関数を定義した上でそれが最小となるパラメータを求めます。

\Omega_{L1}(\bm{w}) = \sum_{m=0}^{M-1}|w_m| \\ J(\bm{w}) = E(\bm{w}) + \lambda\Omega_{L1}(\bm{w}) \\(0 < \lambda < 1)

パラメータの表現力が高いほど $E(\bm{w})$ は小さくなりますが $\Omega_{L1}(\bm{w})$ は大きくなってしまうため、この新たな誤差関数を設定することでパラメータ $w_m$ に極端な値が出ることを抑制することになります。

またL1正則化では、得られたパラメータが疎になりやすいことがわかっています(非常に重要なところですが、これもわかりやすい説明が多く存在しているため詳細はそちらを御覧ください)。

画像のスパース性(JPEG)

※画像処理などを学ばれた方は既習の内容かもしれません

コンピュータで画像を表現するとき、一般的に広く使われているのは

PNG
JPEG
GIF
などだと思いますが、この3つの手法の中でJPEGだけが持つ特徴として、不可逆性が挙げられます。またJPEGは高い圧縮性を誇りますが、そこには不可逆性、そしてスパース性と関係してくるのです。

JPEGによる画像の表現方法

JPEGは、以下の流れで画像の情報を変換しています。

入力(出力)
↕
離散コサイン変換
↕
量子化
↕
符号化

太字の各ステップについて一行で説明します。

離散コサイン変換：画素値を周波数成分に変換
量子化：高周波数成分を0にする

厳密には色々と異なるのですが、だいたいこのような感じです。
符号化については、JPEGではハフマン符号を用いているようですがそれはここではあまり関係ないので触れません。

そしてJPEGの高い圧縮性というのは、量子化のステップでデータを高周波成分を0にすることに起因します。(必ず全ての高周波成分が0になるわけではないです)
「高周波成分を0にしていいのか？」と思われるかもしれませんが、人間の目では高周波成分というのは一般的な画像の認識において高い重要度はありません。

※高周波成分は画素値の差が大きい箇所、色が移り変わるところに該当するため目的によってはJPEGを採用するのは悪手です。

実際に、低周波成分のみの画像データと高周波成分のみの画像を比較してみると以下のようになります。

離散コサイン変換

JPEGでは、画像の周波数成分への分解は離散コサイン変換を利用します。
離散フーリエ変換でもいいじゃないか！と思われそうですがJPEGでは離散コサイン変換を使用します。一番の理由は、圧縮効果を高めることです(エントロピーが小さくなる)。詳しい説明は専門書などでを閲覧していただきたいですが、離散コサイン変換のほうが離散フーリエ変換より効率がいい、ということです。
使用する離散コサイン変換、逆変換は以下の式で表されます。
JPEGでは、8*8の64ピクセルごとに変換を行うため、Nは8となります。

DCT: F(u,v) = C(u)C(v)\sum_{i=0}^{N-1}\sum_{j=0}^{N-1}f_{i,j}cos(\frac{(2i+1)u\pi}{2N})cos(\frac{(2j+1)v\pi}{2N})\\ IDCT: f(i,j) = \sum_{u=0}^{N-1}\sum_{v=0}^{N-1}C(u)C(v)F_{u,v}cos(\frac{(2i+1)u\pi}{2N})cos(\frac{(2j+1)v\pi}{2N})\\ \\ C(p) = \begin{cases} \sqrt{\frac{1}{N}} & (p = 0) \\ \sqrt{\frac{2}{N}} & (p ≠ 0) \end{cases}

$f_{i,j}$ が画素値、 $F_{u,v}$ が周波数成分を表しています。
一見複雑そうな式ですが、二次元方向、横と縦にそれぞれフーリエ級数展開をしているイメージです。

ここで、フーリエ級数を思い出してください。フーリエ級数では周期関数を周波数成分と三角関数との線形和で近似していましたが、今回の画像でも同等です。画像にも基底ベクトルのようなものが存在し、離散コサイン変換で得られた周波数成分との線形和で画像が表現可能となります。以下は授業の発表で使用したスライドからとってきたものです。

量子化

量子化は至ってシンプルです。量子化では量子化テーブルと呼ばれるものを予め用意しておきます。離散コサイン変換によって得られた周波数成分は一つ前の図のように8×8の表形式で与えられますが、それの各要素とその量子化テーブルの各要素との商を整数で丸めたものが量子化後の値です。
JPEGの一般的な規格では、量子化テーブルは以下で規定されているらしいです。

量子化テーブルの高周波成分に対応する値は大きくなっており、高周波成分が0になりやすくなっているのがわかるかと思います。

量子化テーブル

言葉だけでは説明が難しい(面倒くさい)ため、具体的に計算を行いましょう。
仮に画像が以下で与えられたとします。

原画像

DCT係数を計算します。

DCT係数

量子化後の係数は

量子化後のDCT係数
逆量子化は先程と全く逆の手順で、量子化テーブルの値と掛け合わせることになります。すると、DCT係数は以下にようになる。量子化により、高周波成分(右下に行くほど高周波)の多くが0になったことがわかります。

逆量子化後のDCT係数
逆変換を行います。

復元画像
下の画像の左側が元画像、右側がJPEGの表現によって得られた画像です。
若干の誤差はあれど、全体としては元画像をよく表現できています。

左: 原画像右: 復元画像

JPEGにL1正則化を組み込む

JPEGでは、高周波成分をある意味意図的にカットすることによって高い圧縮率を実現していました。
これは、スパースモデリングでの仮定であるデータにはスパース性が存在している、を満たしているとも言えます。

ここで、このアルゴリズムにスパースモデリングの考え方を組み込めないかと考えました。

そこで注目したのが、離散コサイン変換の定義式です。
前述したように、離散コサイン変換では以下が成り立ちます。

f: 原画像
IDCT: 逆離散コサイン変換
F: 周波数成分

f = IDCT ( F )

この式を変形し以下のうような形にすれば、L1正則化による画像の表現方法が獲得できるだろうと考えました。

\bm{f} = \psi \bm{F}\\ \bm{f}_{L1} = argmin_f ||\bm{f} - \psi \bm{F}||_{2}^2 + \lambda |\bm{F}|_1

ここで、 $\bm{f}_{L1}$ は何を意味するのでしょうか。JPEGで得られる表現とは何が異なるのでしょうか。

JPEGの本質は、高周波成分を0にすることです。JPEGの画像は一般的な認識においてはデメリットは少ないのですが、用途によってはJPEGは使えません。高周波成分がカットされた情報は色が移り変わる箇所などを表現するのが不得手なためです。しかし $\bm{f}_{L1}$ は人間の認識などは考慮せず、ただ本質的な情報が抽出されたものと期待されます。

$\phi$ はどのように表現されるでしょうか。
その前に、原画像 $\bm{f}$ と周波数成分 $\bm{F}$ ですが、これはいずれも $1\times 64$ の列ベクトルとなっています。ここでの原画像は8*8のサイズを想定しており、実際は画像を8*8の少画像に分解しそれぞれに処理を行います。 $\bm{f}$ の模式図です。

IDCTの定義は以下の式でした。

IDCT: f(i,j) = \sum_{u=0}^{N-1}\sum_{v=0}^{N-1}C(u)C(v)F_{u,v}cos(\frac{(2i+1)u\pi}{2N})cos(\frac{(2j+1)v\pi}{2N})\\ \\ C(p) = \begin{cases} \sqrt{\frac{1}{N}} & (p = 0) \\ \sqrt{\frac{2}{N}} & (p ≠ 0) \end{cases}

そして、以下の形にするための $\psi$ は、64 $\times$ 64の行列となりますが、要素の値は複雑になってしまい書き下すのが難しいため割愛させていただきます。

\left( \begin{array}{c} f_{0,0} \\ f_{0,1} \\ \vdots \\ f_{1,0} \\ f_{1,1}\\ \vdots \\ f_{N-1,N-1} \end{array} \right) = \psi \left( \begin{array}{c} F_{0,0} \\ F_{0,1} \\ \vdots \\ F_{1,0} \\ F_{1,1}\\ \vdots \\ F_{N-1,N-1} \end{array} \right)

さて、ここまで来たらL1正則化もしくはLASSOを実装するのみですが、長くなったので記事を分けます。

Discussion