Definicja i zbieżność iteracyjnie ważonych najmniejszych kwadratów

Używałem iteracyjnie ponownie ważonych najmniejszych kwadratów (IRLS), aby zminimalizować funkcje następującej formy,

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

gdzie jest liczbą wystąpień , jest dokładnym oszacowaniem, którego chcę, a jest odpowiednią solidną funkcją kary. Powiedzmy, że jest wypukły (choć niekoniecznie ściśle) i na razie różnicowalny. Dobrym przykładem takiego jest funkcja straty Huber . $N$ $x_i \in \mathbb{R}$ $m \in \mathbb{R}$ $\rho$ $\rho$

To, co robiłem, to różnicowanie względem (i manipulowanie) w celu uzyskania, $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)$

i iteracyjnie rozwiązujemy ten problem, ustawiając go na 0 i ustawiając wagi przy iteracji na (uwaga, że postrzegany osobliwość w jest naprawdę zdejmowany osobliwość w ogóle „s może mnie obchodzi). Wtedy otrzymuję $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$

i rozwiązuję, aby uzyskać, . $m(k+1) = \frac{\sum_{i=1}^{N} w_i(k) x_i}{ \sum_{i=1}^{N} w_i(k)}$

Powtarzam ten algorytm stałoprzecinkowy aż do „konwergencji”. Zauważę, że jeśli dojdziesz do stałego punktu, jesteś optymalny, ponieważ twoja pochodna ma wartość 0 i jest to funkcja wypukła.

Mam dwa pytania dotyczące tej procedury:

Czy to standardowy algorytm IRLS? Po przeczytaniu kilku artykułów na ten temat (a były one bardzo rozproszone i niejasne co do tego, czym jest IRLS) jest to najbardziej spójna definicja algorytmu, jaką mogę znaleźć. Mogę publikować artykuły, jeśli ludzie tego chcą, ale tak naprawdę nie chciałem nikogo tutaj uprzedzać. Oczywiście możesz uogólnić tę podstawową technikę na wiele innych problemów związanych z wektorem $x_i$ i argumentami innymi niż $\left|x_i-m{(k)}\right|$ , pod warunkiem, że argument jest normą funkcji afinicznej twoich parametrów. Każda pomoc lub wgląd byłaby w tym świetna.
Wydaje się, że konwergencja działa w praktyce, ale mam na to kilka obaw. Jeszcze nie widziałem tego dowodu. Po kilku prostych symulacjach Matlaba widzę, że jedna iteracja tego nie jest mapowaniem skurczu (wygenerowałem dwa losowe wystąpienia i obliczenia i zobaczył, że czasami jest to więcej niż 1). Również mapowanie zdefiniowane przez kilka kolejnych iteracji nie jest ściśle mapowaniem skurczowym, ale prawdopodobieństwo, że stała Lipschitza będzie wyższa niż 1, staje się bardzo niskie. Czy istnieje więc prawdopodobieństwo mapowania skurczu w prawdopodobieństwie ? Jakiej maszyny użyłem, aby udowodnić, że to się zbiega? Czy to w ogóle się zbiega? $m$ $\frac{\left|m_1(k+1) - m_2(k+1)\right|}{\left|m_1(k)-m_2(k)\right|}$

Wszelkie wskazówki są pomocne.

Edycja: Podoba mi się artykuł na temat IRLS dotyczący rzadkiego odzyskiwania / wykrywania ściskania autorstwa Daubechies i in. 2008 „Iteracyjnie ponownie ważona minimalizacja najmniejszych kwadratów dla rzadkiego odzyskiwania” na arXiv. Ale wydaje się, że koncentruje się głównie na wagach dla problemów niewypukłych. Moja sprawa jest znacznie prostsza.

— Chris A.
źródło

Patrząc na stronę wiki IRWLS , walczę o różnicę między opisaną przez ciebie procedurą a IRWLS (używają po prostu jako szczególnej funkcji ). Czy możesz wyjaśnić, w jaki sposób uważasz, że proponowany algorytm różni się od IRWLS?

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

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

ρ

$\rho$

— user603

Nigdy nie powiedziałem, że było inaczej, a jeśli to sugerowałem, nie chciałem.

— Chris A.

Jeśli chodzi o twoje pierwsze pytanie, należy zdefiniować „standard” lub potwierdzić, że „model kanoniczny” został stopniowo ustanowiony. Jak wskazano w komentarzu, wydaje się, że przynajmniej sposób korzystania z IRWLS jest raczej standardowy.

