SpaceLesson

変調・符号化#172

圧縮センシング — ナイキストを出し抜くスパース信号復元

帯域幅の2倍以上のサンプリング周波数を要求するナイキスト定理に対し、信号が『スパース』であるという追加情報さえあれば、はるかに少ない観測数から元の信号を復元できる。L1ノルム最小化とRIP(制限等長性)を軸に、圧縮センシングの数学的な仕組みを理解する。

前提知識: 情報理論の基礎 — エントロピーと相互情報量で「情報」を測る

圧縮センシングスパース性L1最小化RIP

この回で学ぶこと

PCM/PSK/PMやSDRの回で、私たちはナイキスト・シャノンのサンプリング定理を当然の前提として扱ってきました。帯域幅 BB の信号を完全に復元するには、サンプリング周波数 fsf_s を

fs>2Bf_s > 2B

以上に取らなければならない。これは信号処理の教科書のもっとも基本的な結果であり、これより少ないサンプル数で元の信号を復元することは原理的に不可能だと教わります。実際、信号について何の追加情報も持っていなければその通りです。

しかしこの回で扱う圧縮センシング(compressed sensing, CS)は、「信号について何も知らない」という前提を少しだけ崩すと、この“不可能”が覆るという驚くべき結果を扱います。もし信号が、ある基底で見たときにほとんどの成分がゼロである(これをスパースであると言います)という追加情報を使えるなら、ナイキストレートよりもはるかに少ない数の観測値から、元の信号を完全に復元できる場合があるのです。2004〜2006年にEmmanuel Candès、Justin Romberg、Terence Tao、David Donohoらによって確立されたこの理論は、MRI(磁気共鳴画像)の高速撮像、単一画素カメラ、そして次回扱うスペクトラムセンシング(無線スペクトルの疎な占有状況の高速検出)まで、幅広い応用を生みました。

この回では、「スパース性とは何か」「なぜスパース性があると少ない観測数で復元できるのか」「なぜ本来解くべき組合せ最適化問題(L0ノルム最小化)ではなく、その凸緩和であるL1ノルム最小化を使うのか」を、数式と幾何学的直感の両面から追っていきます。

直感的な導入: なぜ「ほとんどゼロ」なら少ない観測で済むのか

具体例から始めましょう。長さ N=1000N = 1000 のデジタル信号 x∈RNx \in \mathbb{R}^{N} があり、そのうちたった K=5K = 5 個の成分だけが非ゼロで、残り995個はすべて厳密にゼロだと分かっているとします。

もしこの信号のどの5個が非ゼロなのかも分かっているなら、話は簡単です。その5個の値を知るのに必要な観測は、せいぜい5回程度で足ります。問題は、どの位置が非ゼロなのかも分からないという点です。これは「1000個の未知数のうち995個はゼロだが、どれがゼロかは分からない」という組合せ的な不確実性を含んでいます。

素朴に考えると、この不確実性を解消するには N=1000N=1000 個すべてのサンプルを取るしかなさそうに思えます。ところが圧縮センシングの結果は、うまく設計された観測を MM 回行えば(ただし MM は NN よりずっと小さく、KK よりは少し大きい程度でよい)、高い確率で xx を完全に、しかも計算量的に実行可能な方法で復元できることを示します。目安として、

M=O ⁣(Klog⁡NK)M = O\!\left(K \log \frac{N}{K}\right)

程度の観測数で十分であることが理論的に保証されます。K=5K=5, N=1000N=1000 の例なら、MM はおよそ数十程度で足りることが多く、N=1000N=1000 に比べて桁違いに少ない数です。これがナイキスト定理の「帯域幅で決まる最低サンプル数」という発想を覆す、圧縮センシングの核心です。

もちろんこれは魔法ではありません。「スパースである」という強い事前情報を使っているからこそ成立する話であり、その事前情報を復元アルゴリズムにどう組み込むかが、次節以降の数式の主題になります。

数式定式化1: スパース性の定義

信号 x∈RNx \in \mathbb{R}^N が、ある基底(あるいは辞書)行列 Ψ∈RN×N\Psi \in \mathbb{R}^{N \times N} を使って

