5. Arithmetische Algorithmen

5.2 Schnelle Matrix-Exponentiation

Ich zeige Ihnen jetzt eine trickreiche Methode, Fibonacci-Zahlen noch schneller zu berechnen. In Ihrem itertativen Code haben Sie wahrscheinlich etwas in der Art geschrieben:

        newA = b 
        newB = a + b
        a = newA
        b = newB

Mathematisch ausgedrückt: jeder Schleifendurchlauf bildet das Paar $a,b$ auf ein neues Paar ab, nämlich $$ \begin{pmatrix} a \\ b \end{pmatrix} \mapsto \begin{pmatrix} b \\ a + b \end{pmatrix} \ \ . $$ Dies ist eine lineare Abbildung und kann als Matrix dargestellt werden: $$ \begin{pmatrix} a \\ b \end{pmatrix} \mapsto \begin{pmatrix} 0 & 1 \\ 1 & 1 \end{pmatrix} \begin{pmatrix} a \\ b \end{pmatrix} $$ Schreiben wir also $ M := \begin{pmatrix} 0 & 1 \\ 1 & 1 \end{pmatrix} $ und $\mathbf{u}_0 := \begin{pmatrix} 0 \\ 1 \end{pmatrix}$, dann ist der Vektor \begin{pmatrix}a \\ b \end{pmatrix} nach $n$ Schleifendurchläufen zu $M^n \mathbf{u}_0$ geworden. Anders gesagt: $$ M^n \mathbf{u}_0 = \begin{pmatrix} F_{n-1} & F_n \\ F_{n} & F_{n+1} \end{pmatrix} \ . $$ Wir können $F_n$ also berechnen, indem wir $M_n$ berechnen. Wir brauchen also einen Algorithmus für die Matrix-Exponentiation. Einen naiven Algorithmus sehen Sie hier:

function matrixExp (M, n):
    R := die Identitätsmatrix                                
    for i = 1 ... n: 
        R := M * R 
    return R

Übungsaufgabe Implementieren Sie matrixExp(M, n) in Python. Beschränken Sie sich hier auf $2 \times 2$-Matrizen $M$. Schreiben Sie als erstes eine Funktion matrixMult(A,B), die zwei $2 \times 2$-Matrizen multipliziert. Falls Sie es vergessen haben: $$ \begin{pmatrix} a & b \\ c & d \end{pmatrix} \cdot \begin{pmatrix} r & s \\ t & u \end{pmatrix} = \begin{pmatrix} ar + bt & bs + du \\ cr + dt & cs + du \end{pmatrix} \ . $$

Übungsaufgabe Messen Sie empirisch die Laufzeit Ihrer Funktion matrixExp(M,n) für eine feste Matrix $M$, die Sie selbst wählen dürfen, und wachsende Zahlen $n$. Wie verhält sich die Laufzeit in Abhängigkeit von $n$?

Schnelle Exponentiation

Mittels Rekursion können Sie Matrix-Exponentiation viel schneller implementieren. Die Idee ist einfach: wenn $n$ gerade ist, also $n = 2k$, dann können wir $$M \mapsto M^{k} \mapsto (M^{k})^2$$ rechnen. Der erste Schritt ist ein rekursiver Aufruf. Der zweite Schritt multipliziert einfach das Ergebnis mit sich selbst. Wenn $n$ ungerade ist, also $n = 2k + 1$, dann rechnen wir $$ M \mapsto M^{k} \mapsto (M^{k})^2 \mapsto M \cdot (M^{k})^2 \ , $$ also wie im geraden Fall, nur, dass wir noch eine Multiplikation mit $M$ am Ende ausführen. Sie sehen: der Wert von $n$ halbiert sich mit jedem rekursiven Aufruf, sie haben also höchstens $\log_2(n)$ rekursive Aufrufe.

Übungsaufgabe Schreiben Sie in Python eine Funktion fastMatrixExp(M,n), die diese Idee implementiert. Messen Sie die Laufzeit. Setzen Sie die Laufzeit zu $n$ in Bezug. Was geschieht mit der Laufzeit, wenn Sie $n$ verdoppeln?

Übungsaufgabe Wenn Sie fastMatrixExp(M,n) richtig implementiert haben, wird ein Aufruf nur $O(\log n)$ Operationen ausführen, da sich $n$ ja in jedem Schritt halbiert. Warum ist die gemessene Laufzeit dennoch nicht $O(\log n)$?

Um die Laufzeit von fastMatrixExp auch theoretisch zu verstehen, nehmen wir vereinfachend an, dass $n = 2^d$ gilt, $n$ also eine Zweierpotenz ist. Dies ist besonders bequem, weil in diesem Falle die "Zwischenergebnisse" von fastMatrixExp die folgende Form haben:

$$ M \mapsto M^2 \mapsto M^4 \mapsto M^8 \mapsto \cdots \mapsto M^{2^{d-1}} \mapsto M^{2^d} = M^n \ . $$ In jeden Schritt nehmen wir also die vorherige Matrix und quadrieren sie. Wie lange benötigen wir, um eine Matrix zu quadrieren? $$ \begin{pmatrix} u & v \\ w & x \end{pmatrix} \cdot \begin{pmatrix} u & v \\ w & x \end{pmatrix} = \begin{pmatrix} u^2 + vw & uv + vx \\ wu + xw & wv + x^2 \end{pmatrix} \ . $$

Wir brauchen also acht Multiplikationen und vier Additionen. Intuitiv gesehen ist Multiplikation schwieriger als Addition, daher beschränken wir uns darauf, die Komplexität der acht Multiplikationen verstehen zu wollen. Nehmen wir an, die Zahlen $u,v,w,x$ haben Bit-Länge höchstens $k$. Sei $\mult(k)$ die Komplexität der Multiplikation von $k$-Bit-Zahlen. Dann brauchen wir also $O(\mult(k))$ Schritte für die Quadrierung unserer Matrix.

Wenn $M := \begin{pmatrix} 0 & 1 \\ 1 & 1 \end{pmatrix}$, dann ist

$$ M^{2^{d-1}} = M^{n/2} = \begin{pmatrix} F_{n/2-1} & F_{n/2} \\ F_{n/2} & F_{n/2+1} \end{pmatrix} \ . $$

Die vier Einträge in $M^{n/2}$ haben also Bit-Länge $\Theta(\bitsize(F_{n/2})) = \Theta(n/2) = \Theta(n)$, und somit dauert die letzte Matrix-Quadrierung, $ M^{2^{d-1}} \mapsto M^{2^d}$ auch $\Theta(\mult(n))$ viele Schritte. Was ist mit den vorherigen Schritten, also $ M^{2^{i}} \mapsto M^{2^{i+1}}$? Der braucht $\Theta(\mult(2^{i}))$ viele Schritte. Es wird sich herausstellen, dass diese Schritte nicht so sehr ins Gewicht fallen. Formell gesprochen: wenn wir $M^n$ berechnen, dann entfällt ein konstant großer Anteil der Laufzeit auf den letzten Schritt: $M^{n/2} \mapsto M^n$.

Übungsaufgabe Messen Sie empirisch, wie lange es in Python braucht, zwei $n$-stellige Zahlen miteinander zu multiplizieren. Legen Sie eine Tabelle an! Was geschieht mit der Laufzeit, wenn Sie $n$ verdoppeln?

Übungsaufgabe Schreiben Sie eine Funktion fibFME, die $F_n$ wie oben beschrieben mithilfe der Fast Matrix Exponentation berechnet. Messen Sie die Laufzeit und untersuchen Sie, wie sich eine Verdopplung von $n$ auf die Laufzeit auswirkt.