Dynamische Programmierung - Basics

Geschrieben von Sandro Feuz.

Die Mathematik wird oft als die Königin der Wissenschaften bezeichnet. Analog dazu ist die Dynamische Programmierung, aus meiner Sicht, die Prinzessin der Algorithmik. Während die meisten anderen vorgestellten Algorithmen lediglich ein spezifisches Problem lösen (Dijkstra’s Algorithmus findet (nur) kürzeste Wege, Kruskal findet (nur) Minimale Spannende Bäume), ist die Dynamische Programmierung (DP) ungleich mächtiger. DP ist nicht einfach ein einzelner Algorithmus, sondern beschreibt vielmehr eine ganze Klasse von Algorithmen. Eine Art algorithmische Probleme anzugehen und zu lösen.

Ein motivierendes Beispiel

Wenn ich DP in einem Satz zusammenfassen müsste, dann wäre das

Verschwende keine Zeit Dinge doppelt zu berechnen.

Lass uns dies an einem Beispiel anschauen.

\begin{equation*} 0,1,1,2,3,5,8,\dots \end{equation*}

Das ist der Beginn der berühmten Folge von Fibonacci und die \(i\)-te Fibonacci-Zahl \(f_i\) (für \(i=0, 1, \dots\)) kann formal ausgedrückt werden als

\begin{align*} f_0 &= 0\\ f_1 &= 1\\ f_i &= f_{i-2} + f_{i-1} \end{align*}

Lass uns dies in Python implementieren. Die zwei Basisfälle (\(i=0\) und \(i=1\)) können wir in einer if-Bedingung abfangen und der allgemeine Fall kann eins zu eins aus der Definition übernommen werden.

def fib(i):
  if i < 2:
    return i
  else:
    return fib(i-2) + fib(i-1)

Wenn wir versuchen, so die Fibonacci-Zahlen \(f_{5}, f_{10}, f_{50}\) und \(f_{100}\) zu berechnen, merken wir schnell, dass das Programm viel zu langsam ist und für \(f_{100}\) selbst nach Stunden nicht fertig ist. Was ist passiert? Wir verletzen den Grundsatz der Dynamischen Programmierung und verschwenden Zeit damit, das Gleiche wieder und wieder zu berechnen.

Unsere Funktion \(fib\) arbeitet rekursiv und beim Aufruf von \(fib(100)\) wird \(fib(99)\) und \(fib(98)\) aufgerufen. Der \(fib(99)\) Aufruf ruft dann selbst wieder \(fib(98)\) und \(fib(97)\) auf. Und so weiter. Das ist auf dem folgenden Diagramm zusammengefasst.

fib(100)+
        |
        +---fib(99)+
        |          |
        |          +----fib(98)+
        |          |           |
        |          |           +...
        |          |
        |          +----fib(97)+
        |                      |
        |                      +...
        +---fib(98)+
                   |
                   +...

Insbesondere ist zu sehen, dass der \(fib(98)\) Aufruf doppelt passiert, und natürlich wird er auch zwei mal das Gleiche zurückgeben (die 98-te Fibonacci-Zahl). Wir wollen nun schauen, dass wir die ganzen Berechnungen, die beim ersten Aufruf von \(fib(98)\) gemacht werden, nicht über den Haufen werfen müssen. So wäre es doch viel geschickter den Rückgabewert zu speichern und dann beim zweiten \(fib(98)\) Aufruf ohne zusätzliche Berechnungen diesen einfach zurück zu geben. Gehen wir mal davon ausgehen, dass die drei “\(\dots\)”-Aufrufe im Diagramm alle gleich lange dauern (was einigermassen stimmt). Durch das Speichern des \(fib(98)\) Rückgabewertes sparen wir uns einen dieser drei Aufrufe. So reduzieren wir die Laufzeit unseres Programms um ganze 33%, und das lediglich durch das Speichern einer Zahl! Während das Programm vorher 10 Jahre gebraucht hätte, um die hunderste Fibonacci-Zahl zu berechnen, schaffen wir es nun in weniger als 7. Yay!!