x=Ψs,s∈RNx = \Psi s, \qquad s \in \mathbb{R}^N

と表せるとします。Ψ\Psi の列は正規直交基底(たとえばウェーブレット基底やDCT基底)を成すと考えてよく、ss は xx をその基底で展開した係数ベクトルです。xx 自身が時間領域で疎とは限らなくても、適切な基底 Ψ\Psi で見ると係数 ss の大部分がゼロ(または無視できるほど小さい)ということがよくあります。たとえば自然画像はピクセル領域では疎ではありませんが、ウェーブレット領域で見ると係数の大部分がほぼゼロに潰れます(これがJPEG2000のようなウェーブレット圧縮が効く理由でもあります)。

このとき、係数ベクトル ss がK-スパースであるとは、

∥s∥0≜#{i:si≠0}=K,K≪N\|s\|_0 \triangleq \#\{i : s_i \neq 0\} = K, \qquad K \ll N

を満たすことを言います。ここで ∥⋅∥0\|\cdot\|_0 は「L0ノルム」と呼ばれますが、実際には三角不等式を満たさないため数学的な意味でのノルムではなく、単に非ゼロ成分の個数を数える関数です。慣習的にこの呼び方が定着しています。

現実の信号は厳密に KK 個以外がゼロということは稀で、多くの成分が「小さいが完全にはゼロでない」という近似的スパース性を持ちます。しかし議論の骨格を掴むために、この回ではまず厳密なK-スパース性を仮定して理論を組み立て、実務のところで近似的スパース性への拡張に触れます。

数式定式化2: 劣決定な観測系とL0ノルム最小化

次に、信号 xx(あるいは基底変換前の x=Ψsx=\Psi s)を、観測行列 Φ∈RM×N\Phi \in \mathbb{R}^{M \times N} を使って低次元に射影した観測ベクトル

y=Φx∈RM,M<Ny = \Phi x \in \mathbb{R}^{M}, \qquad M < N

を考えます。Φ\Phi の各行は、xx の全成分の重み付き和(内積)を1回測る「観測」を表していると考えてください。ナイキストのサンプリング(等間隔の時間サンプル)は Φ\Phi の特殊な一形態にすぎず、圧縮センシングでは Φ\Phi の行として、たとえばランダムなガウス乱数ベクトルや、ランダムに間引いたフーリエ基底ベクトルを使います。

問題は、y=Φxy=\Phi x という線形方程式系が M<NM < N のとき**劣決定(underdetermined)**であることです。線形代数の基本的な事実として、この方程式を満たす xx は(ある1つの解が存在すれば)無限に存在し、解空間は次元 N−MN-M のアフィン部分空間になります。追加の制約なしにはどの解が「正しい」元の信号かを一意に決めることはできません。

ここでスパース性の仮定を使います。もし元の信号 xx(またはその基底表現 ss)がK-スパースだと分かっているなら、無限にある解の中からもっとも非ゼロ成分の少ない解を探せばよいはずです。これは次の最適化問題として定式化されます。

x^=arg⁡min⁡x∥x∥0s.t.y=Φx\hat{x} = \arg\min_{x} \|x\|_0 \quad \text{s.t.} \quad y = \Phi x

これは理論上もっとも自然な定式化ですが、致命的な欠点があります。この問題を厳密に解くには、原理的には「NN個の成分から k=0,1,2,…k=0,1,2,\dots 個を選んで非ゼロにする、すべての組み合わせを試す」という組合せ探索が必要で、その計算量は NN に対して指数関数的に増大します。これは計算複雑性理論の言葉で言えば NP困難 な問題であり、NN が数百〜数千を超える実用的な信号サイズでは、厳密解を求めることは事実上不可能です。

数式定式化3: L1ノルム最小化への凸緩和

そこでL0ノルムの代わりに、計算可能な凸関数であるL1ノルム

∥x∥1=∑i=1N∣xi∣\|x\|_1 = \sum_{i=1}^{N} |x_i|

を使った、次の最適化問題(basis pursuitと呼ばれます)を考えます。

x^=arg⁡min⁡x∥x∥1s.t.y=Φx\hat{x} = \arg\min_{x} \|x\|_1 \quad \text{s.t.} \quad y = \Phi x

