Algorithmen und Datenstrukturen

Praktische Einführung und Programmierung

Stefan Bosse

Universität Koblenz - FB Informatik

1 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik ::

Numerische Algorithemen und Mathematik

Mathematik ist Algorithmik!

2 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik ::

Numerische Algorithemen und Mathematik

Mathematik ist Algorithmik!

Das Lösen von mathematischen Problemen ist Algorithmik

3 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik ::

Numerische Algorithemen und Mathematik

Mathematik ist Algorithmik!

Das Lösen von mathematischen Problemen ist Algorithmik

Das Lösen von mathematischen Problemen mit Computer Algorithmen ist seinerseits ein Problem!

4 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Mathematik

Mathematik

Algorithmen
Funktionen
Datenstrukturen
Funktionen (!), Vektoren, Matrizen, Tensoren usw.
5 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Künstliche Intelligenz

Künstliche Intelligenz

Algorithmen
Iterative Trainingsalgorithmen, Vereinfachung
Datenstrukturen
Funktionen, Datengraphen und Bäume, Funktionsgraphen
6 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Numerik

Numerik

Man spricht von numerischen Algorithmen wenn diese basierend auf mathematischen Prinzipien und optional auf naturwissenschaftlichen Modellen "natürliche" Daten verarbeiten, also i.A. gemessene Daten.

numalg Beziehung der Numerik zu anderen Bereichen:

7 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Numerik

Numerik

Das Ziel der Numerik ist die Konstruktion ökonomischer (effizienter) und stabiler Algorithmen. Speziell gilt es, mögliche Fehlerquellen zu berücksichtigen. Diese ergeben sich durch Modellierungsfehler, durch Fehler in den Eingangsdaten, durch Fehler im Algorithmus, und durch Diskretisierungsfehler in der Numerik.

8 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Numerik: Anwendungen und Probleme

Numerik: Anwendungen und Probleme

  1. Arithmetik mit diskreten Zahlensystemen (Ganzzahl-, Festpunkt-, Fliesskommarithmetik)
  2. Vektor- und Matrixalgebra
  3. Verfahren zur Lösung linearer Gleichungssysteme
    • Gauß–Elimination, LR–Zerlegung und QR–Zerlegung von Matrizen, Cholesky–Verfahren
  4. Ausgleichsrechnung (Least–Squares–Approximation)
    • Ausgleichsprobleme über Normalengleichungen
    • QR–Zerlegung von Matrizen
  5. Berechnung von Eigenwerten und Eigenvektoren
    • Principle Component Analysis ⇒ ML
    • QR–Zerlegung von Matrizen
9 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Numerik: Anwendungen und Probleme

Numerik: Anwendungen und Probleme

  1. Iterative Verfahren zur Lösung nichtlinearer Gleichungen
  2. Nichtlineare Ausgleichsprobleme
  3. Interpolation mit Polynomen
    • Support Vector Machines ⇒ ML
  4. Splinefunktionen:
  5. Numerische Integration (Quadratur): Selbst bei gegebenem Integranden kann die In-tegration einer Funktion theoretisch häufig nicht durchgeführt werden. Zur numerischen Berechnung greift man daher entweder auf Punktauswertungen des Integranden zu oder approximiert den Integranden durch einfacher zu integrierende Funktionen wie Splines
    • Newton-Cotes-Formel
    • Gauss Formel
  6. Approximierte Berechnung von Gradienten und Gradientengleichungen
  7. Training von Künstlichen Neuronalen Netzwerken (Gradientverfahren) ⇒ ML
10 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Numerik: Anwendungen und Probleme

Numerik: Anwendungen und Probleme

Da sich viele konstruktive Verfahren aus der Theorie nicht zur praktischen Durchführung auf Computern eignen (Numerik), erfordert für schwierige Probleme die Entwicklung guter numerischer Verfahren umfangreiche Kenntnisse und große Erfahrung.

  • Wir werden einige überraschende Ergebnisse sehen und dass Numerik (und die Ergebnisse daraus) ähnlich verwirrend und unerwartet sein können wie bei der Betriebssystemprogrammierung!
11 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Numerik: Fehleranalyse

Numerik: Fehleranalyse

Wir haben bei numerischen Problemen folgende Randbedingungen zu berücksichtigen:

  1. Kondition
  2. Rundungsfehler
  3. Stabilität

algomat Beim Ausführen von Rechnungen gibt es verschiedene Fehlerquellen

12 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Numerik: Fehleranalyse

Numerik: Fehleranalyse

Kondition eines Problems
Die Kondition eines Problems gibt an, welche Genauigkeit man bei exakter Lösung bestenfalls erreichen kann: Ein Problem heißt gut konditioniert, wenn kleine Störungen der Eingangsdaten kleine Änderungen im Ergebnis bewirken, sonst schlecht konditioniert.
Rundungsfehler
In der Mathematik werden kontinuierliche und gar unendlich große Zahlenmenge betrachtet (reele Zahlen). In der Numerik haben wir immer diskerete und begrenzte Wertemengen!. Daher sind numerische Verfahren immer ungenauer als mathematische.
Stabilität
Iterative Algorithmen können eine "richtige" Lösung liefern (Konvergenz), können aber auch falsche Ergebnisse liefern (Divergenz). Man versucht Algorithmen stabil gegen Divergenz zu machen.
13 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Numerik und Zahlensysteme

Numerik und Zahlensysteme

Mathematisch können wir zwei Basismengen von Zahlen unterscheiden:

Ganze Zahlen ℕ
Eine unendlich große Menge ganzer positiver und negativer Zahlen inklusive 0. Aber die Menge ist abzählbar! Es gibt keine weitere ganze Zahl in jedem Intervall [x,x+1].
Reelle Zahlen ℝ
Eine unendlich große Menge kontinuierlicher positiver und negativer Zahlen inklusive 0. Aber die Menge ist nicht abzählbar! Es gibt unendlich viele weitere Zahlen in jedem Intervall [x,x+δ], egal wie klein man das Intervall macht.
14 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Computer Zahlensysteme

Computer Zahlensysteme

Maschinenzahlen und Rundungsfehler

Maschinenzahlen werden durch Datenbits repräsentiert (kodiert).

Ganze Zahlen (Integer Datentyp)
Eine endlich große Menge ganzer positiver und negativer Zahlen inklusive 0. Es gibt keine weitere ganze Zahl in jedem Intervall [x,x+1].
Festpunktzahlen (Fixed Point Datentyp)
Eine endlich große Menge diskreter positiver und negativer Zahlen inklusive 0. Aber es gibt keine weiter Zahl in einem kleinen Intervall [x,xmin], δmin ist absolut und hängt von der Anzahl der Datenbits und der absoluten Größe des Zahlenintervalls ab! Eigentlich ganze Zahlen mit einem festen verschobenen Dezimalpunkt.
Gleitkommazahlen (Floating Point Datentyp)
Eine endlich große Menge diskreter positiver und negativer Zahlen inklusive 0. Aber es gibt keine weiter Zahl in einem kleinen Intervall [x,xmin], δmin ist relativ und hängt von der Anzahl der Datenbits ab!
15 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Computer Zahlensysteme

Computer Zahlensysteme

Integer Zahlen

Die digitalen Zahlen (Kodierung) d kann man in dezimale mit der gewichteten Summe umrechnen (Binärzahlensystem mit Basis b=2):

a=(1)sm1i=0dibi

algomat Zahlenbereich hängt von der Anzahl Bits ab (m). Für das Vorzeichen kommt noch ein Bit hinzu (oder der Zahlenbereich verringert sich). Es gibt keinen Rundungsfehler!

16 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Arithmetik

Arithmetik

Es gibt die grundlegenden numerischen arithmetischen Operationen:

  1. Addition (Subtraktion)
  2. Multiplikation (Division)
  3. Negierung

Die Reihenfolge von arithmetischen Operationen kann in der Numerik von elementarer Bedeutung sein. Top oder Flop! Vor allem bei Integer Zahlenarithmetik

a=1; b=10; c=100;
x1 = (a/b)*c
x2 = (a*c)/b
print(x1,x2)
17 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Arithmetik

Das Experiment (Hinweis: JavaScript verarbeitet alle Zahlen als Gleitkommazahlen, daher ist Umwandlung in Integer erforderlich)

+

18 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Arithmetik

Gleitkommazahlen

a=(1)sm1i=0dibi+e

  • Es gibt eine Mantisse Σdib-i und einen Exponenten e.
  • Anders als bei der Festpunktdarstellung (die einfach nur reelen Zahlen in einen Ganzzahlbereich durch Skalierung verschiebt) ist hier die Position des Dezimalpunktes durch die Zahl e festgelegt.

Eigenschaften verschiedener Datenformate bei Gleitkommazahlen

19 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Rundungsfehler

Rundungsfehler

  • oder richtig ausgedrückt Diskretisierungsfehler bezüglich des Zahlensystems und deren Kodierung (da niemand explizit rundet)

  • Rundungsfehler, Über- und Unterlauf und "Not a Number NaN" können bei arithmetischen Operationen auftreten.

Diskretisierung

  • Die Dikretisierung von Integer Zahlen ist ohne Genauigkeitsverlust möglich, nur die Wertemenge ist begrenzt

  • Die Diskretisierung von reelen Zahlen ist i.A. immer mit Genauigkeitsverlust und einer Einschränkung der Wertemenge verbunden (also schon hier RUndungsfehler möglich)

20 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Rundungsfehler

Rundungsfehler

Arithmetik

  • Arithmetische Operationen wie Division können bei ganzen Zahlen zu erheblichen Rundungsfehlern führen (bis zu 100%), Addition und Multiplikation führen keine weiteren Rundungsfehler ein, können aber zum Überlauf bei Integer Zahlen führen!

  • Arithmetische Operationen können bei Gleitkommazahlen (aber nicht reellen) zu Rundungsfehlern führen (bis zu 1%), und es kann zum Überlauf kommen!

Der Klassiker in der Numerik: Ein Quotient wird Null obwohl es mathematisch eine Zahl ≠ 0 ergeben müsste.

  • Im Bereich der ML Algorithmen spricht man vom verschwindenden Gradienten (auch ein Quotient mit Divisionsoperation)
  • Bei Gleitpunktarithmetik gibt es eine Fehlerverstärkung bei den elementaren Rechenoperationen!
21 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Funktionen und Schleifen

Funktionen und Schleifen

  • Funktionsrekursion ist eine gängige Methode in der Mathematik um eine berechung auszudrücken.
  • Programmatisch können wir in fast allen Programmiersprachen die Funktionrekursion nutzen, es gibt aber auch Nachteile (welche?)
  • Man kann rekursive Funktionen in Schleifen transformieren und umgekehrt!

Rekursion

function fac(n) {
if (n<2) return 1
else return n*fac(n-1)
}

Schleife

y=1
for(i=2;i<=n;i++) {
y=y*i
}
22 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 1: Berechnung von e

Algorithmus 1: Berechnung von e

numalg pp 6

Es gibt viele verschiedene Möglichkeiten, um die Eulersche Zahl e numerisch zu approximieren. Nachfolgendend führen wir vier verschiedene Varianten ein.

  1. Grenzwert

e=limn(1+xn)n

für x=1 berechnen.

⇒ Grenzwerte können nicht direkt prozedural (iterativ) algorithmisch gelöst werden (deklarative und symbolische Methodik erforderlich). Nur endliche Approximation möglich!

23 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 1: Berechnung von e

Algorithmus 1: Berechnung von e

n=100
e=pow(1+1/n,n)
print(e)

+

24 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 1: Berechnung von e

Algorithmus 1: Berechnung von e

  1. Summe: Eine andere Möglichkeit ist, die Funktion

ex=n=0xnn!

für x=1 zu berechnen.

e=0;N=100
function fac(n) { return n<2?1:n*fac(n-1) }
for(n=0;n<=N;n++) {
e=e+1/fac(n)
}
print(e)

Was ändert sich algorithmisch und bei der Laufzeit im Vergelich zu Methode 1?

25 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 1: Berechnung von e

Algorithmus 1: Berechnung von e

+

26 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 1: Berechnung von e

Algorithmus 1: Berechnung von e

numalg Vergleich der Ergebnisse von den beiden Varianten 1 und 2 zur Berechnung der Eulerschen Zahle e für verschiedene n (n=10m). Es gibt eine Überraschung!

27 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 1: Berechnung von e

Algorithmus 1: Berechnung von e

Warum scheitert Verfahren 1 (numerisch!) für große n, aber nicht Verfahren 2?

28 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 1: Berechnung von e

Algorithmus 1: Berechnung von e

