[PYTHON] Sparsame Modellierung: Kapitel 5 Von der exakten Lösung zur ungefähren Lösung

Überblick

Sparse Modeling Implementierte die numerischen Beispiele in Kapitel 5 in Python. Jupyter-Notizbuch, das den Code und die experimentellen Ergebnisse zusammenfasst. Einführung

Kapitel 5 befasst sich mit dem Fall, in dem die beobachteten Daten Rauschen enthalten. Da die beobachteten Daten Rauschen enthalten, wird die Einschränkung des Fehlerterms auf $ \ epsilon $ oder weniger abgeschwächt.

(P_{0}^{\epsilon}): \min_{\mathbf{x}} \|\|\mathbf{x}\|\|_{0} \text{ subject to } \|\|\mathbf{b}-\mathbf{Ax}\|\| _{2} \leq \epsilon

Finden Sie in Orthogonal Matching Tracking (OMP) einfach die Lösung mit der Stoppbedingung $ \ epsilon_ {0} = \ epsilon $ für die iterative Berechnung.

Das Lockern der $ \ ell_ {0} $ -Norm auf die $ \ ell_ {1} $ -Norm ergibt eine Variante von $ (P_ {1}) $, Basisverfolgungsentrauschung (BPDN).

(P_{1}^{\epsilon}): \min_{\mathbf{x}} \|\|\mathbf{x}\|\|_{1} \text{ subject to } \|\|\mathbf{b}-\mathbf{Ax}\|\| _{2} \leq \epsilon

Für einen geeigneten Lagrange-Multiplikator $ \ lambda $ stimmt die Lösung von $ (P_ {1} ^ {\ epsilon}) $ mit der Lösung des folgenden uneingeschränkten Optimierungsproblems überein.

(Q_{1}^{\lambda}): \min_{\mathbf{x}} \lambda \|\|\mathbf{x}\|\|_{1} + \frac{1}{2} \|\|\mathbf{b}-\mathbf{Ax}\|\| _{2}

Eine einfache Möglichkeit, $ (Q_ {1} ^ {\ lambda}) $ zu lösen, sind die iterativ neu gewichteten kleinsten Quadrate (IRLS).

\mathbf{X}=\rm{diag}(\mathbf{\|x\|})Dann$ ||\mathbf{x}||_{1} = \mathbf{x}\mathbf{X}^{-1}\mathbf{x}Es wird. Mit anderen Worten\ell_{1}Norm ist (gewichtet) quadratisch\ell_{2}Es kann als Norm angesehen werden. Aktuelle ungefähre Lösung\mathbf{x_{\rm{k-1}}}Gegeben,\mathbf{X_{\rm{k-1}}}=\rm{diag}(\left|\mathbf{x_{\rm{k-1}}}\right|)$Lösen Sie auf jeden Fall das folgende Problem.

(M_{k}): \min_{\mathbf{x}} \lambda \mathbf{x^{\rm{T}}} \mathbf{X^{\rm{-1}}} \mathbf{x} + \frac{1}{2} \|\| \mathbf{b} - \mathbf{Ax} \|\|_{2}^{2}

Der IRLS-Algorithmus verwendet also die aktuellen Gewichte $ \ mathbf {X_ {\ rm {k-1}}} $, um $ \ mathbf {x_ {k}} $ und dann $ \ mathbf {x_ {k} zu erhalten } $ Wird verwendet, um abwechselnd $ \ mathbf {X_ {\ rm {k}}} $ zu erhalten, und die Lösung von $ (Q_ {1} ^ {\ lambda}) $ wird ungefähr erhalten. .. Kapitel 5 besteht hauptsächlich aus einer Erklärung von IRLS, die ungefähr $ (Q_ {1} ^ {\ lambda}) $ löst, und einem Umriss von LARS (stufenweise Regression des kleinsten Winkels).

IRLS-Algorithmus

Initialisieren

Als $ k = 0 $

(Beliebige) Anfangslösung $ \ mathbf {x_ {\ rm {0}}} = \ mathbf {1} $
Anfangsgewichtsmatrix $ \ mathbf {X_ {\ rm {0}}} = \ mathbf {I} $

Hauptschleife

Setzen Sie $ k \ leftarrow k + 1 $ und führen Sie die folgenden Schritte aus.

Minimales Quadrat mit Regularisierung: Die folgenden linearen Simultangleichungen $ \ left (2 \ lambda \ mathbf {X_ {\ rm {k-1}} ^ {\ rm {-1}}} + \ mathbf {A ^ {\ rm {T}} A} \ right) \ mathbf {x} = \ mathbf {A ^ {\ rm {T}} b} $ gelöst mit der iterativen Methode, ungefähre Lösung $ \ mathbf {x_ { Holen Sie sich \ rm {k}}} $. Das Lehrbuch besagt, dass mehrere Iterationen der konjugierten Gradientenmethode ausreichend sind. In dieser Studie wird scipy.optimize.minimize verwendet.

*Gewichtsaktualisierung:\mathbf{x_{\rm{k}}}Verwenden vonX_{k}(j,j)=\|x_{k}(j)\|+\epsilonAls Gewichtsdiagonalmatrix\mathbf{X}Ist aktualisiert. *Stoppbedingung: Wenn\|\|\mathbf{x_{\rm{k}}}-\mathbf{x_{\rm{k-1}}}\|\|_{2}Wenn der vorgegebene Schwellenwert kleiner ist, wird er beendet, andernfalls wird er wiederholt.

Experiment

Simulation

