Obsah

Úvod

V této práci představíme binární paint shop problem. Jedná se o úlohu, kterou neumíme efektivně řešit (protože je $\NP$ úplná) a ani aproximovat (za předpokladu Unique games conjecture je konstantní aproximace také $\NP$ těžká), takže mimo jiné probíhá aktivní výzkum snažící se najít algoritmus, který je dobrý v průměrném případě (pro náhodný vstup).

Zadání

Zadání binárního paint shop problému (dále těž BPS) je následující:

Počet změn barev řešení považujeme za skóre algoritmu. Relativní skóre pak je poměr skóre a počtu typů aut.

Problém můžeme chápat buď jako optimalizační problém, kde účelová funkce je počet změn v řešení a snažíme se ji minimalizovat, nebo jako rozhodovací problém, kde se ptáme, jestli existuje řešení s nejvýše nějakým zadaným počtem změn.

Cíle práce

Naším cílem bude najít co nejlepší algoritmus řešící BPS a odhadnout u něj střední hodnotu skóre pro náhodný vstup. Dále se pokusíme porovnat již známé algoritmy.

Struktura práce

Nejprve zavedeme notaci potřebnou pro pohodlnou práci s BPS a definujeme aproximační algoritmy. Kapitola 2 obsahuje shrnutí doposud známých algoritmů a jiných výsledků ohledně BPS. V kapitole 3 je představen princip semidefinitního programování, které má uplatnění v algoritmu představeném o kapitole 4. Kapitola 5 se věnuje praktické implementacím algoritmů a naměřeným datům o nich.

Notace

Auta budeme indexovat od $0$ do $2n-1$ a typy od $0$ do $n-1$.

Pomocí $a_{i,0}$ budeme značit pozici prvního auta typu $i$ a pomocí $a_{i,1}$ pozici druhého.

Opačně $t_i$ bude značit typ auta na pozici $i$ a $p_i$ bude značit, jestli se jedná o první nebo druhé auto daného typu (tedy bude nabývat hodnoty $0$ nebo $1$).

Protože někdy bude výhodnější zacházet se znaménky, zavedeme $P_i = -1+2 p_i$, které bude nabývat hodnot $-1$ a $1$. Také rozšíříme notaci o $a_{i, -1} = a_{i,0}$.

Nakonec ještě zavedeme zkratku indexu druhého auta stejného typu jako je auto na pozici $i$: $o_i = a_{t_i,-P_i}$.

Barvy pro nás budou $0$ a $1$. Pokud budeme uvažovat nějaké konkrétní obarvení $c$, tak $c(i)$ bude značit barvu auta na pozici $i$.

Pro algoritmus $\alg$ označíme skóre na vstupu $\alpha$ pomocí $\gamma_{\alg}(\alpha)$. Dále označme všechny vstupy délky $n$ jako $W_n$. Průměrné skóre algoritmu na všech vstupech délky $n$ tedy bude $$\gamma_{\alg}(n) = {\sum_{\alpha \in W_n} \gamma_{\alg} (\alpha) \over |W_n|}.$$ Tuto hodnotu mimo jiné můžeme chápat jako střední hodnotu skóre pro uniformně náhodný vstup délky $n$. K tomu ještě zavedeme notaci pro relativní skóre $\delta_{\alg}(\alpha) = \gamma_{\alg}(\alpha)/n_{\alpha}$ a analogicky průměrné relativní skóre $\delta_{\alg}(n) = \gamma_{\alg}(n)/n$. Nakonec označme $\delta_{\alg} = \lim_{n\to\infty} \delta_{\alg}(n) = \lim_{n\to\infty} \gamma_{\alg}(n)/n$. Protože limita nemusí vždy existovat, zavedeme ještě limes superior a limes inferior: $\delta^+_{\alg} = \limsup_{n\to\infty} \delta_{\alg}(n)$ a $\delta^-_{\alg} = \liminf_{n\to\infty} \delta_{\alg}(n)$. Význam těchto definic bude vysvětlen v následující kapitole.

Dále označme $\gamma(\alpha)$ optimální skóre na vstupu $\alpha$. A následně analogicky definujeme $\gamma(n), \delta(\alpha), \delta(n)$, $\delta$, $\delta^+$ a $\delta^-$ stejně jako u variant s algoritmem, jen skóre daného algoritmu nahradíme za optimální skóre.

Definice aproximačních algoritmů

Protože ne pro všechny problémy známe polynomiální algoritmus, který je schopný je vyřešit, zajímavý výsledek může být, i když se k řešení zvládneme jen v nějakém smyslu alespoň přiblížit. Na to nejprve musíme říct, co pro nás znamená, že nějaké řešení je několikrát horší než jiné. To nám poskytne obecná definice optimalizačního problému.

A nyní již přejdeme k samotným definicím algoritmů blížících se optimu.

Předešlá definice záleží na tom, co považujeme za velikost vstupu, což se pro různé problémy a jejich interpretace liší. Naštěstí nás většinou budou zajímat $c$-aproximace, tedy aproximace, kde $g(n)$ je konstantní funkce rovna $c$.

U maximalizačních problémů je drobný problém v terminologii, protože není shoda na tom, jestli má platit $f(\alg(I)) \ge g(|I|) \cdot \opt(I)$ nebo $f(\alg(I)) \ge {1\over g(|I|)} \cdot \opt(I)$. Naštěstí podle kontextu jde snadno rozhodnout, která definice se používá, protože je potřeba násobit optimum číslem menším rovno jedné (jinak by pro kladné hodnoty účelové funkce nemohl existovat žádný vyhovující algoritmus).

Bohužel někdy asi ani s aproximačními algoritmy nevystačíme a proto zavedeme pravděpodobnostní relaxaci.

Snadným důsledkem definice (a Markovovy nerovnosti) je, že pro každé $\lambda > 1$ pravděpodobnostně $g(n)$-aproximační algoritmus vrátí s pravděpodobností zdola omezenou konstantou $1- {1\over \lambda} \in (0,1)$ řešení s hodnotou nejvýše $\lambda \cdot g(n)\cdot \opt(I)$.

Všimněme si, že na hodnotě $1-{1\over\lambda}$ příliš nezáleží. Opakovaným spouštěním algoritmu a pak vybráním nejlepšího z dodaných řešení zvládneme pravděpodobnost libovolně těsně přiblížit k jedné.

Pokud řešíme složitost algoritmů a úloh, většinou vyžadujeme, aby účelová funkce i rozhodování přípustnosti řešení byly vyčíslitelné v polynomiálním čase. Navíc chceme, aby všechna přípustná řešení měla délku omezenou nějakým polynomem v délce vstupu. Snadno nahlédneme, že za takovýchto podmínek je rozhodovací verze, jestli je optimum alespoň zadané číslo, v $\NP$. Pro takovou úlohu většinou hledáme polynomiální aproximační algoritmus. Binární paint shop vyhovuje všem těmto podmínkám.

Pokud existuje $(1+\varepsilon)$-aproximační algoritmus pro každé $\varepsilon > 0$, říkáme, že máme aproximační schéma.

Doposud známé výsledky

Bonsmaa, Epping a Hochstättler [1] dokázali, že optimalizační verze BPS je $\APX$-těžký problém, což je silnější tvrzení, než že rozhodovací verze je $\NP$-těžká. Za předpokladu $\P \neq \NP$ víme, že kromě polynomiálního algoritmu na rozhodovací verzi nemůže existovat ani polynomiální aproximační schéma na optimalizační verzi.

Snadno nahlédneme, že rozhodovací problém patří do $\NP$. Když nám někdo dá přiřazení barev autům, v lineárním čase jsme schopni ověřit, že počet změn je dostatečně malý a od každého typu a barvy auta je právě jeden výskyt. Tedy rozhodovací problém je $\NP$-úplný.

Dále pak Gupta a spol. [2] ukázali, že za předpokladu Unique games conjecture je NP-těžké i problém libovolně konstantně aproximovat. K tomu využili převod na problém minimálního ne-řezu. To je problém podobný maximálnímu řezu, který bude představen později. Ovšem místo maximalizace počtu hran v řezu minimalizujeme počet hran mimo něj (tedy v rámci partit). Tyto dva problémy mají stejné optimum, ovšem aproximovatelnost se na nich chová drasticky jinak. V momentě, kdy jsou skoro všechny hrany v řezu, tak přidání další hrany do řezu skoro nezmění počet hran v něm. Ovšem odebrání hrany z ne-řezu bude mít velký vliv (procentuálně) na počet hran v ne-řezu.

Vzhledem k předešlým výsledkům je na místě zkoumat různé heuristiky a obecně algoritmy, co jsou nás schopny k optimálnímu řešení alespoň částečně přiblížit.

Jedním z takových algoritmů je hladový algoritmus popsaný Andresem and Hoch-stättlerem [3]:

V daném článku autoři o tomto algoritmu dokázali, když vezmeme uniformně náhodný vstup délky $n$, tak střední hodnota počtu změn ve vygenerovaném řešení je $\sum_{0\le k < n} {2k^2-1 \over 4k^2-1}$.

Hladový algoritmus budeme značit $\algo g$, tedy předešlou hodnotu značíme $\gamma_{\algo g}(n)$. Připomeňme, že pomocí $\delta_{\algo g}(n)$ značíme tuto hodnotu vydělenou velikostí vstupu, tedy $\gamma_{\algo g}(n) / n$.

Pro dostatečně velká $n$ se $\gamma_{\algo g}(n)$ pohybuje zhruba okolo $(1 / 2) n$ (formálně řečeno: platí, že $\delta_{\algo g} = \lim_{n\to\infty} {\gamma_{\algo g}(n) / n} = {1 / 2}$) (viz obrázek ).

U hladového řešení je tedy střední hodnota skóre přes uniformně náhodný vstup přesně vyčíslená. Pro nás je zejména důležité, že jsme ji schopni shora odhadnout. Střední hodnotu totiž můžeme považovat za ukazatel kvality algoritmu (čím menší je, tím se jedná o lepší algoritmus) a tedy horní odhad nám dává záruku kvality algoritmu.

Zajímavé je pro nás zkoumat chování algoritmů na velkých vstupech, tedy dává smysl uvažovat limitu střední hodnoty do nekonečna. Ovšem s narůstajícím počtem aut narůstá i počet potřebných změn, takže samotná limita střední hodnoty moc nedává smysl, protože by byla nekonečná. Místo ní budeme uvažovat $\lim_{n\to\infty} {\gamma_{\alg}(n) / n} = \lim_{n\to\infty}\delta_{\alg}(n)$. O ní víme, že pro libovolný algoritmus (pokud existuje) bude v intervalu $[0, 2]$, protože maximální počet změn musí být mezi $0$ a $2n-1$.

I když nejsme schopní hledat optimum efektivně, pro každý vstup je určitě optimum dobře definovaná hodnota. Můžeme se tedy ptát na otázku, kolik vyjde limita, kde místo skóre algoritmu vložíme optimum pro daný vstup, tedy $\lim_{n\to\infty} \delta(n) = \lim_{n\to\infty} \gamma(n)/n$. Označme hodnotu této limity $\delta$.

Bohužel není zřejmé, že limita skutečně existuje, proto místo limity budeme často uvažovat limes supervisor a inferior. Ty budeme značit přidáním $\pm$ jako horního indexu. Stejnou notaci budeme používat i u limit algoritmů.

Jelikož optimum musí být alespoň tak dobré jako výstup hladového řešení, dostáváme konstruktivní horní odhad $\delta^+ \le 0.5$. Hledáním lepších algoritmů můžeme horní odhad zlepšovat.

Překvapivě však je, že Hančl a kol. [4] ukázali dolní odhad $\delta^- \ge 0.214$ pomocí počítání pravděpodobností pro náhodné obarvení. Samotný fakt, že $\delta^- > 0$ je již poměrně zajímavé tvrzení. To nám říká, že každý algoritmus na průměrném vstupu musí použít aspoň $\Omega(n)$ změn barev. Nemůže tedy například existovat algoritmus, kterému by stačilo jen $\O(n / \log n)$ změn. To také říká, že odhadovat algoritmy pomocí $\lim_{n\to\infty} \gamma_{\alg}/n$, je alespoň co se týče asymptotiky dostačující, protože vystihuje konstantu nejvýznačnějšího členu polynomiální aproximace.

