DP Optimization

Geschrieben von Daniel Rutschmann.

In diesem Script geht es um Techniken, die dazu dienen, um die Laufzeit eines DPs von \(\Theta(n^2)\) auf \(\Theta(n\ log(n))\) oder \(\Theta(n)\) bzw. vom \(\Theta(kn^2)\) auf \(\Theta(kn\ log(n))\) zu reduzieren.

Binäre Suche

Die Einfachste aller Optimierungstechniken verwendet die binäre Suche. Nehmen wir an, wir haben einen DP der folgenden Form:

\begin{equation*} DP[l][r] = min_{l\leq k\leq r} (cost[k] + max(DP[l][k], DP[k][r])) \end{equation*}

Mit \(DP[l][l+1]=0\) und gesucht ist \(DP[0][N]\). Ein naiver Weg diesen zu berechnen, ist, für jeden State über alle mögliche k-Werte zu iterieren und den Besten zu nehmen. Da dies \(\Theta(n)\) Zeit für jeden State braucht, läuft unser DP in \(\Theta(n^3)\). Nun betrachten wir unsere DP-Gleichung ein wenig genauer: Gilt für ein \(k_1\) \(DP[l][k_1] > DP[k_1][r]\), dann gilt für alle \(k_2 > k_1\)

\begin{equation*} max(DP[l][k_2], DP[k_2][r]) \geq DP[l][k_2] > DP[l][k] = max(DP[l][k], DP[k][r]) \end{equation*}

D.h. für wir müssen alle k-Werte grösser als \(k_1\) gar nicht mehr überprüfen, da sie sicher nicht besser als \(k_1\) sind. Analog müssen wir falls \(DP[l][k_3] < DP[k_3][r]\) alle k-Werte kleiner als \(k_3\) nicht mehr überprüfen. Wir können also das optimale k über eine binäre Suche finden, indem wir irgendein k in der Mitte überprüfen und je nachdem welcher DP-Wert im maximum ist, die Linke oder Rechte hälfte nach dem gleichen Prinzip absuchen. Da wir so nur \(\Theta(log(n))\) k-Werte überprüfen, hat sich unsere Laufzeit auf \(\Theta(n^2\ log(n))\) verbessert.

Implementation

vector<vector<long long>(n+1, vector<long long>cache(n+1, -1));
long long DP(int l, int r){
    if(l+1==r) return 0;
    if(l+2==r) return cost[l+1]; // nur eine einzelner möglicher k-Wert
    if(cache[l][r]!=-1) return cache[l][r];
    int a=l+1, b=r-1;
    while(a+1<b){
    int mid = (a+b)/2;
        if(DP[l][mid] > DP[mid][b]){
            b=mid
        } else {
            a=mid;
        }
    }
    cache[l][r] = min(cost[a]+max(DP(l, a), DP(a, r)), cost[b]+max(DP(l, b), DP(b, r)));
    return cache[l][r];
}

Aufgaben:

  • Bottles
  • Costly binary search (Small Subtask)

Convex hull optimization

1-Dimensionaler Fall

Für diese Optimierungstechnik betrachten wir die Aufgabe Good Inflation. Es ist nicht schwierig auf einen DP Ansatz in \(\Theta(n^2)\) zu kommen: \(DP[i] :\) Maximale erreichbare Ballongrösse nach i Minuten Um Sonderfälle, wie z.B. Ballone negativer grösse zu vermeiden, fügen wir am Anfang und Ende jeweis noch ein Pseudoangebot ein (d.h. \(a_0 = 0, d_0=0\ a_{N+1}=0, d_{N+1}=0\)). Die gesuchte Antwort ist hiermit \(DP[N+1]\). Berechnet wird der DP durch Iteration über alle möglichen zuletzt akzeptierten Angebote:

\begin{equation*} DP[i] = a[i] + max_{k<i}(DP[k] - (i-k)*d[k]) \end{equation*}

Da \(N\) jedoch bis zu \(10^6\) gross sein kann, ist ein schnellerer Ansatz gesucht. Schreiben wir hierzu die DP-Formel etwas um:

\begin{align} DP[i] &=a[i]+&max_{k<i}&(DP[k] - (i-k)*d[k]) \\ &=a[i] + &max_{k<i}&(DP[k] - i*d[k]+k*d[k]) \\ &=a[i] + &max_{k<i}&(i*(-d[k])\quad+(DP[k]+k*d[k]))\\ \end{align}