L1ノルムは絶対値の和なので凸関数であり、この問題は等式制約付きの凸最適化問題——実際には線形計画法(LP)に帰着できる問題——です。標準的な数値解法(単体法や内点法)を使えば、NN が数万〜数百万でも実用的な時間で解けます。

圧縮センシング理論の中心的な結果は、適切な条件下では、このL1ノルム最小化の解 x^\hat{x} が、本来解きたかったL0ノルム最小化の解と完全に一致する、というものです。組合せ的に手に負えない問題が、凸最適化という「解ける」問題に置き換えられるにもかかわらず、同じ解が得られる——これが圧縮センシングが実用上の理論として成立する理由です。

観測に雑音 nn が乗る現実的な場合(y=Φx+ny = \Phi x + n、∥n∥2≤ε\|n\|_2 \le \varepsilon)には、等式制約を緩めた

x^=arg⁡min⁡x∥x∥1s.t.∥y−Φx∥2≤ε\hat{x} = \arg\min_{x} \|x\|_1 \quad \text{s.t.} \quad \|y - \Phi x\|_2 \le \varepsilon

(basis pursuit denoising、統計学で言うLASSOと本質的に同じ形)を使います。この場合も、後述のRIP条件のもとで復元誤差が観測雑音のレベル ε\varepsilon に比例して抑えられることが保証されます。

数式定式化4: RIP — L1緩和が正しく機能するための十分条件

L1ノルム最小化がL0ノルム最小化と同じ解を与えるための十分条件として、Candès と Tao が導入したのが**制限等長性(Restricted Isometry Property, RIP)**です。観測行列 Φ\Phi が次数 KK のRIPを、定数 δK∈(0,1)\delta_K \in (0,1) について満たすとは、任意のK-スパースなベクトル xx に対して

(1−δK) ∥x∥22  ≤  ∥Φx∥22  ≤  (1+δK) ∥x∥22(1-\delta_K)\,\|x\|_2^2 \;\le\; \|\Phi x\|_2^2 \;\le\; (1+\delta_K)\,\|x\|_2^2

が成り立つことを言います。この不等式は、「Φ\Phi による射影が、K-スパースなベクトルどうしのユークリッド距離(ノルム)をほぼ保存する」ことを意味します。もし δK\delta_K が小さければ、Φ\Phi はK-スパースな部分空間の中では、ほとんど等長写像(isometry、距離を変えない写像)のように振る舞うということです。これがあれば、異なる2つのK-スパース信号が同じ観測 yy を生む(つまり Φ(x1−x2)=0\Phi(x_1-x_2)=0 となり、区別がつかなくなる)ことを避けられます。

Candès(2008年)は、2K2K 次のRIP定数が

δ2K<2−1≈0.414\delta_{2K} < \sqrt{2} - 1 \approx 0.414

を満たせば、L1ノルム最小化がK-スパース信号を(雑音がなければ)厳密に、(雑音があっても)安定に復元できることを証明しました。これはあくまで十分条件であり、実務ではこれより緩い条件でも復元が成功することがほとんどですが、理論保証の基準として広く引用されます。

重要な実用上の事実は、成分が独立同分布のガウス乱数、あるいは ±1\pm 1 の等確率ランダム変数であるような行列 Φ\Phi は、高い確率でRIPを満たすということです。しかも必要な観測数はおおよそ

M≥C⋅Klog⁡ ⁣(NK)M \ge C \cdot K \log\!\left(\frac{N}{K}\right)

(ある定数 CC)で済むことが確率論的な議論(ジョンソン・リンデンシュトラウスの補題に関連する集中不等式)から導かれます。「ランダムに観測すれば(狙って設計しなくても)高確率でうまくいく」という点が、圧縮センシングが実装しやすい大きな理由の1つです。

幾何学的直感: なぜL1ノルムはスパース解を選ぶのか

なぜL0ではなくL1を使うと、しかも同じスパース解が出てくるのか。これを2次元の簡単な図で直感的に理解しましょう(実際には NN 次元ですが、本質は2〜3次元で見て取れます)。