Komplexitätsanalyse

Wir haben gesehen, dass die hunderste Fibonacci Zahl mit unserem Programm schon nicht mehr terminiert, und darum liegt die Annahme nahe, dass die Laufzeit exponentiell ist. Das ist auch nicht so schwierig zu sehen: Wir merken, dass wir im Code nur die Zahlen \(0\) und \(1\) zurück geben, und in den rekursiven Aufrufen lediglich zusamenzählen. Das heisst, am Schluss wird die Laufzeit des Programms sicher mindestens die Anzahl “Einsen” sein, die zusammengezählt werden. Das ist genau der Wert der \(n\)-te Fibonacci-Zahl selbst. Die \(n\)-te Fibonacci-Zahl ist ungefähr \(1.6^n\), also exponentiell in \(n\).

Fibonacci mit Dynamischer Programmierung

Wir können die obige Idee verallgemeinern. Jedes Mal wenn wir einen Fibonacciwert berechnet haben, speichern wir ihn ab. Bevor wir einen neuen Wert berechnen schauen wir nach ob wir ihn nicht schon kennen. Das Programm wird leicht komplizierter.

gespeichert = {}
def fib(i):
  if i < 2:
    return i
  if i not in gespeichert:
    gespeichert[i] = fib(i-2) + fib(i-1)
  return gespeichert[i]

Wenn wir dieses neue Programm ausführen, sehen wir, dass die hunderste Fibonacci-Zahl statt in 7 Jahren in Millisekunden berechnet wird: für jedes \(i\) von \(1\) bis \(n\) wird nur konstante Arbeit gemacht. Wow, lang lebe Prinzessin DP!

Memoization vs Bottom Up

Wir haben die wichtigste Idee von DP gesehen: Berechne das Gleiche nicht mehrfach, sondern speichere es ab. Oft gibt es zwei “Reihenfolgen”, in der die Werte abgespeichert werden können. Entweder von “unten her” (bottom up) oder von “oben her” (memoization (merke, kein “r”)).

Im DP Fibonacci Code vom vorherigen Abschnitt haben wir die Werte von oben her berechnet. Wir nehmen die schon existierende (langsame) rekursive Formulierung des Problems und speichern die Werte ab, kurz bevor wir sie zurück geben. Deshalb ist der obige Code eine Memoization Lösung für Fibonacci.

Der analoge Code für die bottom up Lösung, welche die Werte “von unten” her berechnet, überlegt sich welche Werte ich von Anfang an schon weiss. Das sind die Verankerungswerte der Rekursion, also \(f_0=0\) und \(f_1=1\). Danach überlegen wir uns: wenn wir jetzt die ersten zwei Werte wissen, können wir den dritten direkt berechnen? Aber ja, denn \(f_2 = f_0+f_1 = 0 + 1 = 1\). Jetzt haben wir die ersten drei Werte, können wir den vierten berechnen? Ja, und so weiter.

In Code könnte das ganze dann wie folgt ausschauen:

gespeichert = {0: 0, 1: 1}
for i in xrange(2, n+1):
  gespeichert[i] =  gespeichert[i-2] + gespeichert[i-1]

Wir verlassen die rekursive Formulierung komplett und füllen die Hilfstabelle direkt von unten her auf.

Meist spielt es keine grosse Rolle, ob wir eine Memoization oder Bottom Up Lösung schreiben, beide gehören zur Dynamischen Programmierung und typischerweise haben beide die gleiche Komplexität (wie auch in unserem Beispiel).

Die Vier Schritte der Dynamischen Programmierung

Verlassen wir nun das Fibonacci-Beispiel. Die Überlegungs-Schritte, die wir ausgeführt haben um zu einer effizienten Lösung zu gelangen, hat Cormen et al [1] etwas formaler wie folgt zusammengefasst

  1. Beschreiben - Die Struktur einer optimalen Lösung beschreiben.
  2. Definieren - Rekursive Formulierung einer optimalen Lösung.
  3. Berechnen - Berechnen des Wertes einer optimalen Lösung (Memoization oder Bottom Up).
  4. Ausgeben - Rekonstruieren der optimalen Lösung.

