Problém

V této práci si představíme binary paintshop problem. Jedná se o úlohu, kterou neumíme efektivně řešit a ani aproximovat (Meunier a Neveu 2012), 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).

V této zápočtové práci bych se rád podíval na některé z těchto algoritmů a pomocí statistických metod pro zvolenou velikost vstupu určil, který z nich je lepší.

Zadání

Nejprve si pojďme představit zadání binary paint shop problému: V řadě je $2n$ aut $n$ různých typů – od každého typu 2. Chtěli bychom od každého typu nabarvit jedno auto červeně a druhé modře. Auto však na barvící linky vjíždí v pořadí, v jakém jsou v řadě. Barvící linka je optimalizovaná na barvení velkého počtu aut jednou barvou. Tedy měnit barvu, kterou se barví, je složitá a drahá záležitost. Chceme tedy nají obarvení aut tak, aby od každého typu bylo jedno červené a jedno modré, přitom počet změn barev v řadě byl co nejmenší. Počet změn barev značíme jako skóre algoritmu.

Programy

Hladové řešení

Jednoduché řešení je následující: Půjdeme postupně po řadě. První auto v řadě obarvíme řekněme červeně. Dále budeme vždy první auto daného typu barvit stejně jako předcházející auto. Druhé auto daného typu obarvíme vždy zbývající barvou.

Snažit se nějak statisticky analyzovat a porovnávat tento algoritmus není moc zajímavé, protože je známo, že střední hodnota počtu změn hladového algoritmu je $\sum_{k=0}^{n-1} {2k^2-1 \over 4k^2-1}$ (Andres a Hochstättler 2011), což je pro velká $n$ zhruba ${1\over 2} n$ (formálně řečeno $\lim_{n\to\infty} {\mathbb{E}_n({\rm g}) \over n} = {1\over 2}$).

Následuje graf skóre $100$ běhů hladového algoritmu pro $n=200$.

Rekurzivní hladové řešení (Andres a Hochstättler 2011)

Rekurzivní hladové řešení postupuje následovně: Z řady aut odstraní první auto a auto stejného typu. Zbytek aut rekurzivně obarví a pak do řady přidá odebranou dvojici tak, aby byl počet změn co nejmenší možný.

Střední hodnotu počtu změn barev tohoto algoritmu lze odhadnout pomocí: ${2\over 5}\ n - {8\over 15} \le \mathbb{E}_n({\rm rg}) \le {2\over 5}\ n + {7\over 10}$.

Hvězdičkové rekurzivní řešení (Šámal et al. 2019)

Při běhu rekurzivního hladového řešení se občas stane, že některé dvojice aut lze obarvit oběma způsoby a počet změn bude stále stejný. Takovéto auta označíme novou barvou *. Přitom budeme dodržovat invariant, že * nikdy není na okrajích a nejsou dvě vedle sebe. Pokud budeme dále přidávat auto poblíž * a jedno nastavení * na barvu bude lepší, tak * přenastavíme na danou barvu.

Existuje domněnka, která tvrdí, že $\mathbb{E}_n({\rm srg}) \le 0.361 n$.

Řešení pomocí semidefinitního programování (Jiří Kalvoda, doposud nepublikováno)

Řešení se částečně podobá aproximačnímu algoritmu na maximální řez pomocí semidefinitního programování (Gärtner a Matousek 2012). V něm si můžeme představit vrcholy jako body $n-1$ rozměrné koule v $n$ dimenzionálním prostoru a semidefinitní programování nám je rozmístí tak, aby součet skalárních součinů bodů spojených hranou byl co možná největší. Pak náhodně vybereme nadrovinu procházející počátkem a vrcholy rozdělíme do partit podle toho na které straně nadroviny je daný bod.

V našem problému také necháme rozmisťovat body po $n-1$ rozměrné sféře. Body nyní budou reprezentovat typy aut. Přesněji řečeno každý bod bude reprezentovat první auto daného typu a druhé auto daného typu bude reprezentovat virtuální bod k němu středově symetrický dle počátku. Když pak tedy rozřízneme body nadrovinou (která žádným bodem neprochází) a autům přidělíme barvy dle toho, na které straně skončí, bude od každého typu právě jedno auto každé barvy. Semidefinitní programování se bude snažit maximalizovat součet skalárních součinů sousedních aut, tedy sousední auta budou blízko, což znamená, že je malá pravděpodobnost, že budou mít různou barvu.

Pro účely této práce vždy vyzkoušíme 10 náhodných řezů koule nadrovinu a z nich vybereme nejlepší možný.

Je známo, že tento algoritmus na každé instanci ve střední hodnotě přes nadroviny vrátí řešení, co je nejhůře o $0.22n$ horší než optimum na daném vstupu. Navíc pomocí hodnoty optimalizační funkce semidefinitního programu pro každou instanci můžeme získat i dolní odhad na skóre optima.

Statistická práce

Na níže uvedeném grafu jsou výsledky všech výše uvedených algoritmů pohromadě (pozor na to, že v grafu je zobrazena jen zajímavá část vertikální osy mezi $0.25 n$ a $0.6n$.

Z grafu můžeme odhadnout, že jako nejlepší algoritmus se jeví semidefinitní programování a po něm hvězdičkové rekurzivní řešení. V následující části bych rád ukázal, že tomu tak skutečně je.

To vypadá jako jednoduchá úloha pro Welchův t-test. Vygeneruji dalších 30 běhů každého algoritmu a na nich spustím:

import zapoctak_past.g as g, zapoctak_past.data_lib as data_lib
import scipy.stats

data = g.load("log-t-test")

rsg          = [ x.score for x in data.pipelines["rsg"] ]
semidef_prog = [ x.score for x in data.pipelines["semidef_prog(10)"] ]

print(scipy.stats.ttest_ind(rsg, semidef_prog, equal_var=False))

TtestResult(statistic=8.852830897345148, pvalue=3.4414910719986256e-12, df=55.35835190768672)

Z tohoto bychom mohli usuzovat, že řešení pomocí semidefinitního programování je na $99\%$ lepší, protože pvalue je dokonce o několik řádů menší než $0.01$.

V čem je tedy problém? Problém t-testu je že nesamplovaná data musí pocházet z normální distribuce a já vůbec netuším, jak bych dokazoval, že naše data se jí alespoň blíží (pokud to vůbec platí).

Pojďme tedy využít trošku více dřevorubecké řešení, kterým to dokážeme.

Budeme se snažit dokázat, že algoritmus pomocí semidefinitního programování je lepší než hvězdičkové rekurzivní řešení. Tedy:
$H_0:$ Hvězdičkové řešení je stejně dobré nebo dokonce lepší.
$H_1:$ Řešení pomocí semidefinitního programování je lepší.

Víme, že rozptyl náhodné veličiny skóre řešení je nejvýše $n^2$, protože se jedná o náhodnou veličinu z intervalu $0$ až $2n$. Když víme, že $\sigma^2(X) \le n^2$, tak $\sigma^2(X_1 + \cdots + X_k) \le kn^2$ (pro nezávislé $X_1,\dots,X_n$, což v našem případě jsou, protože vybíráme vstupy nezávisle). Tedy $\sigma^2(\overline{X_k}) = \sigma^2({X_1 + \cdots + X_k \over k}) \le {k\cdot n^2 \over k^2} = {n^2 \over k}$.

Dle Cebyševovy nerovnosti tedy (pro $t>0$):  $$P\left(\mathbb{E}(X) \ge \overline{X_k} + t {n^2 \over k}\right) \le {1 \over 2 t^2}$$  Což můžeme využít následovně:

Nechť $X$ a $Y$ jsou náhodné veličiny skóre dvou algoritmů (semidefinitního programování a hvězdičkového rekurzivního), kde naměříme $\overline{X_k} < \overline{Y_k}$.  $$P(H_0) = P\left( \mathbb{E}(X) \ge \mathbb{E}(Y) \right) = 1 - P\left( \mathbb{E}(X) < \mathbb{E}(Y) \right) \le 1 - P\left( \mathbb{E}(X) < { \overline{X_k}+\overline{Y_k} \over 2} < \mathbb{E}(Y) \right) \le$$   $$\le 1 - P\left( \mathbb{E}(X) < { \overline{X_k}+\overline{Y_k} \over 2} \right) + 1 - P\left( { \overline{X_k}+\overline{Y_k} \over 2} < \mathbb{E}(Y) \right) \le P\left( \mathbb{E}(X) \ge \overline{X_k} + { \overline{Y_k}-\overline{X_k} \over 2} \right) + P\left( \overline{Y_k} - { \overline{Y_k}-\overline{X_k} \over 2} \ge \mathbb{E}(Y) \right)$$

Nyní použijeme výše uvedenou nerovnost pro $t{n^2\over k} = { \overline{Y_k}-\overline{X_k} \over 2}$, tedy $t = k {\overline{Y_k}-\overline{X_k} \over 2n^2}$. Symetrická nerovnost funguje i pro $Y$ a otočený směr odhadu.

$$P\left( H_0 \right) \le 2 {1 \over 2 \left({k(\overline{Y_k}-\overline{X_k}) \over 2n^2}\right)^2} = {4n^4 \over k^2 \left({\overline{Y_k}-\overline{X_k} }\right)^2}$$

Pokud chceme, abychom věděli, že algoritmus $X$ je lepší než algoritmus $Y$ (pro danou velikost vstupu) s jistotou $99\%$, musí platit, že:

$${4n^4 \over k^2\left({\overline{Y_k}-\overline{X_k} }\right)^2} < 0.01$$

Z výše uvedených grafů můžeme doufat, že $\overline{Y_k}-\overline{X_k}$ bude alespoň zhruba $10$. Pro $k = 500$ by tedy mohla být výsledná hodnota mohla dost malá.

Navrhneme tedy následující pokus: Spustíme oba dva algoritmy na $100\,000$ náhodně vybraných vstupů a spočítáme jejich průměry.

Pokud výsledek výše uvedené formule bude menší než $0.01$, můžeme zamítnout $H_0$ s pravděpodobností $99\%$.

import zapoctak_past.g as g, zapoctak_past.data_lib as data_lib
import scipy.stats


n = 200

data = g.load("log-correct-test")



rsg          = [ x.score for x in data.pipelines["rsg"] ]
semidef_prog = [ x.score for x in data.pipelines["semidef_prog(10)"] ]

k = len(rsg)
assert k == len(semidef_prog)

rsg_avg = sum(rsg) / k
semidef_prog_avg = sum(semidef_prog) / k

print(f"n: {n} k: {k}")
print("rsg:         ", rsg_avg)
print("semidef_prog:", semidef_prog_avg)

assert semidef_prog_avg < rsg_avg
print(4*n**4 / k**2 / (rsg_avg - semidef_prog_avg)**2)

n: 200 k: 100000
rsg:          74.84675
semidef_prog: 64.94531
0.006528046717635241

Nulovou hypotézu se nám tedy podařilo zamítnout a tudíž řešení pomocí semidefinitního programování je lepší.

Všimněte si, že pokud by platilo, že se jedná o normální rozdělení, tak nám stačí mnohem méně dat a máme mnohem větší jistotu.

References