[PYTHON] Modélisation parcimonieuse: Chapitre 5 De la solution exacte à la solution approximative

Aperçu

Modélisation éparse Implémentation des exemples numériques du chapitre 5 en Python. Cahier Jupyter résumant le code et les résultats expérimentaux. introduction

Le chapitre 5 traite du cas où les données observées contiennent du bruit. Puisque les données observées contiennent du bruit, la contrainte du terme d'erreur est affaiblie à $ \ epsilon $ ou moins.

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

Dans le suivi de correspondance orthogonale (OMP), trouvez simplement la solution comme condition d'arrêt $ \ epsilon_ {0} = \ epsilon $ pour le calcul itératif.

Le relâchement de la norme $ \ ell_ {0} $ à la norme $ \ ell_ {1} $ donne une variante de $ (P_ {1}) $, débruitage de poursuite de base (BPDN).

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

Pour un multiplicateur de Lagrange $ \ lambda $ approprié, la solution de $ (P_ {1} ^ {\ epsilon}) $ est cohérente avec la solution du problème d'optimisation sans contrainte ci-dessous.

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

Un moyen simple de résoudre $ (Q_ {1} ^ {\ lambda}) $ est l'itératif repondéré par les moindres carrés (IRLS).

\mathbf{X}=\rm{diag}(\mathbf{\|x\|})Puis$ ||\mathbf{x}||_{1} = \mathbf{x}\mathbf{X}^{-1}\mathbf{x}Il devient. En d'autres termes\ell_{1}La norme est (pondérée) au carré\ell_{2}Cela peut être considéré comme une norme. Solution approximative actuelle\mathbf{x_{\rm{k-1}}}Donné,\mathbf{X_{\rm{k-1}}}=\rm{diag}(\left|\mathbf{x_{\rm{k-1}}}\right|)$Quoi qu'il en soit, résolvez le problème suivant.

(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}

L'algorithme IRLS utilise donc les poids actuels $ \ mathbf {X_ {\ rm {k-1}}} $ pour obtenir $ \ mathbf {x_ {k}} $ puis $ \ mathbf {x_ {k} } $ Est utilisé pour obtenir $ \ mathbf {X_ {\ rm {k}}} $ en alternance, et approximativement $ (Q_ {1} ^ {\ lambda}) $ est obtenu. .. Le chapitre 5 consiste principalement en une explication de IRLS qui résout approximativement $ (Q_ {1} ^ {\ lambda}) $ et un aperçu de LARS (régression du moindre angle par étape).

Algorithme IRLS

Initialisation

Comme $ k = 0 $

Solution initiale (arbitraire) $ \ mathbf {x_ {\ rm {0}}} = \ mathbf {1} $
Matrice de poids initial $ \ mathbf {X_ {\ rm {0}}} = \ mathbf {I} $

Boucle principale

Définissez $ k \ leftarrow k + 1 $ et exécutez les étapes suivantes.

Carré minimum avec régularisation: Les équations linéaires simultanées suivantes $ \ left (2 \ lambda \ mathbf {X_ {\ rm {k-1}} ^ {\ rm {-1}}} + \ mathbf {A ^ {\ rm {T}} A} \ right) \ mathbf {x} = \ mathbf {A ^ {\ rm {T}} b} $ résolu en utilisant la méthode itérative, solution approximative $ \ mathbf {x_ { Récupérez \ rm {k}}} $. Le manuel indique que plusieurs itérations de la méthode du gradient conjugué sont suffisantes. Dans cette étude, scipy.optimize.minimize est utilisé.

*Mise à jour du poids:\mathbf{x_{\rm{k}}}En utilisantX_{k}(j,j)=\|x_{k}(j)\|+\epsilonEn tant que matrice diagonale de poids\mathbf{X}Est mis à jour. *Condition d'arrêt: Si\|\|\mathbf{x_{\rm{k}}}-\mathbf{x_{\rm{k-1}}}\|\|_{2}Si est inférieur au seuil prédéfini, il sort, sinon il se répète.

Expérience

simulation

$ \ mathbf {x} vecteur de 200 dimensions $ $ \ mathbf {b} vecteur $ 100 dimensions $ \ mathbf {A} Matrice dimensionnelle de 100 $ × 200, normalisation de colonne de nombres aléatoires uniformes de $ [-1, 1) $

