Tło matematyczne do PDLP

Ta strona zawiera matematyczne tło dla programistów i zaawansowanych użytkowników PDLP – to rozwiązanie do programowania liniowego i kwadratowego dostępnego w Narzędziach z operatorem LUB. Służy jako odniesienie do części kodu i nie należy go odczytywać samodzielnie. Zainteresowani czytelnicy powinni najpierw zapoznać się z publikacją „Practical Large-Scale Linear Programming using Primal-Dual Hybrid Gradient”, a następnie zapoznać się z kodem, a potem wrócić do tego dokumentu, gdy się w nim odwołuje.

Podstawowe

PDLP uwzględnia następujący problem z programowaniem kwadratowym wypukłym:

$$ \begin{align} \min_x & \, c^Tx + \frac{1}{2}x^TQx \\ \text{s.t.}\; & l^{c} \le Ax \le u^{c} \\ & l^{v} \le x \le u^{v} \end{align} $$

gdzie $A$ jest macierz $m \times n$, a $Q$ jest ukośną, nieujemną macierz $n \times n$ 1. Górne wektory $u^{c}$ i $u^{v}$ mają wpisy w $\mathbb{R} \cup \{ \infty\}$, a dolne wektory $l^{c}$ i $l^{v}$ mają wpisy w $\math$l} \cup \{ \cup \{ -}.

Dual

Niech $a \in \mathbb{R}$. Niech $[a]_+$ oznacza dodatnią część, a $[a]_-$ oznacza jej ujemną część, czyli $a = [a]_+ - [a]_-$. W przypadku zastosowania do wektora części dodatnie i ujemne są obliczane według elementu.

Wartość podwójnym wcześniejszego zadania pierwotnego to ponad $x \in \mathbb{R}^n$, $y \in \mathbb{R}^m$ i $r \in \mathbb{R}^n$. Wektor $y$ zawiera podwójne $multipliers w ograniczeniach liniowych ($l^{c} \le Ax).

$$ \begin{align} \max_{x, y, r} & \, -\frac{1}{2}x^TQx + \left((l^{c})^T[y]_+ - (u^{c})^T[y]_- \right) + \left((l^{v})^T[r]_+ - (u^{v})^T[r]_- \right) \\ \text{s.t.}\; & Qx + c - A^Ty = r \end{align} $$

Gdy $Q = 0$, $x$ można usunąć z dualu, przywracając dwójkę LP.

Podwójne granice zmiennych

Mówimy, że element $y$ spełnia granice podwójnej zmiennej, jeśli wartość parametru $y$-term jest w celu ograniczona, czyli:

$$ y_i \geq 0 \qquad \text{if }u^c_i = \infty, \\ y_i \leq 0 \qquad \text{if }l^c_i = -\infty. $$

Wyprowadź przy użyciu podwójnych sprzężonych

Eliminacje

Ustaw $a \in \mathbb{R} \cup \{-\infty\}$ i $b \in \mathbb{R} \cup \{\infty\}$ z $b \ge a$ i weź pod uwagę interwał $[a, b] \subseteq \ma$bby{R} \cup \{

Niech $\mathcal{I}_{[a, b]} : \mathbb{R} \to \mathbb{R} \cup \{ \infty\}$ będzie funkcją wskaźnika interwału, czyli $\mathcal{I}_{[a, b]}(x)$y, jeśli $x, \bin [

Zdefiniuj $p(y; a, b): \mathbb{R} \to \mathbb{R} \cup \{\infty\}$ jako:

$$ p(y; a, b) = \begin{cases} ay & \text{ if } y < 0 \\ 0 & \text{ if } y = 0 \\ by & \text{ if } y > 0 \end{cases}. $$

Jeśli $a$ lub $b$ mają wartość nieskończoną, korzystaj ze standardowej rozszerzonej arytmetyki rzeczywistej.

Wynik podstawowy: $p(y; a, b) = (\mathcal{I}_{[a, b]})^*(y)$, gdzie $(\cdot)^*$ oznacza sprzężenie wypukłe.

W przypadku wektorów $l \subseteq (\mathbb{R} \cup \{-\infty\})^n$ and $u \subseteq (\mathbb{R} \cup \{\infty\})^n$, funkcja wskaźnika $\mathcal{I}_${[l, u]bb \

Derywacja

Wprowadzamy zmienne pomocnicze $\tilde a \in \mathbb{R}^m$ oraz $\tilde x \in \mathbb{R}^n$, przedstawimy główny problem w następujący sposób:

$$ \begin{align} \min_{x, \tilde x, \tilde a} & \, c^Tx + \frac{1}{2}x^TQx + \mathcal{I}_{[l^c,u^c]}(\tilde a) + \mathcal{I}_{[l^v,u^v]}(\tilde x) \\ \text{s.t.}\; & \tilde a = Ax \\ & \tilde x = x \end{align} $$

Podnosząc w ten sposób ograniczenia równości, otrzymujemy:

$$ \min_{x, \tilde x, \tilde a} \max_{y, r} c^Tx + \frac{1}{2}x^TQx + y^T\tilde a - y^TAx + r^T\tilde x - r^Tx + \mathcal{I}_{[l^c,u^c]}(\tilde a) + \mathcal{I}_{[l^v,u^v]}(\tilde x) $$

Wymiana minimalnej z wartością maksymalną i przegrupowanie:

$$ \max_{y, r} \min_{x, \tilde x, \tilde a} c^Tx + \frac{1}{2}x^TQx - y^TAx - r^Tx + \left( \mathcal{I}_{[l^c,u^c]}(\tilde a) + y^T\tilde a \right) + \left(\mathcal{I}_{[l^v,u^v]}(\tilde x) + r^T\tilde x \right) $$

Wspólna minimalizacja w skali $x$, $\tilde x$ i $\tilde a$ się rozkłada. Dla $x$ widzimy, że minimalizator, o ile istnieje, spełnia warunki $Qx + c – A^Ty = r$, w którym przypadku minimalna wartość wynosi $-\frac{1}{2} x^TQx$. W przypadku $\tilde x$ i $\tilde a$ stosujemy definicję sprzężeń wypukłych z niewielkimi poprawkami dla znaków.

Generuje to podwójne:

$$ \begin{align} \max_{x, y, r} & \, -\frac{1}{2}x^TQx - p(-y, l^c, u^c) - p(-r, l^v, u^v) \\ \text{s.t.}\; & Qx + c - A^Ty = r \end{align} $$

Rozszerzając definicję terminu $p$, otrzymujemy wartość podwójną podaną u góry.

Formacja w kształcie siodła

Gradient hybrydowy podwójny (patrz Chambolle i Pock) jest związany z podstawowym problemem dotyczącym postaci

$$ \begin{align} \min_x f(x) + g(x) + h(Kx) \end{align} $$

które przez sprzężenie podwójnych odpowiada problemowi z punktem siodłowym

$$ \begin{align} \min_x \max_y f(x) + g(x) + y^TKx - h^*(y) \end{align} $$

PDLP powoduje, że problem z programowaniem wypukłym kwadratowym jest tworzony w ten sposób, przez ustawienie:

  • $f(x) = 0 USD
  • $g(x) = c^T x + \frac{1}{2} x^T Q x + \mathcal{I}_{[l^v, u^v]}(x)$
  • $h(a) = \mathcal{I}_{[-u^c,-l^c]}(a)$
  • tys. PLN = -A$

Jak już wspomniano, $h^*(y) = p(y; -u^c,-l^c)$ to liniowa funkcja konweksowa. Zarówno $g$, jak i $h^*$ mogą przyjmować wartości nieskończone, co skutecznie ogranicza domeny odpowiednio $x$ i $y$.

Zwróć uwagę, że operator bliższy $g$ jest obliczany w formie zamkniętej zgodnie z założeniem, że $Q$ jest przekątną. Trzeba pamiętać, że $g$ jest rozdzielany i ma tę właściwość, która obowiązuje w przypadku wszystkich funkcji $f_1, f_2$:

$$ f_2(t) = f_1(t) + \frac{\mu}{2} t^2 \Rightarrow \mathrm{prox}_{\lambda f_2}(t) = \mathrm{prox}_{\frac{\lambda}{1 + \lambda \mu} f_1}\left( \frac{t}{1 + \lambda \mu} \right). $$

Na dowód tego znajdziesz np. twierdzenie 6.13 w artykule Metody Pierwszeństwa w optymalizacji. Powstałe wyrażenie jest

$$ \begin{equation} \mathrm{prox}_{\lambda g}(x) = \mathrm{proj}_{[l^v, u^v]}\left( (I + \lambda Q)^{-1} (x - \lambda c) \right) \end{equation} $$

Zmniejszone koszty, podwójne wartości resztowe i poprawiony podwójny cel

Do obliczania wartości z punktami siodłowymi bezpośrednio działa tylko element $(x,y)$. Obniżone koszty $r$ są niejawne. Aby zwrócić wynik $(x,y,r)$ podczas rozwiązywania formuły w punkcie siodłowym, definiujemy $r$ jako $r = Qx + c – A^Ty$. Skorygowany cel podwójny to obiektywna wartość problemu podwójnego i zawsze wyznacza dolną granicę wartości celu, ale zawsze wynosi -\infty$ przy obniżeniu kosztów.

Zmniejszone koszty i podwójny cel zgłoszony przez PDLP zostały zmienione między narzędziami LUB w wersjach 9.7 i 9.8. Aby sprawdzić, która wersja jest używana, sprawdź komentarz z opisem SolverResult::reduced_costs w primal_dual_hybrid_gradient.h i sprawdź, czy jest w nim wspomniana wersja 9.8.

Wersja 9.7 i wcześniejsze

Aby uzyskać bardziej zrozumiałą wartość podwójną, gdy poprawiony cel podwójny to $-\infty$, raportujemy też cel podwójny, który ignoruje nieskończone terminy w wartości celu. Podwójne reszty to wartości $r$ z terminów nieskończonych w poprawionym celu podwójnym, z 0 w pozostałych komponentach, a zmniejszone koszty zwracane przez PDLP to wartości $r$ z terminów skończonych w poprawionym celu podwójnym, przy 0 w pozostałych komponentach (co daje $r = \mbox{residuals}).

Wersja 9.8 i nowsze

Aby uzyskać bardziej zrozumiałą wartość podwójną, gdy poprawiony cel podwójny to $-\infty$, raportujemy cel podwójny, który zastępuje nieskończone w wartości celu wyrażeniami skończonymi. Jeśli jedna z granic jest ograniczona, zamiast niej zostanie użyta granica nieskończona. W przeciwnym razie granica będzie wynosić zero. Ten wybór zachowuje wklęsłość celu podwójnego, ale niekoniecznie wyznacza dolną granicę wartości celu. Podwójne reszty to wartości $r$ z haseł nieskończonych w poprawionym celu podwójnym. Obniżone koszty zwracane przez PDLP wynoszą $r$.

Traktowanie niektórych zmiennych granic jako nieskończonych

Jeśli w obu wersjach rozwiązaniem funkcja handle_some_primal_gradients_on_finite_bounds_as_residuals ma wartość true (domyślnie), dodatkowe granice zmiennych mogą być traktowane jako nieskończone podczas obliczania celu podwójnego i podwójnych reszt. Jeśli $|x_i - l^v_i| > |x_i|$, $l^v_i$ jest traktowany tak, jakby był nieskończony, i podobnie, jeśli $|x_i - u^v_i| > |x_i|$, $u^v_i$ jest traktowany tak, jakby był nieskończony.

Pamiętaj, że funkcja handle_some_primal_gradients_on_finite_bounds_as_residuals nie ma wpływu na obliczone iteracje. Ma wpływ tylko na podwójny cel i wartości resztowe używane w testach zakończenia i raportowanych statystykach.

Przeskalowanie

Załóżmy, że otrzymano ukośną (kolumnę) macierz przeskalowania $C$ i macierz przeskalowania po przekątnej $R$ z wpisami dodatnimi po przekątnej. Stosując przeskalowanie tak jak w ShardedQuadraticProgram::RescaleQuadraticProgram, otrzymujemy taki przekształcony problem:

$$ \begin{align} \min_{\tilde x} & \, (Cc)^T{\tilde x} + \frac{1}{2}{\tilde x}^T(CQC){\tilde x} \\ \text{s.t.}\; & Rl^{c} \le RAC\tilde x \le Ru^{c} \\ & C^{-1}l^{v} \le \tilde x \le C^{-1}u^{v} \end{align} $$

Rozwiązanie pierwotnego problemu jest przywracane jako $x = C\tilde x$. Jeśli $\tilde y$ to rozwiązanie podwójne, a $\tilde r$ pozwala zmniejszyć koszty przekształconego problemu, $y = R\tilde y$ to rozwiązanie podwójne, a $r = C^{-1}\tilde r$ to zmniejszenie kosztów pierwotnego problemu (derywowanie).

Identyfikacja niewykonalności

Certyfikat pierwotnej niewykonalności to punkt $(y, r) \in \mathbb{R}^m \times \mathbb{R}^n$ satysfakcjonujący:

$$ \begin{equation} \left((l^{c})^T[y]_+ - (u^{c})^T[y]_- \right) + \left((l^{v})^T[r]_+ - (u^{v})^T[r]_- \right) > 0 \\ -A^T y = r . \end{equation} $$

Istnienie takiego punktu wskazuje, że problem pierwotny nie ma rozwiązania.

Podobnie certyfikat podwójnej niewykonalności to punkt $x \w \mathbb{R}^n$, który:

$$ \begin{equation} Qx = 0 \\ c^T x < 0 \\ (Ax)_i \begin{cases} = 0 & \text{if }l^c_i , u^c_i \in \mathbb{R}, \\ \geq 0 & \text{if }l^c_i \in \mathbb{R}, u^c_i = \infty, \\ \leq 0 & \text{if }l^c_i = -\infty, u^c_i \in \mathbb{R}, \end{cases} \\ x_i \begin{cases} = 0 & \text{if }l^v_i , u^v_i \in \mathbb{R}, \\ \geq 0 & \text{if }l^v_i \in \mathbb{R}, u^v_i = \infty, \\ \leq 0 & \text{if }l^v_i = -\infty, u^v_i \in \mathbb{R}. \end{cases} \end{equation} $$

Pamiętaj, że certyfikaty dla programu liniowego można uzyskać po ustawieniu $Q=0$.

  1. Symetrycznie dodatnią macierz celów symetrycznych $S$ można przekształcić do tej postaci, rozliczając $S$ jako $S = R^T R$ (na przykład rozłożenie na czynniki Cholesky'ego). Wprowadzamy dodatkowe zmienne $z$ zdefiniowane przez ograniczenia $R x - z = 0$, co daje $x^T S x = ↩