スパースモデリングを理解したい

最近、青の本*1でスパーズモデリングの勉強を始めた。

最適化の記事同様、ここではあまり式を使わずに全体を俯瞰できるように努める。

まずは種々アルゴリズムの概要を掴み、後の記事で機能する理由について数学的側面を纏めていく予定。

随時更新。

この書籍は、劣決定の連立方程式

$$ \begin{align} \boldsymbol{Ax} = \boldsymbol{b} \quad \boldsymbol{A} \in \mathbb{R}^{n \times m} \end{align} $$

についてスパースな解を得る様々な方法について理論と実践の両方を解説してくれている。どちらかというと理論重きで、応用も画像処理のことがメイン。
「動かしたい」というよりも「理解したい」と思う人が読むべき本だと思う。

l0最小化問題（（P0）問題）
l0からl2への緩和
- FOCUSSアルゴリズム（Focal undetermined system solver）
l0からl1への緩和
- 一般的な緩和法
制約条件の緩和（等式制約から不等式制約へ）
- 解の安定性（最もスパースな場合）
- 追跡アルゴリズムの拡張

l0最小化問題（（P0）問題）

OMP, LS-OMP, WMPなどサポートを選んで解を更新することを繰り返す追跡アルゴリズム（MP）と呼ばれる手法がある。
また、最初にサポートを全て決定してから解を求めるしきい値アルゴリズムと呼ばれる手法もある。 LS-OMPなど逆行列計算が必要なアルゴリズムにおいては、行列が特異行列にならないなど注意が必要。
書籍では性能比較もしているが、検証の範囲ではOMPとLS-OMPはあまり性能差がないように思えた（個人の感想）。
計算が不安定になるくらいならLS-OMPではなく普通のOMPを使えば良いと思う。

直交マッチング追跡（Orthogonal Matching Pursuit; OMP）

残差との内積が最大となる $\boldsymbol{A}$ の列を選んでサポートを追加
以下を解いて暫定解を更新

$$ \begin{align} ||\boldsymbol{b}-\boldsymbol{Ax}||_2^2 \quad \text{subject to Support }{\boldsymbol{x}}=S_{k} \end{align} $$

得られた $x$ を用いて残差を更新

最小二乗-直交マッチング追跡（Least-squares OMP; LS-OMP）

OMPを効率的に解く方法。暫定解と追加されるサポートの候補を以下の式で一度に求める。

$$ \begin{align} \left[\begin{array}{c} \boldsymbol{x}_{S} \\ z \end{array}\right]_{\mathrm{opt}}=\left[\begin{array}{cc} \boldsymbol{A}_{S}^{\top} \boldsymbol{A}_{S} & \boldsymbol{A}_{S}^{\top} \boldsymbol{a}_{i} \\ \boldsymbol{a}_{i}^{\top} \boldsymbol{A}_{S} & ||\boldsymbol{a}_{i}||_{2}^{2} \end{array}\right]^{-1}\left[\begin{array}{c} \boldsymbol{A}_{S}^{\top} \boldsymbol{b} \\ \boldsymbol{a}_{i}^{\top} \boldsymbol{b} \end{array}\right] \end{align} $$

$\boldsymbol{A}_{S}^{\top} \boldsymbol{A}_{S}$ は最初に一度だけ計算しておけばよく、あとは最小の誤差を与える $\boldsymbol{a}_i$ を決定すれば良い。
$\boldsymbol{a}_i$ と $\boldsymbol{A}_{S}$ の列が直交している場合はブロック対角となり、これまでに求めた解は変化せず、要素の追加だけが行われる。

マッチング追跡（MP）

上記よりもより雑な方法。サポート追加時に解の再計算を行わない。 OMPとは異なり前ステップまでに得られた解を更新せずに、新しい要素だけを追加するより単純な方法。そのため精度はOMPより低い。

弱マッチング追跡（WMP）