und wir wollen versuchen sie schrittweise an einem neuen Beispiel durchzumachen.

Maximum Triangular Sum Path

Die Aufgabe, die wir lösen werden, ist wahrscheinlich älter als du! Sie wurde an der Internationalen Informatikolympiade 1994 gestellt [2]:

Gegeben ist eine Zahlenpyramide der Höhe \(n\). Gesucht ist ein maximaler Pfad von der Spitze zur Basis, wobei in jedem Schritt diagonal nach links oder rechts gezogen wird. Der Wert eines Pfades ist die Summer seiner Zahlen. Beispiel:

    7
   3 8
  8 1 0
 2 7 4 5
4 5 2 6 5

Die Korrekte Ausgabe für das Beispiel ist 30, realisiert durch den Pfad 7, 3, 8, 7, 5.

0) Greedy is Wrong

Bevor wir uns auf eine Dynamische Programmierlösung stürzen, sollten wir sicherstellen, dass keine einfachere Lösungsstrategie funktionert. Ein erster Kandidat ist oft “Greedy”, was soviel heisst wie einfach immer die lokal beste Zahl zu nehmen. Im obigen Beispiel würden wir damit den Pfad 7, 8, 1, 7, 5 zurückgeben, der lediglich den Wert 28 hat. Also schlechter als die beste Lösung und somit falsch.

1) Beschreiben - Die Struktur einer optimalen Lösung

Denken wir darüber nach, was eine optimale Lösung charakterisiert. In jedem Schritt gehen wir eine Zeile tiefer und dabei eins nach links oder rechts. Betrachten wir nun den letzten Schritt in der optimalen Lösung.

    .
   . .
  . . .
 2 7 . .
. 5 . . .

Das Feld mit der 5 kann nur erreicht werden, indem man im Schritt davor entweder bei der 2 oder der 7 war. Das gilt für jeden Pfad bis zur 5, also insbesondere auch für einen optimalen Pfad bis zur 5. Wir haben im Beispiel oben gesehen, dass der optimale Pfad bis zur 5 von der 7 her kommt. Was wissen wir dann über den Teilpfad von der Spitze des Dreicks bis zu der 7? Aus allen Pfaden vom Start bis zu der 7 muss dieser Pfad einer sein mit maximaler Summe. Wieso? Es gilt ein “Cut & Paste” - Argument: Wäre dem nicht so, so gäbe es einen Pfad von der Spitze bis zu der 7 der eine höhere Summe hat. Dann könnten wir aber mit diesem Pfad einen insgesamt besseren Pfad bis zur 5 finden, indem man diesen besseren Pfad bis zur 7 geht verlängert. Dies ist aber ein Widerspruch, da wir von einer optimalen Lösung ausgegangen sind.

Mit anderen Worten: eine optimale Lösung bis zur 5 muss immer aus einer optimalen Lösung zu einem seiner zwei möglichen Vorgängerfelder bestehen. Diese Einsicht benutzen wir im nächsten Schritt.

2) Definieren - Rekursive Formulierung einer optimalen Lösung.

Versuchen wir das Problem rekursiv zu formulieren. Sei \(opt(i, j)\) die maximale Summe eines Pfades von der Spitze des Dreiecks bis zum Feld \((i, j)\), dies ist das \(j\)-te Feld in der \(i\)-ten Zeile. Z.B. ist die \(8\) in der dritten Zeile das Feld \((3, 1)\). Die Lösung (der Wert des besten Pfades von der Spitze zu einem beliebigen Feld in der untersten Zeile) wäre dann das Maximum der opt-Felder der letzten Zeile:

\begin{equation*} max_{j=1,\dots,n} opt(n, j). \end{equation*}

Wie kann \(opt(i, j)\) berechnet werden? Dazu überlegen wir uns, wie kann ich auf das Feld \((i, j)\) kommen? Entweder von dem Feld rechts \((i - 1, j)\) oder links \((i - 1, j - 1)\) in der Zeile darüber. Da wir von 1) wissen, dass ein optimaler Weg aus einem optimalen Teilweg zu einem seiner zwei Vorgängerfelder bestehen muss, können wir \(opt(i, j)\) nun direkt angeben: Es ist der bessere der beiden optimale Wegen zu den Vorgängern. Also:

wobei \(T[i][j]\) die Zahl angibt, die in der Pyramide im Feld \((i, j)\) steht. Für das Startfeld \((1,1)\) selbst gibt es nur einen Pfad (der dann natürlich auch der optimale bis dort ist):

\begin{equation*} opt(1, 1) = T[1][1]. \end{equation*}

Diese rekursive Formulierung kann fast direkt in Code übersetzt werden. Hier ist eine Beispielsimplementierung in C++ (wir ändern die Nummerierung der Felder, so dass sie nun bei 0, statt bei 1 beginnt).

std::vector<std::vector<int>> T = {
{     7       },
{   3,  8     },
{  8, 1,  0   },
{ 2, 7, 4, 5  },
{4, 5, 2, 6, 5}};

// Gib den Wert eines optimalen Pfades bis zum Feld in
// Zeile i Spalte j zurück.
int MaximumTriangularSumPathRecursive(int i, int j) {
  // Die Spitze der Pyramide.
  if (i == 0 && j == 0) {
    return T[i][j];
  }
  int left = 0;
  int right = 0;

  // Alle Felder ganz links haben keinen linken Vorgänger.
  if (j > 0) {
    left  = MaximumTriangularSumPathRecursive(i - 1, j - 1);
  }
  // Alle Felder ganz rechts haben keinen rechten Vorgänger.
  if (j < i) {
    right = MaximumTriangularSumPathRecursive(i - 1, j);
  }
  return std::max(left, right) + T[i][j];
}

Aufgabe für den Leser: Was ist die Laufzeit dieser rekursiven Lösung? Findet das Programm noch eine Lösung in absehbarer Zeit bei einem Dreieck mit 100 Zeilen?

3) Berechnen - Berechnen des Wertes einer optimalen Lösung (Memoization oder Bottom Up).

Spoiler Alert: Das Programm wird nicht nützlicher Zeit fertig werden, weil die Laufzeit \(O(2^n)\) (exponentiel) ist, und darum für \(n=100\) die Anzahl Operationen viel zu hoch (eine weitere Aufgabe für den Leser: Wie kommt man auf die Laufzeit von \(O(2^n)\)?)

Was jetzt? Nun, die Lösung ist natürlich wieder Dynamische Programmierung. Wir wählen diesmal direkt eine Bottom Up Lösung, ich möchte aber darauf hinweisen, dass Memoization auch funktionieren würde.

Für die Berechnung verwenden wir Array \(opt[n][n]\) und berechnen Schritt für Schritt die Lösung. Das heisst, wir fangen in der Spitze der Pyramide an (wo es nur genau einen Pfad gibt), und arbeiten uns zeilenweise runter. Dadurch haben wir beim Bearbeiten der i-ten Zeile schon alle Zeilen

\begin{equation*} 0, \dots, i-1 \end{equation*}

bearbeitet, und können darum die rekursive Beschreibung der Lösung von 2) anwenden.

std::vector<std::vector<int>> T = {
{     7       },
{   3,  8     },
{  8, 1,  0   },
{ 2, 7, 4, 5  },
{4, 5, 2, 6, 5}};

// Gib den Wert eines optimalen Pfades von der Spitze bis zur Basis
// zurück.
int MaximumTriangularSumPathBottomUp() {
  const int n = T.size();
  vector<vector<int>> opt(n, std::vector<int>(n));
  opt[0][0] = T[0][0];
  for (int i=1; i < n; i++) {
    for (int j=0; j<=i; j++) {
      if (j == 0) {
        // Linker Rand.
        opt[i][j] = T[i][j] + opt[i-1][j];
      } else if (j == i) {
        // Rechter Rand.
        opt[i][j] = T[i][j] + opt[i-1][j-1];
      } else {
        opt[i][j] = T[i][j] + max(opt[i-1][j], opt[i-1][j-1]);
      }
    }
  }
  // Die Funktion max_element in <algorithm> gibt einen Iterator auf
  // das grösste Element zurück, und der Stern dereferenziert
  // den Iterator (nimmt den wert des Iterators).
  return *std::max_element(opt[n-1].begin(), opt[n-1].end());
}

