Advanced dynamic programming

Geschrieben von Daniel Rutschmann.

2-Dimensionales DP

In diesem Kapitel geht es darum, wie man die Technik des dynamischen Programmierens in höheren Dimensionen verwendet. Eine häufige Anwendung ist es, die Aufgabe auf jeder zusammenhängenden Teilsequenz zu lösen und mit der so gewonnenen Information die Lösung des ursprünglichen Problems zu berechnen. Hierzu einige Typische Aufgaben.

Längste palindromische Teilsequenz

Aufgabe

Gegeben ein String s, finde den längsten, nicht notwendigerweise zusammenhängenden Teilstring, welcher ein Palindrom ist. Ein String ist genau dann ein Palindrom, wenn er von links nach rechts gelesen genau dasselbe ist, wie von rechts nach links gelesen, wie z.B.: Kajak, Rentner oder Radar.

Lösungsansatz

Wie bereits in der Einleitung erwähnt, wollen wir in unserem DP die Aufgabe für jeden zusammenhängenden Teilstring lösen. Wir verwenden also folgenden DP-State:

\(DP[l][r] \colon\) Länge des längsten palindromischen Teilstrings des Teilstrings \([l, r]\). Die gesuchte Antwort ist nun \(DP[0][N-1]\).

Für die DP Berechnung gibt es nun zwei Fälle: * \(s[l] \neq s[r] \colon\) In diesen Fall können wir nicht s[l] und s[r] zusammen verwenden, da der resultierende String kein Palindrom wäre. Also ist

\begin{equation*} DP[l][r] = max(DP[l+1][r], DP[l][r-1]) \end{equation*}
  • \(s[l] = s[r]\) Hier gibt es zusätzlich noch die Möglichkeit, sowohl s[l], als auch s[r] zu verwenden. Also ist
\begin{equation*} DP[l][r] = max(2+DP[l+1][r-1], DP[l+1][r], DP[l][r-1]) \end{equation*}

Zusammen mit den Basisfällen: \(DP[l][l]=0\) und \(DP[l][l+1]=1\) lässt sich direkt eine rekursive Lösung mit Memoisation implementieren:

vector<vector<int> > cache(n, vector<int>(n, -1));
int longestPalindromicSubsequence(int l, int r){
    if(l>r) return 0; //Basisfall: leerer String
    if(l==r) return 1;//ein einzelner char ist ein Palindrom

    if(cache[l][r]!=-1) return cache[l][r];//Wert im cache

    if(s[l]!=s[r]) cache[l][r] = max(longestPalindromicSubsequence(l+1, r), longestPalindromicSubsequence(l, r-1)); //Fall 1
    else cache[l][r] = max(2 + longestPalindromicSubsequence(l+1, r-1), max(longestPalindromicSubsequence(l+1, r), longestPalindromicSubsequence(l, r-1)));//Fall 2

    return cache[l][r];
}

Da es jeweils \(n\) verschiedene Werte für l, bzw. r gibt, gibt es \(\Theta(n^2)\) verschiedene DP-States. Die Berechnung eines States braucht konstante Zeit und aufgrund der Memoisation wird jeder State nur einmal berechnet. Die Laufzeit ist also \(\Theta(n^2)\) Da die Cachetabelle \(n^2\) Einträge hat, brauchen wir \(\Theta(n^2)\) Speicher.

Die DP-Tabelle lässt sich auch iterativ berechnen:

vector<vector<int> > DP(n, vector<int>(n, -1));
for(int l=n-1;l>=0;--l){ //!!Es ist sehr wichtig, diese Schleife von rechts nach links zu machen, da wir sonst auf noch nicht berechntete Felder der DP-Tabelle zugreifen!!
    for(int r=0;r<n;++r){
        if(l>r){ //Basisfall: Leerer String
            DP[l][r]=0;
        } else if(l==r){//Basisfall: Einzelner char
            DP[l][r]=1
        } else if(s[l]!=s[r]){ //Fall 1
            DP[l][r] = max(DP[l+1][r], DP[l][r-1]);
        } else { //Fall 2
            DP[l][r] = max(2+DP[l+1][r-1], max(DP[l+1][r], DP[l][r-1]));
        }
    }
}
return DP[0][n-1];

