1. Introduction

1.1 Rundungsfehler und Fehlerpropagation

Maschinenzahlen (1.1.1–1.1.3)

Absoluter und relativer Fehler

Sei $\tilde{x}\in {\mathbb{R}}^n$ eine Approximation von $x\in {\mathbb{R}}^n$.

• der absolute Fehler ist $\Vert \tilde{x}-x\Vert$
• für $x\ne 0$ ist der relative Fehler $\frac{\Vert \tilde{x}-x\Vert }{\Vert x\Vert }$

Maschinenzahlen

Seien eine Basis $B\in \mathbb{N}$, $B\ge 2$, eine Mantissenlänge $m\in \mathbb{N}$ und ein Exponentenbereich $\{e_{\min },\dots ,e_{\max }\}\subset \mathbb{Z}$ gegeben. Die Menge der Maschinenzahlen ist

$$\mathbb{M}=\{0\}\cup \left\{\pm 0.a_1{a}_2\dots a_m\cdot B^e:a_i\in \{0,\dots ,B-1\},a_1\ne 0,e_{\min }\le e\le e_{\max }\right\}$$

Die Bedingung $a_1\ne 0$ heisst Normalisierung.

Hier ist die Zahlenlinie, mit den darauf liegenden floats. Für float64 ist $B=2$ und $m=53$.

Note: $\mathbb{M}$ ist eine endliche, diskrete Teilmenge von $\mathbb{R}$, die unter $+,-,\cdot ,/$ nicht abgeschlossen ist, also ist Rundung unvermeidbar.

Note (eps): Die Konstante $\text{eps}:=B^{1-m}$ ist der Abstand von $1$ zur nächstgrösseren Maschinenzahl, für float64 also $\text{eps}=2^{-52}\approx 2.2\cdot 10^{-16}$ (np.finfo(np.float64).eps).

Der IEEE-754-Standard garantiert außerdem, dass jede Grundoperation korrekt gerundet wird: Für $x,y\in \mathbb{M}$ und $\circ \in \{+,-,\cdot ,/\}$ gilt $\mathrm{fl}(x\circ y)=(x\circ y)(1+\delta )$ mit $|\delta |\le \text{eps}$. Eine einzelne Operation ist also immer relativ genau.

Auslöschung (1.1.9–1.1.10)

Durch Subtraktion ergibt dich das folgende noise amplifizierende Problem.

Der Verstärkungsfaktor $\frac{|x|+|y|}{|x-y|}$ wird beliebig gross, wenn $x\approx y$.

Der Rundungsfehler der Subtraktion selbst ist dabei vernachlässigbar. Die bereits vorhandenen Fehler werden lediglich relativ zum kleinen Ergebnis sichtbar.

Example: Seien $x$ und $y$ mit relativem Fehler der Grössenordnung $10^{-7}$ behaftet, d.h. ab der 7. Nachkommastelle verfälscht (markiert mit $\ast$):

$$\begin{array}{rl}x& =0.123467\ast \\ y& =0.123456\ast \\ x-y& =0.000011\ast =0.11\ast \cdot 10^{-4}\end{array}$$

Die verfälschten Ziffern rücken nach vorne, und das Ergebnis ist bereits ab der 3. Stelle falsch. Der relative Fehler ist also von $10^{-7}$ auf $10^{-3}$ gewachsen.

Beispiel: Auslöschung beim Differenzquotienten $f^{\prime }(x):=\frac{f(x+h)-f(x)}{h}$

Rein imaginärer Schritt

• Nehmen an Funktion sehr glatt
  • dann muss sie analytisch sein → darstellbar als Potenzreihe
  • Mit Taylor finden wir die Koeffizienten.
• daraus kriegen wir $f^{\prime }(x_0)\approx \text{Im}(f(x_0+ih))/h$
  • die Subtraktion ist weg → keine Auslöschung
  • Die Restterme sind $O(h^2)$, konvergiert also schneller als naive methode

Richardson

• wir starten mit dem symmetrischen Differenzquotient $d(h):=\frac{f(x+h)-f(x-h)}{2h}$
• Taylor-Entwicklungen von $f(x+h)$ und $f(x-h)$ heben sich alle geraden Potenzen von $h$ auf
  • $f^{\prime }(x)-d(h)=c_1{h}^2+c_2{h}^4+c_3{h}^6+\dots \text{mit}{c}_k=-\frac{f^{(2k+1)}(x)}{(2k+1)!}$
  • die Koeffizienten $c_k$ nicht von $h$ abhängen und nur gerade Potenzen auftreten
• $f^{\prime }(x)\approx$ $R_{\ell ,k}:=\frac{4^k{R}_{\ell ,k-1}-R_{\ell -1,k-1}}{4^k-1}$ mit $l$ dem kleinsten $h$ und $k$ der Elimination
• Es gilt dann für den Fehler $f^{\prime }(x)-R_{\ell ,k}=O\left({\left(\frac{h}{2^{\ell }}\right)}^{2k+2}\right)$
• wir haben also nicht die Auslöschung verhindert, die Approximation konvergiert aber so schnell, dass diese zu vernachlässigen ist.

V-Form: Verkleinern von $h$ hilft nur bis $h^{\ast }\sim \sqrt{\text{eps}}$, danach dominiert die Auslöschung.

• nie eine bessere Genauigkeit als $\sqrt{\text{eps}}\approx 10^{-8}$

→ Schranke ist pessimistisch, da sich Rundungsfehler teilweise gegenseitig aufheben, aber die Größenordnung von $h^{\ast }$ stimmt

Auslöschung lässt sich auf zwei Arten bekämpfen:

• schneller Konvergieren (Richardson)
• vermeidet die Subtraktion vollständig (imaginärer Schritt)

1.2 Rechenaufwand

Die asymptotische Komplexität zählt elementare Operationen ($+,-,\cdot ,/$) im Worst Case für $n\to \infty$.
→ kein verlässliches Mass für die tatsächliche Laufzeit

Hängt auf modernen Rechnern von … ab:

• RAM bandwidth
• cache levels
• prozessor kerne, BLAS/CUDA, etc…

1.3 Rechnen mit Matrizen

Durch schlaues erkennen von Strukturen, kann man Matrixmultiplikation oft schneller als $O(n^3)$ machen!

• Struktur ausnutzen: Statt $DA$ mit Diagonalmatrix $D$ in $O(n^3)$ zu berechnen, skaliert man die Zeilen: D.diagonal()[:, np.newaxis] * A in $O(n^2)$.

• Klammerung wählen: $a(b^{\top }x)$ kostet $O(m+n)$, $(a{b}^{\top })x$ dagegen $O(mn)$.

• Niedrigen Rang ausnutzen: $\mathrm{triu}(A{B}^{\top })x$ mit $A,B\in {\mathbb{R}}^{n,p}$ berechnet man spaltenweise über Teilsummen (cumsum) in $O(np)$ statt $O(p{n}^2)$.

• Kronecker-Produkt nie aufstellen: Mit $X=[x_1,\dots ,x_n]$ ist $(A\otimes B)x=\mathrm{vec}(BX{A}^{\top })$, also genügen zwei kleine Matrixprodukte.

• Vektorisieren: Keine expliziten Schleifen über Array-Elemente, und nie np.vectorize verwenden, da es intern doch nur eine Schleife ist.

• Speicher bewusst verwalten:

  • Views statt Kopien verwenden.
  • Arrays einmal vorab mit dem richtigen dtype allozieren, statt sie zu erweitern.
  • Eine echte Kopie erzwingen, wo sie nötig ist.