Wie schnell ist diese Lösung? Wir merken, dass wir für jedes Feld der Pyramide lediglich konstant viel Arbeit haben. Da es \(\frac{n^2}{2}\) Felder gibt, ist die Laufzeit \(O(n^2)\) (also linear in der Eingabegrösse, da wir ja auch \(O(n^2)\) viele Zahlen in der Eingabe haben für eine Pyramide der Höhe \(n\)).

4) Ausgeben - Rekonstruieren der optimalen Lösung.

Wir haben nun eine effiziente Lösung, um den Wert des optimalen Pfades zu finden, und je nach Aufgabenstellung sind wir schon fertig. Aber unter Umständen ist das nicht genug, und wir möchten neben dem Wert eines optimalen Pfads auch noch den Pfad selbst berechnen können. Dazu finden wir zuerst das Ende des optimalen Pfades (wie oben hergeleitet ist es das Maximum der Felder der letzten Zeile). Dann schauen wir rückwärts von welchem Feld aus der optimale Wert erreicht wurde, und gehen so von unten rückwärts bis zur Spitze der Pyramide. Die Lösung sollte in folgendem Code klarer werden:

// T = Wie oben...
// Gib einen optimalen Pfad für T als vector der Länge n zurück, wobei
// path[i] das Verwendete Feld in der i-ten Zeile ist.
std::vector<int> MaximumTriangularSumPathBottomUpWithPath() {

  // Wie oben...

  // Der Pfad besteht aus den Feldern mit Koordinaten
  //   (0, path[0]), (1, path[1]), ..., (n-1, path[n-1])
  std::vector<int> path(n);

  // Finde das Feld in der letzten Zeile, in welcher der optimale Pfad endet.
  // Wie oben gibt std::max_element einen iterator auf dieses Feld zurück,
  // und std::distance gibt dann die Distanz zwischen dem ersten Element
  // und diesem Iterator zurück (also gerade den Spalten Index des Feldes,
  // beginnend bei 0).
  path[n-1] = std::distance(opt[n-1].begin(), std::max_element(opt[n-1].begin(), opt[n-1].end()));

  for (int i=n-2;i>=0;--i) {
    if (path[i+1] == 0) {
      // Feld ganz links hat nur einen möglichen Vorgänger.
      path[i] = 0;
    } else if (path[i+1] == i+1) {
      // Feld ganz rechts hat nur einen möglichen Vorgänger.
      path[i] = i;
    } else {
      if (opt[i+1][path[i+1]-1] > opt[i][path[i+1]]) {
        // Der Wert kommt vom linken Vorgänger.
        path[i] = path[i+1] - 1;
      } else {
        // Der Wert kommt vom linken Vorgänger.
        path[i] = path[i+1];
      }
    }
  }
  return path;
}

Weitere Materialien

Dynamische Programmierung ist ein sehr mächtiges aber auch ein schwieriges Thema, und ich glaube persönlich, dass es am einfachsten zu meistern ist, wenn man selbst möglichst viele Probleme selbst löst (oder wenigstens mindestens paar Stunden versucht zu lösen, bevor man die Lösung anschaut). Deshalb solltest du nach dem Artikel hier am besten mit Beispielen weitermachen. Eine Liste von Beispielaufgaben gibts zum Beispiel hier oder hier nochmals mit Theorie und Lösungen, oder du fragst einfach deinen Lieblings-SOI-Organisator (oder mich [3] :-)) für eine Liste seiner/ihrer interessantesten DP-Tasks.

[1]In dem Buch hier (deutsch) resp. hier (english)
[2]Aufgabe 1 von Tag 1 hier. Eine Version der Aufgabe ist auch auf Project Euler zu finden
[3]sandro.feuz@gmail.com