Jeśli chodzi o twoje drugie pytanie, „mapowanie skurczu w prawdopodobieństwie” może być powiązane (choć nieformalnie) ze zbieżnością „rekurencyjnych algorytmów stochastycznych”. Z tego, co przeczytałem, istnieje ogromna literatura na ten temat głównie w inżynierii. W ekonomii używamy jej trochę, zwłaszcza przełomowych prac Lennarta Ljunga - pierwszą pracą był Ljung (1977) - która pokazała, że zbieżność (lub nie) rekurencyjnego algorytmu stochastycznego może być określona przez stabilność (lub not) powiązanego równania różniczkowego zwyczajnego.

(co zostało ponownie opracowane po owocnej dyskusji z PO w komentarzach)

Konwergencja

Użyję jako odniesienia Sabre Elaydi „An Introduction to Difference Equations”, 2005, 3d ed. Analiza jest uwarunkowana pewną próbką danych, więc są traktowane jako ustalone. $x's$

Warunek pierwszego rzędu minimalizacji funkcji celu, postrzegany jako funkcja rekurencyjna , $m$

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)]} [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]$

ma stały punkt (argmin funkcji celu). Według Twierdzenia 1.13 s. 27–28 Elaydiego, jeśli pierwsza pochodna w odniesieniu do RHS z , oceniona w punkcie stałym , oznacza to , jest mniejsza niż jedność w wartość bezwzględna, wtedy jest asymptotycznie stabilny (AS). Co więcej, w Twierdzeniu 4.3 p.179 rozumiemy, że oznacza to również, że punktem stałym jest jednolicie AS (UAS). „Asymptotycznie stabilny” oznacza, że dla pewnego zakresu wartości wokół stałego punktu sąsiedztwo , niekoniecznie małe, punkt stały jest atrakcyjny $m$ $[1]$ $m^*$ $A'(m^*)$ $m^*$
$(m^* \pm \gamma)$ , a więc jeśli algorytm podaje wartości w tym sąsiedztwie, zbiegnie się. Właściwość „jednolita” oznacza, że granica tego sąsiedztwa, a tym samym jego wielkość, jest niezależna od początkowej wartości algorytmu. Punkt stały staje się globalnie UAS, jeśli . Więc w naszym przypadku, jeśli to udowodnimy $\gamma = \infty$

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

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

udowodniliśmy właściwość UAS, ale bez globalnej konwergencji. Następnie możemy spróbować ustalić, czy sąsiedztwo przyciągania jest w rzeczywistości całymi rozszerzonymi liczbami rzeczywistymi, lub, że konkretna wartość początkowa, jaką stosuje PO, jak wspomniano w komentarzach (i jest to standard w metodologii IRLS), tj. Średnia z próby z „s, , zawsze należy do dzielnicy przyciągania punktu stałego. $x$ $\bar x$

pochodną

