› 圧縮センシング L1 再構成 シミュレーター 戻る
信号処理・スパース表現

圧縮センシング L1 再構成 シミュレーター

N 次元の中で非ゼロ要素が K 個しかない「スパース信号」を、Nyquist より遥かに少ない M 観測から L1 最小化で復元する圧縮センシング（Compressive Sensing）を体験するツールです。信号次元・観測数・観測行列・雑音・再構成アルゴリズムを変えると、Candes-Tao の理論条件・RIP・再構成誤差がリアルタイムで分かります。

パラメータ設定
信号次元 N
復元したい信号 x のサンプル数
スパース数 K
非ゼロ要素の個数（K-sparse の K）
観測数 M
圧縮された観測ベクトル y のサンプル数
観測行列 Φ
RIP を満たす行列を選ぶほど少ない M で再構成できる
雑音標準偏差 σ
観測ノイズ y = Φx + n（n 〜 N(0, σ²)）
再構成アルゴリズム
BP=厳密／IST=反復／OMP=グリーディ
計算結果
—
圧縮率 M/N (%)
—
理論最小 M (Candes-Tao)
—
情報比 M/M_req
—
RIP 定数 δ_2K
—
再構成誤差 ‖x̂−x‖₂
—
計算量 log10(ops)
—
Φx → y → x̂ パイプライン可視化

左：N 次元のスパース信号 x（K 個のスパイク）、中央：M×N 観測行列 Φ、右：M 次元の圧縮観測 y と L1 再構成結果 x̂ の比較。色の濃さは振幅、緑＝再構成成功、赤＝失敗。

元信号 vs 再構成信号
再構成成功率 vs 圧縮率 M/N
理論・主要公式

$$\min_x \|x\|_1 \text{ s.t. } \|y - \Phi x\|_2 \leq \epsilon,\quad M \geq C\,K\,\log(N/K)$$

L1 最小化問題（Basis Pursuit）と Candes-Tao の理論観測数。y=観測ベクトル(M×1)、Φ=観測行列(M×N)、x=スパース信号(N×1)、ε=雑音許容、定数 C ≈ 2〜5。

$$\text{RIP}: (1-\delta_{2K})\|x\|_2^2 \leq \|\Phi x\|_2^2 \leq (1+\delta_{2K})\|x\|_2^2$$

制限等長性（Restricted Isometry Property）。δ_2K < 0.4 程度なら L1 再構成が安定に成立。Gaussian/Bernoulli 行列は M ≥ 5K·log(N/K) で漸近的に満たす。

$$\|\hat{x} - x\|_2 \leq C_1 \sigma \sqrt{M}\,/\,\sqrt{M - C_2 K \log(N/K)}$$

雑音下の再構成誤差の上界。観測数 M が理論最小に近づくほど分母がゼロに近づき、誤差が爆発する。雑音増幅率はアルゴリズムにより BP < IST < OMP。

圧縮センシング (Compressive Sensing) — L1再構成

🙋
「圧縮センシング」って名前はよく聞くんですが、要は普通の「データ圧縮（ZIPみたいな）」と何が違うんですか？
🎓
いい質問だ。ZIP のような圧縮は、まず信号を 全部観測してから 余分を捨てる。これに対して圧縮センシング（CS）は、最初から少ない観測しかしない。N 個サンプルすべき信号を M（≪ N）回しか測らず、後で計算で N 次元を復元する。MRI なら撮影時間が 1/4、Single-Pixel Camera なら 1 ピクセルセンサで画像が撮れる。Candes・Tao・Donoho が 2006 年に「ある条件下では確率1で完全復元できる」と証明したのが革命的だったんだ。
🙋
えっ、観測してないところを「計算で復元」できるんですか？ 普通そんなの足りない方程式（M 本の式で N 未知数）でしょう…
🎓
そう、線形代数的には不定方程式（解が無限にある）。でも、信号が「K-sparse」、つまり N 個中ほとんどがゼロで非ゼロは K 個だけ、という事前情報を入れると一意に決まる。本ツールの「スパース数 K」がそれだね。あとは L1 ノルム最小化、つまり ‖x‖₁ = Σ|xᵢ| が最小になる解を探すだけで、奇跡的に正解が出る。L2 最小化（普通の最小二乗）だと全部の要素に少しずつ値が乗ってしまって失敗するけど、L1 は「角」がスパース解を引き当てるんだ。
🙋
じゃあ観測数 M はどこまで減らせるんですか？ 上の結果で「情報比 M/M_req」が出てますね。
🎓
Candes-Tao が示した魔法の式が M ≥ C·K·log(N/K)。N=500・K=20 なら M_req ≈ 129、N=10⁶・K=10⁴ でも M ≈ 50,000 で足りる。線形ではなく対数で効くのがミソだね。デフォルトの M=100 は M_req=129 を下回るので「情報比 0.78」で再構成保証ナシ判定が出ているはず。M スライダーを 130 以上に上げてみると、verdict が緑に変わるよ。
🙋
観測行列 Φ を Gaussian から Hadamard に変えると RIP δ が悪くなりました。これは何ですか？
🎓
RIP（Restricted Isometry Property、制限等長性）は「Φ を任意の 2K-sparse ベクトルに掛けてもエネルギーがほぼ保存される」という性質。δ_2K < 0.4 程度を満たせば L1 再構成が安定する。Gaussian / Bernoulli random は理論最良に近い RIP を満たすけれど、実応用では物理的に Φ を選べないことが多い。MRI なら DFT 行のランダムサブセット、Single-Pixel Camera なら Hadamard、というふうにハードで決まる。RIP が悪くなった分、観測数 M を上乗せして補うのが実務。
🙋
最後に、BP / IST / OMP の使い分けってどう考えればいいですか？ 計算量 log10(ops) で 1〜2 桁違うんですね。
🎓
Basis Pursuit は L1 最小化を線形計画／SOCP として厳密に解く。CVXPY や MOSEK で書けるけど N³ オーダーで、N=10⁴ あたりが現実的な上限。MRI のように N=512²=262,144 の世界では FISTA（IST の加速版）が定番で、N·M·iter で済む。OMP は K が小さく（数十）既知のときに最速、ただし誤って間違った位置を選ぶと取り返しがつかない弱点がある。研究ではまず BP で「ベスト精度」を測り、本番では FISTA で速度を稼ぐ、というのが王道だね。

よくある質問

違反していません。Nyquist 定理は「任意の帯域制限信号」を完全復元するための必要観測数を述べたものです。圧縮センシングは対象を「K-sparse な（あるいは可圧縮な）信号」というクラスに限定し、その事前情報を利用することで、観測数を M = O(K·log(N/K)) まで減らせることを示しました。Candes-Tao-Donoho の 2006 年の結果で、観測行列が RIP（制限等長性）を満たせば L1 最小化で確率1に近い精度で完全再構成できます。Nyquist が成り立たないのではなく、適用範囲が違うのです。
理論的にはガウシアン random やベルヌーイ ±1 行列が RIP を最良に近い定数で満たし、最少観測数で再構成できます。ただし実応用では物理的に実現可能な行列に制限されます。MRI なら DFT（離散フーリエ）の行をランダムサンプリング、Single-Pixel Camera ならアダマール / Bernoulli が自然な選択です。本ツールでは Gaussian → Bernoulli → DFT → Hadamard の順に RIP 定数 δ_2K が悪化する傾向を再現しています。実装はハードと相談しつつ、δ_2K < 0.4 を目安に設計します。
Basis Pursuit（BP）は L1 最小化を線形計画／二次錐計画として厳密に解くため、最小観測数で最良の精度が得られますが、N³ オーダーの計算量が必要です。IST／FISTA はソフト閾値処理を反復するだけで N·M·iter で済み、大規模問題（MRI 画像 256² や 512² 等）で実用的です。OMP（Orthogonal Matching Pursuit）はグリーディに非ゼロ位置を1つずつ選ぶ手法で、K が小さく既知なら高速ですが、間違った位置を選ぶと取り返しがつかない弱点があります。本ツールの「再構成アルゴリズム」を切り替えると、雑音増幅率と計算量の差が確認できます。
条件次第ですが、実用化されています。MRI の場合、画像は wavelet 領域でほぼスパース（K/N ≈ 5〜10%）なので、k 空間サンプルを 1/4 にしても診断品質の再構成ができることが Siemens／GE／Philips の臨床機で実証されています。Single-Pixel Camera は Rice 大の Baraniuk らが 2008 年に発表し、1 ピクセルセンサだけで N=64×64 画像を M=1638 観測（圧縮率 40%）で復元できました。地震波・無線通信・ゲノミクスでも普及しています。ただし K-sparse という前提が崩れる信号（例：白色雑音）には全く効きません。

実世界での応用