Unser DP versucht also eigentlich die lineare Funktion der Form \((-d[k]) * i + (DP[k]+k*d[k])\) zu maximieren. Wenn wir nun eine Datenstruktur hätten, die folgenden Operationen in \(\Theta(log(n))\) unterstützt, dann könnten wir unseren DP in \(\Theta(n\ log(n))\) berechnen.

  • Füge eine neue lineare Funktion hinzu
  • Finde das Maximum aller bisher eingefügten linearen Funktion für ein bestimmtes x

Eine solche Datenstruktur existiert sogar. Betrachten wir einige linearen Funktionen grafisch:

Wie wir sehen, gibt es manchmal gewisse lineare Funktionen, die für kein x-Wert zu einem Maximum beitragen (orange). Der grün eingefärbte Teil ist jeweils für jeden x-Wert diejenige Funktion, welche dort das Maximum erreicht. Da dieser wie die untere Hälfte einer konvexen Hülle aussieht, nennt man diese Technik auch konvexe Hülle Optimierung. Wie wir sehen, sind die einzelnen Segmente der grünen Teilhülle aufsteigen nach Steigung sortiert, dies ist auch im Allgemeinen der Fall. Wenn wir nun eine Liste der teilweise grünen Funktionen speichern in der Reihenfolge, in welcher sie das Maximum erreichen, für jede Funktion zudem noch speichern, von wo bis wo sie das Maximum ist, dann können wir mit einer binären Suche das Maximum für einen bestimmten x-Wert finden. Auf diese weise können wir also die 2. Operation in \(\Theta(\log(n))\) durchführen.

Wenn wir nun eine neue Funktion hinzufügen, fügen wir sie zuerst an ihren korrekten Platz hinzu und kümmern uns erst danach um die Funktionen, die nun nicht mehr in der Hülle enthalten sind. Um eine Funktion bei einer beliebigen Position einfügen zu können, müssen wir die Funktionen in einem balancierten Baum, wie z.B. einem std::multiset speichern. Nachdem wir die neue Funktion eingefügt haben, schauen wir zuerst, ob sie überhaupt Teil der Hülle ist. Wenn nicht, dann löschen wir sie wider aus der Multiset und sind fertig. Ansonsten kann die Funktion noch andere Funktionen links oder rechts von ihr überdecken. Ist dies der Fall, so löschen wir diese Funktionen, solange bis es keine Überdeckungen mehr gibt. Es könne zwar im schlimmsten Fall \(\Theta(n)\) Funktionen auf einmal gelöscht werden, doch da jede Funktion nur einmal eingefügt und einmal gelöscht wird, ist die amortisierte Laufzeit einer Einfügeoperation \(\Theta(log(n)))\)

Implementation:

const long long xLeftLeftest = -(1LL<<61);//default wert für die linkeste line
const long long qQuery = -(1LL<<60);
/*Spezieller Wert für q,
wird für die binäre Suche im std::multiset verwendet,
da dieses nur binäre Suche mithilfe des < Operator erlaubt.
*/

struct Line{
  long long m, q; //Steigung und y-Achsenschnittpunkt
  mutable long double xLeft;//xCoordinate des Schnittpunktes mit der Line links von dieser
  Line(long long _m, long long _q):m(_m), q(_q), xLeft(xLeftLeftest){}
  bool operator<(const Line& other)const{
    if(q==qQuery){//Speziallfall für die binare Suche
      return m < other.xLeft;
    }
    if(other.q==qQuery){//Speziallfall für die binare Suche
      return xLeft < m;
    }
    return m < other.m;//Normaler < Operator für die sortierung des multisets.
  }
   void recalcXLeft(const Line & pre)const{// Berechnet der Schnittpunkt mit der linken Linie neu
    xLeft = -((long double)pre.q-q ) / (pre.m-m);
  }
};

struct Hull{
  multiset<Line> slopes;// multiset unserer Funktionen.

  bool bad(multiset<Line>::iterator it){//überprüft, ob eine Linie vollständig überdeckt wird
    auto suc = next(it);//Rechte Nachbarlinie
    if(it==slopes.begin()){
      if(suc==slopes.end()) return false;//Nur eine einzelne Linie vorhanden
      return it->m==suc->m && it->q <=suc->q;//Linie am linken Rand
    }
    auto pre = prev(it);//Linke nachbarlinie
    if(suc==slopes.end()){
      return it->m==pre->m && it->q <= pre->q;//Linie am rechten Rand
    }
    return ((long double)it->q - suc->q) / (suc->m - it->m) <= ((long double)pre->q - it->q) / (it->m - pre->m);//Linie in der Mitte
    //return ((__int128_t)it->q - suc->q) * (it->m - pre->m) <= ((__int128_t)pre->q - it->q) * (suc->m - it->m);//Sicherer, geht aber nur, falls ein 128-bit Integer exestiert.

  }
  void insert(Line const & l){//fügt eine neue Linie hinzu
    multiset<Line>::iterator e = slopes.insert(l);
    if(bad(e)){//Neue linie vollständig überdeckt.
      slopes.erase(e);
      return;
    }
    while(next(e)!=slopes.end() && bad(next(e)))//überdeckungen auf der rechten Seite
    slopes.erase(next(e));
    if(next(e)!=slopes.end()){next(e)->recalcXLeft(*e);}//Schnittpunkt neu berechen.
    while(e!=slopes.begin() && bad(prev(e)))//überdeckungen auf der rechten Seite
    slopes.erase(prev(e));
    if(e!=slopes.begin()){//linken Schnittpunkt neu berechnen.
      e->recalcXLeft(*prev(e));
    } else e->xLeft = xLeftLeftest;//Spezialfall am linken Rand
  }

  long long query(long long x){//Frag eine linie ab.
    auto e = slopes.upper_bound(Line(x, qQuery));
    assert(e!=slopes.begin());//Dies wäre eine x Wert kleiner als xLeftLeftest, was wir nicht wollen.
    --e;
    return e->m * x + e->q;
  }
};

2-Dimensionaler Fall

Die konvexe Hülle Optimierung lässt sich auch auf zweidimensonale DPs der Form:

\begin{equation*} DP[i][j] = \mathrm{maxim}_{k<i} (DP[k][j-1]+q[k] + m[k]*i) \end{equation*}

anwenden. Hierzu baut man einfach für jede Koordinate der zweiten Dimension eine eigene konvexe Hülle Datenstruktur.

Deque Optimierung

In bestimmten Spezialfällen läst sich die konvexe Hülle Optimierung in \(\Theta(1)\) amortisierter Laufzeit pro Abfrage implementieren, dies geschieht dann mit einem Deque: Im Spezialfall, dass die Funktionen mit monoton steigenden m werten hinzugefügt werden und die Abfragen monoton steigende x-Werte besitzen ist hier implementiert:

struct Line{
/*  Fast das gleiche struct,
*   mit der Ausnahme,
*   dass wir die linken schnittpunkte nicht speichern müssen,
*   da wir keine binäre Suche druchführen.
*/
  long long m, q;
  Line(long long _m, long long _q):m(_m), q(_q){}
  long long eval(long long x){
    return m*x+q;
  }
};

struct monoHull{
  deque<Line> data;
  bool badBack(Line const & newLine){//Da wir neue Funktionen nur am rechten ende hinzufügen, müssen wir auch nur dort auf überlappungen testen.
    if(data.size()<2) return false;
    Line secondLast = data[(int)data.size()-2];
    Line last = data.back();
    return ((__int128_t)last.q - newLine.q) * (last.m - secondLast.m)  <= ((__int128_t)secondLast.q - last.q) * (newLine.m - last.m);//check x intersection
    //return ((long double)last.q - newLine.q) / (newLine.m - last.m) <= ((long double)secondLast.q - last.q) / (last.m - secondLast.m);//Falls kein 128 bit integer vorhanden ist

  }

  void insert(Line const & l){
    assert(data.empty() || l.m > data.back().m);//steigungen müssen streng monoton wachsend sein!
    while(badBack(l))
      data.pop_back();
    data.push_back(l);//Neue Funktion am ende hinzufügen.
  }
  long long query(long long x){
    if(data.empty()) return -~0ull>>4;
    while(data.size()>1 && data[0].eval(x)<data[1].eval(x))
      data.pop_front();// anstelle der binären Suche wird ein Pointerwalk durchgeführt.
    return data.front().eval(x);
  }
};

Aufgaben:

Allgemeine Deque Optimierung

Diese Deque Technik lässt sich auch verallgemeinern auf DPs, welche keine linearen Funktionen beinhalten, aber dennoch gewisse Eigenschaften besitzen: DP hat die Form:

\begin{equation*} DP[i] = max_{k<i}(DP[k] + cost[k][i]) \end{equation*}

Wobei cost irgendeine Funktion in Abhängigkeit von k und i sein kann. Erfüllt cost die konkave Vierecksungleichung

\begin{equation*} cost[a][d]+cost[b][c]\geq cost[a][c]+cost[b][d]\\ \forall a, b, c, d\quad a<b<c<d \end{equation*}