Dále si představíme několik algoritmů, které se snaží konstruktivně zlepšit horní odhad na $\delta^+$.

Autoři algoritmu, Andres and Hochstättler [3], o něm dokázali, že ${2\over 5}\,n - {8\over 15} \le \gamma_{\algo{rg}}(n) \le {2\over 5}\,n + {7\over 10}$, tedy $\delta_{\algo{rg}} = 2/5 = 0.4$.

Dále se o zlepšení tohoto algoritmu pokusili Hančl a kol. [4]. Ti si všimli, že ve výstupu rekurzivního řešení se občas stane, že přebarvením dvojic typů aut se zlepší skóre. Ukázali, že takových dvojic je vždy ve výstupu lineárně mnoho a z toho pak ukázali horní odhad $\delta^+ \le 0.4 - \varepsilon < 0.4$ (pro $\varepsilon$ zhruba $2\cdot 10^{-6}$).

Hančl a kol. [4] přišli ještě s dalším způsobem, jak vylepšit rekurzivní hladové řešení. Při běhu hladového řešení se občas stane, že některá dvojice aut stejného typu si může prohodit barvy bez změny skóre řešení. Když přidáváme auto do okolí takovéto dvojice, v některých případech můžeme provést toto prohození a tím snížit počet nově vytvořených změn barev. Budeme si tedy udržovat některé z takovýchto dvojic a pokud budeme přidávat auto do okolí evidované dvojice tak, že prohození barev dané dvojice by pomohlo, prohodíme její barvy.

Budeme pracovat nad rozšířenou abecedou barev o znak „$*$“ reprezentující neurčenou barvu. Při rekurzi budeme udržovat invariant, že pro každý typ obě auta buď budou označena $*$ a nebo ani jedno z nich nebude označeno $*$ a pak nutně budou mít různé barvy. Navíc bude platit, že $*$ nikdy není na okrajích a nejsou dvě vedle sebe. Specificky tedy na $*$ budeme přebarvovat dvojici aut v momentě, kdy získá všechny sousedy.

Pro popis algoritmu zavedeme notaci $N(i)$, což bude reprezentovat multimnožinu barev aut sousedících s autem na pozici $0<i<2n-1$, tedy $\{c(i-1), c(i+1)\}$.

Autoři algoritmu o něm vyslovili domněnku, že $\delta_{\algo{rsg}} = {1\over 13} \cdot \sqrt{61} - {3\over 13} \doteq 0.370$. Na obrázku je zobrazen histogram skóre tohoto algoritmu spuštěného na náhodných vstupech a hodnota z předešlé hypotézy.

Zobecněné verze a související problémy

Zobecněním binárního paint shopu je obecný paint shop problém, zavedený Eppingem a spol. [5]. Z něj pochází motivace k binární verzi.

Binární verze problému je tedy speciální případ obecného paint shop problému, kde od každého typu máme dvě auta, máme jen dvě požadované barvy a navíc platí, že pro každou kombinaci typu a barvy chceme právě jedno auto.

Souvisejícím problémem je dělení náhrdelníku zavedený v Alonově článku [6]:

Nejprve si všimněme, že dělení náhrdelníku je speciálním případem obecného paint shop problému. Auta pro nás budou drahokamy. Barva auta bude označovat, který loupežník získá daný drahokam. Minimalizovat počet intervalů je to stejné jako minimalizovat počet hranic mezi nimi (kterých je o jedna méně něž intervalů). Ve vstupu problému dělení náhrdelníku navíc musí platit, pro každý typ má být stejný počet aut obarvený jednotlivými barvami (a tedy počet aut daného typu musí být násobkem $k$). Naopak binární paint shop problém je speciálním případem dělení náhrdelníku pro dva lupiče, kde navíc platí, že všechna $a_i$ jsou 1 (tedy od každého typu jsou na náhrdelníku právě dva drahokamy).

Semidefinitní programování

V této kapitole si představíme princip semidefinitního programování (dále též SDP), jak jej popisují Gärtner a Matoušek [7], a jeho použití na problém maximálního řezu, z něhož vychází algoritmus na binární paint shop.

Nejprve zavedeme a připomeneme notaci důležitou v této kapitole. Nechť $\R^{n\times m}$ značí množinu $n$-řádkových $m$-sloupcových matic složených z reálných čísel. Řádky i sloupce indexujeme od $0$, tedy matice $A\in \R^{n\times m}$ obsahuje prvky $A_{i,j}$ pro všechna $0\le i < n$ a $0 \le j < m$. Nechť $\SYM_n = \{X \in \R^{n\times n} \mid x_{i,j} = x_{j,i} \hbox{\ pro všechna\ } 0 \le i,j < n\}$ je třída všech symetrických matic a nechť $X \circ Y = \sum_{0\le i<n} \sum_{0\le j<m} x_{i,j} y_{i,j}$ značí součet součinu matic po složkách. Nakonec $X \succeq 0$ bude značit skutečnost, že matice $X$ je pozitivně semidefinitní (bude vysvětleno později).

Úloha semidefinitního programování je optimalizační úloha (podobně jako lineární programování) následujícího formátu: $$\progline{maximalizuj}{C \circ X}{}$$ $$\progline{za podmínek}{A_i \circ X = b_i}{0 \le i < m-1}$$ $$\progline{}{X \succeq 0.}{}$$

Vstupem programu tedy je velikost matice $n$, počet podmínek $m \in \N$ a pak pro každou podmínku matice $A_i \in \SYM_n$ a číslo $b_i \in \R$. Výstupem je pak $X\in \SYM_n$ splňující výše uvedené podmínky.

Pro nás bude důležitý následující fakt z lineární algebry:

Pro nás bude důležitá zejména třetí podmínka, protože navíc platí, že ze semidefinitní matice $X$ zvládneme zkonstruovat $Y$ v čase $\O(n^3)$ pomocí tzv. Choleského dekompozice.

Navíc pro libovolnou reálnou matici $Y$ je $Y^{\rm T} Y$ symetrická. Tedy semidefinitní programování můžeme chápat jako optimalizační úlohu na $Y\in \R^{n\times n}$. Pojďme se zamyslet nad tím, co v takovémto pohledu znamenají podmínky a účelová funkce. Matici $Y$ můžeme považovat za $n$ sloupcových vektorů $\vec{y_0}, \dots, \vec{y_{n-1}}$ vedle sebe. Matice $Y^{\rm T}$ pak odpovídá těmto vektorům zapsaných v řádcích pod sebou. Matice $X = Y^{\rm T} Y$ tedy na pozici $i,j$ obsahuje skalární součin vektorů $\vec{y_i}$ a $\vec{y_j}$. Účelová funkce tedy je lineární kombinací skalárních součinů a podmínky odpovídají vynucení rovnosti lineární kombinace skalárních součinů vektorů a konstanty. Speciálně tedy můžeme mít podmínku na délku vektoru: $|\vec{y_i}|^2 = \vec{y_i}^{\rm T}\vec{y_i} = C$. Semidefinitní programování v tomto tvaru budeme nazývat semidefinitní programování v dekomponovaném tvaru.

Maximální řez

Představíme si pravděpodobnostní aproximační algoritmus založený na semidefinitním programování na následující problém:

Problém maximálního řezu (resp. rozhodovací verze, kde se ptáme na existenci řezu alespoň dané velikosti) je NP-úplný [8], proto se u něj zkoumají aproximační algoritmy a pravděpodobnostní aproximační algoritmy.

Nejprve si ukážeme triviální pravděpodobnostní $0.5$-aproximační algoritmus:

Povšimněme si, že triviální algoritmus dosáhl poměrně zajímavého aproximačního poměru a to se ani nedíval na vstup.

Předešlý algoritmus lze derandomizovat, jak popisuje Dimitrakakis [9]. Výsledný algoritmus pak vždy najde řešení obsahující alespoň $1/2$ hran, tedy alespoň $1/2$ optima.

Lepšího aproximačního poměru dosáhneme Goemansovým-Williamsonovým algoritmem, založeným na semidefinitním programování a proto popsaným mimo jiné v již dříve zmíněném úvodu do semidefinitního programování [7].

Naivní implementace je, že si pro každý vrchol $u$ vyrobíme proměnnou $x_u$, která může nabývat hodnot $\pm 1$, která bude říkat, do jaké množiny máme vrchol umístit. Do účelové funkce pak zakódujeme graf. Účelová funkce bude $$\sum_{uv \in E} h(u,v)\cdot\left(- \frac{x_u x_v}{2} + \frac{1}{2}\right) = \frac{F}{2} + \frac{1}{2} \sum_{uv \in E} -h(u,v)x_u x_v ,$$ kde $F=\sum_{uv\in E} h(u,v)$ značí součet hodnot všech hran. Pro každou hranu tedy přičte $h(u,v)$, pokud $x_u$ a $x_v$ mají různá znaménka, a $0$, pokud stejná.

Takovýto optimalizační problém bohužel ani pomocí semidefinitního programování neumíme řešit. Na to, abychom uměli vyjádřit nějaký součin, musíme použít semidefinitního programování v dekomponovaném tvaru. Pak ale místo toho, abychom měli pro každý vektor jen jednu proměnnou, budeme mít pro každý vrchol celý vektor proměnných a místo součinu proměnných budeme mít skalární součin těchto vektorů. Již ale současně neumíme zařídit, aby tyto vektory byly buď $(1,0,\dots,0)$ nebo $(-1,0,\dots,0)$. Můžeme ale přidat podmínku na to, aby délka každého vektoru byla $1$. Tím jsme vyrobili relaxaci výše uvedeného programu, protože za těchto podmínek může vektor nabývat i jiných poloh než $(\pm 1,0,\dots,0)$. Celý algoritmus pak vypadá následovně:

Semidefinitní program v algoritmu si lze představit jako umísťování $n$ bodů na $n-1$ dimenzionální sféru v $\R^n$. Účelová funkce se snaží umístit body vrcholů spojených hranou co nejdále od sebe.

V případě, kdy jsme povolovali pouze hodnoty $\pm 1$, bylo jasné, které vrcholy patří do které množiny. Nyní však výstupem optimalizace jsou vektory a tedy přiřazení množinám není tak jednoznačné. Rádi bychom umístili od sebe vzdálené body do různých množin. Na to můžeme rovnoměrně náhodně zvolit nadrovinu procházející počátkem a rozdělit body do množin podle toho, do které z polorovin určených danou nadrovinou patří. (Předpokládejme, že žádný z bodů neleží na nadrovině. Pravděpodobnost náležení nadrovině je $0$. Pokud se tak stane, je vcelku jedno, do jaké množiny ho vložíme.)

To nám zaručí, že dvojice bodů daleko od sebe budou mít velkou pravděpodobnost toho, že budou v různých poloprostorech. Semidefinitní programování se tedy snaží dostat dvojice bodů odpovídající vrcholům spojených hranou daleko od sebe a podle toho se zvyšuje pravděpodobnost toho, že vrcholy budou v jiných množinách a tedy hrana bude v řezu.

Výběr nadroviny a rozdělování bodů jde implementovat pomocí výběru vektoru $\vec{z}$ kolmého na nadrovinu. Pro každý bod pak stačí spočítat skalární součin se $\vec{z}$ a podle znaménka víme, do kterého poloprostoru patří. Pokud vybereme uniformně náhodný jednotkový vektor, tak jsme uniformně náhodně vybrali nadrovinu. Jednotkový vektor můžeme generovat tak, že vygenerujeme náhodný vektor nezávisle po složkách z normálního rozdělení, a pak ho znormujeme. Povšimněme si, že normalizace ani není potřeba, protože to na znaménku součinů nic nemění.

Nyní pojďme precizněji spočítat pravděpodobnost toho, že se dvojice bodů (vektorů) $\vec{y_u}$, $\vec{y_v}$ rozdělí náhodnou nadrovinou $\rho$ v závislosti na hodnotě skalárního součinu mezi nimi. Můžeme se podívat na rovinu, ve které leží počátek a oba vektory $\vec{y_u}$, $\vec{y_v}$ (viz obrázek ). Průnik náhodné nadroviny s touto rovinou tvoří rovnoměrně náhodně vybranou přímku procházející počátkem.

