Definition und Konvergenz iterativ neu gewichteter kleinster Quadrate

Ich habe iterativ die kleinsten Quadrate (IRLS) neu gewichtet, um Funktionen der folgenden Form zu minimieren:

$J(m) = \sum_{i=1}^{N} \rho \left(\left| x_i - m \right|\right)$

Dabei ist die Anzahl der Instanzen von , die robuste Schätzung, die ich möchte, und ist eine geeignete robuste Straffunktion. Nehmen wir an, es ist konvex (wenn auch nicht unbedingt streng) und vorerst differenzierbar. Ein gutes Beispiel für eine solche ist die Huber-Verlustfunktion . $N$ $x_i \in \mathbb{R}$ $m \in \mathbb{R}$ $\rho$ $\rho$

Was ich getan habe, ist, in Bezug auf differenzieren (und zu manipulieren), um zu erhalten, $J(m)$ $m$

$\frac{dJ}{dm}= \sum_{i=1}^{N} \frac{\rho'\left( \left|x_i-m\right|\right) }{\left|x_i-m\right|} \left( x_i-m \right)$

und iteratives Lösen, indem es gleich 0 gesetzt und die Gewichte bei der Iteration auf (beachte, dass die wahrgenommene Singularität bei wirklich eine entfernbare Singularität in allen , die mich interessieren könnten). Dann erhalte ich, $k$ $w_i(k) = \frac{\rho'\left( \left|x_i-m{(k)}\right|\right) }{\left|x_i-m{(k)}\right|}$ $x_i=m{(k)}$ $\rho$

$\sum_{i=1}^{N} w_i(k) \left( x_i-m{(k+1)} \right)=0$

und ich löse, um zu erhalten, . $m(k+1) = \frac{\sum_{i=1}^{N} w_i(k) x_i}{ \sum_{i=1}^{N} w_i(k)}$

Ich wiederhole diesen Fixpunktalgorithmus bis zur "Konvergenz". Ich werde feststellen, dass Sie optimal sind, wenn Sie an einen festen Punkt gelangen, da Ihre Ableitung 0 ist und es sich um eine konvexe Funktion handelt.

Ich habe zwei Fragen zu diesem Verfahren:

Ist dies der Standard-IRLS-Algorithmus? Nach dem Lesen mehrerer Artikel zum Thema (und sie waren sehr verstreut und vage darüber, was IRLS ist) ist dies die konsistenteste Definition des Algorithmus, den ich finden kann. Ich kann die Papiere posten, wenn die Leute wollen, aber ich wollte eigentlich niemanden hier voreingenommen machen. Natürlich können Sie diese grundlegende Technik auf viele andere Arten von Problemen verallgemeinern, die Vektor und andere Argumente als betreffen Die Angabe des Arguments ist eine Norm einer affinen Funktion Ihrer Parameter. Jede Hilfe oder Einsicht wäre großartig. $x_i$ $\left|x_i-m{(k)}\right|$
Konvergenz scheint in der Praxis zu funktionieren, aber ich habe einige Bedenken. Ich habe noch keinen Beweis dafür zu sehen. Nach einigen einfachen Matlab-Simulationen sehe ich, dass eine Iteration davon keine Kontraktionsabbildung ist (ich habe zwei zufällige Instanzen von generiert und und sah, dass dies gelegentlich größer als 1 ist). Auch das durch mehrere aufeinanderfolgende Iterationen definierte Mapping ist nicht unbedingt ein Kontraktions-Mapping, aber die Wahrscheinlichkeit, dass die Lipschitz-Konstante über 1 liegt, wird sehr gering. Gibt es also eine Vorstellung von einer Kontraktionskartierung der Wahrscheinlichkeit ? Mit welcher Maschinerie würde ich beweisen, dass dies konvergiert? Konvergiert es überhaupt? $m$ $\frac{\left|m_1(k+1) - m_2(k+1)\right|}{\left|m_1(k)-m_2(k)\right|}$

Jede Anleitung ist hilfreich.

Edit: Ich mag das Papier über IRLS für die sparsame Wiederherstellung / Druckmessung von Daubechies et al. 2008 "Iterativ neu gewichtete Minimierung der kleinsten Quadrate für eine sparsame Wiederherstellung" auf dem arXiv. Es scheint sich jedoch hauptsächlich auf Gewichte für nicht konvexe Probleme zu konzentrieren. Mein Fall ist wesentlich einfacher.

estimation least-squares robust irls Chris A.
quelle

Wenn ich mir die Wiki-Seite über IRWLS ansehe, habe ich Probleme mit dem Unterschied zwischen der von Ihnen beschriebenen Prozedur und IRWLS (sie verwenden einfach als ihre spezielle Funktion). Können Sie erklären , auf welche Weise denken Sie der Algorithmus Sie vorschlagen , ist verschieden von IRWLS?

| y_{i} - x x_{i}^{'} β β |^{2}

$|y_i-\pmb x_i'\pmb\beta|^2$

ρ

$\rho$

user603

Ich habe nie gesagt, dass es anders ist, und wenn ich es andeute, wollte ich es nicht.

Chris A.

Bei Ihrer ersten Frage sollte man "Standard" definieren oder anerkennen, dass nach und nach ein "kanonisches Modell" erstellt wurde. Wie ein Kommentar zeigt, scheint es zumindest so, als ob die Art und Weise, wie Sie IRWLS verwenden, eher Standard ist.

Was Ihre zweite Frage betrifft, könnte "Kontraktionsabbildung in der Wahrscheinlichkeit" (jedoch informell) mit der Konvergenz von "rekursiven stochastischen Algorithmen" verbunden sein. Nach allem, was ich lese, gibt es eine riesige Literatur zu diesem Thema, hauptsächlich im Ingenieurwesen. In Economics verwenden wir ein winziges Stück davon, insbesondere die wegweisenden Arbeiten von Lennart Ljung - das erste Papier war Ljung (1977) -, die zeigten, dass die Konvergenz (oder nicht) eines rekursiven stochastischen Algorithmus durch die Stabilität (oder nicht) einer verwandten gewöhnlichen Differentialgleichung.

(Das Folgende wurde nach einer fruchtbaren Diskussion mit dem OP in den Kommentaren überarbeitet.)

Konvergenz

Ich werde Sabre Elaydi "An Introduction to Difference Equations", 2005, 3d ed als Referenz verwenden . Die Analyse hängt von einer bestimmten Datenstichprobe ab, sodass die Werte als fest behandelt werden. $x's$

Die Bedingung erster Ordnung für die Minimierung der Zielfunktion, betrachtet als rekursive Funktion in , $m$

m (k + 1) = \sum_{ich = 1}^{N} v_{ich} [m (k)] x_{ich}, v_{ich} [m (k)] \equiv \frac{w_{ich} [m (k)]}{\sum_{ich = 1}^{N} w_{ich} [m (k)]} [1]

$m(k+1) = \sum_{i=1}^{N} v_i[m(k)] x_i, \;\; v_i[m(k)] \equiv \frac{w_i[m(k)]}{ \sum_{i=1}^{N} w_i[m(k)]} \qquad [1]$

hat einen festen Punkt (das Argmin der Zielfunktion). Nach Theorem 1.13, S. 27-28 von Elaydi ist, wenn die erste Ableitung in Bezug auf der RHS von , die am Fixpunkt ausgewertet wird , , kleiner als Eins in Absolutwert, dann ist asymptotisch stabil (AS). Weiter haben wir nach Satz 4.3 S.179, dass dies auch impliziert, dass der Fixpunkt einheitlich AS (UAS) ist. "Asymptotisch stabil" bedeutet, dass für einen gewissen Wertebereich um den Fixpunkt, eine Nachbarschaft , die nicht notwendigerweise klein ist, der Fixpunkt attraktiv ist $m$ $[1]$ $m^*$ $A'(m^*)$ $m^*$
$(m^* \pm \gamma)$ Wenn der Algorithmus also Werte in dieser Nachbarschaft angibt, konvergiert er. Die Eigenschaft "einheitlich" bedeutet, dass die Grenze dieser Nachbarschaft und damit ihre Größe unabhängig vom Anfangswert des Algorithmus ist. Der Fixpunkt wird global UAS, wenn . Also in unserem Fall, wenn wir das beweisen $\gamma = \infty$

| {EIN}^{'} (m^{*}) | \equiv | \sum_{ich = 1}^{N} \frac{\partial v_{ich} (m^{*})}{\partial m} x_{ich} | < 1 [2]

$|A'(m^*)|\equiv \left|\sum_{i=1}^{N} \frac{\partial v_i(m^*)}{\partial m}x_i\right| <1 \qquad [2]$

Wir haben die UAS-Eigenschaft bewiesen, jedoch ohne globale Konvergenz. Dann können wir entweder versuchen festzustellen, dass die Nachbarschaft der Anziehung tatsächlich die ganzen erweiterten reellen Zahlen sind, oder dass der spezifische Startwert, den das OP verwendet, wie in den Kommentaren erwähnt (und er ist Standard in der IRLS-Methodik), dh der Stichprobenmittelwert von den 's gehört immer zur Nachbarschaft der Anziehungskraft des Fixpunktes. $x$ $\bar x$

Wir berechnen die Ableitung

\frac{\partial v_{ich} (m^{*})}{\partial m} = \frac{\frac{\partial w_{ich} (m^{*})}{\partial m} \sum_{ich = 1}^{N} w_{ich} (m^{*}) - w_{ich} (m^{*}) \sum_{ich = 1}^{N} \frac{\partial w_{ich} (m^{*})}{\partial m}}{{(\sum_{ich = 1}^{N} w_{ich} (m^{*}))}^{2}}

$\frac{\partial v_i(m^*)}{\partial m} = \frac {\frac{\partial w_i(m^*)}{\partial m}\sum_{i=1}^{N} w_i(m^*)-w_i(m^*)\sum_{i=1}^{N}\frac{\partial w_i(m^*)}{\partial m}}{\left(\sum_{i=1}^{N} w_i(m^*)\right)^2}$

= \frac{1}{\sum_{ich = 1}^{N} w_{ich} (m^{*})} \cdot [\frac{\partial w_{ich} (m^{*})}{\partial m} - v_{ich} (m^{*}) \sum_{ich = 1}^{N} \frac{\partial w_{ich} (m^{*})}{\partial m}]

$=\frac 1{\sum_{i=1}^{N} w_i(m^*)}\cdot\left[\frac{\partial w_i(m^*)}{\partial m}-v_i(m^*)\sum_{i=1}^{N}\frac{\partial w_i(m^*)}{\partial m}\right]$ Dann

{EIN}^{'} (m^{*}) = \frac{1}{\sum_{ich = 1}^{N} w_{ich} (m^{*})} \cdot [\sum_{ich = 1}^{N} \frac{\partial w_{ich} (m^{*})}{\partial m} x_{ich} - (\sum_{ich = 1}^{N} \frac{\partial w_{ich} (m^{*})}{\partial m}) \sum_{ich = 1}^{N} v_{ich} (m^{*}) x_{ich}]

$A'(m^*) = \frac 1{\sum_{i=1}^{N} w_i(m^*)}\cdot\left[\sum_{i=1}^{N}\frac{\partial w_i(m^*)}{\partial m}x_i-\left(\sum_{i=1}^{N}\frac{\partial w_i(m^*)}{\partial m}\right)\sum_{i=1}^{N}v_i(m^*)x_i\right]$

= \frac{1}{\sum_{ich = 1}^{N} w_{ich} (m^{*})} \cdot [\sum_{ich = 1}^{N} \frac{\partial w_{ich} (m^{*})}{\partial m} x_{ich} - (\sum_{ich = 1}^{N} \frac{\partial w_{ich} (m^{*})}{\partial m}) m^{*}]

$=\frac 1{\sum_{i=1}^{N} w_i(m^*)}\cdot\left[\sum_{i=1}^{N}\frac{\partial w_i(m^*)}{\partial m}x_i-\left(\sum_{i=1}^{N}\frac{\partial w_i(m^*)}{\partial m}\right)m^*\right]$

und

| {EIN}^{'} (m^{*}) | < 1 \Rightarrow | \sum_{ich = 1}^{N} \frac{\partial w_{ich} (m^{*})}{\partial m} (x_{ich} - m^{*}) | < | \sum_{ich = 1}^{N} w_{ich} (m^{*}) | [3]

$|A'(m^*)| <1 \Rightarrow \left|\sum_{i=1}^{N}\frac{\partial w_i(m^*)}{\partial m}(x_i-m^*)\right| < \left|\sum_{i=1}^{N} w_i(m^*)\right| \qquad [3]$

wir haben

\begin{aligned} \frac{\partial w_{ich} (m^{*})}{\partial m} = & \frac{- ρ^{″} (| x_{ich} - m^{*} |) \cdot \frac{x_{ich} - m^{*}}{| x_{ich} - m^{*} |} | x_{ich} - m^{*} | + \frac{x_{ich} - m^{*}}{| x_{ich} - m^{*} |} ρ^{'} (| x_{ich} - m^{*} |)}{| x_{ich} - m^{*} |^{2}} \\ = \frac{x_{ich} - m^{*}}{| x_{ich} - m^{*} |^{3}} ρ^{'} (| x_{ich} - m^{*} |) - ρ^{″} (| x_{ich} - m^{*} |) \cdot \frac{x_{ich} - m^{*}}{| x_{ich} - m^{*} |^{2}} \\ = \frac{x_{ich} - m^{*}}{| x_{ich} - m^{*} |^{2}} \cdot [\frac{ρ^{'} (| x_{ich} - m^{*} |)}{| x_{ich} - m^{*} |} - ρ^{″} (| x_{ich} - m^{*} |)] \\ = \frac{x_{ich} - m^{*}}{| x_{ich} - m^{*} |^{2}} \cdot [w_{ich} (m^{*}) - ρ^{″} (| x_{ich} - m^{*} |)] \end{aligned}

$\begin{align}\frac{\partial w_i(m^*)}{\partial m} = &\frac{-\rho''(|x_i-m^*|)\cdot \frac {x_i-m^*}{|x_i-m^*|}|x_i-m^*|+\frac {x_i-m^*}{|x_i-m^*|}\rho'(|x_i-m^*|)}{|x_i-m^*|^2} \\ \\ &=\frac {x_i-m^*}{|x_i-m^*|^3}\rho'(|x_i-m^*|) - \rho''(|x_i-m^*|)\cdot \frac {x_i-m^*}{|x_i-m^*|^2} \\ \\ &=\frac {x_i-m^*}{|x_i-m^*|^2}\cdot \left[\frac {\rho'(|x_i-m^*|)}{|x_i-m^*|}-\rho''(|x_i-m^*|)\right]\\ \\ &=\frac {x_i-m^*}{|x_i-m^*|^2}\cdot \left[w_i(m^*)-\rho''(|x_i-m^*|)\right] \end{align}$

Einfügen in wir $[3]$

| \sum_{ich = 1}^{N} \frac{x_{ich} - m^{*}}{| x_{ich} - m^{*} |^{2}} \cdot [w_{ich} (m^{*}) - ρ^{″} (| x_{ich} - m^{*} |)] (x_{ich} - m^{*}) | < | \sum_{ich = 1}^{N} w_{ich} (m^{*}) |

$\left|\sum_{i=1}^{N}\frac {x_i-m^*}{|x_i-m^*|^2}\cdot \left[w_i(m^*)-\rho''(|x_i-m^*|)\right](x_i-m^*)\right| < \left|\sum_{i=1}^{N} w_i(m^*)\right|$

\Rightarrow | \sum_{ich = 1}^{N} w_{ich} (m^{*}) - \sum_{ich = 1}^{N} ρ^{″} (| x_{ich} - m^{*} |) | < | \sum_{ich = 1}^{N} w_{ich} (m^{*}) | [4]

$\Rightarrow \left|\sum_{i=1}^{N}w_i(m^*)-\sum_{i=1}^{N}\rho''(|x_i-m^*|)\right| < \left|\sum_{i=1}^{N} w_i(m^*)\right| \qquad [4]$

Dies ist die Bedingung, die erfüllt sein muss, damit der Festpunkt UAS ist. Da in unserem Fall die Straffunktion konvex ist, sind die Beträge positiv. Bedingung ist also äquivalent zu $[4]$

\sum_{ich = 1}^{N} ρ^{″} (| x_{ich} - m^{*} |) < 2 \sum_{ich = 1}^{N} w_{ich} (m^{*}) [5]

$\sum_{i=1}^{N}\rho''(|x_i-m^*|) < 2\sum_{i=1}^{N}w_i(m^*) \qquad [5]$

Wenn Huberts Verlustfunktion ist, dann haben wir einen quadratischen ( ) und einen linearen ( ) Zweig, $\rho(|x_i-m|)$ $q$ $l$

ρ (| x_{ich} - m |) = {\begin{cases} (1 / 2) | x_{ich} - m |^{2} | x_{ich} - m | \leq δ \\ δ (| x_{ich} - m | - δ / 2) | x_{ich} - m | > δ \end{cases}

$\rho(|x_i-m|)=\cases{ (1/2)|x_i- m|^2 \qquad\;\;\;\; |x_i-m|\leq \delta \\ \\ \delta\big(|x_i-m|-\delta/2\big) \qquad |x_i-m|> \delta}$

und

ρ^{'} (| x_{ich} - m |) = {\begin{cases} | x_{ich} - m | | x_{ich} - m | \leq δ \\ δ | x_{ich} - m | > δ \end{cases}

$\rho'(|x_i-m|)=\cases{ |x_i- m| \qquad |x_i-m|\leq \delta \\ \\ \delta \qquad \qquad \;\;\;\; |x_i-m|> \delta}$

ρ^{″} (| x_{ich} - m |) = {\begin{cases} 1 | x_{ich} - m | \leq δ \\ 0 | x_{ich} - m | > δ \end{cases}

$\rho''(|x_i-m|)=\cases{ 1\qquad |x_i-m|\leq \delta \\ \\ 0 \qquad |x_i-m|> \delta}$

{\begin{cases} w_{ich, q} (m) = 1 | x_{ich} - m | \leq δ \\ w_{ich, l} (m) = \frac{δ}{| x_{ich} - m |} < 1 | x_{ich} - m | > δ \end{cases}

$\cases{ w_{i,q}(m) =1\qquad \qquad \qquad |x_i-m|\leq \delta \\ \\ w_{i,l}(m) =\frac {\delta}{|x_i-m|} <1 \qquad |x_i-m|> \delta}$

Da wir nicht wissen, wie viele der uns in den quadratischen Zweig und wie viele in den linearen, zerlegen wir Bedingung als ( ) $|x_i-m^*|$ $[5]$ $N_q + N_l = N$

\sum_{ich = 1}^{N_{q}} ρ_{q}^{″} + \sum_{ich = 1}^{N_{l}} ρ_{l}^{″} < 2 [\sum_{ich = 1}^{N_{q}} w_{ich, q} + \sum_{ich = 1}^{N_{l}} w_{ich, l}]

$\sum_{i=1}^{N_q}\rho_q''+\sum_{i=1}^{N_l}\rho_l'' < 2\left[\sum_{i=1}^{N_q}w_{i,q} +\sum_{i=1}^{N_l}w_{i,l}\right]$

\Rightarrow N_{q} + 0 < 2 [N_{q} + \sum_{ich = 1}^{N_{l}} w_{ich, l}] \Rightarrow 0 < N_{q} + 2 \sum_{ich = 1}^{N_{l}} w_{ich, l}

$\Rightarrow N_q + 0 < 2\left[N_q +\sum_{i=1}^{N_l}w_{i,l}\right] \Rightarrow 0 < N_q+2\sum_{i=1}^{N_l}w_{i,l}$

was gilt. Für die Huber-Verlustfunktion ist der Fixpunkt des Algorithmus also unabhängig von den -Werten gleichmäßig asymptotisch stabil . Wir stellen fest, dass die erste Ableitung für jedes kleiner als die absolute Einheit ist , nicht nur für den festen Punkt. $x$ $m$

Was wir jetzt tun sollten, ist entweder zu beweisen, dass die UAS-Eigenschaft auch global ist, oder dass, wenn dann zur Nachbarschaft der Anziehungskraft von . $m(0) = \bar x$ $m(0)$ $m^*$

Alecos Papadopoulos
quelle

Danke für die Antwort. Geben Sie mir etwas Zeit, um diese Antwort zu analysieren.

Chris A.

Bestimmt. Immerhin wartete die Frage 20 Monate.

Alecos Papadopoulos

Ja, ich wurde an das Problem erinnert und beschloss, ein Kopfgeld zu zahlen. :)

Chris A.

Ich Glückspilz. Ich war vor 20 Monaten nicht dort - ich hätte diese Frage aufgegriffen, Kopfgeld oder nicht.

Alecos Papadopoulos

Vielen Dank für diese Antwort. Es sieht so aus, als hättest du dir das Kopfgeld verdient. Übrigens, Ihre Indizierung auf die Ableitung von wrt ist merkwürdig. Konnten die Summierungen in der zweiten Zeile dieses Dokuments keine andere Variable wie ?

v_{i}

$v_i$

m

$m$

j

$j$

Chris A.

Definition und Konvergenz iterativ neu gewichteter kleinster Quadrate

Antworten: