Programowanie liniowe z ograniczeniami macierzowymi

Przegląd

Możesz wypróbować wariant metody mnożników alternatywnych kierunków (ADMM), która okazała się zaskakująco szybko problemów typu lasso . Strategia polega na sformułowaniu problemu za pomocą rozszerzonego Lagrangianu, a następnie przejściu gradientem na podwójny problem. Jest to szczególnie miłe w przypadku tego konkretnego problemu z regulacją ponieważ niepłaska część każdej iteracji metody ma dokładne rozwiązanie, które można po prostu ocenić element po elemencie, podczas gdy część gładka wymaga rozwiązania układu liniowego. $l_1$ $l^1$

W tym poście my

uzyskaj ogólną formułę ADMM w celu uogólnienia problemu,
wyłapać podproblemy dla każdej iteracji ADMM i specjalizować je w swojej sytuacji, a następnie
zbadania wynikowego układ liniowy, który musi być rozwiązany w każdej iteracji oraz opracowania szybkiego solver (lub przygotowującym) w oparciu o precomputing dekompozycji wartości własnej (lub niskich przybliżeń rangi ich) dla i . $M^TM$ $YY^T$
podsumuj kilkoma uwagami końcowymi

Większość wielkich pomysłów tutaj omówiono w następującym znakomitym artykule przeglądowym,

Boyd, Stephen i in. „Optymalizacja rozproszona i uczenie statystyczne za pomocą metody mnożników zmiennego kierunku”. Podstawy i trendy® w uczeniu maszynowym 3.1 (2011): 1-122. http://www.stanford.edu/~boyd/papers/pdf/admm_distr_stats.pdf

Zanim przejdę do szczegółów, chcę zauważyć, że jest to odpowiedź na metodę / algorytm, a nie praktyczna odpowiedź na istniejący kod - jeśli chcesz skorzystać z tej metody, musisz rzucić własną implementację.

Preparat ADMM

Ogólnie załóżmy, że chcesz rozwiązać

\begin{aligned} min_{x} & \sum_{i} | x_{i} | \\ s.t. & A x = b \end{aligned} .

$\begin{array}{rl} \min_{x} & \sum_{i} |x_i|\\ \textrm{s.t.} & Ax = b \end{array}.$

Problem z oryginalnym postem należy do tej kategorii po odpowiedniej wektoryzacji. (jest to tylko w zasadzie - zobaczymy, że wektoryzacja nie musi być wykonywana w praktyce)

Zamiast tego możesz rozwiązać równoważny problem, który ma Lagrangian

\begin{aligned} min_{x, z} & \sum_{i} | x_{i} | + \frac{α}{2} | | x - z | |^{2} + \frac{β}{2} | | A z - b | |^{2} \\ s.t. & A z = b \\ & & x = z, \end{aligned}

$\begin{array}{rl} \min_{x,z} & \sum_{i} |x_i| + \frac{\alpha}{2}||x-z||^2 + \frac{\beta}{2}||Az-b||^2 \\ \textrm{s.t.} & Az = b \\ \textrm{&} & x = z, \end{array}$

\begin{aligned} L (x, z, λ, γ) = & \sum_{i} | x_{i} | + \frac{α}{2} | | x - z | |^{2} + \frac{β}{2} | | A z - b | |^{2} + λ^{T} (A z - b) + γ^{T} (x - z) \\ = & \sum_{i} | x_{i} | + \frac{α}{2} | | x - z + \frac{1}{α} γ | |^{2} + \frac{β}{2} | | A z - b + \frac{1}{β} λ | |^{2} \\ + \frac{α}{2} | | \frac{1}{α} γ | |^{2} + \frac{β}{2} | | \frac{1}{β} λ | |^{2} . \end{aligned}

$\begin{align} L(x,z,\lambda,\gamma) =& \sum_{i} |x_i| + \frac{\alpha}{2}||x-z||^2 + \frac{\beta}{2}||Az-b||^2 + \lambda^T(Az-b) + \gamma^T(x-z) \\ =& \sum_{i} |x_i| + \frac{\alpha}{2}||x-z + \frac{1}{\alpha}\gamma||^2 + \frac{\beta}{2}||Az-b + \frac{1}{\beta}\lambda||^2 \\ &+ \frac{\alpha}{2}||\frac{1}{\alpha}\gamma||^2 + \frac{\beta}{2}||\frac{1}{\beta}\lambda||^2. \end{align}$

Metoda mnożników w naprzemiennym kierunku rozwiązuje podwójny problem, poprzez wznoszenie gradientu na podwójnych zmiennych, z wyjątkiem niedokładne naprzemienne wyświetlanie podwójnych podproblemów. Tj. się iterację

max_{λ, γ} min_{x, z} L (x, z, λ, γ),

$\max_{\lambda,\gamma} \min_{x,z} L(x,z,\lambda,\gamma),$

\begin{aligned} x^{k + 1} & = {a r g m i n}_{x} L (x, z^{k}, λ^{k}, γ^{k}) \\ z^{k + 1} & = {a r g m i n}_{z} L (x^{k + 1}, z, λ^{k}, γ^{k}) \\ γ^{k + 1} & = γ^{k} + α (x^{k + 1} - z^{k + 1}) \\ λ^{k + 1} & = λ^{k} + β (A z^{k + 1} - b) . \end{aligned}

$\begin{align} x^{k+1} &= \mathrm{argmin}_x L(x,z^k,\lambda^k,\gamma^k) \\ z^{k+1} &= \mathrm{argmin}_z L(x^{k+1},z,\lambda^k,\gamma^k) \\ \gamma^{k+1} &= \gamma^k + \alpha(x^{k+1}-z^{k+1}) \\ \lambda^{k+1} &= \lambda^k + \beta(Az^{k+1}-b). \end{align}$

W pewnych łagodnych warunkach dla parametrów i (wyjaśnionych w powiązanym powyżej dokumencie Boyda i Parikha) metoda ADMM zbiegnie się z prawdziwym rozwiązaniem. Współczynnik zbieżności jest liniowy, ponieważ stanowi rdzeń metody wynurzania gradientowego. Często można go przyspieszyć, aby był superliniowy, 1) zmieniając parametry i miarę przechodzenia na podstawie heurystyki lub 2) stosując przyspieszenie Nesterowa. Uwagi na temat zmiany parametrów kary znajdują się w artykule z ankietą Boyda, a na temat korzystania z przyspieszenia Niestierowa z ADMM patrz następujący artykuł: $\alpha$ $\beta$ $\alpha$ $\beta$

Goldstein, Tom, Brendan O'Donoghue i Simon Setzer. „Szybkie metody optymalizacji naprzemiennego kierunku”. Raport CAM (2012): 12–35. ftp://ftp.math.ucla.edu/pub/camreport/cam12-35.pdf

Jednak nawet jeśli ogólny wskaźnik zbieżności jest tylko liniowy, w przypadku problemów z zaobserwowano metodę bardzo szybkiego znajdowania wzoru rzadkości, a następnie zbiegania się wolniej na dokładnych wartościach. Ponieważ znalezienie wzoru rzadkości jest najtrudniejszą częścią, jest to bardzo przypadkowe! Dokładne powody, dla których wydają się być obszarem obecnych badań. Wszyscy widzą, że wzór rzadkości zbiega się szybko, ale wydaje się, że nikt nie wie dokładnie, dlaczego tak się dzieje. Jakiś czas temu zapytałem o to Boyda i Parikha, a Parikh pomyślał, że można to wyjaśnić, interpretując metodę w kontekście systemów sterowania. Inne heurystyczne wyjaśnienie tego zjawiska znajduje się w załączniku do następującego artykułu: $l^1$

Goldstein, Tom i Stanley Osher. „Podzielona metoda Bregmana dla problemów regulowanych przez L1”. SIAM Journal on Imaging Sciences 2.2 (2009): 323-343. ftp://ftp.math.ucla.edu/pub/camreport/cam08-29.pdf

Oczywiście teraz leży trudność w rozwiązywaniu i aktualizacji podproblemów do sytuacji konkretnego. Ponieważ Lagrangian jest kwadratowy w , podproblem aktualizacji wymaga po prostu rozwiązania układu liniowego. subproblem wydaje się trudniejsze, ponieważ jest nondifferentiable, ale okazuje się, że istnieje dokładny wzór na rozwiązanie, które można zastosować element po elemencie! Teraz omawiamy te podproblemy bardziej szczegółowo i określamy je dla problemu w oryginalnym poście. $x$ $z$ $z$ $z$ $x$

Konfiguracja podproblemu aktualizacji (system liniowy) $z$

Do aktualizacji mamy $z$

{a r g m i n}_{z} L (x_{k}, z, λ_{k}, γ_{k}) = {a r g m i n}_{z} \frac{α}{2} | | x - z + \frac{1}{α} γ | |^{2} + \frac{β}{2} | | A z - b + \frac{1}{β} λ | |^{2} .

$\mathrm{argmin}_z L(x_k,z,\lambda_k,\gamma_k) = \mathrm{argmin}_z \frac{\alpha}{2}||x-z + \frac{1}{\alpha}\gamma||^2 + \frac{\beta}{2}||Az-b + \frac{1}{\beta}\lambda||^2.$

Specjalizuje się w tym problemie:

\begin{aligned} {a r g m i n}_{Z_{J}, Z_{B}} & \frac{α}{2} | | J^{k + 1} - Z_{J} + \frac{1}{α} Γ_{J} | |_{F r o}^{2} + \frac{α}{2} | | B^{k + 1} - Z_{B} + \frac{1}{α} Γ_{B} | |_{F r o}^{2} \\ + \frac{β}{2} | | M Z_{J} + Z_{B} Y - X + \frac{1}{α} Λ | |_{F r o}^{2}, \end{aligned}

$\begin{align} \mathrm{argmin}_{Z_J,Z_B} &\frac{\alpha}{2}||J^{k+1}-Z_J + \frac{1}{\alpha}\Gamma_J||_{Fro}^2 + \frac{\alpha}{2}||B^{k+1}-Z_B + \frac{1}{\alpha}\Gamma_B||_{Fro}^2 \\ &+\frac{\beta}{2}||MZ_J + Z_BY - X + \frac{1}{\alpha}\Lambda||^2_{Fro}, \end{align}$

gdzie oznacza normę Frobenius (elementwise ). Jest to kwadratowy problem minimalizacji, w którym warunki optymalności pierwszego rzędu można znaleźć, biorąc częściowe pochodne celu w odniesieniu do i i ustawiając je na zero. To jest, $||\cdot||Fro$ $l_2$ $Z_J$ $Z_B$

\begin{aligned} 0 & = - \frac{α}{2} (J^{k + 1} - Z_{J} + \frac{1}{α} Γ_{J}) + \frac{β}{2} M^{T} (M Z_{J} + Z_{B} Y - X + \frac{1}{β} Λ), \\ 0 & = - \frac{α}{2} (B^{k + 1} - Z_{B} + \frac{1}{α} Γ_{B}) + \frac{β}{2} (M Z_{J} + Z_{B} Y - X + \frac{1}{β} Λ) Y^{T} . \end{aligned}

$\begin{align} 0 &= -\frac{\alpha}{2}(J^{k+1} - Z_J + \frac{1}{\alpha}\Gamma_J) + \frac{\beta}{2}M^T(MZ_J + Z_BY - X + \frac{1}{\beta}\Lambda), \\ 0 &= -\frac{\alpha}{2}(B^{k+1} - Z_B + \frac{1}{\alpha}\Gamma_B) + \frac{\beta}{2}(MZ_J + Z_BY - X + \frac{1}{\beta}\Lambda)Y^T. \end{align}$

Jak zauważono w komentarzach oryginalnego plakatu Justina Solomona, ten system dla jest symetryczny, więc gradient sprzężony jest idealną metodą bez matrycy. W dalszej części bardziej szczegółowo omówiono ten system i sposób jego rozwiązania / warunku wstępnego. $Z_J,Z_B$

Rozwiązywanie problemu aktualizacja podproblemu (analityczne rozwiązanie progowania) $x$

Teraz do podproblemu , $x$

{a r g m i n}_{x} L (x, z^{k}, λ^{k}, γ^{k}) = {a r g m i n}_{x} \sum_{i} | x_{i} | + \frac{α}{2} | | x - z^{k} + \frac{1}{α} γ^{k} | |^{2}

$\mathrm{argmin}_x L(x,z^k,\lambda^k,\gamma^k) = \mathrm{argmin}_x \sum_{i} |x_i| + \frac{\alpha}{2}||x-z^k + \frac{1}{\alpha}\gamma^k||^2$

Pierwszą rzeczą do zobaczenia jest to, że suma może być podzielona element po elemencie,

\sum_{i} | x_{i} | + \frac{α}{2} | | x - z^{k} + \frac{1}{α} γ^{k} | |^{2} = \sum_{i} | x_{i} | + \frac{α}{2} \sum_{i} (x_{i} - z_{i}^{k} + \frac{1}{α} γ_{i}^{k})^{2},

$\sum_{i} |x_i| + \frac{\alpha}{2}||x-z^k + \frac{1}{\alpha}\gamma^k||^2 = \sum_{i} |x_i| + \frac{\alpha}{2}\sum_i (x_i-z_i^k + \frac{1}{\alpha}\gamma_i^k)^2,$

Możemy więc rozwiązać problem optymalizacji element po elemencie równolegle, uzyskując

x_{i}^{k + 1} = {a r g m i n}_{x_{i}} | x_{i} | + \frac{α}{2} (x_{i} - z_{i}^{k} + \frac{1}{α} γ_{i}^{k})^{2} .

$x_i^{k+1} = \mathrm{argmin}_{x_i} |x_i| + \frac{\alpha}{2}(x_i-z_i^k + \frac{1}{\alpha}\gamma_i^k)^2.$

Ogólna postać tego równania to

min_{s} | s | + \frac{α}{2} (s - t)^{2} .

$\min_s |s| + \frac{\alpha}{2}(s-t)^2.$

Funkcja wartości bezwzględnej próbuje przyciągnąć optymalny punkt w kierunku , podczas gdy człon kwadratowy próbuje przyciągnąć optymalny punkt w kierunku . prawdziwe rozwiązanie leży więc gdzieś na segmencie między nimi, ze zwiększeniem zmierzającym do przyciągnięcia optymalnego punktu w kierunku , a zmniejszeniem przyciągającym optymalny punkt w kierunku . $s=0$ $s=t$ $[0,t)$ $\alpha$ $t$ $\alpha$ $0$

Jest to funkcja wypukła, ale nie można jej odróżnić od zera. Warunkiem punktu minimalizującego jest to, że sub pochodna celu w tym punkcie zawiera zero. Kwadratowy składnik ma pochodną , a funkcja wartości bezwzględnej ma pochodną dla , pochodną o ustalonej wartości jako przedział gdy , i pochodną dla . W ten sposób otrzymujemy pochodną dla ogólnej funkcji celu, $\alpha(s-t)$ $-1$ $s < 0$ $[-1,1]$ $s=0$ $1$ $s > 0$

\partial_{s} (| s | + \frac{α}{2} (s - t)^{2}) = {\begin{cases} 1 + α (s - t) & s > 0 \\ [- 1, 1] + α t, & s = 0, \\ - 1 + α (s - t), & s < 0. \end{cases}

$\partial_s \left(|s| + \frac{\alpha}{2}(s-t)^2\right) = \begin{cases} 1 + \alpha (s-t)\, & s > 0 \\ [-1,1] + \alpha t, & s = 0, \\ -1 + \alpha (s-t), & s < 0. \end{cases}$

Z tego widzimy, że pochodna celu przy zawiera wtedy i tylko wtedy, gdy , w którym to przypadku jest minimalizatorem. Z drugiej strony, jeśli nie jest minimalizatorem, wówczas możemy ustawić pochodną jednowartościową równą zero i rozwiązać dla minimalizatora. Wykonanie tego daje $s=0$ $0$ $|t| \le \frac{1}{\alpha}$ $s=0$ $s=0$

{a r g m i n}_{s} | s | + \frac{α}{2} (s - t)^{2} = {\begin{cases} t - \frac{1}{α}, & t > \frac{1}{α}, \\ 0, & | t | \leq \frac{1}{α}, \\ t + \frac{1}{α}, & t < - \frac{1}{α} \end{cases}

$\mathrm{argmin}_s |s| + \frac{\alpha}{2}(s-t)^2 = \begin{cases} t - \frac{1}{\alpha}, & t > \frac{1}{\alpha}, \\ 0, & |t| \le \frac{1}{\alpha}, \\ t + \frac{1}{\alpha}, & t < -\frac{1}{\alpha} \end{cases}$

Specjalizując ten wynik ponownie w prawdziwym problemie, który próbujemy rozwiązać w pierwotnym pytaniu, gdzie daje, Aktualizacja dla to po prostu $t = Z_{ij}^k - \frac{1}{\alpha}\Gamma_{ij}^k$

J_{i j}^{k + 1} = {\begin{cases} Z_{i j}^{k} - \frac{1}{α} Γ_{i j}^{k} - \frac{1}{α}, & Z_{i j}^{k} - \frac{1}{α} Γ_{i j}^{k} > \frac{1}{α}, \\ 0, & | Z_{i j}^{k} - \frac{1}{α} Γ_{i j}^{k} | \leq \frac{1}{α}, \\ Z_{i j}^{k} - \frac{1}{α} Γ_{i j}^{k} + \frac{1}{α}, & Z_{i j}^{k} - \frac{1}{α} Γ_{i j}^{k} < - \frac{1}{α} . \end{cases}

$J_{ij}^{k+1} = \begin{cases} Z_{ij}^k - \frac{1}{\alpha}\Gamma_{ij}^k - \frac{1}{\alpha}, & Z_{ij}^k - \frac{1}{\alpha}\Gamma_{ij}^k > \frac{1}{\alpha}, \\ 0, & |Z_{ij}^k - \frac{1}{\alpha}\Gamma_{ij}^k| \le \frac{1}{\alpha}, \\ Z_{ij}^k - \frac{1}{\alpha}\Gamma_{ij}^k + \frac{1}{\alpha}, & Z_{ij}^k - \frac{1}{\alpha}\Gamma_{ij}^k < -\frac{1}{\alpha}. \end{cases}$

B

$B$

B^{k + 1} = Z_{B} - \frac{1}{α} Γ_{B},

$B^{k+1} = Z_B - \frac{1}{\alpha}\Gamma_B,$

jak zauważono w oryginalnym plakacie Justina Solomona w komentarzach. Ogólnie rzecz biorąc, wykonanie aktualizacji dla wymaga tylko zapętlenia wpisów macierzy i oceny powyższych wzorów dla każdego wpisu. $J,B$

Uzupełnienie schura dla systemu $Z_J,Z_B$

Najdroższym krokiem iteracji jest rozwiązanie systemu,

\begin{aligned} 0 & = - \frac{α}{2} (J^{k + 1} - Z_{J} + \frac{1}{α} Γ_{J}) + \frac{β}{2} M^{T} (M Z_{J} + Z_{B} Y - X + \frac{1}{β} Λ), \\ 0 & = - \frac{α}{2} (B^{k + 1} - Z_{B} + \frac{1}{α} Γ_{B}) + \frac{β}{2} (M Z_{J} + Z_{B} Y - X + \frac{1}{β} Λ) Y^{T} . \end{aligned}

W tym celu warto trochę wysiłku skonstruować dobrego solvera / warunku wstępnego dla tego systemu. W tej sekcji robimy to poprzez wektoryzację , utworzenie dopełnienia Schur , wykonanie pewnych manipulacji produktem Krnoecker, a następnie cofnięcie wektoryzacji. Powstały układ dopełniacza Schura jest nieznacznie zmodyfikowanym równaniem Sylvester .

Poniższe dane dotyczące wektoryzacji i produktów Kronecker są absolutnie kluczowe:

$\mathrm{vec}(ABC) = (C^T \otimes A)\mathrm{vec}(B),$
$(A \otimes B)(C \otimes D) = AC \otimes BD$ ,
$(A \otimes B)^{-1} = A^{-1} \otimes B^{-1}$ , oraz
$(A \otimes B)^T = A^T \otimes B^T$ .

Tożsamości te zachowują się zawsze, gdy rozmiary macierzy i odwracalność są takie, że każda strona równania jest prawidłowym wyrażeniem.

Wektoryzowana forma systemu to:

(α I + β [\begin{matrix} I \otimes M^{T} M & (Y \otimes M)^{T} \\ Y \otimes M & Y Y^{T} \otimes I \end{matrix}]) [\begin{matrix} v e c (Z_{J}) \\ v e c (Z_{B}) \end{matrix}] = [\begin{matrix} v e c (α J + β M^{T} X + Γ_{J} - M^{T} Λ) \\ v e c (α B + β X Y^{T} + Γ_{B} - Λ Y^{T}) \end{matrix}],

$\left(\alpha I +\beta\begin{bmatrix}I \otimes M^TM & (Y \otimes M)^T \\ Y \otimes M & YY^T \otimes I\end{bmatrix}\right)\begin{bmatrix}\mathrm{vec}(Z_J) \\ \mathrm{vec}(Z_B)\end{bmatrix} = \begin{bmatrix}\mathrm{vec}(\alpha J + \beta M^TX + \Gamma_J - M^T\Lambda) \\ \mathrm{vec}(\alpha B + \beta XY^T + \Gamma_B - \Lambda Y^T)\end{bmatrix},$

lub

[\begin{matrix} I \otimes (α I + β M^{T} M) & β (Y \otimes M)^{T} \\ β Y \otimes M & (α I + β Y Y^{T}) \otimes I \end{matrix}] [\begin{matrix} v e c (Z_{J}) \\ v e c (Z_{B}) \end{matrix}] = [\begin{matrix} v e c (F) \\ v e c (G) \end{matrix}],

$\begin{bmatrix}I \otimes (\alpha I + \beta M^TM) & \beta (Y \otimes M)^T \\ \beta Y \otimes M & (\alpha I + \beta YY^T) \otimes I\end{bmatrix} \begin{bmatrix}\mathrm{vec}(Z_J) \\ \mathrm{vec}(Z_B)\end{bmatrix} = \begin{bmatrix}\mathrm{vec}(F) \\ \mathrm{vec}(G)\end{bmatrix},$

gdzie i są skróconą notacją dla prawej strony. Teraz przeprowadzamy eliminację bloku-gaussa / eliminację Schura, aby wyeliminować lewy dolny blok matrycy w procesie kondensacji produktów Kroneckera. To jest, $F$ $G$

[\begin{matrix} I \otimes (α I + β M^{T} M) & β (Y \otimes M)^{T} \\ 0 & (α I + β Y Y^{T}) \otimes I - β^{2} Y Y^{T} \otimes M (α I + β M^{T} M)^{- 1} M^{T} \end{matrix}] \dots \cdot [\begin{matrix} v e c (Z_{J}) \\ v e c (Z_{B}) \end{matrix}] = [\begin{matrix} v e c (F) \\ v e c (G) - β Y \otimes M (α I + β M^{T} M)^{- 1} v e c (F) \end{matrix}] .

$\begin{bmatrix}I \otimes (\alpha I + \beta M^TM) & \beta (Y \otimes M)^T \\ 0 & (\alpha I + \beta YY^T) \otimes I - \beta^2 YY^T \otimes M(\alpha I + \beta M^TM)^{-1} M^T\end{bmatrix} \dots \\ \cdot \begin{bmatrix}\mathrm{vec}(Z_J) \\ \mathrm{vec}(Z_B)\end{bmatrix} = \begin{bmatrix}\mathrm{vec}(F) \\ \mathrm{vec}(G) - \beta Y \otimes M(\alpha I + \beta M^TM)^{-1}\mathrm{vec}(F)\end{bmatrix}.$

Bez wektorowania dwa równania, które musimy rozwiązać kolejno, to:

$Z_{B} (α I + β Y Y^{T}) - (β M (α I + β M^{T} M)^{- 1} M^{T}) Z_{B} (β Y Y^{T}) \dots = G - β M (α I + β M^{T} M)^{- 1} F Y^{T}$ $Z_B (\alpha I + \beta YY^T) - (\beta M (\alpha I + \beta M^TM)^{-1} M^T)Z_B(\beta YY^T) \dots \\ = G - \beta M (\alpha I + \beta M^TM)^{-1} F Y^T$
$(α I + β M^{T} M) Z_{J} = F - β M^{T} Z_{B} Y .$ $(\alpha I + \beta M^TM) Z_J = F - \beta M^T Z_B Y.$

Rozwiązanie układu dopełniacza Schur, gdy są kwadratowe, mają wysoką rangę $Y,M$

W tej sekcji rozwiązujemy układ dopełniacza Schura dla (równanie 1. powyżej), stosując wstępnie obliczone pełne SVD macierzy i stosując zmodyfikowaną wersję algorytmu Bartelsa-Stewarta dla Sylvester równanie. Algorytm został nieco zmodyfikowany w stosunku do standardowej wersji, aby uwzględnić dodatkowy w drugim członie, co sprawia, że nie jest to równanie Sylvester. Po za pomocą pierwszego równania, można łatwo znaleźć z drugiego równania. Drugie równanie jest łatwe do rozwiązania za pomocą dowolnej metody. $Z_B$ $YY^T, MM^T, M^TM$ $\beta YY^T$ $Z_B$ $Z_J$

Ta metoda wymaga wstępnego obliczenia dwóch pełnych plików SVD przed rozpoczęciem procesu ADMM, ale szybko można ją zastosować w rzeczywistych iteracjach ADMM. Ponieważ metoda dotyczy pełnych SVD macierzy więzów, właściwe jest, gdy są one zbliżone do kwadratu i wysokiej rangi. Możliwa jest również bardziej skomplikowana metoda wykorzystująca SVD niskiej rangi, ale została przedstawiona w dalszej części.

Metoda rozwija się w następujący sposób. Niech oznaczają Precomputed pełne dekompozycje wartości pojedyncze, i skroplenia prawą stronę za . Wtedy pierwszym równaniem jest: Mnożenie przez czynniki ortogonalne, aby wyczyścić lewą i prawą stronę i ustawić nowy tymczasowy nieznany , to dalej staje się,

Q D Q^{T} = Y Y^{T}, W Σ W^{T} = M M^{T}, V T V^{T} = M^{T} M

$Q D Q^T = YY^T, \\ W\Sigma W^T = MM^T, \\ VTV^T = M^TM$

H

$H$

Z_{B} Q (α I + D) Q^{T} - W β Σ (α I + Σ)^{- 1} Σ W^{T} Z_{B} Q D Q^{T} = H .

$Z_B Q (\alpha I + D) Q^T - W \beta \Sigma (\alpha I + \Sigma)^{-1}\Sigma W^T Z_B Q D Q^T = H.$

A = W^{T} Z_{B} Q

$A = W^T Z_B Q$

A (α I + D) - β Σ (α I + Σ)^{- 1} Σ A D = W H Q^{T} .

$A (\alpha I + D) - \beta \Sigma (\alpha I + \Sigma)^{-1}\Sigma A D = W H Q^T.$

Teraz możemy znaleźć , rozwiązując układ diagonalny , $A$

((α I + D) \otimes I + D \otimes β Σ (α I + Σ)^{- 1} Σ) v e c (A) = v e c (W H Q^{T}) .

$\left((\alpha I + D) \otimes I + D \otimes \beta \Sigma (\alpha I + \Sigma)^{-1}\Sigma \right)\mathrm{vec}(A) = \mathrm{vec}(W H Q^T).$

Po znalezieniu obliczamy i znając rozwiązujemy drugie równanie powyżej dla , co jest trywialne, ponieważ mamy już rozkład wartości własnej dla . $A$ $Z_B = W A Q^T$ $Z_B$ $Z_J$ $M^TM$

Koszt początkowy polega na obliczeniu dwóch symetrycznych dodatnich określonych wartości własnych i , a następnie koszt za iterację dla pełnego rozwiązania jest zdominowany przez garstkę mnożenia macierzy-macierzy, która jest tego samego rzędu wielkość jako wykonanie 1 podrzędności CG. Jeśli początkowe dekompozycje wartości własnych są zbyt kosztowne, można je obliczyć niedokładnie, na przykład przez wcześniejsze zakończenie iteracji Lanczosa i utrzymanie największych wektorów własnych. Następnie metodę tę można zastosować jako dobry warunek wstępny dla CG, a nie jako bezpośredni solver. $M^TM$ $YY^T$

Metoda rozwiązania, gdy są bardzo prostokątne lub mają przybliżenie niskiej rangi $M,Y$

Teraz zwracamy uwagę na rozwiązywanie lub kondycjonowanie gdy albo a) macierze wejściowe są bardzo prostokątne - co oznacza, że mają znacznie więcej wierszy niż kolumn lub odwrotnie - lub b) mają przybliżenie niskiej rangi. Poniższe wyprowadzenie obejmuje szerokie zastosowanie formuły Woodbury'ego, uzupełnienia Schura i innych podobnych manipulacji. $Z_J,Z_B$ $M,Y$

Zaczynamy od naszego systemu uzupełniającego Schur,

(α I + β Y Y^{T}) \otimes I - β^{2} Y Y^{T} \otimes M (α I + β M^{T} M)^{- 1} M^{T} .

$(\alpha I + \beta YY^T) \otimes I - \beta^2 YY^T \otimes M(\alpha I + \beta M^TM)^{-1} M^T.$

Kilka manipulacji przekształca ten system w bardziej symetryczną formę,

(α I + β I \otimes M M^{T} + β Y Y^{T} \otimes I) v e c (Z_{B}) = (I \otimes (I + \frac{β}{α} M M^{T})) v e c (H) .

$(\alpha I + \beta I \otimes MM^T + \beta YY^T \otimes I)\mathrm{vec}(Z_B) = \left(I \otimes (I + \frac{\beta}{\alpha}MM^T)\right)\mathrm{vec}(H).$

Teraz przybliżamy niskie rangi. Niech będzie przybliżonym zmniejszeniem SVD lub przybliżeniem niskiego stopnia i ( jest symbolem zastępczym i nie jest używany). Podstawienie ich do naszego układu daje następującą macierz odwrotną, którą chcemy zastosować,

Q D^{1 / 2} Q_{2}^{T} = Y W Σ^{1 / 2} V^{T} = M

$Q D^{1/2} Q_2^T = Y \\ W \Sigma^{1/2} V^T = M$

Y

$Y$

M

$M$

Q_{2}

$Q_2$

(α I + β I \otimes W Σ W^{T} + β Y Y^{T} \otimes I)^{- 1} .

$(\alpha I + \beta I \otimes W \Sigma W^T + \beta YY^T \otimes I)^{-1}.$

Ponieważ macierz, którą mamy do odwrócenia, jest aktualizacją tożsamości niskiej rangi, logiczną strategią jest próba użycia formuły Woodbury'ego

(A + U C U^{T})^{- 1} = A^{- 1} - A^{- 1} U (C^{- 1} + U^{T} A^{- 1} U)^{- 1} U^{T} A^{- 1} .

$(A + UCU^T)^{-1} = A^{-1} - A^{-1}U(C^{-1}+U^TA^{-1}U)^{-1}U^TA^{-1}.$

Jednak należy zachować ostrożność, ponieważ elementy niskiej rangi i nie są ortogonalne. Tak więc, aby zastosować formułę Woodbury, zbieramy zarówno aktualizacje niskiej rangi w jedną dużą aktualizację. Zrób to i stosując formułę Woodbury, $I \otimes W$ $Y \otimes I$

{(\frac{1}{α} I + β [\begin{matrix} I \otimes W & Q \otimes I \end{matrix}] [\begin{matrix} I \otimes Σ \\ D \otimes Y \end{matrix}] [\begin{matrix} I \otimes Σ^{T} \\ Q^{T} \otimes I \end{matrix}])}^{- 1} = α I - \frac{β}{α^{2}} [\begin{matrix} I \otimes W & Q \otimes I \end{matrix}] {[\begin{matrix} I \otimes (Σ^{- 1} + \frac{β}{α} I) & \frac{β}{α} Q \otimes W^{T} \\ \frac{β}{α} Q^{T} \otimes W & (D^{- 1} + \frac{β}{α} I) \otimes Y \end{matrix}]}^{- 1} [\begin{matrix} I \otimes Σ^{T} \\ Q^{T} \otimes I \end{matrix}] .

$\left(\frac{1}{\alpha} I + \beta \begin{bmatrix}I\otimes W & Q \otimes I\end{bmatrix}\begin{bmatrix}I \otimes \Sigma & \\ & D \otimes Y\end{bmatrix}\begin{bmatrix}I \otimes \Sigma^T \\ Q^T \otimes I\end{bmatrix}\right)^{-1} \\ = \alpha I - \frac{\beta}{\alpha^2}\begin{bmatrix}I\otimes W & Q \otimes I\end{bmatrix}\begin{bmatrix}I \otimes (\Sigma^{-1}+\frac{\beta}{\alpha}I) & \frac{\beta}{\alpha}Q \otimes W^T\\ \frac{\beta}{\alpha}Q^T\otimes W & (D^{-1} + \frac{\beta}{\alpha}I) \otimes Y\end{bmatrix}^{-1}\begin{bmatrix}I \otimes \Sigma^T \\ Q^T \otimes I\end{bmatrix}.$

Odwrotność rdzenia można obliczyć według odwrotnej formuły odwrotnej 2x2,

{[\begin{matrix} A & B \\ B^{T} & C \end{matrix}]}^{- 1} = [\begin{matrix} (A - B C^{- 1} B^{T})^{- 1} & - A^{- 1} B (C - B^{T} A^{- 1} B)^{- 1} \\ - C^{- 1} B^{T} (A - B C^{- 1} B^{T})^{- 1} & (C - B^{T} A^{- 1} B)^{- 1} \end{matrix}] .

$\begin{bmatrix}A & B \\ B^T & C\end{bmatrix}^{-1} = \begin{bmatrix}(A-BC^{-1}B^T)^{-1} & -A^{-1}B(C-B^TA^{-1}B)^{-1} \\ -C^{-1}B^T(A-BC^{-1}B^T)^{-1} & (C-B^TA^{-1}B)^{-1}\end{bmatrix}.$

Ten post jest już wystarczająco długi, więc oszczędzę obszernych szczegółów obliczeń, ale efekt końcowy jest taki, że podłączenie niezbędnych podmacierzy do odwrotnego bloku i pomnożenie wszystkiego przez to daje następującą jawną formę dla ogólnego odwrotności,

(α I + β I \otimes M M^{T} + β Y Y^{T} \otimes I)^{- 1} = \frac{1}{α} I - \frac{β}{α^{2}} (t_{11} + s_{11} + t_{12} + s_{12} + t_{21} + s_{21} + t_{22} + s_{22}),

$(\alpha I + \beta I \otimes MM^T + \beta YY^T \otimes I)^{-1} = \frac{1}{\alpha} I - \frac{\beta}{\alpha^2}(t_{11} + s_{11} + t_{12} + s_{12} + t_{21} + s_{21} + t_{22} + s_{22}),$

gdzie

\begin{aligned} t_{11} & = \frac{α}{β} I \otimes W l^{- 1} W^{T} \\ s_{11} & = (Q \otimes W l^{- 1}) D_{11} (Q^{T} \otimes l^{- 1} W^{T}) \\ t_{12} & = - \frac{α}{β} Q h^{- 1} Q^{T} \otimes W l^{- 1} W^{T} \\ s_{12} & = - (Q h^{- 1} \otimes W l^{- 1}) D_{22} (h^{- 1} Q^{T} \otimes W^{T}) \\ t_{21} & = t_{12} \\ s_{21} & = - (Q h^{- 1} \otimes W) D_{22} (h^{- 1} Q^{T} \otimes l^{- 1} W^{T}) \\ t_{22} & = \frac{α}{β} Q h^{- 1} Q^{T} \otimes I \\ s_{22} & = (Q h^{- 1} \otimes W) D_{22} (h^{- 1} Q^{T} \otimes W^{T}) \\ D_{11} & = \frac{α}{β} {(h \otimes I - I \otimes l^{- 1})}^{- 1} \\ D_{22} & = \frac{α}{β} {(I \otimes l - h^{- 1} \otimes I)}^{- 1} \\ l & = \frac{α}{β} Σ^{- 1} + I \\ h & = \frac{α}{β} D^{- 1} + I . \end{aligned}

$\begin{align} t_{11} &= \frac{\alpha}{\beta}I \otimes W l^{-1} W^T \\ s_{11} &= (Q \otimes W l^{-1})D_{11}(Q^T \otimes l^{-1}W^T) \\ t_{12} &= -\frac{\alpha}{\beta} Q h^{-1} Q^T \otimes W l^{-1} W^T \\ s_{12} &= -(Q h^{-1} \otimes W l^{-1})D_{22}(h^{-1} Q^T \otimes W^T) \\ t_{21} &= t_{12} \\ s_{21} &= -(Q h^{-1} \otimes W)D_{22}(h^{-1} Q^T \otimes l^{-1} W^T) \\ t_{22} &= \frac{\alpha}{\beta}Q h^{-1} Q^T \otimes I \\ s_{22} &= (Q h^{-1} \otimes W)D_{22}(h^{-1}Q^T \otimes W^T) \\ D_{11} &= \frac{\alpha}{\beta}\left(h \otimes I - I \otimes l^{-1} \right)^{-1} \\ D_{22} &= \frac{\alpha}{\beta}\left(I \otimes l - h^{-1} \otimes I \right)^{-1} \\ l &= \frac{\alpha}{\beta} \Sigma^{-1} + I \\ h &= \frac{\alpha}{\beta} D^{-1} + I. \end{align}$

W tej formie możemy zastosować odwrotność i znaleźć po semestrze przez 8 lewych i prawych kanapek mnożenia macierzy. Ogólna formuła stosowania sumy produktów Kronecker to: $Z_B$

((A_{1} \otimes B_{1}) + (A_{2} \otimes B_{2}) + \dots) v e c (C) = v e c (B_{1}^{T} C A_{1} + B_{2}^{T} C A_{2} + \dots) .

$\left((A_1 \otimes B_1) + (A_2 \otimes B_2) + \dots\right)\mathrm{vec}(C) = \mathrm{vec}(B_1^T C A_1 + B_2^T C A_2 + \dots ).$

Zauważ, że wszystkie jawne odwrotności, na których się skończyliśmy, są ukośne, więc nie ma nic do „rozwiązania”.

Liniowy kod solvera

W Matlabie zaimplementowałem powyższe dwa . Wydaje się, że działa dobrze. Kod solvera jest tutaj. $z_J,Z_B$

https://github.com/NickAlger/MeshADMM/blob/master/zkronsolve.m

Skrypt testowy do sprawdzania działania solverów jest tutaj. Pokazuje również przykładowo, jak wywołać kod solvera.

https://github.com/NickAlger/MeshADMM/blob/master/test_zkronsolve.m

Uwagi końcowe

Metody typu ADMM dobrze nadają się do takich problemów, ale konieczne byłoby wdrożenie własnej implementacji. Ogólna struktura metody jest dość prosta, więc implementacja nie jest zbyt trudna w czymś takim jak MATLAB.

Brakującym elementem tego postu, który należałoby określić, aby w pełni zdefiniować metodę problemu, jest wybór parametrów karnych . Na szczęście metoda jest ogólnie dość solidna, o ile wartości parametrów nie są szalone. W artykule Boyda i Parikha znajduje się sekcja o parametrach karnych, podobnie jak odniesienia do nich, ale po prostu eksperymentowałbym z tymi parametrami, aż uzyskasz rozsądne wskaźniki zbieżności. $\alpha,\beta$

Przedstawione strategie solvera są wysoce skuteczne, jeśli macierze ograniczeń są albo a) gęste, kwadratowe i mają wysoką rangę, albo b) mają dobre przybliżenie niskiej rangi. Inną użyteczną Solver które mogą być temat pracy przyszłości będzie Solver zoptymalizowane w następnym przypadku - matryca ograniczenie jest niewielkie i squareish i wysoki stopień, ale istnieje dobra prekondycjonerze do . Byłoby tak, gdyby na przykład był zdyskretowanym Laplacianinem. $Z_J,Z_B$ $M$ $\alpha I + MM^T$ $M$

— Nick Alger
źródło

Wdrażanie tego teraz! Aby sprawdzić, rozwiązanie macierzy dla i powinno być określone symetrycznie / dodatnio, ponieważ pochodzi z najmniejszych kwadratów, prawda? To empirycznie wydaje się prawdą :-). Czy więc CG jest lepszą opcją niż GMRES?

Z_{B}

$Z_B$

Z_{J}

$Z_J$

— Justin Solomon

Myślę też, że aktualizacja B jest nieprawidłowa? Pracuję nad tym bardziej szczegółowo, ale przypominam, że B nie pojawia się w mojej funkcji energetycznej (brak term), więc nie jestem pewien, czy powinien przyjmować wartości tylko w Czy myślę o tym źle? Dzięki!

| B |

$|B|$

\pm (1 - 1 / α) .

$\pm (1-1/\alpha).$

— Justin Solomon

[raczej, ]

B = Z_{B} - Γ_{B} / α

$B = Z_B-\Gamma_B/\alpha$

— Justin Solomon

Niesamowity! Po wprowadzeniu własnych formuł dla i (prawdopodobnie bliskie / równoważne z tym, co opublikowałeś, ale coś nie działało), jest to znacznie lepsze niż metoda IRLS. Dzięki!

J

$J$

B

$B$