\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}}{{(\sum_{i = 1}^{N} w_{i} (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_{i = 1}^{N} w_{i} (m^{*})} \cdot [\frac{\partial w_{i} (m^{*})}{\partial m} - v_{i} (m^{*}) \sum_{i = 1}^{N} \frac{\partial w_{i} (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]$ Następnie

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

$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_{i = 1}^{N} w_{i} (m^{*})} \cdot [\sum_{i = 1}^{N} \frac{\partial w_{i} (m^{*})}{\partial m} x_{i} - (\sum_{i = 1}^{N} \frac{\partial w_{i} (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]$

| A^{'} (m^{*}) | < 1 \Rightarrow | \sum_{i = 1}^{N} \frac{\partial w_{i} (m^{*})}{\partial m} (x_{i} - m^{*}) | < | \sum_{i = 1}^{N} w_{i} (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]$

mamy

\begin{aligned} \frac{\partial w_{i} (m^{*})}{\partial m} = & \frac{- ρ^{″} (| x_{i} - m^{*} |) \cdot \frac{x_{i} - m^{*}}{| x_{i} - m^{*} |} | x_{i} - m^{*} | + \frac{x_{i} - m^{*}}{| x_{i} - m^{*} |} ρ^{'} (| x_{i} - m^{*} |)}{| x_{i} - m^{*} |^{2}} \\ = \frac{x_{i} - m^{*}}{| x_{i} - m^{*} |^{3}} ρ^{'} (| x_{i} - m^{*} |) - ρ^{″} (| x_{i} - m^{*} |) \cdot \frac{x_{i} - m^{*}}{| x_{i} - m^{*} |^{2}} \\ = \frac{x_{i} - m^{*}}{| x_{i} - m^{*} |^{2}} \cdot [\frac{ρ^{'} (| x_{i} - m^{*} |)}{| x_{i} - m^{*} |} - ρ^{″} (| x_{i} - m^{*} |)] \\ = \frac{x_{i} - m^{*}}{| x_{i} - m^{*} |^{2}} \cdot [w_{i} (m^{*}) - ρ^{″} (| x_{i} - 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}$

Wstawiamy to do mamy $[3]$

| \sum_{i = 1}^{N} \frac{x_{i} - m^{*}}{| x_{i} - m^{*} |^{2}} \cdot [w_{i} (m^{*}) - ρ^{″} (| x_{i} - m^{*} |)] (x_{i} - m^{*}) | < | \sum_{i = 1}^{N} w_{i} (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_{i = 1}^{N} w_{i} (m^{*}) - \sum_{i = 1}^{N} ρ^{″} (| x_{i} - m^{*} |) | < | \sum_{i = 1}^{N} w_{i} (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]$

Jest to warunek, który musi być spełniony, aby punktem stałym był UAS. Ponieważ w naszym przypadku funkcja kary jest wypukła, zaangażowane kwoty są dodatnie. Zatem warunek jest równoważny $[4]$

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

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

Jeśli jest funkcją straty Huberta, to mamy gałąź kwadratową ( ) i liniową ( ), $\rho(|x_i-m|)$ $q$ $l$

ρ (| x_{i} - m |) = {\begin{cases} (1 / 2) | x_{i} - m |^{2} | x_{i} - m | \leq δ \\ δ (| x_{i} - m | - δ / 2) | x_{i} - 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}$

ρ^{'} (| x_{i} - m |) = {\begin{cases} | x_{i} - m | | x_{i} - m | \leq δ \\ δ | x_{i} - 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_{i} - m |) = {\begin{cases} 1 | x_{i} - m | \leq δ \\ 0 | x_{i} - m | > δ \end{cases}

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

{\begin{cases} w_{i, q} (m) = 1 | x_{i} - m | \leq δ \\ w_{i, l} (m) = \frac{δ}{| x_{i} - m |} < 1 | x_{i} - 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}$

Ponieważ nie wiemy, ile zumieść nas w gałęzi kwadratowej, a ile w liniowej, rozkładamy warunek jako ( ) $|x_i-m^*|$ $[5]$ $N_q + N_l = N$

\sum_{i = 1}^{N_{q}} ρ_{q}^{″} + \sum_{i = 1}^{N_{l}} ρ_{l}^{″} < 2 [\sum_{i = 1}^{N_{q}} w_{i, q} + \sum_{i = 1}^{N_{l}} w_{i, 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_{i = 1}^{N_{l}} w_{i, l}] \Rightarrow 0 < N_{q} + 2 \sum_{i = 1}^{N_{l}} w_{i, 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}$

który trzyma. Zatem dla funkcji utraty Hubera stały punkt algorytmu jest jednakowo asymptotycznie stabilny, niezależnie od . Zauważmy, że pierwsza pochodna jest mniejsza niż jedność w wartości bezwzględnej dla dowolnego , nie tylko punktu stałego. $x$ $m$

To, co powinniśmy teraz zrobić, to albo udowodnić, że właściwość UAS jest również globalna, albo, że jeśli to należy do sąsiedztwa przyciągania . $m(0) = \bar x$ $m(0)$ $m^*$

— Alecos Papadopoulos
źródło

Dzięki za odpowiedzi. Daj mi trochę czasu na przeanalizowanie tej odpowiedzi.

— Chris A.

Na pewno. W końcu pytanie czekało 20 miesięcy.

— Alecos Papadopoulos

Tak, przypomniano mi o problemie i postanowiłem wystawić nagrodę. :)

— Chris A.

Mam szczęście. Nie było mnie tam 20 miesięcy temu - podjąłbym to pytanie, nagrodę czy nie.

— Alecos Papadopoulos

Dziękuję bardzo za tę odpowiedź. Wygląda na to, że jak dotąd zasłużyłeś na nagrodę. BTW, twoje indeksowanie na pochodnej wrt jest notacyjnie dziwne. Czy w podsumowaniach w drugim wierszu tej funkcji nie można użyć innej zmiennej, takiej jak ?

v_{i}

$v_i$

m

$m$

j

$j$

— Chris A.