Zajímá nás tedy, jaká je pravděpodobnost toho, že jedno z ramen nadroviny bude uvnitř konvexního úhlu mezi $\vec{y_u}$ a $\vec{y_v}$. Velikost tohoto úhlu označme $\alpha$. Víme, že $\vec{y_u}^{\rm\, T}\vec{y_v} = 1 \cdot 1 \cdot \cos \alpha$. Náhodnou přímku můžeme vygenerovat jako náhodný úhel $\beta$ mezi $0$ a $2\pi$ s tím, že přímka pak povede směrem $\beta$ a $(\beta + \pi) \mod 2\pi$. Nikdy však uvnitř konvexního úhlu nebudou obě ramena přímky, tedy pravděpodobnost, že alespoň jedno bude uvnitř a tedy body budou oddělené je:

$$\frac{2 \cdot \alpha}{2\pi} = \frac{\alpha}{\pi} = \frac{\arccos \vec{y_u}^{\rm\,T}\vec{y_v}}{\pi}$$

U každé hrany $uv$ optimalizujeme $-h(u,v) \vec{y_u}^{\rm T}\vec{y_v}$, což je ekvivalentní optimalizování $h(u,v)\cdot\left(\frac{1}{2} - \frac{\vec{y_u}^{\rm T}\vec{y_v}}{2}\right)$.

Nyní se podívejme na poměr mezi pravděpodobností toho, že hrana bude v řezu, a účelovou funkcí před vynásobením hodnotou hrany jakožto funkci proměnné $x = \vec{y_u}^{\rm\,T} \vec{y_v}$ v intervalu $\left[ -1, 1 \right]$:

$$\frac{\arccos \vec{y_u}^{\rm\,T} \vec{y_v}}{\pi} / \frac{1 - \vec{y_u}^{\rm T}\vec{y_v}}{2} = \frac{\arccos x}{\pi} / \frac{1 - x}{2}$$

Označme tuto hodnotu $c = 0.8785$.

Tedy když je účelová funkce dané hrany $h(u,v)r$, tak pravděpodobnost výběru hrany do řezu bude alespoň $cr$. Sečtením přes všechny hrany (z linearity středních hodnot) tedy získáme, že střední hodnota součtu hran v řešení je alespoň $cR$, kde $R = \sum_{uv \in E} h(u,v)\cdot\left(\frac 12 - \frac{\vec{y_u}^{\rm T}\vec{y_v}}2\right)$, tedy jedno z ekvivalentních vyjádření účelové funkce.

Pokud tedy najdeme dobré řešení semidefinitního programu, ve střední hodnotě pak najdeme i dobré řešení maximálního řezu. Dále však musí platit, že optimální řešení semidefinitního programu je alespoň součet hran optimálního řešení maximálního řezu, když za účelovou funkci považujeme $\sum_{uv \in E} h(u,v)\cdot\left(\frac 12 - \frac{\vec{y_u}^{\rm T}\vec{y_v}}2\right)$. Jedno z možných řešení totiž snadno sestrojíme z maximálního řezu: Vrcholům jedné množiny přidělíme vektory $(1,0,\dots,0)$ a druhé $(-1,0,\dots,0)$. Účelová funkce pak přičte $h(u,v)$ za každou hranu v řezu a $0$ jinak.

Bohužel na rozdíl od lineárního programování, u semidefinitního programování nejsme schopní v polynomiálním čase najít optimální řešení. Naštěstí však platí, že v polynomiálním čase za určitých podmínek jsme schopní se mu libovolně přiblížit. Přesněji řečeno, pro každé $\varepsilon > 0$ jsme schopní najít řešení s účelovou funkcí nejdále $\varepsilon$ od optima. Ovšem musí platit, že všechna přípustná řešení jsou dostatečně malá. Přesněji řečeno, musí platit, že Frobeniova norma všech přípustných matic je omezena konstantou s polynomiální délkou zápisu. Tohoto výsledku jde dosáhnout pomocí elipsoidové metody, která má nejlepší teoretické výsledky. V praxi se však často používají jiné algoritmy. Všechny naše programy budou splňovat podmínky na řešení, protože každá z hodnot hledané matice nutně bude ležet v intervalu $[-1, 1]$.

Řešení semidefinitním programováním

Řešení vychází z Goemansova-Williamsonova algoritmu na maximální řez. Pro každé auto nám bude semidefinitní program umísťovat bod (někdy také chápán jako vektor) na jednotkovou sféru. Poté náhodně rozřízneme nadrovinou sféru na dvě poloviny a podle toho, do které poloviny auto patří, zvolíme jeho barvu.

Rozdílné barvy aut $i, j$ stejného typu zařídíme tak, že vynutíme $\vec{y_i} = -\vec{y_j}$. Na to nám stačí jediná podmínka – říct, že jejich skalární součin je $-1$.

Účelovou funkcí pak řekneme, že sousední auta mají preferovat stejnou barvu, tedy jejich body na sféře mají být blízko sebe, což znamená, že skalární součin má být co největší.

Pro každou dvojici do účelové funkce přičteme číslo mezi $-1$ a $1$, kde $1$ značí, že vektory jsou stejné, a $-1$, že jsou protilehlé. Účelovou funkci můžeme ekvivalentně zapsat jako $$\hbox{minimalizuj}\qquad \sum_{0\le i < 2n-1} \frac{1}{2} - \frac{\vec{y_i}^{\rm T} \vec{y_{i+1}}}{2} = \frac{2n-1}2 - \frac12 \sum_{0 \le i < 2n-1} \vec{y_i}^{\rm T} \vec{y_{i+1}}$$

Nyní tedy pro každou dvojici sousedních aut přičítáme číslo mezi $0$ a $1$, kde $0$ přičteme pro stejné vektory, mezi nimiž nikdy nebude změna barvy a $1$ přičteme v případě opačných vektorů, mezi nimiž se nutně barva změní. Účelovou funkci se snažíme minimalizovat.