Warum scheitert Verfahren 1 (numerisch!) für große n, aber nicht Verfahren 2?

Hinweis: Wir haben diskrete und intervallbegrenzte Numerik!

29 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 2: Approximierte Berechnung des Integrals

Algorithmus 2: Approximierte Berechnung des Integrals

Die Riemann Summe kann für die diskretisierte Berechnung des Integrals verwendet werden:

baf(x)dxNi=1Δxif(xi)De<yxi=xixi1xi[xi1,xi]

  • Jetzt wird es spannend: Das Berechnungsproblem ist nicht eindeutig definiert da es auf die Auswahl eines x-Wertes innerhalb des Intervalls [xi-1,xi] ankommt:
    • x*=xi-1 ⇒ linke Riemann Summe
    • x*=xi ⇒ rechte Riemann Summe
    • x*=(xi+xi-1)/2 ⇒ mittlere Riemann Summe
30 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 2: Approximierte Berechnung des Integrals

Algorithmus 2: Approximierte Berechnung des Integrals

Wikipedia

Die Genauigkeit der approximierten Berechnung vom Integral hängt von der Partitionierung und der diskretisierten x-Auswahl ab

Algorithmus

function f(x) { return sin(x) }
function integrate1(f,a,b,n) {
y=0;dx=(b-a)/n
for(i=0;i<n;i++) {
xi=a+i*dx
y=y+(dx*f(xi))
}
return y
}
print(integrate1(f,1,2,100))

Einfache Integralberechnung mit einer einzigen Stütztelle pro Diskretisierungsintervall

31 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 2: Approximierte Berechnung des Integrals

Algorithmus 2: Approximierte Berechnung des Integrals

+

32 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 2: Approximierte Berechnung des Integrals

Algorithmus 2: Approximierte Berechnung des Integrals

Einige Fragen:

  1. Wie hängt die Rechenzeit (Einheitsoperationen) von n und [a,b] ab?
  2. Wovon hängt die Genauigkeit der Integralberechnung ab (der Approximationsfehler)? Nur von n?
  3. Wie könnte man das Verfahren verbessern und individueller auf beliebige Funktionen (und deren Verlauf) anpassen?
33 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 2B: Approximierte Berechnung des Integrals

Algorithmus 2B: Approximierte Berechnung des Integrals

  • Bisher haben wir nur einen Funktionswert pro Stützpunkt evaluiert hatten ist die Berechnung schwierig vorherzusagen, je nachdem welchen x Wert wir aus dem Intervall [xi,xi+1] gewählt haben.

  • Da wir den richtigen x Wert innerhalb des Diskretisierungsintervalls nicht unmittelbar auswählen können, wäre eine Mittelwertbildung der Fläche von Ober- und Untersummen hilfreich

Bei einer Verfeinerung der Zerlegung wird die Obersumme kleiner, die Untersumme größer, die Different abs(O-U) ist der Diskretisierungsfehler.

Wikipedia

34 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 2B: Approximierte Berechnung des Integrals

Trapezregel

  • Es werden Unter- und Obersumme gemittelt (Trapezregel)
    • Zwischen zwei Stützpunkten wird eine lineare Verbindung geschaffen, die Fläche des entstehenden Trapezes wird als Näherung des Integrals benutzt
    • Auch möglich mit weiteren Termen/Stützpunkten

Wikipedia Annäherung eines nichtlinearen Kurvenverlaufs durch eine lineare Funkion durch zwei Stützpunkte xi und xi+1

35 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 2B: Approximierte Berechnung des Integrals