Hier ist es noch einfacher einzusehen, dass wir \(\Theta(n^2)\) Zeit und \(\Theta(n^2)\) Speicher brauchen.

Zusatzaufgabe: Implementiere den oben beschriebenen DP mit \(O(n)\) Speicher.

Matrix Kettenmultiplikation

Ein in vielen Anwendungsbereichen auftretendes Problem ist die Multiplikation von Matrizen. Die Matrixmuliplikation ist assoziativ, d.h. für 3 Matrizen \(A,\ B,\ C\) gilt: \(A\times (B\times C) = (A\times B)\times C\) D.h. das Resultat ist unabhängig von der Klammerung. Die zur Multiplikation benötigte Zeit kann jedoch mit der Klammerung stark variieren, da der naive Algorithmus zur Multiplikation einer \(x\times y\) Matrix mit einer \(y\times z\) Matrix \(x*y*z\) elementare Multiplikationen benötigt. Ist \(A\) eine \(5\times 10\), \(B\) eine \(10\times 1\) und C eine \(1\times 100\) Matrix, so braucht die Klammerung \(A \times (B \times C)\) ganze \(6000\) Multiplikationen, während die Klammerung \((A \times B) \times C\) nur \(550\) Multiplikationen benötig, also deutlich weniger. Wenn wir also einen effizienten Weg finden, die optimale Klammerung zu berechnen, kann viel Zeit gespart werden. Da exponentiell viele Klammerungen möglich sind, dauert es viel zu lange, wenn wir jede mögliche ausprobieren.

Aufgabe:

Gegeben n Matrizen \(A_0, A_1, ..., A_{n-1}\), wobei \(A_i\) Grösse \(s[i] \times s[i+1]\) hat, berechne die minimale benötigte Anzahl von Multiplikationen zur berechnung von \(A_0 \times A_1 \times ...\times A_{n-1}\).

Lösungsansatz

Bei jeder Klammerung gibt es eine zuletzt ausgeführte Multiplikation. Zuvor wurden schon die \(k\) Matrizen links miteinander und die \((n-k)\) Matrizen rechts miteinander auf über eine optimale Klammerung multipliziert.

\begin{equation*} (.. k \ Matrizen ..) \times (.. (n-k)\ Matrizen ..) \end{equation*}

Die Aufgabe hat also optimale Teilstruktur. Die letzte Multiplikation braucht nun \(s[0]*s[k]*s[n]\) Zeit. Da wir nicht wissen, welche Multiplikation in der optimalen Klammerung zuletzt ausgeführt wurde, bleibt uns nichts Besseres übrig, als alle möglichen letzten Multiplikationen auszuprobieren.

Folgender DP-State ist als naheliegend: \(DP[l][r]\colon\) Minimale Anzahl elementarer Multiplikationen zur multiplikation der Matrizen \(A_l, A_{l+1}, ..., A_{r-1}\) Die gesuchte Antwort ist nun \(DP[0][n]\). Berechnet wird der DP durch Iteration über alle möglichen letzten Multiplikationen:

\begin{equation*} DP[l][r] = \min_{l\ < \ k\ < \ r} (DP[l][k] + DP[k][r] + s[l] * s[k] * s[r]) \end{equation*}

Mit dem Basisfall \(DP[l][l+1]=0\)

Rekursive Implementation mit Memoisation:

vector<long long> s(n+1);
vector<vector<long long> > cache(n+1, vector<long long>(n+1, -1));
const long long INF = numeric_limits<long long>::max();

long long matrixChainMultiplication(int l, int r){
    if(l==r+1) return 0;//Basisfall: Nur eine einzelne Matrix
    if(cache[l][r]!=-1) return cache[l][r];
    long long bestVal = INF;
    for(int k=l+1, k<r;++k){
        bestVal = min(bestVal, matrixChainMultiplication(l, k) + matrixChainMultiplication(k, r) + s[l]*s[k]*s[r]);
    }
    return cache[l][r] = bestVal;
}