$ \ mathbf {x} $ nombre d'éléments non nuls $ k = 4 $, et la valeur des éléments non nuls est une distribution uniforme dans la plage de $ [-2, -1] \ cup [1, 2] $ Extrayez iid de.
Calculez $ \ mathbf {b} = \ mathbf {Ax} + \ mathbf {e} $ en utilisant le bruit additif $ \ mathbf {e} $ extrait de la distribution gaussienne avec un écart type $ \ sigma = 0,1 $. Faire.

L'un des problèmes à résoudre est de savoir comment déterminer la valeur de $ \ lambda $. Ici, selon la règle empirique, $ \ lambda $ est une valeur proche du rapport $ \ sigma $ et de l'écart type des éléments non nuls. Comme l'écart type des éléments non nuls est d'environ 2, la valeur proche de $ \ sigma / 2 = 0,05 $ est définie comme $ \ lambda $.

résultat

\lambdaLa solution obtenue en modifiant la valeur de$ \frac{|| \hat{\mathbf{x_{\rm{\lambda}}}} - \mathbf{x_{\rm{0}}} ||}{||\mathbf{x_{\rm{0}}}||} $Évalué en. Le nombre de répétitions IRLS était de 100.

La ligne brisée est\left\|\log\left(\|\|\mathbf{A\hat{x_{\rm{\lambda}}}}-\mathbf{b}\|\|^{2}/\left(n\sigma^{2}\right)\right)\right\|Indique la valeur de\lambdaRésiduel à quelle valeur\|\|\mathbf{A\hat{x_{\rm{\lambda}}}}-\mathbf{b}\|\|Est la puissance du bruitn\sigma^{2}Il montre si ce sera à peu près le même que. c'est\lambdaC'est une autre règle empirique qui détermine la valeur de.

Trois solutions lors de l'utilisation de la meilleure valeur $ \ lambda $ et lors de l'utilisation de valeurs de plus en plus petites

Par rapport à la vraie solution $ \ mathbf {x_ {\ rm {0}}} $ (cercle rouge), en utilisant la valeur optimale $ \ lambda $, $ \ mathbf {x_ {\ rm {0}} } La position de l'élément non nul de $ est mieux restaurée. Si la valeur de $ \ lambda $ est petite, la solution est dense, et si elle est grande, la solution est clairsemée.

Afficher les solutions pour tous les $ \ lambda $ à la fois dans un graphique

Chaque graphique montre la valeur d'un élément de $ \ mathbf {x} $ pour $ \ lambda $. La ligne brisée est la valeur de $ \ mathbf {x_ {\ rm {0}}} $.

*Optimale\lambdaValeur de fonction pour le nombre d'itérations IRLS lorsque(\lambda\|\|\mathbf{x}\|\|_{1}+\frac{1}{2}\|\|\mathbf{b}-\mathbf{Ax}\|\|)

Algorithme LARS

Je ne pouvais pas le comprendre, donc je n'ai expérimenté qu'avec Lasso Lars de sklearn.

Expérience

simulation

$ \ mathbf {x} $ 50 dimensions vecteur $ \ mathbf {b} vecteur $ 30 dimensions $ \ mathbf {A} $ 30 × 50 matrice dimensionnelle, normalisation de colonne de nombres aléatoires uniformes de $ [-1, 1) $

En supposant que le nombre d'éléments non nuls k de $ \ mathbf {x} $, la valeur de l'élément non nul est iid extraite de la distribution uniforme dans l'intervalle $ [-2, -1] \ cup [1, 2] $. ..
Calculez $ \ mathbf {b} = \ mathbf {Ax} + \ mathbf {e} $ en utilisant le bruit additif $ \ mathbf {e} $ extrait de la distribution gaussienne avec un écart type $ \ sigma = 0,1 $. Faire.
Simulez 200 fois chacun avec $ k = 1, ..., 15 $
Comparez OMP et LARS

indice

$ L_2 $ norme estimée de $ \ mathbf {\ hat {x}} $ et vraie $ \ mathbf {x} $ *Estimé\mathbf{\hat{x}}Et vrai\mathbf{x}Distance entre les supports$dist(\hat{S},S)=\frac{max\\{|\hat{S}|, |S| \\} - |\hat{S}\cap S|}{max\\{|\hat{S}|,|S|\\}}$
Prenez la moyenne de 200 simulations.

résultat

Considération

Si k est petit, OMP a une erreur plus petite. Lorsque k est grand, LARS a une erreur plus petite.
La restauration du support semble plus précise dans OMP. Il s'inverse lorsque K augmente.