サポートに追加する次の要素の選択を準最適解で満足する方法。MPよりも雑な方法（ただしメリットもちゃんとある）
最適な洗濯から $t \in (0, 1)$ 倍だけ離れているいずれかの列を選択する。要は候補列と残差の最大の内積を探すのではなく、tだけ小さい閾値を超え得た最初の列を採用することにする（その分１ループあたりの計算が速くなる）。
探索範囲を限定することで、探索が高速になる可能性がある（全体として高速になるという保証はない）

しきい値アルゴリズム

OMPアルゴリズムを簡単にしたもの。最初の射影で内積が大きい順にサポートを $k$ 列選んでしまう。つまり、最初に $\boldsymbol{A}$ の全ての列について

$$ \begin{align} \min _{z_{j}} ||\boldsymbol{a}_{j} z_{j} - \boldsymbol{b}||_{2}^{2} \end{align} $$

を計算して、ソートしておけば良い。ただ、閾値は自分で決めるので、非ゼロ要素の個数が既知であると想定していることに注意。または、許容誤差をクリアするまで選ぶサポートの数を増やしていく方法もある。

l0からl2への緩和

$l_0$ 最小化問題は離散最適化なので扱いにくい。そこで、 $l_0$ 最小化問題を $l_p \quad (p \in (0, 1$ )]最小化問題など、滑らかな関数の問題に緩和する手法や $l_1$ 最小化問題に置き換える方法がある。これにより一般的な最適化問題として取り扱えるようになる。ただし、大域最適解へ収束すると言う保証はなくなる。

FOCUSSアルゴリズム（Focal undetermined system solver）

反復重み付け最小二乗法（Itereative-reweighted-least-squares; IRLS）を用いて、 $l_p$ ノルム（ $p \in (0, 1$ ]）を重み付き $l_2$ ノルムで表現する方法。現在の近似解を $\boldsymbol{x}_k$ と書くことにして、以下の問題を反復して解く。

$$ \begin{align} \left(M_{k}\right): \min_{\boldsymbol{x}} || \boldsymbol{X}_{k-1}^{+} \boldsymbol{x} ||_{2}^{2} \quad \text { subject to } \quad \boldsymbol{b}=\boldsymbol{Ax} \end{align} $$

ここで、

$$ \begin{align} \boldsymbol{X}_{k-1} = \mathrm{diag} \left( |\boldsymbol{x}|^q \right) \end{align} $$

とすると

$$ \begin{align} \boldsymbol{X}_{k-1}^{-1} \boldsymbol{x} &= \left( \frac{x_1}{\left| x_1^q \right|}, \ldots, \frac{x_m}{\left| x_m^q \right|} \right)^{\top} \\ || \boldsymbol{X}_{k-1}^{-1} \boldsymbol{x} ||_2^2 &= ||\boldsymbol{x}||_{2-2q}^{2-2q} \end{align} $$

であり、 $p=2-2q$ と考えれば $||\boldsymbol{x}||_{p}^{p}$ を模したことになる。つまり、 $q$ の値によっていろいろな $l_p$ ノルムを $l_2$ ノルム最小化問題に置き換えることができる。特に、 $q=1$ とすれば $l_0$ の置き換えとなる。IRLSの「反復重み付け」は、この $l_2$ ノルムにおける $\boldsymbol{X}_{k-1}^{-1}$ に由来してるのだと思う（個人の感想）。実際には一般化逆行列 $\boldsymbol{X}_{k-1}^{+}$ を用いる点に注意。また、そのまま用いると対角成分にゼロが現れるので注意。計算に影響するのは選ばれたサポートに関する対角成分と暫定解の内積なので、計算上は対角ゼロ要素の部分は逆行列をとってもゼロとして仕舞えばよい。