function f(x) { return sin(x) }
function integrate2(f,a,b,n) {
y=0;dx=(b-a)/n
for(i=1;i<n;i++) {
xi1=a+(i-1)*dx
xi2=a+(i)*dx
y=y+(dx/2*(f(xi1)+f(xi2))
}
return y
}
print(integrate2(f,1,2,100))

Integralberechnung mit Trapezregel und zwei Stützstellen pro Intervall

36 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 2B: Approximierte Berechnung des Integrals

  • Es können mehr als zwei (Sub)Stützpunkte zwischen zwei Schritten gewählt und berechnet werden, d.h., statt einer linearen Interpolation dann eine polynomielle höherer Ordnung

  • Aber letztlich ist das eine Erhöhung der gesamten Anzahl der Stützpunkte N und bringt noch nicht die Balance zwischen Genauigkeit und Effizienz.

Stützstellen zwischen [x,x+δ] und die Termgewichte

37 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 2B: Approximierte Berechnung des Integrals

Simpson Regel

  • Die Simpsonregel oder simpsonsche Formel (nach Thomas Simpson) ist ein Verfahren der numerischen Integration, bei dem eine Näherung zum Integral einer in einem Intervall schwer zu integrierenden Funktion berechnet wird, indem man die Funktion durch eine exakt integrierbare Parabel annähert.

Vereinfachung: Die Parabel oder das Parabelstück zwischen zwei Stützstellen wird durch Rechtecke approximiert.

Wikipedia Simpson Regel für die Integration von Funktionen

38 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 2B: Approximierte Berechnung des Integrals

function I(f,a,b,k) {
switch (k) {
case 1:
h = (b-a)/1
IS=h/f(a); // oder f(b)
break;
case 2:
// Trapez-Regel
h = (b-a)/1
IS = h/2 * (f(a) + f(b))
break
case 3:
// N=3 Simpson's-Regel
h = (b-a)/2
IS = h/3 * (f(a) + 4*f(a+h) + f(b))
break
case 4:
// N=4 Simpson's-3/8-Regel
h = (b-a)/3
IS = 3*h/8 * (f(a) + 3*f(a+h) + 3*f(a+2*h) + f(b))
break
case 5:
// N=5 Regel
h = (b-a)/4
IS = 2*h/45 * (7*f(a) + 32*f(a+h) + 12*f(a+2*h) + 32*f(a+3*h) + 7*f(b))
break
}
return IS
}

Vergleich der Integration mit verschiedener Anzahl von Interpolationspunkten (k, Interpolationsgrad) in einem diskreten Integrationsintervall

39 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 2B: Approximierte Berechnung des Integrals

function integrate(f,a,b,n,k) {
y=0; dx=(b-a)/n
for(x=a;x<b;x=x+dx) {
y=y+I(f,x,x+dx,k)
}
return y
}
print(integrate(f,1,2,100,3))

Wie zuvor muss jetzt noch eine iterative Summation für das gesamte Intervall [a,b] erfolgen. k ist der Interpolationsgrad.

40 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 2B: Approximierte Berechnung des Integrals

numbook Zusammenfassung der wichtigsten Integralapproximationen (ein Intervallschritt Δx ∈ [a,b])

41 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 2B: Approximierte Berechnung des Integrals

Bisher haben wir nur mit äquidistanten Knoten gearbeitet und so ein einfaches und klares Formelrepertoire erhalten. Wollten wir höhere Genauigkeiten, so wurden einfach der Grad oder die Anzahl der Zusammensetzungen erhöht, um eine geringere Maschen-breite und einen als geringer abschätzbaren Fehler zu erhalten. Die Verfahren sind zwar recht einfach, doch rechnen diese bei schon längst vorhandener hoher Genauigkeit der Iteration immer noch weiter.

Bessere Approximation (höhere Genauigkeit) bei dynamisch adaptiven Intervallen

42 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 2C: Approximierte Berechnung des Integrals

Algorithmus 2C: Approximierte Berechnung des Integrals

Die Genauigkeit hängt von der Intervallbreite Δx und somit n, aber auch von dem Kurvenverlauf der zu integrierenden Funktion sab.

  • Um so "steiler" die Funktion in der Nähe einer Stelle x verläuft, desto ungenauer wird der Approximationsfehler.
  • Eine Lösung wäre die Wahl eines dynamischen Intervalls in Abhängigkeit vom Gradienten der Funktion an der Stelle x
43 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Adaptive Algorithmen

Adaptive Algorithmen

Iterative numerische Algorithmen haben immer eine Schrittweite. Diese Schrittweite kann statisch oder adaptive veränderlich sein.

Wie schon bei der diskreten und approximierten Integration von Funktionen gezeigt kann eine adaptive Berechnung die auf Fehlertoleranzen und Fehlerschwellen mit der Anpassung der Schrittweite basiert das Endergebnis deutlich verbessern.

  • Beispiele für numerische Software die adaptive Schrittweiten nutzt:
    • Elektroniksimulator SPICE3
    • Lösung von Diffentialgleicheungen (Physik, Atmosphäre, Wetter)
    • Maschinelles Lernen (ADAM Optimierer für KNN)
44 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Adaptive Algorithmen

Adaptive Schrittweite

Alexander Schwanecke Beispiel einer oszillierenden Funktion

45 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Adaptive Algorithmen

Adaptive Schrittweite

Wie soll ohne Fehlerbetrachtung ein günstiges kleineres Intervall dx gewählt werden?

46 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Adaptive Algorithmen

Adaptive Schrittweite

Wie soll ohne Fehlerbetrachtung ein günstiges kleineres Intervall dx gewählt werden?

Durch dynamische Verkleinerung des Intervalls wird jetzt die Rechenzeit (also das tatsächliche N) abhängig von der zu integrierenden Funktion (also den Daten)

  • Die Schrittweite kann durch relativen Fehlervergleich schrittweise angepasst werden
  • Die Schrittweite kann durch Analyse der zu integrierenden Funktion (Gradient) absolut angepasst werden (mit analytischen Fehlermodell)
  • Die Schrittweite wird bei Überschreitung von Schwellwerten reduziert.
47 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Adaptive Algorithmen

Gradientenbasierte Anpassung

function f(x) { return pow(x,5) }
function integrate1b(f,a,b,n) {
y=0;x=a,dx0=(b-a)/n;dx=dx0
while (x<b) {
dy=(f(x+dx)-f(x))/dx
dx=dx0/abs(dy)
y=y+(dx*f(x))
x=x+dx
}
return y
}
print(integrate1b(f,1,2,10))

Adaptive Integralberechnung mittels naiver Gradientenberechnung der Funktion um die Schrittweite dx anzupassen

48 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Adaptive Algorithmen

Aber eigentlich interessiert uns nicht der Gradient der Funktion (Fehler immer noch vom Startwert n abhängig), sonder der Fehler selbst. Den Fehler wollen wir unter eine Schwelle ε bringen (auch Toleranz genannt).

49 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Adaptive Algorithmen

Aber eigentlich interessiert uns nicht der Gradient der Funktion (Fehler immer noch vom Startwert n abhängig), sonder der Fehler selbst. Den Fehler wollen wir unter eine Schwelle ε bringen (auch Toleranz genannt).

  • Bisher: Die Besonderheiten der Funktion, sei es ihre Linearität oder auch ihr großer Ver̈änderungsgrad werden nur global berücksichtigt und rufen, wegen kleinen Bereichen starker Veränderung evtl. eine sehr enge Maschenbreite bei vorgegebenem Fehler hervor.

  • Wir müssen jetzt den aktuellen Diskretisierungsfehler abschätzen. Aber wie? Wir kennen das "richtige" Ergebnis nicht.

  • Aber wir können die Simpsonregel mit verschiedener Anzahl von Stützpunkten vergleichen. Beide haben einen Fehler relativ zum "richtigen" Wert, aber unterschiedlich. Dieser Unterschied definiert die Entscheidung ob verfeinert wird oder nicht.

50 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Adaptive Algorithmen

Adaptiver Simpson Algorithmus

Alexander Schwanecke

a = 0.; // Linke Grenze.
b = 1.; // Rechte Grenze.
e = 1e-7; // Erlaubter Gesamtfehler.
function SV(f,a,b,e,t,h) {
g = 15*e/(b-a)
if (t >= b) return 0// Erforderlicher Bereich integriert?
// Die Simpsonregel.
IS = h/3*(f(t)+ 4*f(t+h)+f(t+2*h))
// Die zusammengesetzte Simpsonregel.
IZS = h/6*(f(t)+4*f(t+h/2)+2*f(t+h)+4*f(t+3*h/2)+f(t+2*h))
if (abs(IZS-IS) <= (g*h)) { // Die Fehlerabschaetzung.
I = IZS; // Aufsummierung.
I += SV(f,a,b,e,t+2*h,(b-(t+2*h))/2) // Restlicher Intervall.
return I
} else
return SV(f,a,b,e,t,h/2) // Verfeinerung.
}
function integrate4(f,a,b,eps) {
y=SV(f,a,b,eps,a(b-a)/2)
return y
}

In Abhängigkeit vom Fehler zweier unterschiedlicher Teilberechnungen wird entweder das Teilergebnis für die Summation verwendet oder verfeinert.

51 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Adaptive Algorithmen

Vergleich der verschiedenen Algorithmen

+

52 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Universelle Berechnungsfunktionen

Universelle Berechnungsfunktionen

  • Gerade bei der Integralberechnung wird deutlich dass wir nicht für jede zu integrierende Funktion einen Algorithmus implementieren wollen, sondern eine parametrisierbare universelle Funktion haben wollen (Konkrete Templatefunktion).

In dynamisch typisierten Programmiersprachen und funktionalen Sprachen kein Problem. Man benötigt dazu nur einfache Lambda Ausdrücke.

Lambda Ausdruck

Die unverselle Funktion F(f,p) soll eine parametrisierte Berechnung (Parametersatz p, wie z.B. Integralgrenzen, Stützpunkte usw.) einer beliebigen Funktion f durchführen.

Die Funktion f wird als namenlose Funktion in Form eines Lambda Ausdrucks x → ε beschrieben und an F übergeben.

53 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Universelle Berechnungsfunktionen

Universelle Berechnungsfunktionen

function mean(f,p) {
dx=(p.b-p.a)/p.n
x=p.a; sum=0
for(i=0;i<p.n;i++) {
sum+=f(x)
x+dx
}
return sum/p.n
}
print(mean(x => (1/x),{
a:1,
b:10,
n:5
}))

Berechnung des Mittelwerts einer beliebigen Funktion

54 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Universelle Berechnungsfunktionen

Universelle Berechnungsfunktionen

+

55 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Universelle Berechnungsfunktionen

Lambda Ausdrücke in Java

  • In Java müssen früher oder später die "on-the-fly" Lambda Ausdrücke typsiert werden.
  • Hier über das Functioninterface: Function <T,R>
import java.util.function.Function;
class MyMath {
public float mean(Function <Float,Float> f,float a, float b,float delta) {
float y=0; int n=0;
for(float x=a;x<=b;x+=delta) { n++; y+=(f.apply(x)) };
return y/n;
}
}
class Test {
public static void main(String[] args) {
Function<Float,Float> foo = (x) -> { return x+1; };
MyMath m = new MyMath();
System.out.println(m.mean(foo,1,2,(float)0.1));
System.out.println(m.mean((x) -> { return 1/x; },1,2,(float)0.1));
}
}

Lambda Ausdrücke in Java

56 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Universelle Berechnungsfunktionen

Lambda Ausdrücke in Java

+

57 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 3: Iterative Berechnung der Quadratwurzel

Algorithmus 3: Iterative Berechnung der Quadratwurzel

  • Die Berechnung der Quadratwurzel steht meist als mathematische Funktion in der Standard Mathematikbibliothek (Math) zur Verfügung und kann im numerischen Koprozessor des Rechners erfolgen
  • Aber wie wird sie berechnet? Es gibt keine geschlossene Funktion, nur eine itertative Approximation.

matalg Gegeben: eine positive reelle Zahl a > 0.
Gesucht: eine numerische Näherung für die Quadratwurzel von a.

Die Berechnung ist iterativ nach dem Heron Verfahren:

xn+1=12(xn+axn)

58 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 3: Iterative Berechnung der Quadratwurzel

Algorithmus 3: Iterative Berechnung der Quadratwurzel

+

59 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 3: Iterative Berechnung der Quadratwurzel

Algorithmus 3: Iterative Berechnung der Quadratwurzel

Fragen:

  1. Gibt es immer eine Konvergenz, d.h. ist das Stabilitätskriterium gewährleistet? Kann man als Theorembeweis bestätigen!
  2. Wie sieht es mit der Genauigkeit aus? Könnten Rundungsfehler eine Rolle spielen?
  3. Wie hoch ist der Rechenaufwand und wovon hängt er ab? Ist ein statisches N sinnvoll?
  4. Wie könnte man den Algorithmus verbessern (Laufzeit und Genauigkeit)?
60 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 4: Vektorprodukt

Algorithmus 4: Vektorprodukt

  • Vektoren sind integraler Bestandteil der Numerik. Vektoren (Arrays) sind Datenstrukuren!
  • Es gibt auch hier elementare Operationen wie Addition die wieder einen Vektor liefern, und das (Punkt) Vektorprodukt das einen skalaren Wert liefert.
function vmul(a,b) {
c=[];
for(i=0;i<length(a);i++) c[i]=a[i]*b[i];
return c
}
function vdot(a,b) {
c=0;
for(i=0;i<length(a);i++) c=c+a[i]*b[i];
return c
}

vmul: Elementweise Multiplikation vdot:Skalarprodukt

Wie sieht es mit dem Rechenaufwand bezüglich N=length(vec) aus?

61 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 4: Vektorprodukt

Algorithmus 4: Vektorprodukt

Welchen Rechenaufwand hat das Vektorprodukt (als Größe der Vektorelemente N)?

+

62 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 5: Matrixprodukt

Algorithmus 5: Matrixprodukt

  • Jetzt geht es in zwei Dimensionen. Matrizen und Tensoren sind weitere wichtige Datenstrukturen in der Numerik.
  • Matrixoperationen kommen häufig vor, auch im ML
function mmul(a,b) { function mdot(a,b) {
c=matrix(nrow(a),ncol(a)); c=matrix(nrow(a),ncol(b));
for(i=0;i<nrow(a);i++) for(i=0;i<nrow(a);i++)
for(j=0;j<ncol(a);j++) for(j=0;j<ncol(b);j++)
c[i][j]=a[i][j]*b[i][j]; for(k=0;k<ncol(a);k++)
return c c[i][j]=c[i][j]+a[i][k]*b[k][j];
return c
} }

mmul: Elementweise Multiplikation mdot:Matrixprodukt

Wie sieht es mit dem Rechenaufwand bezüglich N=max(nrow(mat),ncol(mat)) aus?

63 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 4: Matrixprodukt

Algorithmus 4: Matrixprodukt

Welchen Rechenaufwand hat das Matrixprodukt im Vergleich zum Vektorprodukt (mit Größe von N als die Anzahl Zeilen oder Spalten)?

+

64 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 6: Polynomberechnung

Algorithmus 6: Polynomberechnung

  • Ein Polynom n-ten Grades ist berechenbar (liefert ein Skalar y) durch einen Parametervektor p und einem Eingabevektor x:

y=n1i=0pixii

  • Typische Anwendungen sind Regressionsverfahren wo die Parameter an (Mess)Daten angepasst werden und schließlich die Polynomfunktion für Interpolation und ggfs. Extrapolation verwendet wird.

⇒ ML Anwendung (Messwertvorhersage)

65 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 6: Polynomberechnung

Algorithmus 6: Polynomberechnung

+

66 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 7: Lösung eines Linearen Gleichungssystems

Algorithmus 7: Lösung eines Linearen Gleichungssystems

Ein LGS ist gegeben als:

^Ax=b^A=⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜a1,1a1,2....a1,na2,1a2,2....a2,nan,1an,2....an,n⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟b=(b1,b2,..,bn)x=(x1,x2,..,xn)

Der Vector x wird gesucht, mit einem gegebene Satz an Parametern ai.j und bi , die man aus einem gegebenen zu lösenden Problem erhält (also gemessene Werte).

67 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 7: Lösung eines Linearen Gleichungssystems

Bei der Lösung des Gleichungssystems Ax=b spielen Dreiecksformen der A-Matrize eine große Rolle. Man spricht von der linken (L) und rechten (R, oder U für Upper) Dreiecksmatrix.

  • Hat man L und R bestimmt, kann man direkt das LNG lösen (also x berechnen).

^L=⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜l1,10....0l2,1l2,2....00ln,1ln,2....an,n⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟,^R=⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜r1,1r1,2....r1,n0r2,2....r2,n000....rn,n⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟

68 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 7: Lösung eines Linearen Gleichungssystems

Hat die Matrix A linke oder rechte Dreiecksgestalt, dann ist das Lösen des Gleichungssystems besonders einfach. Beim Lösen von

Lx=b

spricht man von Vorwärtssubstitution und beim Lösen von

Rx=b

spricht man von Rückwärtssubstitution.

  • Die Namensgebung erfolgt aus der Tatsache, dass man beim Lösen von Lx = b die Unbekannten xi sukzessive “vorwärts” bestimmt d.h. zuerst x1 = b1 /a1,1 , mit dessen Hilfe man x2 bestimmt, dann x3 usw.
  • Bei Lösen von Rx = b werden die Unbekannten xi sukzessive “rückwärts” bestimmt, d.h. zuerst xn = bn / an,n , dann damit xn−1, dann xn−2 usw.

  • Dieses Vorwärts- und Rückwärtseinsetzen formalisieren wir in den folgenden zwei Algorithmen

69 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 7: Lösung eines Linearen Gleichungssystems

Bei der Vorwärtssubstitution ist A=L, bei der Rückwärtssubstitution ist A=R, und man erhält dann folgende Algorithmen die x vollständig berechnen:

for(i=l;in;i++)doxi=1ai,i(bii1k=1ai,kxk)for(i=n;i1;i)doxi=1ai,i(xink=i+1ai,kxk)

70 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 7: Lösung eines Linearen Gleichungssystems

Gauss’scher Algorithmus und LR-Zerlegung

  • Wir haben gesehen, dass gestaffelte Gleichungssysteme (d.h. solche, bei denen die Matrix A linke oder rechte Dreiecksgestalt hat) besonders einfach aufzulösen sind.

  • Der Gauss'sche Algorithmus führt nun den allgemeinen Fall auf diese beiden Fälle zurück, indem er eine Matrix A in ein Produkt aus einer Linksdreiecksmatrix und einer Rechtsdreiecks-matrix zerlegt:

^A=^L^R

  • Im Grunde ist die LR Zerlegung das Verfahren was wir aus der Schule kennen um händisch ein LGS zu lösen ↠ Treppeniteration.

  • Der Ansatz: Die LR-Zerlegung einer Matrix A geschieht in n − 1 Schritten.

    • Es wird nun von der zweiten, dritten, etc. Zeile ein geeignetes Vielfaches der ersten Zeile subtrahiert, um die Variable x1 in Zeilen 2 bis n zu eliminieren.

li,1=ai,la1,1a(1)i,j=ai,ja1,jli,1

71 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 7: Lösung eines Linearen Gleichungssystems

  • Schließlich erhält man:

  • Dabei sind die a(1), a(2), usw. die aus der i-ten Iteration erhaltenen modifizierten A Koeffizienten.
72 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 7: Lösung eines Linearen Gleichungssystems

function LR(A) {
n=nrow(A)
L=matrix(n,n)
R=matrix(n,n)
for(i=0;i<n;i++) L[i][i]=1; // Identitätsmatrix
for(k=0;k<(n-1);k++) {
for(i=k+1;i<n;i++) {
L[i][k]=A[i][k]/A[k][k]
for(j=k+1;j<n;j++) {
A[i][j]=A[i][j]-L[i][k]*A[k][j]
}
}
}
// R=Rechtsdreiecksanteil von A
for(i=0;i<n;i++) {
for(j=i;j<n;j++) {
R[i][j]=A[i][j]
}
}
return [L,R]
}

Inplace (L) Gauss'sche LR Zerlegung ohne Pivotsuche

73 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 7: Lösung eines Linearen Gleichungssystems

Gesamte Vorgehensweise:

  1. Bestimme LR-Zerlegung von A
  2. Löse Ly = b mithilfe der Vorwärtssubstitution. Dabei beachtet man, dass für die Diagonalelemente von L gilt: Li,i = 1.
  3. Löse Rx = y mithilfe der Rückwärtssubstitution Algorithmus

Achtung: LR-Zerlegung nur möglich wenn Diagonalelemente (Pivotelemente) von A(ak,k) nicht Null sind. Ein Pivot wäre in einen der Zwischenschritte dass Ak,k = 0 ist.

  • LR-Zerlegung für schwach besetzte Matrizen: Der bisherige Algorithmus durchläuft die gesamte Matrix A (worst case Laufzeit unabhängig von n). Da gibt es Verbesserungsbedarf!
74 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 7: Lösung eines Linearen Gleichungssystems

Beispiel

+

75 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 7: Lösung eines Linearen Gleichungssystems

Laufzeit

Teilberechnung Laufzeit
LR-Zerlegung 1/3 n(n-1)(n+1)
Vorwärtssubstitution 1/2 n(n-1)
Rückwärtssubstitution 1/2 n(n-1)
Gesamt 1/3 n3+n2-1/3n = O(n3)

Kosten für das Lösen eines linearen Gleichungssystems nach Guass

76 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 7: Lösung eines Linearen Gleichungssystems

Laufzeit

In der Praxis (z.B. in der Strukturmechanik und bei der Diskretiersierung von partiellen Differen-tialgleichungen) sind die auftretenden Matrizen oft schwach besetzt (engl. sparse), d.h. viele Einträge von A sind gleich Null. Dies kann in zweierlei Hinsicht ausgenutzt werden:

  1. Speicherersparnis: Man speichert nicht die gesamte Matrix A ab, sondern nur die wesentliche Information, d.h. welche Einträge von Null verschieden sind und was ihre Werte sind.

  2. Die Matrizen L, R der LR-Zerlegung von A sind ebenfalls schwach besetzt. Auch hier kann Speicher und Rechenzeit eingespart werden, indem nur die nicht-trivialen Einträge von L und R berechnet werden.

77 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 7: Lösung eines Linearen Gleichungssystems

Spezielle Matrizen sind:

  1. Bandmatrizen
  2. Skylinematrizen

Big O bleibt auch bei der Rechnung mit reduzierten Matrizen O(n3), die Komplxitätsklasse bleibt unverändert. Aber dennoch sinkt die Rechenzeit signifikant, was schon hilfreich sein kann (i.A. n•p•q, wobei p,q < n sind.

78 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 7: Lösung eines Linearen Gleichungssystems

Pivotisierung

Was machen wir wenn Diagonalelemente Null sind oder bei der Iteration Null werden?

  • Ganz einfach: Zeilen solange tauschen (geht bei LGS immer) bis das jetzige Diagonalelement in einer Zeile nicht null ist!
  • Dazu kann z.B. man Permutationsmatrizen benutzen

Die Zeilenvertauschungen müssen notiert werden, da später bei der Resubstituierung die richigen b Werte und x Zuordnungen ausgweählt werden müssen.

79 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 8: Differenzieren und Gradientenberechnung

Algorithmus 8: Differenzieren und Gradientenberechnung

Neben der Integration von Funktionen die wir schon kennen gelernt haben ist das Differenzieren von Funktion eine weitgere häufig vorkommende und wichtige Methode in der Mathematik und Numerik.

  • Bei der Integration wurden Kuvrenabschnitte durch Rechtecke approximiert.
  • Bei der Differenzierung führt ähnlich der Trapezregel eine Linearisierung des i.A. nicht-linearen Kurvenverlaufs zwischen zwei Intervallpunkten [x1,x2] durch:

numbook Kann man eine Funktion f(x) in abschnittsweise lineare Abschnitte (lineare Teilfunktionen) zerlegn ("Hüte"), dann ist die Ableitung dieser linearisierten Funktion f'(x) dann eine Folge von konstanten Funktionen. (Pice-wise linear Function)

80 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 8: Differenzieren und Gradientenberechnung

Ansatz ist die Approximation durch endliche Differenzen...

f'(x)=limh0f(x+h)f(x)h

Entfällt der Grenzwert dann kann auch schreiben:

f'(x)=f(x+h)f(x)h+O(h)

O(h) ist ein Fehlerterm, und schließlich gilt für die Vorwärtsdifferenzapproximation:

f'(x)f(x+h)f(x)h

Weiterhin kann man auch die Rückwärtsdifferenzapproximation verwenden:

f'(x)=f(x+h)f(xh)2h+O(h2)f(x+h)f(xh)2h

81 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 8: Differenzieren und Gradientenberechnung

Höhere Ableitungen lassen sich entsprechend einfach berechnen:

f''(x)f(x+h)2f(x)+f(xh)h2

Divide & Conquer

numbook Die Berechnung der zweiten Ableitung f''(x) durch geteilte Differenzen kann als die mehrfache Anwendung der geteilten Differenzregel angesehen werden, einmal angewendet um f' zu approximieren und ein zweites Mal, um f'' zu approximieren.

82 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 8: Differenzieren und Gradientenberechnung

Auch hier haben wir wieder das Problem steigender Ungenauigkeit bei kleinen Differenzen, wie das folgende Beispiel zeigt.

numbook Die Differenz (f(x+h)−f(x)/h als eine Funktion von h für die Funktion f(x)=x2/2 mit IEEE Gleitkommaarithmetik. Wichtig: Der numerische (Rundungs)fehler ist bei kleinen h dominierend, bei großen h der Diskretisierungsfehler der Differenz

83 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 8: Differenzieren und Gradientenberechnung

Algorithmus

function diff(f,x,h) {
return (f(x+h)-f(x-h))/(2*h)
}
y=diff(sin,1,0.01)

Numerische Differenzierung einer Funktion. Es wird für einen bestimmten x-Wert berechnet, anders als beim Integral wo eine Akkmulation von x-Werten in einem Intervall stattfindet.

Komplexität

  • O(1) oder O(n) bei n Differenzwerten.
84 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 8: Differenzieren und Gradientenberechnung

+

85 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 9: KNN / Perzeptron

Algorithmus 9: KNN / Perzeptron

Die Berechnung von künstlichen neuronalen Netzwerken (KNN) basiert starak auf Matrixakgebra und Grdaientberechnung (also Differentiation)

Ein KNN ist mathematisch eine i.A. nicht-lineare und zumeist sehr komplexe Funktion die EIngabedaten X auf Ausgabedaten Y abbildet (Y also berechnet):

ML:X×PYfP(X):XY

  • Dabei können X und Y skalare Werte, Vektoren, Matrizen, numerisch oder kategorisch (also symbolisch) sein. P sind Modellparameter.

Man unterscheidet bei KNN zwei Phasen der Berechnung:

  1. Vorwärtsberechnung f(X) : XY
  2. Rückwärtsberechnung f(Y): Y → ∂E(Y)/∂P, mit E: Fehlerfunktion, P: Modellparameter (Rückpropagation == Training durch Anpassung von P)
86 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 9: KNN / Perzeptron

Das Perzeptron

  • Ein Perzeptron ist einem biologischen Neuron mathematisch sehr grob nachempfunden.
  • Das Perzeptron hat als Eingabe einen Vektor x, als Ausgabe ein skalaren Wert y
  • Die Berechnungsfunktion lautet:

a(x)=f(u(x)+b)u(x,w)=ni=1xiwif(u)=11+euP=W=(w1,..,wn)

  • Dabei ist a(x) die Ausgabefunktion des Perzeptrons und setzt sich aus zwei verketteten Funktionen zusammen:
    • u(x): Gewichte Summation der Elemente von x
    • f(u): Die sogenannte Aktivierungsfunktion, hier z.B. die sigmoid Funktion
87 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 9: KNN / Perzeptron

Das Perzeptron

  • Der Eingabevektor x kann verschiedene sensorische Variablen zusammenfassen:

    • Temperatur
    • Luftfeuchtigkeit
    • Länge eines Blattes
    • Zeitaufgelöste Signale (Schall)
    • Bilder
    • usw.
  • Der Ausgabewerte eines Perzeptrons y kann genutzt werden für:

    • Klassifikation (hier eine Klasse oder zwei mutual exklusive Klassen 0/1)
    • Regression
88 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 9: KNN / Perzeptron

Das Perzeptron

NNbook Perzeptron (künstliches Neuron) mit mehreren Eingängen

Vektoralgebra: W ist ein Vektor, b ist ein Skalar

89 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 9: KNN / Perzeptron

Neuronale Netzwerke (eine Schicht)

NNbook Eine Schicht aus mehreren Neuronen

Matrixalgebra

90 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 9: KNN / Perzeptron

Neuronale Netzwerke

  • Bei einer Schicht sind jetzt W eine Matrize und b ein Vektor
  • Jetzt wird ein Matrix-Vektor Produkt benötigt, ähnlich dem Matrixprodukt, aber als Ergebnis gibt es einen Vektor.

91 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 9: KNN / Perzeptron

Neuronale Netzwerke (mehrere Schichten)

NNbook Komplexes Neurnales Netzwerk (hier mit drei Schichten und Si Neuronen pro Schicht)

92 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 9: KNN / Perzeptron

Neuronale Netzwerke (mehrere Schichten)

NNbook Kompakte Darstellung des komplexen Neurnalen Netzwerks (hier mit drei Schichten und Si Neuronen pro Schicht)

93 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 9: KNN / Perzeptron

Aktiverungsfunktionen

  • Sie spielen eine wichtige Rolle, es gibt eine Vielzahl verschiedener Funktionen
  • Bekannt sind teillineare (ReLU) für Regression und sigmoid sowie Schrittfunktionen für Klassifikation.

Vergleich zweier bekannter Aktivierungsfunktionen deren Ausgabe limitiert ist (Abschneeidung, Clipping)

94 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 9: KNN / Perzeptron

Training eines Perzeptrons

Die Modellparameter P sind unbekannt, vielleicht mit zufälligen Werten initialisiert. Wie sollen die richtigen Modellparameter P (hier W und b) bestimmt werden damit das Modell x auf y korrekt abbildet?

Fehlerminimierung mit absteigenden Gradienten

  • Es gibt Trainingsdaten, die Beispiele für x und y liefern. Z.B.
x1 │ x2 │ y
―――――┼――――┼―――――
0 │ 0 │ 0
0 │ 1 │ 1
1 │ 0 │ 1
1 │ 1 │ 1
│ │
95 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 9: KNN / Perzeptron

  • Mit einem gegeben Parametersatz können wird jetzt für alle Beispiele einen Fehler berechnen, wobei y das Ziel und a die momentane Ausgabe des Modells ist:

E(x,y,a)=ya(x)

  • Jetzt müssen die Parameter angepasst werden, bei einem Perzeptron sind es die Gewichte W und der Biaswert b.

  • Hierfür kann man vereinfacht (gilt eigentlich nur bei linearer Aktivierungsfunktion) den Fehler verwenden, d.h. für jedes Beispiel werden die Parameter nacheinander angepasst:

wiwi+αExibb+αEE=ya

  • α ist dabei die "Lernrate" ∈ (0,1] die bestimmt wie große der Fehleranteil bei der Anpassung der Parameter ist.
    • Zu kleines α: Langsame Anpassung des Modells
    • Zu großes α: Sprünge und Divergenz (keine Konvergenz in Richtung Fehler 0)
96 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 9: KNN / Perzeptron

+

97 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 9: KNN / Perzeptron

Bisher wurde eine lineare Aktivierungsfunktion angenommen. Stückweise ist die sigmoid Funktion auch linear, daher "geht so".

  • Verbesserte Anpassung der Parameter mit dem Ziel der Fehlerminimierung (absteigender Gradient) mit Differenzierung der Perzeptron Funktion (bzw. deren Fehler) nach den einzelnen Parametern:

wiwi+ΔwiΔwi=αEwiE(xj)=yja(xj)

  • In der Neuronfunktion a stecken ja zwei Funktionen drin. Die Ableitung der Summenfunktion u ist konstant, bleibt nur noch die Aktivierungsfunktion f, so dass dann gilt:

Δwi=αEf'(u)xi

98 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 9: KNN / Perzeptron

Ableitung von einer Funktion? Geschlossen analytisch oder approximativ numerisch?

  • Es gibt verschiedene gängige Aktivierungsfunktionen wie z.B. die sigmoid Funktion.

  • U.A. die sigmoid Funktion hat eine Eigenschaft die es erlaubt die Ableitung f' wieder mit f auszudrücken ("Selbstähnlichkeit")!

f(u)=11+euf'(u)=f(u)(1f(u))

  • Numerisch: Auch angenehm da die sigmoid Funktion einen begrenzter Gradienten besitzt, nämlich [0,0.25] da der Wertebereich von der sigmoid Funktion auf (0,1) begrenzt ist!
    • Apprximationsfehler sind auch bei größeren h Intervall (aber h < 1) klein!
    • Numerische Fehler sind nicht zu erwarten (also Divergenz)
99 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 9: KNN / Perzeptron

+

100 / 101

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Zusammenfassung

Zusammenfassung

  1. Numerik ist Algorithmik.
  2. Datenstrukturen sind hier:
    • Funktionen (Werte erster Ordnung, Lambda Ausdrücke!)
    • Vektoren
    • Matrizen
  3. Operationen (Funktionen) auf den Daten sind hier:
    • Iterative approximative Berechnungen von Zahlen und Funktionen mit Summen
    • Integration von Funktionen
    • Differenzierung von Funktionen (Gradienten)
    • Regression
    • Vektor- und Matrixalgebra
    • Lösung linearer Gleichungssysteme (V/M Algebra)
    • Maschinelles Lernen (V/M Algebra, Gradient)
  4. Korrektheit == Konvergenz + Genauigkeit abhängig von Verfahren und Daten!
  5. Effizienz == Abhängig von Verfahren und Daten, i.A. polynomielle Laufzeitklasse O(n|n2|n3)!
101 / 101