— Justin Solomon

Dobre wieści. Tak miło widzieć, kiedy wkład tutaj prowadzi do prawdziwych rezultatów.

— Michael Grant,

Czy stać cię na te SVD, o których wspomniał Geoffrey Irving? Jeśli możesz, rozważę iteracyjnie przeważone podejście metodą najmniejszych kwadratów (IRLS) . To podejście rozwiązałoby problemy z formą gdzie jest macierzą wag.

\begin{array}{ll} minimize & \sum_{i j} W_{i j} J_{i j}^{2} \\ subject to & M J + B Y = X \end{array}

$\begin{array}{ll}\text{minimize}&\sum_{ij} W_{ij}J_{ij}^2\\\text{subject to}&MJ+BY=X\end{array}$

W

$W$

Iteracje zaczynają się od jako macierzy wszystkich; daje to optymalną . Iteracje następują od gdzie jest małą stałą, która zapobiega dzieleniu przez zero. Nie jestem całkowicie pewien kryteriów konwergencji, ale być może link do Wikipedii, który przedstawiłem powyżej, może zawierać odniesienia. $W^{(0)}$ $J^{(0)}$

W_{i j}^{(k + 1)} = {| max {J_{i j}^{(k)}, ϵ} |}^{- 1}

$W_{ij}^{(k+1)}=\left|\max\{J_{ij}^{(k)},\epsilon\}\right|^{-1}$