$l_2$ 最小化問題に置き換えたので、ラグランジュの未定乗数法で解くことができる。得られた解 $\boldsymbol{x}$ を使って $\boldsymbol{X}_{k}$ を更新し、再び同様の手順で計算を行うことを繰り返す。 $\boldsymbol{X}_{k}$ 対角成分が0（つまり近似解 $\boldsymbol{x}_{k}$ にゼロ要素が現れたら）だと以降ずっとゼロになるので、初期解は全要素を非ゼロにする必要がある点に注意。

l0からl1への緩和

$(P_0)$ 問題を $(P_1)$ 問題に置き換える手法としては、基底追跡（Basis pursuit; BP）がある。 $l_0$ では $\boldsymbol{x}$ の要素の値は影響しなかったが、 $l_p$ では影響するので $\boldsymbol{A}$ の列を正規化しておきましょうという話である（当然列ごとの重みは変えられるが）。

計算時間はBPよりOMPが圧倒的に短くて済むらしいが、厳密解との $l_2$ 誤差、サポート間距離などの評価尺度ではIRLS（FOCUSSアルゴリズム）やBPの方が断然良い。

個人的には次元削減したいのか、解が欲しいのかで使い分ければ良いのだろうと思った。

一般的な緩和法

$l_0$ を単純に $l_1$ に置き換えてしまう方法。気をつけるべきは $l_0$ ではxの要素の値は影響しなかったが、 $l_0$ 以外では影響する点である。

つまり、行列 $\boldsymbol{A}$ の列を正規化しておくなど、スケーリングを行っておく必要がある。

制約条件の緩和（等式制約から不等式制約へ）

等式制約条件 $\boldsymbol{b}=\boldsymbol{Ax}$ が厳しすぎるので、許容誤差 $\epsilon > 0$ を導入して $(P_0)$ 問題を以下のように緩和する。

$$ \begin{align} \left(P_{0}^{\epsilon}\right): \min _{\boldsymbol{x}}||\boldsymbol{x}||_{0} \text { subject to }||\boldsymbol{b}-\boldsymbol{A} \boldsymbol{x}||_{2} \leq \epsilon \end{align} $$

一意性と最適性の保証がなくなるので、 $\epsilon$ に応じた解の安定性について考える必要がある。

解の安定性（最もスパースな場合）

以下の問題設定で、近似解 $\boldsymbol{x}_{0}^{\epsilon}$ を求めたいとする。

$$ \begin{align} \boldsymbol{x}_{0}^{\epsilon}=\operatorname{argmin}_{\boldsymbol{x}} || \boldsymbol{x} ||_{0} \quad \text { subject to } \quad || \boldsymbol{b}-\boldsymbol{A} \boldsymbol{x} ||_{2} \leq \epsilon \end{align} $$

この式において、 $\epsilon$ はノイズと解釈できる。すなわち、ノイズの大きさだけ解が不安定になる。非ゼロ要素が１つだけしかない最もスパースな解について考えると、一定の範囲までは、選ばれるサポートは一意だが、ノイズが大きくなるとサポートの一意性すらなくなる。詳細は書籍p96, 図5.2を参照すること。

追跡アルゴリズムの拡張

先述した追跡アルゴリズムを、 $(P_0)$ 問題を解くことから $(P_{0}^{\epsilon})$ を解く問題へと緩和する。

OMP

停止条件を $\epsilon_0 = \epsilon$ とすれば

$$ \begin{align} ||\boldsymbol{b} - \boldsymbol{Ax}||_2 \le \epsilon \end{align} $$

が満たされるまで解ベクトルの非ゼロ要素を追加していけば良い。

BPDN (Basis pursuit denoising; BPDN)

OMPの場合とほとんど同様だが、 $l_0$ ノルムを $l_1$ ノルムへ置き換えて以下の問題を解く。
（行列 $\boldsymbol{A}$ の各列は正規化されているものとする。）

$$ \begin{align} \left(P_{1}^{\epsilon}\right): \min _{\boldsymbol{x}}||\boldsymbol{x}||_{1} \text { subject to }||\boldsymbol{b}-\boldsymbol{Ax}||_{2} \leq \epsilon \end{align} $$