Nyní by se hodilo odhadnout střední hodnotu počtu změn barev stejně jako u maximálního řezu. Jeho důkaz vycházel z pozorování ohledně poměru pravděpodobností výběru a účelové funkce, přesněji řečeno z toho, že tento poměr umíme zdola odhadnout. Ovšem v našem případě místo maximalizace děláme minimalizaci účelové funkce a pravděpodobnosti. Tedy abychom mohli tvrdit, že pravděpodobnost je menší než nějaký násobek účelové funkce, potřebovali bychom horní odhad poměru. Ovšem tento poměr je neomezený (viz obrázky a ). V okolí $x=1$ jsou obě funkce poblíž $0$, ovšem pravděpodobnost se k nule blíží mnohem strměji. Tedy pro dvojici vektorů poblíž sebe je skalární součin skoro 1, ovšem pravděpodobnost oddělení je libovolně krát větší než vzdálenost součinu od jedné.

O binárním paint shop problému je navíc známé, že je za předpokladu Unique game conjecture a $\P\neq \NP$ je polynomiálně neaproximovatelný s konstantním faktorem [2], takže nemožnost výše uvedeného postupu by nás ani neměla zaskočit, protože v případě, že by šlo udělat odhad tímto způsobem, získali bychom pravděpodobnostní aproximační algoritmus.

Místo toho se podíváme na rozdíl pravděpodobnosti a účelové funkce.

Označme konstantu z předchozího lemmatu jako $d = 0.1053$.

Na rozdíl od jiných řešení binárního paint shop problému, zde se jedná o střední hodnotu přes náhodné rozhodnutí algoritmu, nikoliv přes náhodnou distribuci na vstupech. Tedy když spustíme algoritmus na konkrétní instanci, dostaneme řešení a také dolní odhad na optimalní řešení dané instance (ten jednak můžeme odvodit z hodnoty řešení ale také přímo z účelové funkce semidefinitního programu).

Implementační zjednodušení

Předchozí semidefinitni program ještě lze nepatrně zjednodušit se zachováním všech potřebných vlastností. Místo toho, abychom měli vektor pro každé auto, uděláme si jen vektor pro každou dvojici aut stejného typu. Ten bude odpovídat řekněme prvnímu autu z dvojice. Všimneme si, že účelovou funkci zvládneme stále vyjádřit. Kdykoliv jsme používali vektor druhého auta, stačí použít mínus vektor prvního auta. Tedy stačí některé členy jednou až dvakrát vynásobit hodnotou $-1$.

Výsledný program vypadá takto:

$$\progline{maximalizuj}{\sum_{0\le i < 2n-1} \vec{y_{t_i}}^{\rm T} \vec{y_{t_{i+1}}} P_{i} P_{i+1}}{}$$ $$\progline{za podmínek}{|y_i| = 1}{0 \le i < n}$$

Výhodou je, že je potřeba poloviční počet proměnných, matice je čtvrtinová a odebrali jsme $n$ podmínek, proto v implementační části práce využijeme této formy programu.

Měření řešení BPS

Součástí práce je implementace algoritmů řešících Binární paint shop problém. Každý z nich byl následně spuštěn pro různé velikosti, pokaždé na $1\,000$ nezávisle náhodně vybraných vstupech s počtem typů aut 10, 20, 50, 100, 200, 400, 566, 800, 1131, 1600, 2263 a 3200. (Hodnoty, co nejsou hezká čísla jsou zhruba $\sqrt{2}$-násobky předešlé hodnoty, tedy v polovině mezi okolními hodnotami na logaritmické stupnici.) Jedna z implementací $\algo{sdp}$ – pomocí sage vyžaduje velké množství paměti a proto byla spuštěna jen na vstupech do velikosti $566$. Celý test běžel několik dní na čtveřici stojů, na každém dva programy současně, každý z nich běžel jednovlákonvě a využíval nejvýše 16 GB operační paměti. U $\algo{sdp}$ řešení bylo použito 10 náhodných řezů a vždy byl vybraný nejlepší z nich. Gitový repozitář s implementací i naměřenými daty je k dispozici na https://gitlab.kam.mff.cuni.cz/jirikalvoda/binary-paint-shop-problem. Na následujících stranách jsou zpracovaná různá naměřená data.

Pro reprodukovatelnost byly vstupy generovány deterministicky ze seedů. Se stejným seedem se tedy vždy vygeneruje stejný vstup a ideálně i stejné řešení. Pro různé algoritmy a velikosti vstupů byly použité různé seedy.

Praktické řešení semidef. programů

U většiny algoritmů byla implementace poměrně přímočará. Ovšem u řešení semidefinitního programování je většina složitosti algoritmu schovaná právě v řešení semidefinitních programů, což už svojí složitostí nepatří mezi algoritmy, které bychom chtěli (re)implementovat. Proto je žádoucí se spolehnout na funkčnost již existujících implementací. Bohužel kvalita dostupných řešičů semidefinitních programů je poměrně nízká.

SDPA-C

Jedním z použitých programů byl SDPA-C [10]. Jedná se o knihovnu v C++ založenou na metodě vnitřních bodů v primárním a duálním problému. Dle autorů projektu [10]:

Pro účely práce byla použita verze SDPA-C ze dne 2023-06-21 (soubor sdpa-c.7.3.8.src.20180125.tar.gz). Kompilace je poměrně komplikovaná. Jednak vyžaduje velké množství nainstalovaného software (například překladač Fortranu), ale také poskytnutý Makefile není plně kompatibilní s moderními verzemi překladačů.

Bohužel i po úspěšné kompilaci stále není vyhráno. Dodaný program přeložený aktuálním kompilátorem bohužel nefunguje.

Při spuštění i na malých vstupech program rychle spadne s chybou Segmentation Fault ve funkci Newton::compute_bMatgVec_dense_threads_SDP.

Při následujícím ladění programu bylo zjištěno, že se program zacyklí v dané funkci ve smyčce, při které inkrementuje hodnotu proměnné Column_Number až do doby, kdy přeteče, což už přímo vede k pádu programu. Pohledem na kód to vypadá, že toto by se nemělo dít. Program už mnohokrát měl ukončit smyčku a doběhnout na konec funkce, ovšem po přidání dostatečného počtu debugovacích výpisů je vidět, že program cyklus opustí, provede příkazy mezi cyklem a koncem funkce a pak se zase objeví uprostřed cyklu.