Iterative Implementation

const long long INF = numeric_limits<long long>::max();
vector<long long> s(n+1);
vector<vector<long long> > DP(n+1, vector<long long>(n+1, INF));

for(int l=n-1;l>=0;--l){
    for(int r=l+1;r<=n;++r){
        if(l+1==r){//Basisfall: nur eine einzige Matrix;
            DP[l][r]=0;
        } else {
            for(int k=l+1;k<r;++k){
                DP[l][r]=min(DP[l][r], DP[l][k] + DP[k][r] + s[l]*s[k]*s[r]);
            }
        }
    }
}
return DP[0][n];

Da \(\Theta(n^2)\) DP-States vorhanden sind und wir für jeden State \(\Theta(n)\) Zeit benötigen, läuft dieser Algorithmus in \(\Theta(n^3)\) Zeit und \(\Theta(n^2)\) Speicher.

Zusatzaufgabe: Implementiere den oben beschriebenen DP in \(O(n)\) Speicher.

Challenge: Löse das Matrix Kettenmultiplikationsproblem in \(O(n\ log\ n)\)

Längste gemeinsame Teilsequenz

In verschiedenen Bereichen der Biologie müssen verschiedene DNA-Stränge verglichen werden. Ein DNA-Strang kann über einen String ausgedrückt werden, wobei die einzelnen Zeichen für die 4 Basen (A, C, G, T) stehen. Eine Möglichkeit, Zwei DNA Stränge zu vergleichen, um zu ermitteln, wie ähnlich diese sind, ist die Berechnung der längsten gemeinsamen Teilsequenz der Stränge.

Aufgabe

Gegeben ein String \(A\) = \(a_0,a_1,...,a_{n-1}\) der Länge n und ein String \(B\) = \(a_b, b_1,...,b_{m-1}\) der Länge m, berechne die Länge des längstmöglichen Strings \(C\) = \(c_0,c_1,...,c_{k-1}\), der sowohl eine Teilsequenz von A, als auch von B ist.

Lösungsansatz

Wir betrachten \(a_{n-1}\) und \(b_{m-1}\), das letzte Zeichen von \(A\) bzw. \(B\) . * Ist \(a_{n-1}\neq b_{m-1}\), so ist C die längste gemeinsame Teilsequenz von \(a_0,...,a_{n-2}\) und B (wenn \(c_{k-1} \neq a_{n-1}\)), oder von A und \(b_0,...,b_{m-2}\) (wenn \(c_{k-1} \neq b_{m-1}\)) * Ansonsten gilt \(a_{n-1} = b_{m-1}\) und es besteht noch die Mmöglichkeit, dass \(c_{k-1} = a_{n-1}\) und \(c_0,...,c_{k-2}\) die längstmögliche gemeinsame Teilsequenz von \(a_0,...,a_{n-2}\) und \(b_0,...,b_{m-2}\) ist. Die Aufgabe besitzt also ebenfalls optimale Teilstruktur und wir kommen zu folgendem Ansatz mit dynamischer Programmierung:

Sei \(DP[x][y]\) die Länge der längsten gemeinsamen Teilsequenz der ersten \(x\) Zeichen von \(A\) und der ersten \(y\) Zeichen von \(B\).

\begin{equation*} DP[x][y] = \begin{cases} max(DP[x-1][y], DP[x][y-1]) & a_x\neq b_y \\ max(1+DP[x-1][y-1], DP[x-1][y], DP[x][y-1])& a_x= b_y\\ \end{cases} \end{equation*}

Ausserdem ist \(DP[x][0]=DP[0][y]=0\) Die gesuchte Antwort ist nun \(DP[n][m]\).

Rekursive Implementation:

string A, B;
vector<vector<int> > cache(n+1, vector<int>(m+1, -1));
int longestCommonSubsequence(int x, int y){
    if(x==0 || y==0) return 0; //Basisfall: Ein String ist leer
    if(cache[x][y]!=-1)return cache[x][y];
    if(A[x]!=B[y]){// Fall 1
        return cache[x][y] = max(longestCommonSubsequence(x-1, y), longestCommonSubsequence(x, y-1));
    } else {// Fall 2
        return cache[x][y] = max(1 + longestCommonSubsequence(x-1, y-1), max(longestCommonSubsequence(x-1, y), longestCommonSubsequence(x, y-1)));
    }
}

Iterative Implementation:

string A, B;
vector<vector<int> > DP(n+1, vector<int>(m+1, -1));
for(int i=0;i<=n;++i){
    for(int j=0;j<=m;++j){
        if(i==0 || j==0){//Basisfall: Ein string ist leer
            DP[i][j]=0;
        } else if(A[i]==B[i]){// Fall 1
            DP[i][j]=max(DP[i-1][j], DP[i][j-1]);
        } else {// Fall 2
            DP[i][j]=max(1+DP[i-1][j-1], max(DP[i-1][j], DP[i][j-1]));
        }
    }
}
return DP[n][m];

Ähnlich wie beim Problem der längsten palindromischen Teilsequenz berechnen wir \(\Theta(nm)\) verschiedene DP-States mit jeweils konstantem Aufwand. Wir brauchen also \(\Theta(nm)\) Zeit und Speicher.

Zusatzaufgabe: Implementiere diesen DP in \(O(n+m)\) Speicher.

Zusatzaufgabe 2: Augmentiere den DP, sodass nicht nur die Länge der längsten gemeinsamen Teilsequenz berechnet wird, sondern diese Sequenz am Schluss auch ausgegeben wird. Schaffst du es in \(O(n\cdot m)\) Zeit und \(O(n+m)\) Speicher?

Zusätzliche Aufgaben:

  • Editierdistanz
  • 2D-Knapsack
  • Optimal binary search trees in \(\Theta(n^3)\)

Bitmask DP

Eine andere Möglichkeit, einen mehrdimensionalen DP-State zu verwenden, ist, dass man nebst der momentanen Position auch noch eine Bitmask speichert, die aussagt, welche Teilmenge man schon besucht hat. Diese Teilmenge kann eine Teilmenge von Knoten eines Graphen sein, oder aber eine Teilmenge der Bits einer Zahl. Eine Bitmask ist im Wesentlichen eine Ganzahl, wobei das i-te Element genau dann in der Bitmask enthalten ist, falls das i-te bit dieser Zahl gesetzt ist. Wenn unsere Bitmask also 0001001101 ist (in Binärdarstellung), dann sind die Elemente 0, 2, 3, 6 in der Teilmenge enthalten.

Traveling salesman problem

Maus Stofl will eine Städtereise durchführen. Dafür hat sie bereits einige Städte ausgesucht, die sie alle unbedingt besuchen will. Da Stofl gerne möglichst viel Zeit für die Besichtigung der Städte haben will, will sie die Reisezeit minimieren. Da es bis zu \(n!\) verschiedene Reiserouten gibt, gelingt es Stofl nicht, von selbst die optimale Route zu finden. Du musst ihr also helfen!

Aufgabe

Gegeben einen gewichteten Graphen G, sowie einen Startknoten \(s\) und einen Zielknoten \(t\), finde den kürzesten Pfad von \(s\) nach \(t\), der jeden Knoten von G genau einmal besucht. Ein Pfad, welcher jeden Knoten genau einmal besucht, heisst auch Hamiltonpfad.

Lösungsansatz

Unser DP sollte das Problem auf einem Teilpfad lösen. Um weiter zu rechnen, ist die Reihenfolge der Knoten im Teilpfad jedoch nicht weiter wichtig. Es reicht aus, wenn wir wissen, welche Knoten im Teilpfad enthalten sind und welches der letzte Knoten ist. Also verwenden wir folgenden DP-State: \(DP[curPos][mask] :\) Länge des kürzesten Pfades von \(s\) nach \(curPos\), welcher jeden Knoten in \(mask\) genau einmal besucht. Zur Berechnung des DP betrachten wir alle möglichen Vorletzten Knoten und nehmen davon den Besten:

\begin{equation*} DP[curPos][mask] = min _{last\in mask,\ last\neq curPos} (DP[last][mask \setminus \{curPos\}] + cost[last][curPos]) \end{equation*}

Wobei am Anfang \(DP[s][2^s] = 0\) und die gesuchte Antwort \(DP[t][2^n-1]\) ist.

Rekursive Implementation

vector<vector<int> > cost(n, vector<int>(n));
vector<vector<long long> > cache(n, vector<long long>(1<<n, -1));
const long long INF = numeric_limits<long long>::max()/2;
long long TSP(int curPos, int mask){
    if(curPos == s && mask == (1<<s))return 0;//Basisfall beim Startknoten
    if(cache[curPos][mask]!=-1){
        return cache[curPos][mask];
    }
    long long bestVal = INF;
    int prevMask = mask^(1<<curPos);//Bitmask der vorherigen Knoten
    for(int last=0;last<n;++last){
        if(prevMask & (1<<last)){// Iteriere über alle möglichen Vrgägngerknoten.
            bestVal = min(bestVal, TSP(last, prevMask) + cost[last][curPos]);
        }
    }
    return cache[curPos][mask] = bestVal;
}

Da wir \(\Theta (n\ 2^n)\) States berechne mit \(\Theta(n)\) Aufwand für jeden State, brauchen wir \(\Theta(n^2 2^n)\) Zeit und \(\Theta(n\ 2^n)\) Speicher.

Iterative Implementation

vector<vector<int> > cost(n, vector<int>(n)); //Adjazenzmatrix der Kantengewichte
const long long INF = numerick_limits<long long>::max()/2;//Division durch 2, um overflows bei der addition zu vermeiden.
vector<vector<long long> > DP(n, vector<long long>(1<<n, INF));
DP[s][1<<s]=0;
for(int iteration=2;iteration!=n;++iteration){//in jeder ieteration berechnen wir die DP-Feldner, welche genau iteration Elemente in der bitmask haben.
    for(int mask=0;mask < (1<<n);++mask){
        if(__builtin_popcount(mask) != iteration) continue; //Falls die Bitmask nicht genau iteration Elemente enthält, wollen wir den State nicht jetzt berechen.
        for(int curPos=0;curPos<n;++curPos){
            if(!(mask & (1<<curPos)))continue; //Die momentane Position sollte näturlich in der Bitmask enthalten sein
            int prevMask = mask^(1<<curPos);//Bitset der vorherigen Elemente
            for(int last=0;last<n;++last){
                if(!(prevMask & (1<<last))) continue;//auch die vorherige Position sollte in der Bitmask enthalten sein.
                DP[curPos][mask] = min(DP[curPos][mask], DP[last][prevMask] + cost[last][curPos]);
            }
        }
    }
}
return DP[t][(1<<n)-1];

Man könnte aufgrund der 4 verschachtelten for-Schleifen meinen, dass die Laufzeit \(\Theta(n^3 2^n)\) beträgt. Aufgrund der Abbruchbedingung auf Zeile 7 Wird jedoch jeder der \(\Theta(n^2 2^n)\) States nur einmal berechnet. Da wir für jeden State über alle potenziellen Vorgängerknoten iterieren, brauchen wir für jeden State \(\Theta(n)\) Zeit Gesamthaft also \(\Theta(n^2 2^n)\) Zeit und \(\Theta(n\ 2^n)\) Speicher. Dies ist immer noch exponentiell viel, aber doch deutlich besser, als die \(\Theta(n!)\) Zeit, die ein Bruteforce Algorithmus bräuchte. Ausserdem ist es nicht möglich, dieses Problem in polynomieller Laufzeit zu lösen, sofern nicht (P=NP) gilt.