N=2N=2、観測は1本の線形制約 y=Φxy = \Phi x(M=1M=1)だとします。このとき解の候補は、2次元平面上の1本の直線(アフィン部分空間)上に並びます。この直線上で ∥x∥1\|x\|_1 が最小になる点を探す、というのが最適化問題の幾何学的な意味です。

∥x∥1≤r\|x\|_1 \le r が定める領域(L1球)は、2次元では原点を中心とした菱形(正方形を45度回転させた形)になります。原点からこの菱形を徐々に膨らませていき(半径 rr を大きくしていき)、直線に初めて接触する瞬間を考えます。菱形には角があり、その角はちょうど座標軸上——つまり x1=0x_1=0 または x2=0x_2=0 というスパースな点——に位置しています。直線の傾きがよほど菱形の辺と平行でない限り、直線は菱形の角で接触する可能性が非常に高く、その接点はスパースな解になります。

一方、もし同じことをL2ノルム(∥x∥2≤r\|x\|_2 \le r、原点中心の円)でやったらどうなるでしょうか。円には角がなく、なめらかな曲線しかありません。直線が円に接する点は、よほど特殊な位置関係でない限り、x1,x2x_1, x_2 の両方が非ゼロの一般的な点になります。つまりL2ノルム最小化は、スパースな解を選ぶ理由が何もないのです(これが従来の最小二乗法や逆フィルタが、劣決定系に対してスパースでない「なめらかにばらけた」解を返す理由でもあります)。

この直感は次元が上がっても本質的に変わりません。L1球(超立方体の対角切断のような多面体)は、座標軸に近い低次元の面・稜・頂点を持ち、それらはすべてスパースな点に対応します。アフィン部分空間(観測制約が定める解空間)がこれらの低次元な面のどこかで接触する確率は、球のように滑らかな凸体に接触する場合よりもはるかに高い——これがL1ノルム最小化がスパース解を「誘導する」幾何学的なからくりです。

実務での使われ方

  • スペクトラムセンシング(次回予告): 広い周波数帯域を常時監視する無線受信機にとって、実際に電波が使われている(占有されている)周波数帯はごく一部にすぎず、周波数領域で見ればスパースです。次回学ぶコグニティブ無線では、この性質を使って、ナイキストレートでのフルバンドサンプリングを行わずに、圧縮センシングによって疎な周波数占有状況を高速に検出する研究が進んでいます。IEEE 802.22やDARPAのRF関連プログラムなど、スペクトル利用効率化の文脈で活発に研究されている分野です。
  • MRI(磁気共鳴画像)の高速撮像: 2007年のLustigらの研究(Sparse MRI)は、医用画像がウェーブレット領域や有限差分領域でスパースであることを利用し、通常のナイキスト条件を満たさない疎なk空間(周波数領域)サンプリングから画像を再構成する手法を確立しました。これはスキャン時間の大幅な短縮(小児や心臓の撮像で息止め時間を短くできるなど)に直結し、圧縮センシングがもっとも成功した実用分野の1つとされています。
  • 電波天文学の画像再構成: VLA(Very Large Array)やALMA、事象の地平線望遠鏡(EHT)のような電波干渉計は、地球上に離散配置されたアンテナ群でしか天球面の空間周波数(視野のフーリエ成分)をサンプリングできず、これは本質的に劣決定な逆問題です。古典的なCLEAN法(1974年、Högbomによる)は、天体画像が輝点の集まりとしてスパースであるという仮定を利用する点で圧縮センシングと発想を共有しており、近年ではL1正則化を明示的に使った再構成アルゴリズムも実装されています(EHTによるブラックホールシャドウ画像の再構成でも、スパース性を利用した正則化手法が比較検討されました)。
  • 深宇宙探査機での将来的な応用可能性: 探査機側の観測データ(スペクトルデータ、特定波長のみに情報が集中する分光データなど)が疎な表現を持つ場合、探査機の限られた電力・帯域幅の制約下で、圧縮センシング的な「疎な観測を行い地上で復元する」というアプローチはデータ圧縮・伝送の効率化につながる可能性があります。ただし現時点(2020年代)では、CCSDS標準の画像圧縮(CCSDS 122.0-B、ウェーブレット変換ベース)のような確立された圧縮方式が主流であり、圧縮センシングを深宇宙リンクの主要な符号化方式として採用したミッションはまだ限定的です。将来、センサ側の演算資源が限られる小型探査機・分散センサネットワークのようなミッション形態において、「圧縮してから送る」のではなく「最初から少ない観測しか取らない」という圧縮センシングの発想が活きる可能性が、研究レベルで議論されています。