Problém byl v tom, že funkce vracející void* dle normy [11] bodu 6.6.3 musí skončit pomocí return a nemůže jen tak opustit funkci dojitím na její konec. Pokud se tak stane, jedná se o nedefinované chování, nikoliv jen o nedefinovanou hodnotu navrácenou z funkce. Kdyby se jednalo jen o nedefinovanou hodnotu, tak se nic zásadního neděje. Funkce může vrátit cokoliv, ale to je nám jedno, protože výsledek se ignoruje.

Ovšem jelikož se jedná o nedefinované chování, překladač může program přeložit tak, že v takovémto případě udělá prakticky cokoliv. A v tomto případě optimalizace překladače přeložili program tak, že v případě opuštění funkce se vykonávání programu vrátí zpět dovnitř a pokračuje se v iterování smyčky. To je důsledkem toho, jakým způsobem se snaží překladač optimalizovat kód. Překladače nemusí dodržovat naprogramovanou strukturu programu. Tedy není problém přesunout část za smyčkou do ní na místo, kde se smyčka opouští. Když překladač ví, že v korektním programu by se tudy nemělo opouštět funkci, nemusí se překladač obtěžovat umísťováním explicitního návratu z funkce na dané místo.

Naštěstí oprava tohoto problému je přímočará – stačí do takovýchto funkcí explicitně dopsat return NULL;.

Toto bohužel není jediný problém s SDPA-C. Na některých vstupech knihovna dojde do stavu, kdy nemá semidefinitní matici a výpočet spadne s následující chybou:

CHOLMOD warning: matrix not positive definite.
file: ../Supernodal/t_cholmod_super_numeric.c line: 911
getMinEigenValue:: cannot decomposition ::
line 193 in sdpa_linear.cpp

Tedy během výpočtu program něco spočte špatně a jako mezivýsledek mu vyjde matice, co není pozitivně semidefinitní, se kterou již dále nejde pokračovat. Příčinu této chyby se bohužel nepodařilo odhalit, protože na rozdíl od předchozího problému vůbec není jasné, kde začít hledat. Detekce nesemidefinitnosti matice totiž může být v úplně jiné části kódu, než kde nastane chyba.

Problém se projevuje například na vstupu s $n=5$ a seedem $19$.

Další problém (nebo ten stejný) se projevuje zejména u větších vstupů. Uprostřed algoritmu se také objeví chybová hláška o nesemidefinitnosti matice. Následně program udělá mnoho iterací s nekonečnou účelovou funkcí, dokud se nepřekročí maximální počet iterací, a pak je vrácena matice ze samých $\pm \infty$, což dokonce porušuje zadané podmínky. Toto se stává například na seedu $18$ pro $n = 150$.

Vzhledem k výše uvedeným problémům dává smysl se porozhlédnout po alternativách.

Sage

SageMath je dle oficiálních stránek [12] open-source software pro matematické výpočty založený na rozšířené syntaxi pythonu, takže je poměrně snadné ho využívat i jako programovací jazyk. Pro účely práce byla použita verze sage 10.1 ze dne 2023-08-20. Sage mimo jiné obsahuje rozhraní pro řešení semidefinitních programů. Rozhraní je navrženo tak, aby bylo schopné pracovat s více různými backendy pro řešení semidefinitních programů. Bohužel z dokumentace [13] nebylo zřejmé, jaké backendy sage podporuje. Při zavolání default_sdp_solver("") se objeví chybová hláška, která tvrdí, že backend by mělo jít nastavit na CVXOPT (což je výchozí hodnota) nebo Matrix, případně uživatelem dodanou třídu. Ovšem při přepnutí na Matrix spuštění řešení spadne na NotImplementedError, tedy zdá se, že ani tento backend není implementovaný (alespoň ve verzi, co mám k dispozici) a tedy jediná možnost je CVXOPT.

Drobným problémem je, že dle dokumentace [13] rozhraní sage vyžaduje formulování vstupu jako primárního semidefinitního programu, ovšem náš algoritmus je založený na řešení duálního semidefinitního programu. Mezi nimi však lze snadno přecházet pomocí triviální úpravy, ovšem výsledný kód pak vypadá velice neintuitivně.

Hodnoty skalárních součinů

V odhadu střední hodnoty počtu změn barev v $\algo{sdp}$ řešení jsme využívali toho, že funkce rozdílu pravděpodobnosti řezu a účelové funkce (funkce $f$) je nejvýše $d = 0.1053$. Ovšem tato funkce nenabývá takto vysokých hodnot zdaleka všude, na polovině definičního oboru je dokonce záporná (viz obrázek ). Kdyby se tedy alespoň část optimalizovaných skalárních součinů vyskytovala mimo oblast, kde $f$ nabývá vysokých hodnot, šel by odhad zlepšit.

Bohužel dle naměřených dat to vypadá, že hodnoty skalárních součinů se koncentrují pouze poblíž maxima $f$. Viz obrázek . Tento výsledek intuitivně dává smysl, protože jenom takto může semidefinitní program dosáhnout lepší účelové funkce než počet změn barev v optimálním řešení. Ovšem náš jediný dolní odhad je hodnota semidefinitního programu a tedy nejsme schopní více přiblížit dolní a horní odhad k sobě.

Z grafu je vidět, že skoro nikdy semidefinitní program neumístí sousední body do protilehlých částí. Drobnou výjimku tvoří hodnota skalárního součinu okolo $-1$, kterých semidefinitní program za $100$ běhů vyrobil zhruba $100$. K tomuto mohl být donucen existencí dvojic aut stejného typu hned vedle sebe, kdy nemá jinou možnost než mezi nimi mít skalární součin $-1$.

Ve 100 vstupech by tedy sousedních aut stejného typu mělo být ve střední hodnotě zhruba 100, což z části vysvětluje pozorovanou anomálii. Není však vyloučeno, že hodnoty skalárního součinu $-1$ nebo blízké program dosáhne i jinak.

Dimenze SDP

Semidefinitní program rozmísťuje vektory do $n-1$ dimenzionální sféry. Z naměřených dat však vychází, že semidefinitní program většinou generuje řešení, které má mnohem menší dimenzi. Přesněji řečeno, pro řešení často platí, že všechny vektory v něm mají několik prvních souřadnic velké hodnoty a ve zbylých souřadnicích mají hodnoty blízké nule. Tedy kdybychom vektory promítli na méně dimenzionální sféru (vzali místo nich nejbližší bod na ní), tak se účelová funkce moc nezmění.