ϵ

$\epsilon$

Możesz również rozważyć wygładzoną metodę pierwszego rzędu. TFOCS, którego jestem współautorem, może sobie z tym poradzić, używając swojego „wygładzonego stożkowego podwójnego” solvera (SCD), ale nie będzie tak łatwy w użyciu.

Jeśli chcesz wypróbować macierzową metodę punktu wewnętrznego, przeczytaj pracę Jacka Gondzio.

EDYCJA: hmm, może się zdarzyć, że IRLS nie będzie mógł użyć SVD do obliczenia rozwiązań. Jeśli tak, wrócę do jednej z innych opcji.

— Michael Grant
źródło

Nie jestem pewien, czy mógłbym tutaj użyć SVD, ale IRLS to świetny pomysł, niezależnie od tego! Prędkość nie jest tak ważna, jak pamięć, i zawstydzająco użyłem IRLS do pokrewnych badań kilka miesięcy temu i zadziałało świetnie (kopiąc siebie za to, że nie próbowałem tego wcześniej!). Nawet bez SVD dla IRLS powinno być to możliwe przy użyciu solwera liniowego, takiego jak CG, który nie wymaga pełnego systemu. W rzeczywistości CG można prawdopodobnie zatrzymać za pomocą dość luźnych ograniczeń przed dostosowaniem jak sugerujesz. Patrząc również na podejście ADMM, ale mam z tym mniej doświadczenia.

W_{i j}

$W_{ij}$

— Justin Solomon

Tak, ADMM też byłoby świetne. Właściwie napisałem sekcję sugerującą całkowite wyeliminowanie Y, ale później zobaczyłem, że nie było kwadratowe.

M

$M$

— Michael Grant

Wdrożono strategię IRLS - jest zbieżna, ale liczbowo nie radzi sobie zbyt dobrze, ponieważ układ liniowy, który musi rozwiązać, jest źle uwarunkowany dzięki szerokiemu zakresowi ; za pomocą GMRES do rozwiązania systemu. Spróbuję ADMM w następnej kolejności!

w

$w$

— Justin Solomon