演習問題

  1. N=100N=100 の信号のうち K=4K=4 個の成分だけが非ゼロだとします。目安の式 M=O(Klog⁡(N/K))M = O(K\log(N/K)) を使って、対数の底を自然対数として概算した場合の観測数 MM のオーダーを求め、これが N=100N=100 に比べてどの程度小さいか議論してください。
  2. L0ノルム最小化が「なぜNP困難なのか」を、N=20N=20、K=3K=3 の場合に「非ゼロ位置の組み合わせ数」を実際に計算する((203)\binom{20}{3})ことで具体的に確認してください。また NN が2倍、KK が同じままだったら、この組み合わせ数はどう変化しますか。
  3. RIP条件 (1−δK)∥x∥22≤∥Φx∥22≤(1+δK)∥x∥22(1-\delta_K)\|x\|_2^2 \le \|\Phi x\|_2^2 \le (1+\delta_K)\|x\|_2^2 において、もし δK=0\delta_K = 0 だったら Φ\Phi はK-スパースなベクトルに対してどのような性質を持つ写像になるか(ヒント: 等号成立時の意味を考えてください)。また δK\delta_K が1に近づくと、この不等式が保証する内容はどう弱くなっていくか説明してください。
  4. 本文の幾何学的直感の節を踏まえ、もしL1ノルムの代わりに ∥x∥p=(∑i∣xi∣p)1/p\|x\|_p = (\sum_i |x_i|^p)^{1/p}(0<p<10 < p < 1)という「劣加法的」なノルムもどきを使ったらどうなるか考えてみましょう(このような p<1p<1 の球は原点方向に凹んだ、より角の鋭い形になります)。スパース性を誘導する能力という観点で、p<1p<1 はL1(p=1p=1)よりも有利になりそうか、それでも実務でL1が好まれる理由は何だと考えられるか、自分の言葉で論じてください。

まとめと次回予告

圧縮センシングは、「信号がスパースである」という追加情報を使うことで、ナイキスト定理が要求するサンプル数よりもはるかに少ない観測数 M=O(Klog⁡(N/K))M = O(K\log(N/K)) から信号を復元できることを示す理論です。本来解きたい問題(L0ノルム最小化)は組合せ的でNP困難ですが、その凸緩和であるL1ノルム最小化(basis pursuit)が、RIPと呼ばれる十分条件のもとで同じ解を与えることが保証されます。この「なぜL1がスパース解を選ぶのか」という問いには、L1球の角がスパースな点に対応するという明快な幾何学的直感があります。

次回は、この圧縮センシングの発想を無線通信の現場に持ち込んだコグニティブ無線とスペクトラムセンシングを扱います。広大な周波数帯域のうち実際に使われているのはごく一部という「周波数占有のスパース性」を、今回学んだ理論がどう実際の高速スペクトル検出アルゴリズムに落とし込まれるのかを見ていきます。

参考文献

  • E. J. Candès, J. Romberg, T. Tao, “Robust Uncertainty Principles: Exact Signal Reconstruction from Highly Incomplete Frequency Information,” IEEE Transactions on Information Theory, vol. 52, no. 2, 2006
  • D. L. Donoho, “Compressed Sensing,” IEEE Transactions on Information Theory, vol. 52, no. 4, 2006
  • E. J. Candès, M. B. Wakin, “An Introduction To Compressive Sampling,” IEEE Signal Processing Magazine, vol. 25, no. 2, 2008
  • E. J. Candès, “The Restricted Isometry Property and Its Implications for Compressed Sensing,” Comptes Rendus Mathematique, vol. 346, 2008
  • M. Lustig, D. Donoho, J. M. Pauly, “Sparse MRI: The Application of Compressed Sensing for Rapid MR Imaging,” Magnetic Resonance in Medicine, vol. 58, no. 6, 2007
  • CCSDS 122.0-B, Image Data Compression