Nevíme o žádné hypotéze, proč program generuje řešení s malou dimenzí. Taktéž není jasné, jestli či jak by toho šlo využit pro získání lepšího řešení či dolního odhadu. Z naměřených dat zobrazených na obrázcích až na následujících stranách však vychází, že průměrná dimenze je mnohem menší než $n$, ovšem nejspíše není omezena žádnou konstantou, tedy s rostoucím $n$ také roste, ale mnohem pomaleji. Hypotéza zní, že průměrná dimenze leží v $\Theta(\log n)$.

Malá dimenze má ale zásadní výhodu pro vizualizaci řešení, protože trojrozměrný prostor (tedy sféru dimenze 2) si zvládne člověk snadno představit. Z čehož plyne, že zhruba polovinu řešení pro $n=50$ jsme schopni zakreslit jen s drobným zkreslením, tak, jak je ukázáno na obrázku .

Časová složitost algoritmu $\algo{sdp}$

Víme, že existuje implementace algoritmu $\algo{sdp}$ s polynomiální časové složitosti (za použití elipsoidové metody). Použité řešiče však využívají jiných metod řešení semidefinitních programů a ani jeden z nich moc neslibuje, jakou přesně složitost má. Proto je na místě změřit, jak rychlé řešení jsou.

Testování probíhalo na stroji s procesorem AMD Ryzen 5 7600. Programům bylo poskytnuto $16\,{\rm GB}$ operační paměti a běh byl omezen na jedno jádro procesoru. Během výpočtu na stroji neběželo nic kromě výpočtu a základních funkcí operačního systému.

Dle výše uvedeného grafu můžeme regresí v logaritmickém grafu usuzovat, že časová složitost obou algoritmů je v $\O(n^4) \cap \Omega(n^3)$. Implementace pomocí SDPA-C má menší multiplikativní konstantu a navíc mnohem rychlejší čas startu, který se zejména projevuje na malých instancích.

Skóre

Nejzajímavější je pro nás naměřené (relativní) skóre algoritmů.

Mezi uvažovanými řešeními jsou dvě implementace $\algo{sdp}$. Implementace pomocí sage bohužel vyžaduje velké množství paměti a proto nebyla testovaná na tak velkých vstupech. O implementaci pomocí SDPA-C víme, že občas chybuje. Naštěstí chyby nenastávají tak často, aby to moc ovlivňovalo výsledné skóre. Pokud algoritmus chyboval a vrátil nějaké řešení, je uvažováno skóre daného řešení. Pokud algoritmus žádné řešení nevrátil, uvažuje se místo toho hodnota $n$, což je jistě horní odhad na snadno dosažitelné řešení.

Na výše uvedeném grafu a tabulkách na následujících stranách si můžeme všimnout, že s rostoucím $n$ se u všech měřených algoritmů zmenšuje rozptyl relativního skóre.

Z grafu vidíme, že pro dostatečně velká $n$ je naměřené relativní skóre algoritmu $\algo{sdp}$ menší než $0.34$. Z toho můžeme tedy usuzovat, že $\delta_{\algo{sdp}}(n) \le 0.34$ pro $n\in\{100, 200, 400, 566, 800, 1131, 1600, 2263\}$. Z toho pak můžeme vyslovit hypotézu, že $\delta_{\algo{sdp}}(n) \le 0.34$ pro $n \ge 100$ a tedy i $\delta_{\algo{sdp}}^+ \le 0.34$.

Z naměřených dat také můžeme usuzovat, že $\algo{sdp}$ je lepší než libovolný z jiných představených algoritmů.

Dolní odhad

Výstup z každého algoritmu nám dává horní odhad na optimum daného vstupu. Ovšem semidefinitní programování nám navíc dává i dolní odhad na optimum, protože víme, že optimum jednoho z ekvivalentních vyjádření účelové funkce je vždy menši rovno optimu BPS. Použité řešiče semidefinitních programů nám navíc dají i horní odhad na optimum SDP (který se od nalezeného řešení může nepatrně lišit).

Takto vygenerovaný dolní odhad nám bohužel nepřináší žádnou zajímavou informaci o chování na náhodném vstupu. Připomeňme, že Hančl a kol. [4] dokázali, že $\delta^- \ge 0.214$. Průměrný dolní odhad na testovaných vstupech je pro dostatečně velký vstup dle obrázku menší než $0.214$, tedy nedostáváme ani žádnou novou hypotézu na dolní odhad.

Nicméně stále je zajímavé, že pro konkrétní instanci umíme určit nějaký netriviální dolní odhad.

Závěr

V této práci jsme představili algoritmus ${\bf sdp}$ na BPS založený na semidefinitním programování. Bohužel se nám nepodařilo dokázat žádný netriviální odhad na $\delta_{\algo{sdp}}^+$. Nicméně dle naměřených dat můžeme soudit, že $\delta_{\algo{sdp}}^+$ se pohybuje okolo 0.34. Místo toho jsme však dokázali, že pro libovolný vstup bude střední hodnota (přes náhodná čísla generovaná algoritmem nikoliv přes vstup) skóre řešení nejhůře $0.212 n$ od optima, tedy, že platí $$\forall \alpha:\qquad\E[\delta_{\algo{sdp}}(\alpha)] \le \delta(\alpha) + 0.212 n.$$

Toto řešení jsme dále dokonce dvakrát implementovali s využitím různých implementací řešení semidefinitních programů, které jsme tímto i porovnali. Z naměřených dat jsme jednak odhadli střední hodnotu skóre algoritmu přes náhodný vstup. A také jsme si všimli, že dimenze řešení semidefinitního programování je poměrně malá, což jsme zformulovali jako hypotézu.

Stále však zůstává otevřená otázka, kolik přesně je $\delta^+$ a $\delta^-$ (a případně zda se rovnají) a jaké nejlepší $\delta^+_{\alg}$ je možno dosáhnout polynomiálním algoritmem $\alg$.

Seznam použité literatury

Seznam obrázků