医用イメージング（MRI 高速化）：圧縮センシングの最も成功した実応用です。MRI 画像は wavelet 領域でほぼスパースであるため、k 空間サンプル（観測）を 1/2〜1/4 に間引いて撮影しても、L1 最小化（CS-MRI）で診断レベルの画像が復元できます。臨床機（Siemens Compressed Sensing GRASP、GE HyperSense、Philips Compressed SENSE）に搭載され、心臓 MRI の息止め時間短縮、小児 MRI の鎮静時間短縮に直接貢献しています。CT でも線量削減を目的に同様の研究が進行中です。

Single-Pixel Camera と特殊イメージング：Rice 大の Baraniuk らが 2008 年に発表した Single-Pixel Camera は、DMD（デジタルマイクロミラー）で Hadamard / Bernoulli パターンの光を選択し、1 個のフォトダイオードだけで画像を撮影します。CCD アレイが作りにくいテラヘルツ・赤外・X 線領域で威力を発揮し、近赤外イメージングや量子イメージングで実用化が進んでいます。

無線通信・レーダー・地震信号：5G ミリ波の sparse channel estimation、cognitive radio の広帯域スペクトル監視、地中レーダーや地震波のスパース反射係数推定など、信号が物理的に「少数の到来波」で構成される問題で標準的に使われます。CS により AD コンバータのサンプリングレートを 1/10 に下げる sub-Nyquist receiver が実装されています。

ゲノミクス・group testing：COVID-19 検査でも話題になった「集団検査（pooled testing）」は本質的に圧縮センシングです。N 人の感染状態（K 人陽性、K ≪ N）を、M 個の混合検体で推定する問題は CS の枠組みで定式化でき、検査キット数を桁違いに減らせます。1 細胞 RNA-seq の遺伝子発現推定でも同様の応用があります。

よくある誤解と注意点

最大の落とし穴は、「信号が本当に K-sparse か」を確認しないことです。CS の理論はあくまで「K-sparse あるいは K-圧縮可能」という前提に依存しています。実世界の信号がそのままスパースなことは稀で、wavelet・DCT・curvelet などの適切な変換領域でスパースになるかを最初に検証する必要があります。画像なら wavelet で K/N ≈ 5%、音声なら DCT で K/N ≈ 10% 程度が目安。完全にホワイトノイズに近い信号には CS は全く効きません。「適切な基底（dictionary）を選ぶ」ことが応用での 9 割の仕事です。

次に、雑音耐性を過信する誤解。本ツールの再構成誤差式 ‖x̂−x‖₂ ≤ C·σ·√M / √(M−M_req) からも分かるように、観測数 M が理論最小 M_req に近づくほど分母が小さくなり、誤差は爆発的に増えます。論文上の「M = M_req で完全復元」は雑音ゼロの理想条件です。実応用では M ≥ 1.5〜3 倍の M_req を取り、雑音マージンを確保するのが鉄則。SNR が低い系（MRI なら 20dB 以下）では特に M を多めに取る必要があります。

最後に、計算量を見積もらずに大規模問題に BP を当てないこと。Basis Pursuit は美しい厳密解ですが、N³ オーダーで N=10,000 でも 10¹² 演算、画像 512²=262,144 なら 10¹⁶ 演算と現実離れします。実応用では FISTA・ADMM・SPGL1 のような first-order 法、または ISTA をニューラルネット展開した LISTA 等を使います。研究プロトタイプではまず BP で「達成可能な精度の上限」を測り、本番実装では速度と精度のバランスを取る、という二段階のアプローチが標準です。

使い方ガイド

  1. 信号次元 N（例：1024）とスパース性 K（例：32）を設定し、観測数 M を変動させる
  2. 観測行列をガウスランダム行列またはフーリエ部分行列から選択し、Candes-Tao条件 M ≥ C·K·log(N/K) との比較を確認
  3. 雑音レベルを 0～0.1 の範囲で調整し、BP（基底追従）、IST（反復軟閾値）、OMP（直交マッチング追従）の3アルゴリズムで再構成誤差を比較
  4. RIP定数 δ_{2K} が理論値 0.4 以下であることを確認し、復元成功条件を検証

具体的な計算例

N=512 次元信号、K=16 スパース、M=80 観測の場合：圧縮率 M/N=15.6%、Candes-Tao理論最小値 M_min=98（K·log(N/K)≈140）。ガウスランダム観測行列で RIP δ_{2K}≈0.32、白色ガウス雑音σ=0.01 添加時に BP で再構成誤差 ‖x̂−x‖₂≈0.015、IST で 0.018、OMP で 0.022 となり、BP の優位性を実証できます。

実務での注意点