$ \ mathbf {x} $ 200-dimensionaler Vektor $ \ mathbf {b} $ 100-dimensionaler Vektor $ \ mathbf {A} $ 100 × 200 dimensionale Matrix, Spaltennormalisierung einheitlicher Zufallszahlen von $ [-1, 1) $

$ \ mathbf {x} $ Anzahl der Nicht-Null-Elemente $ k = 4 $, und der Wert der Nicht-Null-Elemente ist eine gleichmäßige Verteilung im Bereich von $ [-2, -1] \ cup [1, 2] $ Extrahiere iid aus.
Berechnen Sie $ \ mathbf {b} = \ mathbf {Ax} + \ mathbf {e} $ mit dem additiven Rauschen $ \ mathbf {e} $, das aus der Gaußschen Verteilung mit der Standardabweichung $ \ sigma = 0.1 $ extrahiert wurde. Machen.

Eines der zu lösenden Probleme ist die Bestimmung des Werts von $ \ lambda $. Hier ist nach der empirischen Regel $ \ lambda $ ein Wert, der nahe am Verhältnis von $ \ sigma $ und der Standardabweichung von Nicht-Null-Elementen liegt. Da die Standardabweichung von Nicht-Null-Elementen ungefähr 2 beträgt, wird der Wert in der Nähe von $ \ sigma / 2 = 0,05 $ als $ \ lambda $ festgelegt.

Ergebnis

\lambdaDie Lösung, die beim Ändern des Wertes von erhalten wird$ \frac{|| \hat{\mathbf{x_{\rm{\lambda}}}} - \mathbf{x_{\rm{0}}} ||}{||\mathbf{x_{\rm{0}}}||} $Bewertet in. Die Anzahl der IRLS-Wiederholungen betrug 100.

Die gestrichelte Linie ist\left\|\log\left(\|\|\mathbf{A\hat{x_{\rm{\lambda}}}}-\mathbf{b}\|\|^{2}/\left(n\sigma^{2}\right)\right)\right\|Zeigt den Wert von an\lambdaRest zu welchem Wert\|\|\mathbf{A\hat{x_{\rm{\lambda}}}}-\mathbf{b}\|\|Ist die Kraft des Rauschensn\sigma^{2}Es zeigt, ob es ungefähr das gleiche sein wird wie. das ist\lambdaDies ist eine weitere empirische Regel, die den Wert von bestimmt.

Drei Lösungen bei Verwendung des besten $ \ lambda $ -Werts und bei Verwendung größerer und kleinerer Werte

Im Vergleich zur wahren Lösung $ \ mathbf {x_ {\ rm {0}}} $ (roter Kreis), $ \ mathbf {x_ {\ rm {0}} bei Verwendung des optimalen Wertes $ \ lambda $ } Die Position des Nicht-Null-Elements von $ wird am besten wiederhergestellt. Wenn der Wert von $ \ lambda $ klein ist, ist die Lösung dicht, und wenn er groß ist, ist die Lösung spärlich.

Zeigen Sie Lösungen für alle $ \ lambda $ gleichzeitig in einem Diagramm an

Jedes Diagramm zeigt den Wert eines Elements von $ \ mathbf {x} $ für $ \ lambda $. Die gestrichelte Linie ist der Wert von $ \ mathbf {x_ {\ rm {0}}} $.

*Optimal\lambdaFunktionswert für die Anzahl der IRLS-Iterationen bei(\lambda\|\|\mathbf{x}\|\|_{1}+\frac{1}{2}\|\|\mathbf{b}-\mathbf{Ax}\|\|)

LARS-Algorithmus

Ich konnte es nicht verstehen, also habe ich nur mit sklearns Lasso Lars experimentiert.

Experiment

Simulation

$ \ mathbf {x} $ 50-dimensionaler Vektor $ \ mathbf {b} $ 30-dimensionaler Vektor $ \ mathbf {A} $ 30 × 50 dimensionale Matrix, Spaltennormalisierung einheitlicher Zufallszahlen von $ [-1, 1) $

Unter der Annahme, dass die Anzahl der Nicht-Null-Elemente k von $ \ mathbf {x} $ ist, wird der Wert des Nicht-Null-Elements aus der Gleichverteilung im Bereich von $ [-2, -1] \ cup [1, 2] $ extrahiert. ..
Berechnen Sie $ \ mathbf {b} = \ mathbf {Ax} + \ mathbf {e} $ mit dem additiven Rauschen $ \ mathbf {e} $, das aus der Gaußschen Verteilung mit der Standardabweichung $ \ sigma = 0.1 $ extrahiert wurde. Machen.
Simulieren Sie jeweils 200 Mal mit $ k = 1, ..., 15 $
Vergleichen Sie OMP und LARS

Index

$ L_2 $ Norm des geschätzten $ \ mathbf {\ hat {x}} $ und des wahren $ \ mathbf {x} $ *Geschätzt\mathbf{\hat{x}}Und wahr\mathbf{x}Abstand zwischen den Stützen$dist(\hat{S},S)=\frac{max\\{|\hat{S}|, |S| \\} - |\hat{S}\cap S|}{max\\{|\hat{S}|,|S|\\}}$
Nehmen Sie den Durchschnitt von 200 Simulationen.

Ergebnis

Erwägung

Wenn k klein ist, hat OMP einen kleineren Fehler. Wenn k groß ist, hat LARS einen kleineren Fehler.
Die Wiederherstellung des Supports scheint in OMP genauer zu sein. Es kehrt sich um, wenn K zunimmt.