この問題は、ラグランジュ乗数 $\lambda$ を用いて以下の制約なし最適化を解くことと同じである。

$$ \begin{align} \left(Q_{1}^{\lambda}\right): \min _{\boldsymbol{x}} \lambda||\boldsymbol{x}||_{1}+\frac{1}{2}||\boldsymbol{b}-\boldsymbol{A} \boldsymbol{x}||_{2}^{2} \end{align} $$

IRLSによる解法

$\left(Q_{1}^{\lambda}\right)$ を解く単純な手法として反復再重み付け最小二乗法（IRLS）がある。FOCUSSアルゴリズムで用いている手法だが、枠組みが同じだけで計算の仕方は異なるので注意。

$\boldsymbol{X}=\mathrm{diag}\left(|x|\right)$ とすると、 $||\boldsymbol{x}||_1 = \boldsymbol{x}^{\top}\boldsymbol{X}^{-1}\boldsymbol{x}$ とと書ける。これは $l_1$ ノルムを重み付き $l_2$ ノルムで表現したことになる。これを利用すると $\left(Q_{1}^{\lambda}\right)$ は以下のように置き換えられる。

$$ \begin{align} \left(M_{k}\right): \min _{\boldsymbol{x}} \lambda \boldsymbol{x}^{\top} \boldsymbol{X}_{k-1}^{-1} \boldsymbol{x}+\frac{1}{2}||\boldsymbol{b}-\boldsymbol{Ax}||_{2}^{2} \end{align} $$

二次最適化になっており、標準的な線形代数を用いた反復計算を用いて解くことができる。具体的には、この目的関数を $\lambda(x)$ で微分して0とおくと

$$ \begin{align} \left(2 \lambda \boldsymbol{X}_{k-1}^{-1}+\boldsymbol{A}^{\top} \boldsymbol{A}\right) \boldsymbol{x}=\boldsymbol{A}^{\top} \boldsymbol{b} \end{align} $$

が得られるので、各ステップで $\boldsymbol{x}$ について数値最適化を行い、得た近似解で重み行列 $\boldsymbol{X}$ を更新することを繰り返すだけである。解の変化がしきい値以下になったら計算を終了する。 $\lambda$ は経験的に決定される。

LARSアルゴリズム

以下の目的関数を最小化することを考える。

$$ \begin{align} f(\boldsymbol{x})=\lambda||\boldsymbol{x}||_{1}+\frac{1}{2}||\boldsymbol{b}-\boldsymbol{A} \boldsymbol{x}||_{2}^{2} \end{align} $$

この式の劣微分（劣勾配集合）