so könne wir eine deque aller potentieller Werte \(a_0, a_1, ..., a_l, \ a_0 <a_1<...<a_l\) speichern. Aufgrund der Vierecksungleichung soll momentan \(a_{x-1}\) ein besserer Kandidat sein als \(a_x\). (Da sonst \(a_{x-1}\) nie mehr besser sein wird als \(a_x\) und wir \(a_{x-1}\) somit wegwerfen könten.) Nun berechnen wir (allenfalls mit einer binären Suche) für jedes \(x\) den zeitpunkt \(i_x\), zu dem die Wahl von \(a_x\) besser sein wird, als die Wahl von \(a_{x-1}\). Wenn wir nun ein neuer potentieller wert \(a_{l+1}=i\) einfügen, berechne wir \(i_{l+1}\) und entfernen \(a_l\) solange \(i_{l+1} < i_l\) gilt. (Da zu dem Zeitpunkt, bei dem \(a_l\) besser als \(a_{l-1}\) wäre \(a_{l+1}\) bereits besser als \(a_l\) ist!)

Aufgaben: Einsenden auf Yandex

  • IOI 2016 Aliens (60 Punkte)
  • IOI 2015 Teams (sehr Schwierig)

Divide & Conquer

Die im letzten Kapitel erwähnten Optimierungstechniken sind nicht ganz einfach zu implementieren. Für 2-Dimensionale DP-Probleme gibt es eine einfachere Technik zur berechnung von DPs, dessen Kostenfunktion die konkave Vierecksungleichung erfüllt. Unser DP hat folgende Form.

\begin{equation*} DP[i][j] = max_{k<i}(DP[k][j-1] + cost[k][i]) \end{equation*}

Normalerweise müssen wir für jeden i-Wert über alle mögliche k-Werte iterieren. Doch wenn wir den DP in einer geschickteren Reihenfolge berechnen, können wir erheblich Zeit einsparen: Hierzu berechnen wir zuerst \(DP[N/2][j]\). Seit \(k_{N/2}\) der optimale k-Wert bei dieser Berechnung. Aufgrund der Vierecksungleichung wissen wir nun, dass für alle \(i<N/2\) nur noch die k-Werte in \([0, k_{n/2}]\) besser als \(k_{N/2}\) sein können. Analog können für alle \(i>N/2\) nur noch die k-Wert in \([k_{N/2}, N]\) besser als \(k_{N/2}\) sein. Wenn wir diese Technik rekursiv anwenden, dann erhalten wir eine Laufzeit von \(\Theta(n\ log(n))\), da auf jeder der \(\Theta (log(n))\) Rekursionsebenen jeder der k-Werte höchstens zweimal in betracht gezogen wird.

Implementation

const long long INF = ~0ull>>2;
vector<vector<long long> > DP(N+1, vector<long long>(M+1, INF));
/// Berechnet den DP für alle i in [leftI, rightI] mit mögliche k-Werten in [leftK, rightK]
void calcRec(int const &j, int leftI, int rightI, int leftK, int rightK){
    if(leftI>rightI) return; // Keine zu berechnenden Werte mehr übrig
    long long bestVal = INF;
    int bestPos = -1;
    int i = (leftI+rightI)/2;// Berechen den DP für die Position in der Mitte.
    for(int k=leftK;k<=min(i-1, rightK);++k){
        if(DP[k][j-1] + cost[k][i] < bestVal){
            bestVal = DP[k][j-1] + cost[k][i];
            bestPos = k;
        }
    }
    DP[i][j] = bestVal;
    calcRec(j, leftI, i-1, leftK, bestPos);
    calcRec(j, i+1, rightI, bestPos, rightK);
}
long long calcDP(){
    DP[0][0]=0;
    for(int j=1;j<M;++j){
        calcRec(j, 0, N, 0, N);
    }
    return DP[N][M];
}

Aufgaben

  • Ciel and Gondolas Achtung, riesiger Input! Input als String einlesen.
  • IOI 2016 Aliens (60 Punkte)

Minimum anstelle von Maximum

Die beschriebenen Techniken funktionieren allesamt auch für die Minimierung eines DPs. Dieser Fall läst sich auf dem Maximumsfall reduzieren, indem du statt \(DP[i]\) einfacht \(-DP[i]\) berechnest, und \(-cost[i][k]\) anstelle von \(cost[i][k]\) verwendest.

Für den Fall, dass die Kostenfunktion die konvexe und nicht die konkave Vierecksungleichung erfüllen sollte, kann die Optimierung trotzdem angewendet werden, es muss dazu lediglich das Deque umgedreht werden. (Sodass danach \(a_0 > a_1 > ... > a_l\)). Für die Divide & Conquer Optimierung müssen nur die letzten zwei Zeilen zu

calcRec(j, leftI, i-1, bestPos, rightK);
calcRec(j, i+1, rightI, leftK, bestPos);

angepasst werden.