Kapitel 2: Algebraische Algorithmen
(Effiziente Algorithmen, WS 2019) Gerhard Woeginger
WS 2019, RWTH
Algebraische Algorithmen
Multiplikation von ganzen Zahlen Multiplikation von Matrizen Berechnung von inversen Matrizen Multiplikation von Polynomen
Zum Aufw¨ armen:
Multiplikation von ganzen Zahlen
Multiplikation von ganzen Zahlen: Schulmethode
In der Grundschule haben wir gelernt,
wie man zwei positive ganze Zahlenx undy mit einander multipliziert:
1 1 9 5 8 3 3 ∗ 2 8 5 1
2 3 9 1 6 6 6
9 5 6 6 6 6 4
5 9 7 9 1 6 5
1 1 9 5 8 3 3
3 4 0 9 3 1 9 8 8 3
Diese Schulmethode verwendetΘ(n2)Operationen.
(Wir betrachten hier die Bit-Komplexit¨at von Algorithmen.) Frage
Geht das irgendwie besser? Schneller?
Multiplikation von ganzen Zahlen: Schulmethode
In der Grundschule haben wir gelernt,
wie man zwei positive ganze Zahlenx undy mit einander multipliziert:
1 1 9 5 8 3 3 ∗ 2 8 5 1
2 3 9 1 6 6 6
9 5 6 6 6 6 4
5 9 7 9 1 6 5
1 1 9 5 8 3 3
3 4 0 9 3 1 9 8 8 3
Diese Schulmethode verwendetΘ(n2)Operationen.
(Wir betrachten hier die Bit-Komplexit¨at von Algorithmen.) Frage
Geht das irgendwie besser? Schneller?
Multiplikation von ganzen Zahlen: Eine erste Idee
Wir teilen die Ziffernfolgen der beidenn-stelligen Zahlen
x=10n/2x1+x0 undy=10n/2y1+y0in zwei gleich lange Teile auf.
x: x1 x0
y: y1 y0
Divide and Conquer
Wir berechnen rekursiv die vier Produktex0y0,x0y1, x1y0, x1y1. Wir geben10nx1y1+ (x0y1+x1y0)10n/2+x0y0aus.
Ergo: T(n) =4T(n/2) + Θ(n)
Multiplikation von ganzen Zahlen: Die zweite Idee
Idee von Anatoly Alexeevitch Karatsuba (1962) Statt den vier Produktenx0y0,x0y1,x1y0, x1y1,
berechnen wir nur die drei Produkte x0y0, x1y1, und (x0+x1)(y0+y1)
Damit k¨onnen wir den gemischten Term in der Form x0y1+x1y0 = (x0+x1)(y0+y1)−x0y0−x1y1
schreiben, und das gew¨unschte Produkt wie folgt berechnen:
xy = 10nx1y1+ ((x0+x1)(y0+y1)−x0y0−x1y1)10n/2+x0y0
Multiplikation von ganzen Zahlen: Die zweite Idee
Idee von Anatoly Alexeevitch Karatsuba (1962) Statt den vier Produktenx0y0,x0y1,x1y0, x1y1,
berechnen wir nur die drei Produkte x0y0, x1y1, und (x0+x1)(y0+y1) Damit k¨onnen wir den gemischten Term in der Form
x0y1+x1y0 = (x0+x1)(y0+y1)−x0y0−x1y1
schreiben, und das gew¨unschte Produkt wie folgt berechnen:
xy = 10nx1y1+ ((x0+x1)(y0+y1)−x0y0−x1y1)10n/2+x0y0
Multiplikation von ganzen Zahlen: Resultat
Ergo:T(n) =3T(n/2) + Θ(n)
Satz (Anatoly Karatsuba, 1962)
Das Produkt von zwein-stelligen Zahlen
kann mitΘ(nlog23)vielen Bit-Operationen berechnet werden.
Anmerkung: log23≈1.585
Anmerkungen
Satz (David Harvey & Joris van der Hoeven, M¨arz 2019) Das Produkt von zwein-stelligen Zahlen
kann mitΘ(nlogn)vielen Bit-Operationen berechnet werden.
Spektrum der Wissenschaft (April 2019) https://www.spektrum.de/news/
die-schnellste-art-zu-multiplizieren/1638472
Matrix-Multiplikation
Matrix-Multiplikation: Schulmethode
Problem: Matrix-Multiplikation
Eingabe:Zwein×nMatrizenAund B Gesucht:Das ProduktC :=AB
Trivialer Algorithmus: Berechne EintragCij alsPn
k=1aikbkj
Satz
Das Produkt von zwein×nMatrizen
kann in kubischer ZeitO(n3)berechnet werden.
Anmerkung: Wir z¨ahlen die Anzahl der Multiplikationen.
Frage
Geht das irgendwie besser? Schneller?
Multiplikation von 2 × 2 Matrizen (1)
Mit acht Multiplikationen und vier Additionen:
a b c d
! w x y z
!
= aw+by ax+bz cw+dy cx+dz
!
Multiplikation von 2 × 2 Matrizen (2)
Mit sieben Multiplikationen und achtzehn Additionen:
m1 = (b−d)(y+z) m2 = (a+d)(w +z) m3 = (a−c)(w+x) m4 = (a+b)z m5 = a(x−z) m6 = d(y−w) m7 = (c+d)w
aw+by ax+bz cw+dy cx+dz
!
= m1+m2−m4+m6 m4+m5
m +m m −m +m −m
!
Multiplikation von 2 × 2 Matrizen (3)
Mit sechs Multiplikationen und 8.128 Additionen:
Nein, das geht nicht.
Bei sieben Multiplikationen ist bereits Schluss.
Shmuel Winograd hat 1971 gezeigt, dass man zur Multiplikation von 2×2Matrizen mindestens sieben Multiplikationen braucht (auch wenn man noch so viele Additionen und Subtraktionen zur Verf¨ugung hat).
Multiplikation von 2 × 2 Matrizen (3)
Mit sechs Multiplikationen und 8.128 Additionen:
Nein, das geht nicht.
Bei sieben Multiplikationen ist bereits Schluss.
Shmuel Winograd hat 1971 gezeigt, dass man zur Multiplikation von 2×2Matrizen mindestens sieben Multiplikationen braucht (auch wenn man noch so viele Additionen und Subtraktionen zur Verf¨ugung hat).
Schnelle Matrix-Multiplikation nach Strassen (1)
Idee von Volker Strassen (1969)
Konstruiere einen Divide-and-Conquer Algorithmus f¨ur die Multiplikation vonn×nMatrizen.
Verwende die sieben Multiplikationen und achtzehn Additionen (aus dem Multiplikationsschema f¨ur2×2Matrizen) im Conquer Schritt.
Schnelle Matrix-Multiplikation nach Strassen (2)
Im Divide Schritt unterteilen wir die beiden Matrizenn×nMatrizen Aund B jeweils in viern/2×n/2Matrizen:
A= A11 A12 A21 A22
!
B = B11 B12 B21 B22
!
Im Conquer Schritt berechnen wir durch zehn Matrix-Additionen und durch sieben rekursive Aufrufe die sieben Matrix-Produkte
M1, . . . ,M7f¨ur diesen/2×n/2MatrizenA11,A21, . . . ,B22. Wir f¨ugen die Ergebnisse durch acht Matrix-Additionen zur Ergebnismatrix zusammen:
M1+M2−M4+M6 M4+M5
M6+M7 M2−M3+M5−M7
!
Schnelle Matrix-Multiplikation nach Strassen (3)
Ergo:T(n) =7T(n/2) + Θ(n2)
Satz (Strassen, 1969)
Das Produkt von zwein×nMatrizen
kann mitO(nlog27)Integer-Multiplikationen berechnet werden.
Anmerkung: log27≈2.807
Anmerkungen
Der konstante Faktor (in der O-Notation versteckt) von Strassen’s Algorithmus ist gr¨osser als der konstante Faktor in derO(n3) Schulmethode.
In der Praxis: Schulmethode f¨ur kleinen, Strassen f¨ur grossen.
Der Cross-over Punkt liegt meistens um n=20.
In der Praxis: Wenn Matrizen d¨unn besetzt sind (= mit sehr vielen Nullen), dann gibt es schnellere Spezialalgorithmen aus der Numerik
Don Coppersmith und Shmuel Winograd (1990) haben den Exponenten von Strassen’s2.807auf2.376verbessert
Virginia Vassilevska Williams (2011) hat den Exponenten weiter auf 2.373verbessert
Als untere Schranke f¨ur die Multiplikation von zwein×nMatrizen kennen wir nur die SchrankeΩ(n2).
Anmerkungen
Der konstante Faktor (in der O-Notation versteckt) von Strassen’s Algorithmus ist gr¨osser als der konstante Faktor in derO(n3) Schulmethode.
In der Praxis: Schulmethode f¨ur kleinen, Strassen f¨ur grossen.
Der Cross-over Punkt liegt meistens um n=20.
In der Praxis: Wenn Matrizen d¨unn besetzt sind (= mit sehr vielen Nullen), dann gibt es schnellere Spezialalgorithmen aus der Numerik Don Coppersmith und Shmuel Winograd (1990) haben den
Exponenten von Strassen’s2.807auf2.376verbessert
Virginia Vassilevska Williams (2011) hat den Exponenten weiter auf 2.373verbessert
Als untere Schranke f¨ur die Multiplikation von zwein×nMatrizen kennen wir nur die SchrankeΩ(n2).
Verifikation von
Matrix-Multiplikationen
Verifikation von Matrix-Multiplikationen
Problem: Matrix-Multiplikation Verifikation Eingabe:Drein×nMatrizenA, B,C Frage:GiltAB =C?
Trivialer L¨osungsansatz: Berechne das ProduktAB
Vergleiche dien2 Eintr¨age mit den Eintr¨agen inC Frage
Geht das irgendwie besser? Schneller? Einfacher?
Verifikation von Matrix-Multiplikationen
Problem: Matrix-Multiplikation Verifikation Eingabe:Drein×nMatrizenA, B,C Frage:GiltAB =C?
Trivialer L¨osungsansatz:
Berechne das ProduktAB
Vergleiche dien2 Eintr¨age mit den Eintr¨agen inC Frage
Geht das irgendwie besser? Schneller? Einfacher?
Ein probabilistischer Ansatz
Idee von Rusins Martins Freivalds (1977) W¨ahle zuf¨allig einen Vektorx ∈ {0,1}n. Wenn ABx =Cx, dann return “AB=C”.
WennABx 6=Cx, dann return “AB 6=C”.
Anmerkungen:
Wenn Freivalds Ungleichheit behauptet, dann gilt wirklich AB 6=C. Wenn Freivalds Gleichheit behauptet, dann gilt aber nicht
notwendigerweiseAB =C.
Beispiel: Falls der zuf¨allig gew¨ahlte Vektorx=0ist, dann behauptet Freivalds Gleichheit v¨ollig unabh¨angig vonA, B,C
Ein probabilistischer Ansatz
Idee von Rusins Martins Freivalds (1977) W¨ahle zuf¨allig einen Vektorx ∈ {0,1}n. Wenn ABx =Cx, dann return “AB=C”.
WennABx 6=Cx, dann return “AB 6=C”.
Anmerkungen:
Wenn Freivalds Ungleichheit behauptet, dann gilt wirklich AB 6=C. Wenn Freivalds Gleichheit behauptet, dann gilt aber nicht
notwendigerweiseAB =C.
Beispiel: Falls der zuf¨allig gew¨ahlte Vektorx=0ist, dann behauptet Freivalds Gleichheit v¨ollig unabh¨angig vonA, B,C
Fehleranalyse (1)
Satz
WennD einen×nMatrix mitD 6=0ist undx ein zuf¨allig gew¨ahlter Vektor in {0,1}n,
dann ist die Wahrscheinlichkeit vonDx =0h¨ochstens1/2.
Beweisskizze:
W¨ahle Indizesk und`, sodassdk`6=0in MatrixD gilt.
Setzey :=Dx und betrachte die k-te Komponenteyk iny. Dann giltyk = dk1x1+dk2x2+· · ·+dknxn = dk`x`+R Zufallsexperiment: Wir setzen der Reihe nach die Komponenten xi
miti 6=`zuf¨allig und unabh¨angig von einander auf0/1, und erhalten so die entsprechenden RestsummeR.
Dann ist mindestens einer der beiden Werte R+dk` mit x` = +1 und R mitx`=0ungleich0.
Daher gilt yk 6=0mit Wahrscheinlichkeit mindestens1/2.
Fehleranalyse (2)
Satz
Wenn der Freivalds Algorithmus behauptet,
dassAB 6=C, dann stimmt das auf jeden Fall;
dassAB =C, dann stimmt das mit Wahrscheinlichkeit≥1/2.
Beweis:
Wende Satz von vorhergehender Seite auf MatrixD :=C−AB an
Wiederholt man den Freivalds Algorithmus100-mal (unabh¨angig) so sinkt die Fehlerwahrscheinlichkeit von1/2auf2−100≈10−30
Laufzeit
Satz
Der Freivalds Algorithmus kann inO(n2)Zeit implementiert werden.
Beweis: Das ProduktABx berechnet man alsA(Bx).
Ubung ¨
Ein viel einfacherer Ansatz (“Ingenieursmethode”) W¨ahle zuf¨allig Zeilenindexi und Spaltenindexj Berechne Eintrag[AB]ij im ProduktAB Wenn [AB]ij =Cij, dann return “AB =C”.
Wenn[AB]ij 6=Cij, dann return “AB 6=C”.
Lineare LaufzeitO(n)
Frage: Wie gross ist die Fehlerwahrscheinlichkeit?
(Extremfall: Nur ein einziger Eintrag inC ist falsch)
Frage: Wie oft muss man diesen einfachen Algorithmus wiederholen, damit die Fehlerwahrscheinlichkeit unter 1/2sinkt?
Berechnung von inversen Matrizen
Inverse Matrizen
Definition
Dieinverse MatrixA−1 zu einer quadratischenn×nMatrixA erf¨ullt die Gleichungen A−1A=AA−1=I.
Eine MatrixAmitdet(A) =0heisstsingul¨ar und besitzt keineinverse Matrix.
Eine MatrixAmitdet(A)6=0heisstregul¨arund besitzt eine eindeutigeinverse Matrix.
F¨ur jede regul¨are MatrixAgilt: det(A)·det(A−1) =1
Das Gauss-Jordan Eliminationsverfahren (Jiu Zhang Suanshu, 100 v.Chr.) berechnet die inverse Matrix inO(n3)Zeit.
Frage
Geht das irgendwie besser? Schneller?
Inverse Matrizen
Definition
Dieinverse MatrixA−1 zu einer quadratischenn×nMatrixA erf¨ullt die Gleichungen A−1A=AA−1=I.
Eine MatrixAmitdet(A) =0heisstsingul¨ar und besitzt keineinverse Matrix.
Eine MatrixAmitdet(A)6=0heisstregul¨arund besitzt eine eindeutigeinverse Matrix.
F¨ur jede regul¨are MatrixAgilt: det(A)·det(A−1) =1
Das Gauss-Jordan Eliminationsverfahren (Jiu Zhang Suanshu, 100 v.Chr.) berechnet die inverse Matrix inO(n3)Zeit.
Frage
Geht das irgendwie besser? Schneller?
Inversion → Multiplikation
F¨ur zwein×nMatrizenAundB betrachten wir die3n×3nMatrix M =
In A 0 0 In B 0 0 In
mit M−1=
In −A AB 0 In −B 0 0 In
.
Satz
Wenn die inverse Matrix einern×nMatrix inO(nα)Zeit berechnet werden kann (mit fixemα≥2),
dann kann auch das Produkt zweier beliebiger n×nMatrizen inO(nα)Zeit berechnet werden.
Inversion → Multiplikation
F¨ur zwein×nMatrizenAundB betrachten wir die3n×3nMatrix M =
In A 0 0 In B 0 0 In
mit M−1=
In −A AB 0 In −B 0 0 In
.
Satz
Wenn die inverse Matrix einern×nMatrix inO(nα)Zeit berechnet werden kann (mit fixemα≥2),
dann kann auch das Produkt zweier beliebiger n×nMatrizen inO(nα)Zeit berechnet werden.
Inversion → Multiplikation
F¨ur zwein×nMatrizenAundB betrachten wir die3n×3nMatrix M =
In A 0 0 In B 0 0 In
mit M−1=
In −A AB 0 In −B 0 0 In
.
Satz
Wenn die inverse Matrix einern×nMatrix inO(nα)Zeit berechnet werden kann (mit fixemα≥2),
dann kann auch das Produkt zweier beliebiger n×nMatrizen inO(nα)Zeit berechnet werden.
Multiplikation → Inversion
Unser n¨achstes Ziel ist es nun,
aus einemO(nα)Zeit Algorithmus f¨ur Matrix-Multiplikation einenO(nα)Zeit Algorithmus f¨ur Matrix-Inversion zu bauen.
Satz
Wenn das Produkt zweier beliebigern×nMatrizen inO(nα)Zeit berechnet werden kann (mit fixemα≥2),
dann kann auch die inverse Matrix von einern×nMatrix inO(nα)Zeit berechnet werden.
Gegeben sei also eine regul¨aren×nMatrixA.
Wir wollen die inversen×nMatrixA−1f¨urAfinden.
Multiplikation → Inversion: Erste Vereinfachung
Erste Vereinfachung
Wir k¨onnen/d¨urfen/werden annehmen,
dass die Dimension n=2q der Eingabematrix eine Zweierpotenzist.
A 0 0 Ik
!−1
= A−1 0 0 Ik
!
Multiplikation → Inversion: Zweite Vereinfachung
Zweite Vereinfachung
Wir k¨onnen/d¨urfen/werden annehmen,
dass die EingabematrixAsymmetrischundpositiv-definitist.
1. Berechne Hilfsmatrix B=ATA.
Diese HilfsmatrixB ist symmetrisch und positiv-definit.
2. BerechneB−1= (ATA)−1.
3. Dann giltA−1= (ATA)−1AT =B−1AT.
Schritt 1 kostet eine Matrix-Transposition und eine Matrix-Multiplikation.
Schritt 2 ist die Inversion einer symmetrischen + positiv-definiten Matrix.
Schritt 3 kostet eine Matrix-Multiplikation.
Multiplikation → Inversion: Algorithmus (1)
Wir verwenden Divide & Conquer und unterteilen dien×nMatrixA in vier n2×n2 Untermatrizen:
A = B CT C D
!
Dann gilt:
A−1 = B−1+B−1CTS−1CB−1 −B−1CTS−1
−S−1CB−1 S−1
!
wobei die MatrixS dasSchur Komplement vonB inC ist:
S = D−CB−1CT
Multiplikation → Inversion: Algorithmus (2)
1. Bestimme die Untermatrizen B, C, CT,D.
2. BerechneB−1 rekursiv.
3. Berechne das ProduktW =CB−1und transponiereWT =B−1CT. 4. Berechne das ProduktX =WCT und die Differenz S=D−X.
(Comment: Es giltX =CB−1CT undS=D−CB−1CT) 5. BerechneS−1rekursiv.
6. Berechne das Produkt Y =S−1W und transponiereYT =WTS−1. (Comment: Y =S−1CB−1undYT =B−1CTS−1)
7. Berechne das ProduktZ =WTY. (Comment: Z=B−1CTS−1CB−1)
Multiplikation → Inversion: Algorithmus (3)
Man pr¨uft leicht nach, dass die gesuchte Inverse wie folgt aussieht:
A−1 = B−1+B−1CTS−1CB−1 −B−1CTS−1
−S−1CB−1 S−1
!
= B−1+Z −YT
−Y S−1
!
In anderen Worten:
Aus den Hilfsmatrizen, die in den sieben Schritten berechnet werden, l¨asst sich leicht die inverse MatrixA−1 zusammensetzen.
Multiplikation → Inversion: Algorithmus (4)
In den sieben Schritten
wird zweimal eine n2×n2 Matrix invertiert;
werden vier Produkte von n2×n2 Matrizen berechnet;
werden verschiedene n2×n2 Matrizen transponiert, zu einander addiert, von einander subtrahiert.
F¨ur die Zeitkomplexit¨atT(n)des Algorithmus gilt daher T(n) ≤ 2T(n/2) +4·O(nα) +O(n2)
≤ 2T(n/2) +O(nα)
≤ O(nα)
Inverse Matrizen: Zusammenfassung
Satz
Multiplikation vonn×nMatrizen und Inversion vonn×n Matrizen sind zwei exakt gleich schwere Probleme.
Jeder Algorithmus f¨ur das eine Problem ¨ubersetzt sich in einen Algorithmus f¨ur das andere Problem mit der selben asymptotischen Worst Case Zeitkomplexit¨at.
Multiplikation von Polynomen
Polynom-Multiplikation: Schulmethode
In der Mittelschule haben wir gelernt,
wie man zwei PolynomeA(x)undB(x)mit einander multipliziert:
(x3+5x2−2x+1)·(3x3−x2+x+2)
= 3x6 +15x5 −6x4 +3x3
−x5 −5x4 +2x3 −x2 x4 +5x3 −2x2 +x
+2x3 +10x2 −4x +2
= 3x6 +14x5 −10x4 +12x3 +7x2 −3x +2 Die Schulmethode verwendetΘ(n2)Operationen,
um zwei Polynome vom Gradnmit einander zu multiplizieren.
Polynom-Multiplikation: Illustration
a
0b
0a
1b
0a
2b
0a
3b
0a
4b
0a
5b
0a
6b
0a
7b
0a
0b
1a
1b
1a
2b
1a
3b
1a
4b
1a
5b
1a
6b
1a
7b
1a
0b
2a
1b
2a
2b
2a
3b
2a
4b
2a
5b
2a
6b
2a
7b
2a
0b
3a
1b
3a
2b
3a
3b
3a
4b
3a
5b
3a
6b
3a
7b
3a
0b
4a
1b
4a
2b
4a
3b
4a
4b
4a
5b
4a
6b
4a
7b
4a
0b
5a
1b
5a
2b
5a
3b
5a
4b
5a
5b
5a
6b
5a
7b
5a
0b
6a
1b
6a
2b
6a
3b
6a
4b
6a
5b
6a
6b
6a
7b
6a
0b
7a
1b
7a
2b
7a
3b
7a
4b
7a
5b
7a
6b
7a
7b
7c
6= a
0b
6+ a
1b
5+ a
2b
4+ a
3b
3+ a
4b
2+ a
5b
1+ a
6b
0Polynom-Multiplikation: Genaue Problemstellung
A(x) = an−1xn−1+an−2xn−2+· · ·+a2x2+a1x+a0 B(x) = bn−1xn−1+bn−2xn−2+· · ·+b2x2+b1x+b0 C(x) = A(x)B(x)
= c2n−2x2n−2+c2n−3x2n−3+· · ·+c2x2+c1x+c0
Problem: Polynom-Multiplikation
Eingabe:Ganze Zahlena0, . . . ,an−1 undb0, . . . ,bn−1 Gesucht:Zahlenc0, . . . ,c2n−2 mitck =Pk
i=0aibk−i Anmerkung: Wir nehmen an, dassai=bi=0f¨uri ≥n
Darstellung von Polynomen (1)
Koeffizienten-Darstellung
Das PolynomA(x) =an−1xn−1+· · ·+a1x+a0wird durch die Folgea0, . . . ,an−1 der Koeffizienten spezifiziert.
Punkt-Wert-Darstellung
Das PolynomA(x) =an−1xn−1+· · ·+a1x+a0wird durch nPunkte(x0,y0), . . . ,(xn−1,yn−1)spezifiziert, wobeiyi=A(xi)f¨ur0≤i ≤n−1gilt.
Darstellung von Polynomen (2)
Zur Erinnerung
F¨urnpaarweise verschiedene Zahlen x0, . . . ,xn−1, und f¨urnPunkte(x0,y0), . . . ,(xn−1,yn−1)
existiert genau ein PolynomA(x) =an−1xn−1+· · ·+a1x+a0 mityi =A(xi)f¨ur0≤i ≤n−1.
Beweis f¨ur “es existiert h¨ochstens ein PolynomA(x)”:
Wenn ein Polynom vom Gradn−1mehr alsn−1Nullstellen hat, so ist es identisch gleich 0
Beweis f¨ur “es existiert mindestens ein PolynomA(x)”:
Lagrange Interpolation A(x) =
n−1
Xyk
Q
j6=k(x−xj) Q
j6=k(xk−xj)
Darstellung von Polynomen (3): ¨ Ubersetzung/Hin
Koeffizienten→Punkt-Wert
Gegeben: Koeffizienten-Darstellunga0, . . . ,an−1eines PolynomsA(x);
npaarweise verschiedene Zahlen x0, . . . ,xn−1
Gesucht: Punkt-Wert-Darstellung vonA(x)f¨ur St¨utzstellenx0, . . . ,xn−1
Jeder Wertyi=A(xi)kann inO(n)Zeit berechnet werden: 1 y= a[n-1];
2 for i= n-2 downto 0 do 3 y= x*y + a[i] 4
5 return y; Gesamtzeit:O(n2)
Darstellung von Polynomen (3): ¨ Ubersetzung/Hin
Koeffizienten→Punkt-Wert
Gegeben: Koeffizienten-Darstellunga0, . . . ,an−1eines PolynomsA(x);
npaarweise verschiedene Zahlen x0, . . . ,xn−1
Gesucht: Punkt-Wert-Darstellung vonA(x)f¨ur St¨utzstellenx0, . . . ,xn−1
Jeder Wertyi=A(xi)kann inO(n)Zeit berechnet werden:
1 y= a[n-1];
2 for i= n-2 downto 0 do 3 y= x*y + a[i]
4
5 return y;
Darstellung von Polynomen (4): ¨ Ubersetzung/R¨ uck
Punkt-Wert→Koeffizienten
Gegeben: Punkt-Wert-Darstellung(x0,y0), . . . ,(xn−1,yn−1)f¨urA(x);
Gesucht: Koeffizienten-Darstellunga0, . . . ,an−1vonA(x)
Die Lagrange’sche Interpolationsformel kann inO(n2)Zeit ausgewertet werden
Zusammenfassend:
Koeffizienten-Darstellung und Punkt-Wert-Darstellung k¨onnen in quadratischer ZeitO(n2)in einander ¨ubergef¨uhrt werden
Zur¨ uck zur Polynom-Multiplikation
Unser Hauptziel: Polynom-Multiplikation in sub-quadratischer Zeit Problem: Polynom-Multiplikation (in Koeffizienten-Darstellung) Eingabe:Ganze Zahlena0, . . . ,an−1 undb0, . . . ,bn−1
Gesucht:Zahlenc0, . . . ,c2n−2 mitck =Pk
i=0aibk−i
Nebenproblem: Polynom-Multiplikation in sub-quadratischer Zeit Problem: Polynom-Multiplikation (in Punkt-Wert-Darstellung) Eingabe:Punkt-Wert-Darstellung(x0,y0), . . . ,(xn−1,yn−1)f¨urA(x)
Punkt-Wert-Darstellung(x0,y00), . . . ,(xn−1,yn−10 )f¨urB(x) Gesucht:Punkt-Wert-Darstellung(x0,y000), . . . ,(xn−1,yn−100 )f¨urA(x)B(x) Einfach in linearer ZeitO(n)zu l¨osen
Arbeitsplan
Schritt 1:
Ubersetze die beiden Polynome¨ A(x)undB(x)
von Koeffizienten-Darstellung nach Punkt-Wert-Darstellung Schritt 2:
MultipliziereA(x)undB(x)in Punkt-Wert-Darstellung Schritt 3:
Ubersetze das in Schritt 2 berechnete Produkt¨ A(x)B(x) von Punkt-Wert-Darstellung nach Koeffizienten-Darstellung
Intermezzo:
Rechnen mit komplexen Zahlen
Komplexe Zahlen
Eine komplexe Zahlz kann dargestellt werden:
Algebraisch durch Zerlegung in Realteil und Imagin¨arteil:z=a+i b Polar durch Radius r und Winkelφ: z =rcosφ+i rsinφ
Exponentiell: z=r·eiφ
Rechenoperationen mitz =r ·eiφ undz0=r0·eiφ0 Multiplikation: (r·eiφ) (r0·eiφ0) = rr0·ei(φ+φ0) Potenzierung: zn = rn·ei nφ
Einheitswurzeln (1)
Dien-ten Einheitswurzelnω0n, ω1n, . . . , ωn−1n
sind dien(komplexen) L¨osungen der Gleichung ωn=1.
1 1
i
i
!80D!88
!81
!82
!83
!48
!85
!86
!87
1 2πi/n k 2kπi/n
Einheitswurzeln (2)
ωn=e2πi/nist dien-te Haupt-Einheitswurzel
Alle Einheitswurzeln sind Potenzen der Haupt-Einheitswurzel ωn
Es gelten die ¨ublichen Rechenregelnωknωn` =ωk+`n
Die n-ten Einheitswurzeln bilden eine multiplikative Gruppe, die zur zyklischen GruppeZn (Restklassengruppe modulon) isomorph ist Wenn nZweierpotenz, dann erzeugt jedes Elementωkn mit ungeradem Exponentenk die gesamte Gruppe.
Einheitswurzeln (2)
ωn=e2πi/nist dien-te Haupt-Einheitswurzel
Alle Einheitswurzeln sind Potenzen der Haupt-Einheitswurzel ωn
Es gelten die ¨ublichen Rechenregelnωknωn` =ωk+`n
Die n-ten Einheitswurzeln bilden eine multiplikative Gruppe, die zur zyklischen GruppeZn (Restklassengruppe modulon) isomorph ist Wenn nZweierpotenz, dann erzeugt jedes Elementωkn mit ungeradem Exponentenk die gesamte Gruppe.
Einheitswurzeln (3)
Halbierungslemma
F¨ur eine Zweierpotenz nfallen die Quadrate dern-ten
Einheitswurzeln mit den(n/2)-ten Einheitswurzeln zusammen.
Summenlemma
Jeden-te Einheitswurzelωnk 6=1erf¨ullt die Gleichung (ωkn)n−1+ (ωnk)n−2+· · ·+ (ωnk)2+ (ωnk) +1 = 0
Arbeitsplan: Schritt 1
Schritt 1: Zielsetzung
Zur Erinnerung:
Schritt 1:
Ubersetze die beiden Polynome¨ A(x)undB(x)
von Koeffizienten-Darstellung nach Punkt-Wert-Darstellung
Von jetzt an nehmen wir an, dassn=2q eine Zweierpotenz ist.
Die Punkt-Wert-Darstellung vonA(x)und B(x)werden wir mit den n-ten Einheitswurzelnhx0, . . . ,xn−1i=hω0n, . . . , ωn−1n i als St¨utzstellen bestimmen.
Wir beschreiben das Verfahren nur f¨urA(x).
Das PolynomB(x)wird analog behandelt.
Wir verwenden Divide & Conquer.
Divide & Conquer (1)
Divide & Conquer Ansatz:
Aeven(x) = a0+a2x+a4x2+a6x3+· · ·+an−2xn/2−1 Aodd(x) = a1+a3x+a5x2+a7x3+· · ·+an−1xn/2−1
Koeffizienten-Darstellung:
Aeven(x) : ha0, a2, a4, a6, . . . , an−2i Aodd(x) : ha1, a3, a5, a7, . . . , an−1i
A(x) = Aeven(x2) +x·Aodd(x2)
Divide & Conquer (2)
Wir beobachten, dass die Zahlen(ωn0)2,(ωn1)2, . . . ,(ωnn−1)2 exakt mit den(n/2)-ten Einheitswurzeln zusammen fallen
Wir bestimmen rekursiv die Punkt-Wert-Darstellung vonAeven(x) f¨ur dien/2St¨utzstellen(ω0n)2,(ω1n)2, . . . ,(ωn−1n )2
Wir bestimmen rekursiv die Punkt-Wert-Darstellung vonAodd(x) f¨ur dien/2St¨utzstellen(ω0n)2,(ω1n)2, . . . ,(ωn−1n )2
Wir berechnen aus diesen beiden Punkt-Wert-Darstellungen die Punkt-Wert-Darstellung vonA(x)f¨ur dienSt¨utzstellen ωn0, . . . , ωnn−1 mit Hilfe vonA(x) = Aeven(x2) +x·Aodd(x2) Ergo: T(n) =2T(n/2) + Θ(n)
Schritt 1: Zusammenfassung
Ergo: T(n) =2T(n/2) + Θ(n) Ergo: T(n)∈O(nlogn)
Satz
Schritt 1 mit denn-ten Einheitswurzelnhx0, . . . ,xn−1i=hω0n, . . . , ωn−1n i als St¨utzstellen kann inO(nlogn)Zeit durchgef¨uhrt werden.
A moment of reflection
Was haben wir in Schritt 1 eigentlich getan? (1)
1 1 1 1 1 · · · 1
1 ω ω2 ω3 ω4 · · · ωn−1
1 ω2 ω4 ω6 ω8 · · · ω2(n−1)
1 ω3 ω6 ω9 ω12 · · · ω3(n−1)
1 ω4 ω8 ω12 ω16 · · · ω4(n−1)
... ... ... ... ... ... ...
... ..
. ..
. ..
. ..
. ..
. ..
. 1 ωn−1 ω2(n−1) ω3(n−1) ω4(n−1) · · · ω(n−1)(n−1)
a0
a1
a2
a3
a4
... ... an−1
=
y0
y1
y2
y3
y4
... ... yn−1
Wir haben den Koeffizientenvektoramit einer MatrixV multipliziert,
Was haben wir in Schritt 1 eigentlich getan? (2)
Wir haben den Koeffizientenvektora mit einer MatrixV multipliziert, und als Resultat den Vektory erhalten.
Die MatrixV =V(ω)mit Eintrag ωrs in Zeiler und Spaltes ist eine so-genannteVandermondeMatrix.
Der Vektory =VawirdDiskrete Fourier Transformierte (DFT)des Vektorsa bez¨uglich der Einheitswurzelωgenannt.
Analog f¨ur Einheitswurzelnα:=ωk mit ungeradem Exponentenk: Die Vandermonde MatrixV =V(α)hat den Eintragαrs in Zeile r und Spaltes stehen.
Der Vektory =VawirdDiskrete Fourier Transformierte (DFT)des Vektorsa bez¨uglich der Einheitswurzelαgenannt.
(Analog f¨ur alle anderen Einheitswurzeln. Aber das verwenden wir nicht.)
Was haben wir in Schritt 1 eigentlich getan? (2)
Wir haben den Koeffizientenvektora mit einer MatrixV multipliziert, und als Resultat den Vektory erhalten.
Die MatrixV =V(ω)mit Eintrag ωrs in Zeiler und Spaltes ist eine so-genannteVandermondeMatrix.
Der Vektory =VawirdDiskrete Fourier Transformierte (DFT)des Vektorsa bez¨uglich der Einheitswurzelωgenannt.
Analog f¨ur Einheitswurzelnα:=ωk mit ungeradem Exponentenk: Die Vandermonde MatrixV =V(α)hat den Eintragαrs in Zeile r und Spaltes stehen.
Der Vektory =VawirdDiskrete Fourier Transformierte (DFT)des Vektorsa bez¨uglich der Einheitswurzelαgenannt.
Was haben wir in Schritt 1 eigentlich getan? (3)
Satz
F¨ur jede Zweierpotenznund
f¨ur jede Einheitswurzelα:=ωk mit ungeradem Exponentenk, kann die DFT f¨ur einenn-dimensionalen Vektorabez¨uglich der Einheitswurzelα inO(nlogn)Zeit berechnet werden.
Arbeitsplan: Schritt 3
Schritt 3: Zielsetzung
Zur Erinnerung:
Schritt 3:
Ubersetze das in Schritt 2 berechnete Produkt¨ A(x)B(x) von Punkt-Wert-Darstellung nach Koeffizienten-Darstellung
Unser Startpunkt ist die Punkt-Wert-Darstellung
eines Polynoms C(x)annSt¨utzstellen(x0,y0), . . . ,(xn−1,yn−1), wobeixk =ωnk f¨ur0≤k ≤n−1
Wir suchen die Koeffizienten-Darstellung c0, . . . ,cn−1von C(x)
Illustration (1)
1 1 1 1 1 · · · 1
1 ω ω2 ω3 ω4 · · · ωn−1
1 ω2 ω4 ω6 ω8 · · · ω2(n−1)
1 ω3 ω6 ω9 ω12 · · · ω3(n−1)
1 ω4 ω8 ω12 ω16 · · · ω4(n−1)
... ... ... ... ... ... ...
... ... ... ... ... ... ...
1 ωn−1 ω2(n−1) ω3(n−1) ω4(n−1) · · · ω(n−1)(n−1)
c0
c1
c2
c3
c4
... ... cn−1
=
y0
y1
y2
y3
y4
... ... yn−1
Lineares GleichungssystemVc=y (und wir suchen den Vektorc)
Illustration (2)
c0
c1
c2
c3
c4
... ... cn−1
=
1 1 1 1 1 · · · 1
1 ω ω2 ω3 ω4 · · · ωn−1
1 ω2 ω4 ω6 ω8 · · · ω2(n−1)
1 ω3 ω6 ω9 ω12 · · · ω3(n−1)
1 ω4 ω8 ω12 ω16 · · · ω4(n−1)
... ... ... ... ... ... ...
... ... ... ... ... ... ...
1 ωn−1 ω2(n−1) ω3(n−1) ω4(n−1) · · · ω(n−1)(n−1)
−1
y0
y1
y2
y3
y4
... ... yn−1
Lineares Gleichungssystemc=V−1y (und wir suchen den Vektorc)
Das lineare Gleichungssystem (1)
Es seiω=e2πi/n dien-te Haupt-Einheitswurzel.
Es seiV = (vj,k)dien×nMatrix mitvj,k =ωjk. Es seieny0, . . . ,yn−1 die vorgegebenen Funktionswerte
des PolynomsC(x)an den St¨utzstellenω0, . . . , ωn−1
Unser Ziel: L¨ose das GleichungssystemVc =y nach dem Vektorc auf
Das lineare Gleichungssystem (2)
Satz
Die inverse MatrixV−1= (wj,k)ist gegeben durchwj,k =ω−jk/n.
Beweis: Wir verifizieren, dassV−1V =I gilt.
Der Eintrag in derr-ten Zeile unds-ten Spalte vonV−1V ist [V−1V]r,s =
n−1
X
k=0
(ω−rk/n) (ωks)
= 1 n
n−1
X
k=0
(ωs−r)k
F¨urr =s istωs−r =1, und daherPn−1
k=0(ωs−r)k =n F¨urr 6=s istωs−r 6=1, und daherPn−1
k=0(ωs−r)k =0
Das lineare Gleichungssystem (3)
Altes Ziel:
L¨ose das GleichungssystemVc =y nach dem Vektorc auf
Neues Ziel:
Berechne den Vektorc=V−1y, wobei[V−1]j,k =ω−jk/ngilt.
Neues Ziel, besser formuliert:
Berechne die Diskrete Fourier Transformierte (DFT)d des Vektorsy bez¨uglich der Einheitswurzelω−1.
Der gesuchte Vektor ist dannc =d/n.
Ergo: Schritt 3 kann inO(nlogn)Zeit erledigt werden
Das lineare Gleichungssystem (3)
Altes Ziel:
L¨ose das GleichungssystemVc =y nach dem Vektorc auf
Neues Ziel:
Berechne den Vektorc=V−1y, wobei[V−1]j,k =ω−jk/ngilt.
Neues Ziel, besser formuliert:
Berechne die Diskrete Fourier Transformierte (DFT)d des Vektorsy bez¨uglich der Einheitswurzelω−1.
Der gesuchte Vektor ist dannc =d/n.
Ergo: Schritt 3 kann inO(nlogn)Zeit erledigt werden
Polynom-Multiplikation:
Zusammenfassung
Arbeitsplan
Schritt 1:
Ubersetze die beiden Polynome¨ A(x)undB(x)
von Koeffizienten-Darstellung nach Punkt-Wert-Darstellung Laufzeit:O(nlogn)
Schritt 2:
MultipliziereA(x)undB(x)in Punkt-Wert-Darstellung Laufzeit:O(n)
Schritt 3:
Ubersetze das in Schritt 2 berechnete Produkt¨ A(x)B(x) von Punkt-Wert-Darstellung nach Koeffizienten-Darstellung Laufzeit:O(nlogn)
Hauptresultat
Satz
Das Produkt zweier Polynome vom Gradn−1in
Koeffizienten-Darstellung kann inO(nlogn)Zeit berechnet werden.
Anmerkungen:
Da das ErgebnispolynomC(x)den Grad2n−2hat, m¨ussen wir ganz am Anfang die Koeffizientenvektoren a0, . . . ,an−1und
b0, . . . ,bn−1durch Hinzuf¨ugen von Komponenten mit Wert0auf die Dimension 2n−2bringen.
Durch Hinzuf¨ugen von noch mehr Komponenten mit Wert0machen wir die Dimension zu einer Zweierpotenz
Anmerkungen
Die Diskrete Fourier Transform wurde von James William Cooley und John Wilder Tukey 1965 algorithmisch untersucht.
(“An algorithm for the machine calculation of complex Fourier series”, Mathematics of Computation 11, pp 297–301)
Der Numeriker Carl David Tolm´e Runge verwendete bereits 1903 die Reduktion der n-dimensionalen DFT auf zwei(n/2)-dimensionale DFTs
Der rekursive Algorithmus wurde bereits von Carl Friedrich Gauss im Jahre 1805 benutzt, um die Flugbahnen der Asteroiden Pallas und Juno zu interpolieren. (“Theoria interpolationis methodo nova tractata”, verfasst in Neu-Latein)
Die FFT (Fast Fourier Transformation) wurde als einer der Top 10 Algorithms of the 20th Centurygelistet (Algorithms with the greatest influence on the development and practice of science and engineering in the 20th century)
Die Liste der Top 10 Algorithmen
Metropolis Algorithm for Monte Carlo Simplex Method for Linear Programming Krylov Subspace Iteration Methods
The Decompositional Approach to Matrix Computations The Fortran Optimizing Compiler
QR Algorithm for Computing Eigenvalues Quicksort Algorithm for Sorting
Fast Fourier Transform Integer Relation Detection Fast Multipole Method