$$ \begin{align} \partial f(\boldsymbol{x})=\left\{\boldsymbol{A}^{\top}(\boldsymbol{A} \boldsymbol{x}-\boldsymbol{b})+\lambda \boldsymbol{z}\right\}, \quad \forall \boldsymbol{z}=\left\{\begin{array}{cc} +1 & x[i]>0 \\ {[-1,+1]} & x[i]=0 \\ -1 & x[i]<0 \end{array}\right. \label{lars_subderiv} \end{align} $$

が0となる $\boldsymbol{x}$ と $\boldsymbol{z}$ を同時に求める。なお、この時の $z$ を決める処理をソフトしきい値処理として知られている（ハードしきい値処理は0にするか非ゼロにするかの二値化処理のこと）。

以下の話は、解パス図において $\lambda$ を大きい方から小さい方へ動かすことを想像するとわかりやすい。

$\lambda \rightarrow \infty$ において、 $\boldsymbol{x}_{\lambda}=\boldsymbol{0}$ である。また、このとき

$$ \begin{align} \operatorname{argmin}_{\boldsymbol{x}} f(\boldsymbol{x}) = \boldsymbol{0} \end{align} $$

より、

$$ \begin{align} \partial f(\boldsymbol{x}) &= \boldsymbol{0} \\ -A_{\top}\boldsymbol{b} + \lambda \boldsymbol{z}_{\lambda} &= \boldsymbol{0} \\ \boldsymbol{z}_{\lambda} &= A^{\top}\boldsymbol{b} / \lambda \label{z_inf} \end{align} $$

が導かれる。すなわち、 $\lambda \ge ||A^{\top}\boldsymbol{b}||_{\infty}$ ならば、 $\boldsymbol{z} \le 1$ であることが保証され、式\eqref{lars_subderiv}より最適解は $x[i]=0$ のままであることがわかる。

さらに $\lambda$ を小さくしていき、

$$ \begin{align} \lambda = ||A^{\top} \boldsymbol{b}||_{\infty} \end{align} $$

となったとき、これよりさらに $\lambda$ をさらに小さくしようとすると（すなわち $A^{\top} \boldsymbol{b}$ の要素の最大値よりも小さくしようとすると）、 $\boldsymbol{x} = \boldsymbol{0}$ のままでは式\eqref{lars_subderiv}の $\boldsymbol{z}$ の要請を満たさない。そこで、

$$ \begin{align} z_{\lambda}[i] = \operatorname{sign} (x_{\lambda}[i]) \end{align} $$

を導入して、改めて $\partial f(\boldsymbol{x}) = \boldsymbol{0}$ を考えると

$$ \begin{align} 0=& \boldsymbol{a}_{i}^{\boldsymbol{T}}(\boldsymbol{a} x-\boldsymbol{b})+\lambda z_{\lambda}[i] \\ & \Rightarrow x_{\lambda}[i]=\frac{\boldsymbol{a}_{i}^{\top} \boldsymbol{b}-\lambda z_{\lambda}[i]}{\boldsymbol{a}_{i}^{\top} \boldsymbol{a}_{i}}=\frac{\boldsymbol{a}_{i}^{\top} \boldsymbol{b}-\lambda \operatorname{sign}\left(x_{\lambda}[i]\right)}{\boldsymbol{a}_{i}^{\top} \boldsymbol{a}_{i}} \end{align} $$

このとき、 $\lambda \rightarrow \infty$ における $\boldsymbol{z}$ の条件式\eqref{z_inf}は

$$ \begin{align} A^{\top}(A\boldsymbol{x} - \boldsymbol{b}) + \lambda \boldsymbol{z} &= 0 \\ \boldsymbol{z} &= \frac{A^{\top}(\boldsymbol{b} - A\boldsymbol{x})}{\lambda} \end{align} $$

という式に置き換えられている。

さらに $\lambda$ の値を減らしていくと、現在の解は $\partial f(\boldsymbol{x}) = \boldsymbol{0}$ より、

$$ \begin{align} A^{\top}(A\boldsymbol{x} - \boldsymbol{b}) + \lambda \boldsymbol{z} &= 0 \\ \boldsymbol{x} = (A^{\top}A)^{-1} (A^{\top}\boldsymbol{b} - \lambda\boldsymbol{z}) \end{align} $$

と書け、サポート $S$ を使って表せば

$$ \begin{align} \left\{\begin{array}{c} \boldsymbol{x}^{S}_{\lambda} = (A_{S}^{\top}A_{S})^{-1} (A_{S}^{\top}\boldsymbol{b} - \lambda\boldsymbol{z}_{\lambda}^{S}) \\ \boldsymbol{z}_{\lambda}^{S} = \operatorname{sign} (\boldsymbol{x}_{\lambda}^{S}) \\ \end{array}\right. \end{align} $$

と書ける。

$\lambda$ を小さくしていく過程で、 $\boldsymbol{x}_{\lambda}$ のサポート外の $\boldsymbol{z}_{\lambda}$ が $[-1, +1$ ]の範囲から出ることがある。この場合はこの要素をサポートに追加し、解を更新する。逆に、解の非ゼロだった要素が0になってしまうことがある。この場合は該当するサポートを除いて解を更新する。