\documentclass[11pt,ngerman]{article} \usepackage[ngerman]{babel} \usepackage{fontspec} \usepackage[a4paper, bindingoffset=0.0in, left=25mm, right=25mm, top=25mm, bottom=25mm, footskip=10mm]{geometry} \usepackage{enumitem} \usepackage{epsfig} \usepackage{epsf} \usepackage{amsmath,amssymb,amsfonts,amsthm,bm} \usepackage{mathtools} \usepackage{mathrsfs} \usepackage{array} \usepackage{tabularx} \usepackage{multirow} \usepackage{longtable} \usepackage{pstricks} \usepackage{lipsum} \usepackage{ulem} \usepackage{blindtext} \usepackage{setspace} \usepackage{graphicx} \usepackage{algorithm} \usepackage{algpseudocode} \usepackage{booktabs} \usepackage{geometry} \usepackage{float} \setromanfont{Times New Roman} \setsansfont{Arial} % --- Geometrie und Layout --- \usepackage{titlesec} \usepackage{fancyhdr} \pagestyle{fancy} \fancyhf{} \fancyhead[L]{\leftmark} \fancyfoot[C]{\thepage} % --- Farben und Grafik --- \usepackage{xcolor} \definecolor{tumblue}{RGB}{0, 101, 165} \definecolor{alertred}{RGB}{180, 30, 30} \definecolor{successgreen}{RGB}{30, 130, 30} \definecolor{warmorange}{RGB}{200, 120, 0} \usepackage{tikz} \usetikzlibrary{arrows.meta, angles, quotes, calc} % --- Hyperlinks (nach geometry/titlesec/fancyhdr geladen) --- \usepackage{hyperref} \hypersetup{ colorlinks=true, linkcolor=tumblue, urlcolor=tumblue, citecolor=tumblue } % --- Benutzerbefehle --- \newcommand{\N}{\mathbb{N}} \newcommand{\Z}{\mathbb{Z}} \newcommand{\Q}{\mathbb{Q}} \newcommand{\R}{\mathbb{R}} \newcommand{\C}{\mathbb{C}} \renewcommand{\vec}[1]{\mathbf{#1}} \newcommand{\E}{\mathbb{E}} \renewcommand{\i}{\mathrm{i}} \newcommand{\cond}{\operatorname{cond}} \newcommand{\rank}{\operatorname{rank}} \newcommand{\tr}{\operatorname{tr}} \newcommand{\diag}{\operatorname{diag}} \newcommand{\sign}{\operatorname{sign}} \newcommand{\argmin}{\operatorname{arg\,min}} \newcommand{\diff}{\mathrm{d}} \newcommand{\re}{\operatorname{Re}} \newcommand{\im}{\operatorname{Im}} \newcommand{\argg}{\operatorname{arg}} \titleformat{\section}{\Large\bfseries\color{tumblue}}{\thesection}{1em}{} \titleformat{\subsection}{\large\bfseries\color{tumblue!80!black}}{\thesubsection}{1em}{} \titleformat{\subsubsection}{\normalsize\bfseries\color{tumblue!60!black}}{\thesubsubsection}{1em}{} % --- Boxen für wichtige Hinweise --- \usepackage{tcolorbox} \tcbuselibrary{skins, breakable} \newtcolorbox{warningbox}[1][]{ colback=alertred!5!white, colframe=alertred, fonttitle=\bfseries, title=#1, breakable } \newtcolorbox{tipbox}[1][]{ colback=tumblue!5!white, colframe=tumblue, fonttitle=\bfseries, title=#1, breakable } \newtcolorbox{successbox}[1][]{ colback=successgreen!5!white, colframe=successgreen, fonttitle=\bfseries, title=#1, breakable } \newtcolorbox{orangebox}[1][]{ colback=warmorange!5!white, colframe=warmorange, fonttitle=\bfseries, title=#1, breakable } \begin{document} \title{\vspace{-2cm}\textbf{\color{tumblue} Mathematik, Numerische Methoden \\ und Signalverarbeitung}\\[0.5cm] \large Ein Skript für Studierende der Fakultät für Chemie der TUM} \author{Ilgis Ibragimov \& Elena Ibragimova \& KI\footnote{Korrektur- und Übersetzungsassistent}} \date{September 2026 \\ DOI: 10.5281/zenodo.23002397} \maketitle \thispagestyle{empty} \vfill \begin{center} \textit{Gewidmet unseren Töchtern Maria und Alevtina. \\ Möge die Mathematik für euch nicht eine Sammlung von Formeln werden, \\ sondern die Sprache, in der das Universum geschrieben ist.} \end{center} \newpage \tableofcontents \newpage \section{Vorwort} Liebe Töchter, wir wollten so sehr, dass auf eurem Lebensweg alles klar und einfach ist. Wir wissen, dass Mathematik oft erklärt wird, indem man den Menschen absichtlich mit Terminologie und komplizierten Formeln verwirrt --- so dass man den Wald vor lauter Bäumen nicht sieht. Dieses Skript ist unser Versuch, einen Überblick über die mathematischen und rechnerischen Ideen zu geben, die für einen modernen Chemiker, der mit experimentellen Daten, Modellierung und Berechnungen arbeitet, besonders nützlich sind, damit du ein \textit{klares Bild} vom modernen Stand des mathematischen Wissens hast. Nicht allen Wissens --- das ist in einem Buch unmöglich --- sondern genau dessen, was in Chemie, Biochemie, Spektroskopie, molekularer Modellierung und Datenverarbeitung Anwendung findet. Wir haben gemeinsam einen langen Weg zurückgelegt: \begin{itemize} \item Von den Grundlagen --- Zahlenmengen, komplexe Zahlen, Vektoren und Matrizen, \item über die lineare Algebra --- SVD, Kondition, Methoden zur Lösung von Systemen, \item zu nichtlinearen Problemen --- Optimierung, Gradienten, Newton-Verfahren, \item durch die Signalverarbeitung --- Fourier, Prony, Compressed Sensing, \item zu numerischen Methoden --- Integration, finite Elemente, Basisfunktionen, \item und schließlich --- zum maschinellen Lernen, zu neuronalen Netzen und zur modernen Rechentechnik. \end{itemize} All das ist nicht nur eine Sammlung zusammenhangloser Themen. Es ist eine \textbf{einheitliche Sprache}, in der die moderne Wissenschaft spricht. Und wenn du diese Sprache beherrschst --- werden sich vor dir Türen öffnen, von denen du jetzt nicht einmal ahnst. \subsection{Wie man dieses Buch liest} Obwohl der Text so geschrieben ist, dass er möglichst einfach lesbar ist, ohne detaillierte mathematische Beweise, können manche Stellen zu abstrus formuliert sein. Das ist normal --- Mathematik kommt nicht auf einmal. Deshalb ist der Quelltext dieses Buches in \LaTeX{} beigefügt. Wenn etwas nicht ganz klar ist: \begin{enumerate} \item Finde den entsprechenden Text im Quelltext, \item Kopiere ihn in einen Chat (Qwen, DeepSeek, Kimi, ChatGPT, GROK, Gemini), \item Bitte ihn, dies verständlicher zu erklären, oder stelle die Frage, die du hier nicht verstehst. \end{enumerate} Mit Hilfe dieser Chats kann man auch den ganzen Text ins Englische oder Deutsche übersetzen --- dazu muss man nur bitten, unter Beibehaltung der \LaTeX{}-Auszeichnung zu übersetzen, und den erhaltenen Text mit dem Befehl \texttt{xelatex NumMathForChemists.tex} in eine PDF-Datei kompilieren. \begin{tipbox}[Tipp] Versuche nicht, alles auf einmal zu lesen. Mathematik ist wie Musik: man muss sie ``spielen'', also Aufgaben lösen, Code ausprobieren, experimentieren. Wenn ein Kapitel schwierig erscheint --- komm eine Woche später darauf zurück, und du wirst überrascht sein, wie viel einfacher es geworden ist. \end{tipbox} % \subsection{Wozu braucht ein Chemiker Mathematik?} % % Bevor wir anfangen, Zahlen, Matrizen, Operatoren und Fourier-Transformationen zu untersuchen, ist es nützlich, die natürlichste Frage zu beantworten: % \textit{Wozu braucht ein Chemiker das alles überhaupt?} % Die Antwort ist sehr einfach: weil die moderne Chemie immer häufiger nicht direkt mit Substanzen arbeitet, sondern mit \textbf{Daten, Modellen und Berechnungen}. % % Stellen wir uns ein ganz gewöhnliches Experiment vor. Wir bringen eine Probe in ein Gerät und erhalten ein gewisses Signal. Zum Beispiel kann ein Spektrometer mehrere Tausend Zahlen zurückgeben. Was wollen wir mit diesen Zahlen tun? % % Wir wollen Rauschen entfernen, Peaks finden, ihre Position und Intensität bestimmen, das erhaltene Spektrum mit bekannten Spektren vergleichen, die Konzentration eines Stoffes bestimmen oder die Parameter eines Moleküls rekonstruieren. Das heißt, das Experiment kann in sehr vereinfachter Form dargestellt werden als: % % $$ % \boxed{ % \text{Substanz} % \longrightarrow % \text{Messung} % \longrightarrow % \text{Daten} % \longrightarrow % \text{mathematisches Modell} % \longrightarrow % \text{chemische Schlussfolgerung} % } % $$ % % Und fast jeder Übergang in dieser Kette erfordert Mathematik. Deshalb kann der gesamte weitere Text als eine Reise entlang einer Kette aufgefasst werden: % % $$ % \boxed{ % \begin{array}{c} % \text{Zahlen}\\ % \downarrow\\ % \text{Vektoren und Matrizen}\\ % \downarrow\\ % \text{Operatoren}\\ % \downarrow\\ % \text{lineare und nichtlineare Probleme}\\ % \downarrow\\ % \text{Statistik und Unsicherheit}\\ % \downarrow\\ % \text{Signale und Fourier}\\ % \downarrow\\ % \text{numerische Methoden}\\ % \downarrow\\ % \text{Kompression und Informationsextraktion}\\ % \downarrow\\ % \text{maschinelles Lernen}\\ % \downarrow\\ % \text{moderne Computerchemie} % \end{array} % } % $$ % Es ist nicht notwendig, sich alle Formeln zu merken. Es ist viel wichtiger, allmählich zu lernen, die mathematische Struktur einer Aufgabe zu erkennen. % % % \begin{itemize} % \item Wenn eine Menge von Zahlen vor dir erscheint, frage: \textit{Was für ein mathematisches Objekt ist das?} % \item Wenn viele Messungen erscheinen: \textit{Kann man sie als Vektor oder Matrix darstellen?} % \item Wenn es unbekannte Parameter gibt: \textit{Wie baut man ein Modell und findet diese Parameter?} % \item Wenn es Rauschen gibt: \textit{Wie trennt man Information von zufälligem Fehler?} % \item Wenn es zu viele Daten gibt: \textit{Kann man verborgene Struktur finden oder die Dimension reduzieren?} % \end{itemize} % % % Genau diese Denkweise, und nicht das Auswendiglernen einer großen Zahl von Formeln, ist das Hauptziel dieses Skripts. % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Einleitung: Bezeichnungen und die Natur der Zahlen} \subsection{Die Hierarchie der Zahlenmengen} Bevor wir anfangen, über komplizierte Dinge zu sprechen, lass uns über die Sprache einig werden. In der Mathematik gruppieren wir Zahlen in Mengen, die bestimmte Eigenschaften haben. Wir verwenden die folgenden Standardbezeichnungen: \begin{itemize} \item $\N = \{1, 2, 3, \dots\}$ --- \textbf{natürliche Zahlen}. Sie werden zum Zählen verwendet. \item $\Z = \{\dots, -2, -1, 0, 1, 2, \dots\}$ --- \textbf{ganze Zahlen}. Sie fügen die Null und negative Zahlen hinzu und erlauben es, Schulden und Vermögen, Temperaturen unter Null usw. zu beschreiben. \item $\Q$ --- \textbf{rationale Zahlen}. Jede Zahl, die als Bruch $\frac{p}{q}$ dargestellt werden kann, wobei $p \in \Z$, $q \in \Z \setminus \{0\}$. \item $\R$ --- \textbf{reelle Zahlen}. Sie umfassen alle rationalen sowie irrationale Zahlen (zum Beispiel $\sqrt{2}, \pi, e$), die nicht als gewöhnlicher Bruch dargestellt werden können. Sie füllen die gesamte Zahlengerade kontinuierlich aus. \item $\C$ --- \textbf{komplexe Zahlen}. Über sie sprechen wir weiter unten ausführlich. \end{itemize} \textit{Verschachtelung der Mengen:} $\N \subset \Z \subset \Q \subset \R \subset \C$. \subsection{Zahlen im Computer} Wir sind daran gewöhnt, dass Zahlen in der Mathematik unendlich genau sind. Aber wenn wir zur \textbf{numerischen Mathematik} und zum Programmieren übergehen, ist der Computer gezwungen, Zahlen in Speicherzellen fester Größe (in Bits und Bytes) zu speichern. \begin{itemize} \item \textbf{Ganzzahlige Typen} (int8, int16, int32, int64) speichern exakte Werte aus $\Z$, sind aber im Bereich begrenzt (zum Beispiel kann int32 Zahlen von $-2^{31}$ bis $2^{31}-1$ speichern). \item \textbf{Reelle Typen} (float32, float64 / double) speichern Zahlen aus $\R$ im Exponentialformat (Mantisse und Exponent). Dadurch entstehen Rundungsfehler: das computerliche $0.1 + 0.2$ ist nicht immer exakt gleich $0.3$. Das zu verstehen ist für numerische Methoden von entscheidender Bedeutung! \end{itemize} \section{Wenn eine Zahl nicht genügt: Komplexe Zahlen} Unsere Welt ist so komplex, dass nicht alles durch eine einzige Zahl beschrieben werden kann. In Physik und Chemie müssen wir oft mit Entitäten arbeiten, die \textit{zwei} Zahlen zu ihrer Beschreibung benötigen (zum Beispiel Amplitude und Phase einer Welle oder Real- und Imaginärteil einer Wellenfunktion). \subsection{Woher kommen sie?} Der einfachste und historischste Weg, ``Zahlenpaare'' kennenzulernen, ist der Versuch, die Gleichung zu lösen: \[ x^2 + 1 = 0 \quad \Rightarrow \quad x^2 = -1 \] In der Menge der reellen Zahlen $\R$ ist das Quadrat jeder Zahl nichtnegativ. Es gibt keine Wurzel aus $-1$. Aber die Mathematiker (und nach ihnen die Physiker) sagten: ``Was, wenn wir einfach eine Zahl \textit{erfinden}, die mit sich selbst multipliziert $-1$ ergibt?''. So entstand die \textbf{imaginäre Einheit} $\i$. Wir ``ziehen nicht die Wurzel'' aus $-1$, wir \textit{definieren} eine neue Entität: \[ \i^2 = -1 \] \textit{Wichtige Bemerkung: Die Regel $\sqrt{a}\sqrt{b} = \sqrt{ab}$ gilt nur für $a,b \ge 0$. Deshalb schreiben wir nicht $\sqrt{-1}\sqrt{-1} = \sqrt{(-1)(-1)} = 1$. Wir akzeptieren einfach $\i^2 = -1$ als Tatsache.} \subsection{Algebraische Form und komplexe Ebene} Jede komplexe Zahl $z \in \C$ kann in der Form geschrieben werden: \[ z = x + \i y \] wobei $x = \re(z)$ der \textbf{Realteil} und $y = \im(z)$ der \textbf{Imaginärteil} ist. Beide Teile $x, y \in \R$. Geometrisch ist eine komplexe Zahl ein Punkt in der \textbf{komplexen Ebene}. \begin{itemize} \item Auf der horizontalen Achse (der $\re$-Achse) wird der Realteil $x$ aufgetragen. \item Auf der vertikalen Achse (der $\im$-Achse) wird der Imaginärteil $y$ aufgetragen. \end{itemize} Eine komplexe Zahl kann auch als \textbf{Vektor} betrachtet werden, der vom Ursprung $(0,0)$ zum Punkt $(x,y)$ geht. \subsection{Arithmetik komplexer Zahlen} Die \textbf{Addition} ist intuitiv klar. Wenn wir zwei Zahlen $z_1 = x_1 + \i y_1$ und $z_2 = x_2 + \i y_2$ haben, addieren wir einfach ihre Real- und Imaginärteile getrennt: \[ z_1 + z_2 = (x_1 + x_2) + \i(y_1 + y_2) \] Geometrisch funktioniert dies wie die \textbf{Dreiecksregel} für Vektoren: Wir zeichnen zwei Vektoren vom Ursprung, verschieben den Anfang des zweiten Vektors an das Ende des ersten, und die Summe ist der Vektor vom Anfang des ersten bis zum Ende des zweiten. Die \textbf{Multiplikation} in algebraischer Form sieht etwas komplizierter aus. Wir lösen die Klammern wie in der gewöhnlichen Algebra auf und denken daran, dass $\i^2 = -1$: \begin{align*} z_1 \cdot z_2 &= (x_1 + \i y_1)(x_2 + \i y_2) \\ &= x_1 x_2 + \i x_1 y_2 + \i y_1 x_2 + \i^2 y_1 y_2 \\ &= (x_1 x_2 - y_1 y_2) + \i(x_1 y_2 + y_1 x_2) \end{align*} Aber um die \textit{wahre Magie} der Multiplikation komplexer Zahlen zu verstehen, müssen wir zu einer anderen Schreibweise übergehen. \subsection{Trigonometrische und Exponentialform} Statt der Koordinaten $(x,y)$ kann ein Punkt in der Ebene durch Polarkoordinaten angegeben werden: den Abstand vom Ursprung (den Betrag) $r$ und den Winkel $\varphi$ relativ zur positiven Richtung der $\re$-Achse. Der Zusammenhang mit der algebraischen Form ist aus der Trigonometrie offensichtlich: \[ x = r \cos \varphi, \quad y = r \sin \varphi \] Dann nimmt die komplexe Zahl die \textbf{trigonometrische Form} an: \[ z = r(\cos \varphi + \i \sin \varphi) \] \subsubsection{Eulers Formel und die Exponentialform} Hier betritt die größte Formel der Mathematik die Bühne --- \textbf{Eulers Formel}: \[ e^{\i\varphi} = \cos \varphi + \i \sin \varphi \] Dank ihr kann jede komplexe Zahl in der unglaublich bequemen \textbf{Exponentialform} geschrieben werden: \[ z = r e^{\i\varphi} \] \subsubsection{Die geometrische Bedeutung der Multiplikation} Betrachten wir, was bei der Multiplikation zweier Zahlen in Exponentialform geschieht: \[ z_1 = r_1 e^{\i\varphi_1}, \quad z_2 = r_2 e^{\i\varphi_2} \] \[ z_1 \cdot z_2 = (r_1 e^{\i\varphi_1}) \cdot (r_2 e^{\i\varphi_2}) = (r_1 r_2) e^{\i(\varphi_1 + \varphi_2)} \] \textbf{Die Schlussfolgerung, die man sich für immer merken muss:} Die Multiplikation komplexer Zahlen ist die \textbf{Multiplikation ihrer Längen (Beträge)} und die \textbf{Addition ihrer Winkel (Argumente)}! Geometrisch: Die Multiplikation mit einer komplexen Zahl ist eine Drehung des Vektors um den Winkel $\varphi$ und seine Streckung/Stauchung um den Faktor $r$. \subsection{Komplexe Exponentialfunktionen und Logarithmen} \subsubsection{Eigenschaften der Exponentialfunktion} Da sich $e^{\i\varphi}$ wie eine gewöhnliche Exponentialfunktion verhält, gelten für sie alle Standardregeln: \[ e^{z_1} \cdot e^{z_2} = e^{z_1 + z_2}, \quad \frac{e^{z_1}}{e^{z_2}} = e^{z_1 - z_2}, \quad (e^z)^n = e^{nz} \] Für eine komplexe Zahl $z = x + \i y$ zerfällt die Exponentialfunktion in Real- und Imaginärteil: \[ e^z = e^{x + \i y} = e^x \cdot e^{\i y} = e^x (\cos y + \i \sin y) \] Der Betrag dieser Zahl ist $e^x$, und das Argument ist $y$. \subsubsection{Der komplexe Logarithmus} Der Logarithmus ist die zur Exponentialfunktion inverse Operation. Wenn $e^w = z$, dann ist $w = \ln z$. Sei $z = r e^{\i\varphi}$ und $w = u + \i v$. Dann: \[ e^{u + \i v} = r e^{\i\varphi} \quad \Rightarrow \quad e^u e^{\i v} = r e^{\i\varphi} \] Durch Gleichsetzen von Beträgen und Argumenten erhalten wir: \[ e^u = r \quad \Rightarrow \quad u = \ln r \] \[ v = \varphi + 2\pi k, \quad k \in \Z \] Somit ist der \textbf{komplexe Logarithmus} mehrwertig: \[ \ln z = \ln|z| + \i(\argg z + 2\pi k) \] Der Hauptwert des Logarithmus (bei $k=0$ und $\argg z \in (-\pi, \pi]$) wird mit $\operatorname{Ln} z$ bezeichnet. Die Mehrwertigkeit entsteht dadurch, dass eine Drehung um $2\pi$ uns zum selben Punkt in der komplexen Ebene zurückbringt. \subsection{Spickzettel: Trigonometrie und hyperbolische Funktionen} Hyperbolische Funktionen erschrecken Chemiestudierende oft, aber tatsächlich sind sie nahe Verwandte der gewöhnlichen Sinus- und Kosinusfunktionen, die einfach in der komplexen Ebene ``leben''. \subsubsection{Definitionen über die Exponentialfunktion} \begin{align*} \cos z &= \frac{e^{\i z} + e^{-\i z}}{2}, & \sin z &= \frac{e^{\i z} - e^{-\i z}}{2\i} \\ \cosh z &= \frac{e^{z} + e^{-z}}{2}, & \sinh z &= \frac{e^{z} - e^{-z}}{2} \end{align*} \subsubsection{Zusammenhang zwischen trigonometrischen und hyperbolischen Funktionen} Wenn wir $\i z$ anstelle von $z$ in die Definitionen einsetzen, sehen wir eine erstaunliche Symmetrie (Osbornsche Regeln): \begin{align*} \cos(\i z) &= \cosh z, & \cosh(\i z) &= \cos z \\ \sin(\i z) &= \i\sinh z, & \sinh(\i z) &= \i\sin z \end{align*} \textit{Eselsbrücke:} Beim Übergang von der Trigonometrie zur Hyperbolik wird das Argument mit $\i$ multipliziert, und vor Sinus/Kosinus können imaginäre Einheiten auftreten. Hyperbolische Funktionen beschreiben nicht Schwingungen (wie der Sinus), sondern exponentielles Wachstum/Abklingen (wie $e^x$), was für gedämpfte Prozesse von entscheidender Bedeutung ist. \subsection{Wozu braucht man das in der realen Chemie und Physik?} Es mag scheinen, dass imaginäre Zahlen nur eine schöne mathematische Abstraktion sind. Aber unsere Welt ist so eingerichtet, dass wir \textit{ohne} den imaginären Raum reale Dinge nicht beschreiben könnten. \subsubsection{Quantenmechanik und Orbitale} Die Schrödinger-Gleichung, die das Verhalten von Elektronen in Atomen beschreibt, enthält die imaginäre Einheit $\i$ explizit: \[ \i\hbar \frac{\partial \Psi}{\partial t} = \hat{H}\Psi \] Die Wellenfunktion $\Psi$ ist komplexwertig. Für den Grundzustand des Wasserstoffatoms (1s-Orbital) ist die Wellenfunktion reell, und man kann sie im realen Raum zeichnen. Sobald das Elektron jedoch in einen angeregten Zustand übergeht (zum Beispiel ein 2p-Orbital mit magnetischer Quantenzahl $m \neq 0$), enthält seine Wellenfunktion den Phasenfaktor $e^{\i m\varphi}$. Ohne komplexe Zahlen ist es einfach unmöglich, den Drehimpuls des Elektrons und die Feinstruktur der Spektren zu beschreiben. \subsubsection{Kernspinresonanz (NMR) und Fourier-Analyse} Du wirst viel mit NMR-Spektroskopie arbeiten. Wenn du eine Probe in ein Magnetfeld bringst und einen Radiofrequenzimpuls gibst, beginnen die Kerne zu präzedieren. Der Detektor registriert ein mit der Zeit abklingendes Signal --- dies nennt man \textbf{FID} (Free Induction Decay). Das Signal von einem Kerntyp sieht wie eine gedämpfte Schwingung aus: \[ S(t) = A e^{-t/T_2} \cos(\omega t) \] Wenn die Probe viele verschiedene Kerne enthält, verwandelt sich das Signal in einen Brei aus vielen überlagerten Kosinusfunktionen. Wie kann man herausfinden, welche Frequenzen $\omega$ darin verborgen sind? Hier kommt die \textbf{Fourier-Transformation} zur Hilfe. Aber mit Kosinusfunktionen zu arbeiten ist schwierig. Es ist viel einfacher, zu komplexen Exponentialfunktionen überzugehen, indem man Eulers Formel verwendet: \[ \cos(\omega t) = \frac{e^{\i\omega t} + e^{-\i\omega t}}{2} \] In NMR-Spektrometern verwendet man Quadraturdetektion, die es erlaubt, direkt ein komplexes Signal zu messen: \[ S_{complex}(t) = A e^{-t/T_2} e^{\i\omega t} \] Wenn wir auf dieses Signal die schnelle Fourier-Transformation (FFT) anwenden, geht die Zeit $t$ in die Frequenz $\omega$ über. Eine lange oszillierende Funktion im Zeitbereich verwandelt sich in einen \textbf{schmalen Peak (eine Lorentz-Kurve)} im Frequenzbereich! Jeder Kerntyp entspricht seiner eigenen Frequenz $\omega$, und im resultierenden NMR-Spektrum sehen wir einzelne Peaks. Die gesamte moderne Signalverarbeitung (von der MRT in der Medizin bis zum Audio-MP3) funktioniert genau dank der Magie der komplexen Zahlen und der Euler-Exponentialfunktionen. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Wenn eine Messung nicht genügt: Vektoren und Normen} \subsection{Von der Ebene zum N-dimensionalen Raum} Nun, wenn wir komplexe Zahlen haben (im Wesentlichen Zahlenpaare), dann gibt es sicher auch etwas Komplizierteres? Quaternionen? Oktaven? Ja, tatsächlich haben Mathematiker Zahlen von der Ebene auf höhere Dimensionen verallgemeinert. Aber jenseits der komplexen Zahlen und der Quaternionen (die sich übrigens in der Robotik und der 3D-Grafik zur Beschreibung von Drehungen in unserem dreidimensionalen Raum hervorragend eingebürgert haben) werden spezielle ``Zahlentypen'' praktisch nicht verwendet. Stattdessen gruppieren Menschen einfach Zahlen in Mengen und nennen sie \textbf{N-dimensionale Vektoren}. Das ist einfach der N-dimensionale Raum. $N$ kann zwei sein (dann ist es eine Ebene), drei (unser gewöhnlicher physikalischer Raum) oder sehr, sehr groß. Solche Vektoren bezeichnen wir fett, zum Beispiel $\vec{v} \in \R^N$. Das ist ein Vektor: \[ \vec{v} = (v_1, v_2, \dots, v_N)^T \] \textit{Hier bedeutet $T$ Transposition, um eine Spalte in eine Textzeile zu schreiben.} Wir sind frei, uns auszudenken, was wir mit diesem Raum machen, aber das Erste, was wir brauchen, ist zu verstehen, wie man die ``Größe'' oder ``Länge'' eines solchen Vektors misst. \subsection{Wozu brauchen wir Vektoren? Beispiele aus der Chemie} Auf den ersten Blick ist ein Vektor einfach eine Menge von Zahlen. Aber woher kommen die Zahlen? Angenommen, es sind Daten von einem Chromatographen oder einem NMR-Spektrometer. Stell dir vor, wir wollen herausfinden, welcher Peak hier der ``hellste'' ist. Das ist doch einfach das Maximum unter diesen Zahlen, oder? Und nun nimm ein Chromatogramm von irgendeinem komplexen Schmutz und ein Chromatogramm einer einzigen reinen Substanz. Wir injizieren beide in das Chromatograph in gleichem molarem Volumen (wir wollen vereinbaren, dass der Detektor alle Moleküle gleich empfindet). \begin{itemize} \item Im Fall der reinen Substanz erhalten wir einen schönen, sehr scharfen und hohen Peak. \item Im Fall des Schmutzes erhalten wir einen Berg von Peaks, vielleicht sogar zu einem durchgehenden flachen Buckel verschmolzen. \end{itemize} Aber interessant ist: Die \textbf{Fläche} (integrierte Intensität) beider Chromatogramme wird gleich sein! Der Peak der reinen Substanz wird im Vergleich zum Berg sehr hoch sein, aber die Summe aller Antworten ist gleich der Stoffmenge. Aus diesem Problem ergeben sich auf natürliche Weise drei verschiedene Möglichkeiten, unseren Datenvektor zu bewerten: \begin{enumerate} \item \textbf{Maximum.} Uns interessiert die maximale Konzentration (oder die maximale Absorption, wenn der Peak ``nach unten schaut''). Wir nehmen den größten Wert. \item \textbf{Summe.} Uns interessiert die Gesamtmenge der Substanz. Wir addieren alle Werte (meist im Betrag, da die Grundlinie nicht ideal sein muss). \item \textbf{Wurzel aus der Summe der Quadrate.} Und wozu braucht man das? Stell dir vor, wir beschreiben ein Molekül nicht durch ein Chromatogramm, sondern durch drei Parameter: \textit{Polarität, Molmasse, Siedepunkt}. Das ist ein Vektor in $\R^3$. Um zu verstehen, wie \textit{ähnlich} zwei Moleküle einander sind (zum Beispiel für das Wirkstoffdesign), können wir nicht einfach die Differenzen der Parameter addieren oder die maximale Differenz nehmen. Wir brauchen genau den \textbf{euklidischen Abstand} --- den kürzesten Weg in diesem ``chemischen Raum''. Oder ein anderes Beispiel: In der Physik ist die quadratische Norm eines Vektors (eines Spektrums) proportional zur \textbf{Gesamtenergie} dieses Signals. \end{enumerate} \subsection{Die mathematische Sprache der Normen} Alle diese drei intuitiven Begriffe heißen in der Mathematik \textbf{Normen von Vektoren} und werden durch doppelte senkrechte Striche $\| \vec{v} \|$ bezeichnet. \begin{itemize} \item \textbf{Maximumnorm ($L_\infty$):} \[ \| \vec{v} \|_\infty = \max_{1 \le k \le N} |v_k| \] \item \textbf{Manhattan-Abstand / Summe der Beträge ($L_1$):} \[ \| \vec{v} \|_1 = \sum_{k=1}^N |v_k| \] \item \textbf{Euklidische Norm / Länge des Vektors ($L_2$):} \[ \| \vec{v} \|_2 = \sqrt{\sum_{k=1}^N |v_k|^2} \] \end{itemize} Aber Mathematiker lieben Ordnung, deshalb haben sie all das in einer allgemeinen Formel für die \textbf{$L_p$-Norm} zusammengefasst: \[ \| \vec{v} \|_p = \left( \sum_{k=1}^N |v_k|^p \right)^{1/p} \] In dieser Formel ist $L_1$ der Fall $p=1$. $L_2$ ist $p=2$. Und was wird $p$ für das Maximum? Ja, genau, $p \to \infty$. Wenn man den Grenzwert dieser Formel für $p \to \infty$ bildet, erhält man mathematisch streng genau das maximale Element des Vektors. \subsection{Der kontinuierliche Fall: Funktionen als Vektoren} Darüber hinaus kann unser Vektor ganz kontinuierlich sein. Dann ist er tatsächlich eine \textbf{Funktion} $f(x)$. Vorerst merken wir uns nur, dass man eine Funktion als ``unendlich-dimensionalen Vektor'' betrachten kann, dessen Werte an jedem Punkt $x$ gegeben sind. Und für Funktionen funktionieren die Normen genauso, nur ersetzen wir die Summe durch ein Integral: \[ \| f \|_p = \left( \int |f(x)|^p \, dx \right)^{1/p} \] Wir werden manchmal mit einem Array (einem Vektor), manchmal mit einer Funktion operieren, und später, im Kurs über Signalverarbeitung, werden wir sehen, dass dies zwei Seiten derselben Medaille sind. \subsection{Minkowskis Ungleichung} Es gibt einen entscheidend wichtigen Punkt, den man sich merken muss. Wir werden oft $\| \vec{x} + \vec{y} \|$ mit $\| \vec{x} \|$ und $\| \vec{y} \|$ vergleichen. Intuitiv ist klar, dass, wenn wir von Punkt A nach Punkt B und dann von B nach C gehen, der Gesamtweg nicht kürzer sein wird, als wenn wir direkt von A nach C gegangen wären. Diese geometrische Eigenschaft (die Länge einer Dreiecksseite ist nicht größer als die Summe der beiden anderen) wird für beliebige $L_p$-Normen in \textbf{Minkowskis Ungleichung} formalisiert: \[ \| \vec{x} + \vec{y} \|_p \le \| \vec{x} \|_p + \| \vec{y} \|_p \] \textit{Wie verwendet man das?} Es erlaubt uns, Fehler abzuschätzen. Wenn $\vec{x}$ das wahre Signal und $\vec{y}$ das Rauschen (Messfehler) ist, dann garantiert Minkowskis Ungleichung, dass die Norm (Größe) unseres gemessenen Signals $(\vec{x}+\vec{y})$ die Summe der Norm des wahren Signals und der Norm des Rauschens nicht überschreitet. \textit{(Einen strengen und schönen Beweis dieser Ungleichung findet man in \href{https://de.wikipedia.org/wiki/Minkowski-Ungleichung}{Wikipedia}, um dieses Skript nicht zu überladen.)} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Zwei Dimensionen und mehr: Matrizen und Tensoren} \subsection{Von Vektoren zu Matrizen} Da wir einen Vektor $\vec{a} \in \R^N$ haben (eine eindimensionale Entität, das diskrete Analogon der Funktion $f(x)$), ist es logisch zu fragen: Was, wenn wir eine Funktion zweier Variablen $f(x,y)$ haben? Was ist ihr diskretes Analogon? Und für $f(x,y,z)$? Ja, zweidimensionale Objekte im diskreten Raum sind \textbf{Matrizen}. Und wie im Fall der komplexen Zahlen ist hier alles sehr, sehr gut untersucht. Und alles, was drei oder mehr Dimensionen hat, nennt man üblicherweise \textbf{Tensoren}. Übrigens, wenn du jemals auf alte wissenschaftliche Literatur über Chemometrie oder Datenanalyse stößt, wirst du dich über den Zoo der Namen für Tensoren wundern. Sie werden \textit{multilinear, multidimensional, procrustes, multimodal, three-way, multi-way arrays} genannt. Erschrick nicht: Hinter all diesen schönen Wörtern verbirgt sich einfach ein mehrdimensionales Zahlenarray. Vorerst bleiben wir bei zweidimensionalen Objekten: der Funktion $f(x,y)$ und der Matrix $\mathbf{A} \in \R^{N \times M}$. \subsection{Komplexe Vektoren und Matrizen} Warum habe ich bisher nur $\R$ geschrieben? Sowohl ein Vektor als auch eine Matrix können natürlich komplex sein! Wir können die Räume $\C^N$ und $\C^{N \times M}$ betrachten. Hier ist alles einfach: Statt einer reellen Zahl in jeder Zelle des Vektors oder der Matrix steht ein komplexes Paar. Aber wie berechnet man die Norm einer komplexen Zahl? Genauso! Der Betrag einer komplexen Zahl $z = x + \i y$ wird über die konjugierte Zahl $\bar{z} = x - \i y$ berechnet: \[ |z| = \sqrt{x^2 + y^2} = \sqrt{z \cdot \bar{z}} \] Entsprechend sieht die $L_2$-Norm eines komplexen Vektors $\vec{z} \in \C^N$ so aus: \[ \| \vec{z} \|_2 = \sqrt{\sum_{k=1}^N |z_k|^2} = \sqrt{\sum_{k=1}^N z_k \bar{z}_k} \] In Matrixschreibweise wird dies über die hermitesche Konjugation (Transposition + komplexe Konjugation) $\vec{z}^H$ geschrieben: \[ \| \vec{z} \|_2 = \sqrt{\vec{z}^H \vec{z}} \] \subsection{Normen von Matrizen} Und für Matrizen können wir ebenfalls Normen definieren! Aber hier gibt es eine wichtige Nuance. Die meisten einfachen Normen für Matrizen werden so berechnet: Wir ``schneiden'' die Matrix gedanklich in Zeilen, ordnen alle ihre $N \times M$ Zellen zu einem langen Spaltenvektor an und nehmen die Norm dieses Vektors. Die beliebteste dieser Normen ist die \textbf{Frobenius-Norm} (oder euklidische Norm für Matrizen), die das Analogon der $L_2$-Norm ist: \[ \| \mathbf{A} \|_F = \sqrt{\sum_{i=1}^N \sum_{j=1}^M |a_{ij}|^2} \] Für sie gibt es, wie für die gewöhnliche $L_2$, viele wichtige und schöne Eigenschaften, die wir bald betrachten werden. Jedoch gibt es neben solchen ``elementweisen'' Normen in der linearen Algebra spezielle \textbf{Operatornormen (induzierte Normen)}. Sie beantworten die Frage: ``Wie stark kann eine Matrix einen Vektor strecken, wenn wir sie mit ihm multiplizieren?''. Aber das ist eine Geschichte für das nächste Kapitel, in dem wir uns eingehend mit linearen Operatoren beschäftigen. %%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%% \section{Operatoren: wenn die Mathematik beginnt zu handeln} \subsection{Was ist ein Operator?} Wir sind zu einem sehr wichtigen Konzept gekommen, das überall in der Mathematik und in der Chemie vorkommt. Das ist der \textbf{Operator}. Ein Operator ist eine gewisse \textit{Wirkung} auf ein Objekt. Aber hier wird es sofort irgendwie verwirrend: eine Wirkung auf was? Auf eine Zahl? Auf einen Vektor? Auf eine Funktion? Betrachten wir gemeinsam mehrere Beispiele, und alles wird sich klären. \subsection{Eine Matrix als Operator} Angenommen, wir haben eine quadratische Matrix $\mathbf{A} \in \C^{N \times N}$ (ich werde im Folgenden fast immer $\C$ verwenden, da $\R$ nur eine Teilmenge von $\C$ ist und alles, was für reelle Zahlen funktioniert, auch für komplexe funktioniert). Betrachten wir die Wirkung $\mathbf{A}\vec{b}$, also die Multiplikation einer Matrix mit einem Vektor. Das Ergebnis ist wieder ein Vektor aus $\C^N$. Die Matrix \textit{transformiert} einen Vektor in einen anderen. Dies ist das einfachste Beispiel eines linearen Operators. \subsection{Analogie zum Funktionaloperator} Und nun --- das Interessanteste. Das Analogon eines Matrixoperators in der Welt der Funktionen ist der \textbf{Integraloperator}: \begin{equation} (\hat{K}g)(x) = \int_a^b K(x,y)\, g(y)\, \diff y \label{IntegralOperator} \end{equation} Hier ist $K(x,y)$ der \textit{Kern} des Operators (das Analogon der Matrix $\mathbf{A}$), $g(y)$ die Eingabefunktion (das Analogon des Vektors $\vec{b}$), und das Ergebnis ist eine neue Funktion von $x$. Wo, interessant, kommt eine solche Abstraktion im realen Leben vor? Es stellt sich heraus, überall! Zum Beispiel wirkt in der Quantenmechanik der Hamilton-Operator $\hat{H}$ auf die Wellenfunktion $\Psi$, und dies ist gerade ein Integro-Differentialoperator. In der Spektroskopie wird die Antwort eines Instruments auf ein Eingangssignal oft durch ein Faltungsintegral beschrieben --- auch das ist ein Operator. \subsection{Zurück zur Schule: ein Gleichungssystem} Aber fangen wir mit dem Einfachsten an. In der Schule haben wir bereits Gleichungssysteme gelöst, zum Beispiel: \[ \begin{cases} 2x - 3y = -4 \\ 3x + 2y = 9 \end{cases} \] Du denkst wahrscheinlich: ``Na und? Wir konnten das doch durch Einsetzen oder Additionsverfahren lösen''. Ja, konnten wir. Aber schreiben wir es in Matrixform $\mathbf{A}\vec{x} = \vec{b}$ um: \[ \underbrace{ \begin{pmatrix} 2 & -3 \\ 3 & 2 \end{pmatrix} }_{\mathbf{A}} \underbrace{ \begin{pmatrix} x \\ y \end{pmatrix} }_{\vec{x}} = \underbrace{ \begin{pmatrix} -4 \\ 9 \end{pmatrix} }_{\vec{b}} \] Du wirst mich berechtigterweise fragen: ``Warum so verkomplizieren?''. Und ich antworte: Nein, wir verkomplizieren nicht --- wir \textit{bringen Ordnung} in unser Wissen. Wir isolieren die \textit{Struktur} des Problems. Und nun können wir sagen: Wenn wir einen \textbf{linearen Operator} $\mathbf{A}$ haben, dann wirkt er auf den Vektor $\vec{x}$, und das Ergebnis ist der Vektor $\vec{b}$. \subsection{Erinnerung: grundlegende Operationen} Bevor wir weitergehen, halten wir ausdrücklich die Operationen fest, die wir ständig verwenden werden. \textbf{Ein Vektor als Spaltenmatrix.} Ein Vektor $\vec{x} \in \C^N$ ist tatsächlich eine Matrix der Größe $N \times 1$. Deshalb gelten alle Regeln der Matrixmultiplikation automatisch auch für Vektoren. \textbf{Transposition.} Die Operation $\mathbf{A}^T$ vertauscht Zeilen und Spalten: $(\mathbf{A}^T)_{ij} = A_{ji}$. Für komplexe Matrizen verwendet man häufiger die \textbf{hermitesche Konjugation} $\mathbf{A}^H = \overline{\mathbf{A}^T}$ (Transposition + komplexe Konjugation aller Elemente). \textbf{Skalarprodukt.} Für zwei Vektoren $\vec{x}, \vec{y} \in \C^N$ ist das Skalarprodukt definiert als: \[ \langle \vec{x}, \vec{y} \rangle = \vec{x}^H \vec{y} = \sum_{k=1}^N \bar{x}_k y_k \] Für reelle Vektoren ist dies einfach $\vec{x}^T \vec{y} = \sum x_k y_k$. \textbf{Die wichtigste Eigenschaft:} Das Skalarprodukt eines Vektors mit sich selbst ist gleich dem \textbf{Quadrat seiner euklidischen Norm}: \[ \langle \vec{x}, \vec{x} \rangle = \|\vec{x}\|_2^2 = \sum_{k=1}^N |x_k|^2 \] Dies ist die Brücke zwischen Algebra und Geometrie: Die Länge eines Vektors ist die Wurzel aus seinem ``skalaren Quadrat''. \textbf{Multiplikation einer Matrix mit einem Vektor.} Wenn $\mathbf{A} \in \C^{M \times N}$ und $\vec{x} \in \C^N$, dann wird das Ergebnis $\vec{y} = \mathbf{A}\vec{x} \in \C^M$ nach der Regel berechnet: \[ y_i = \sum_{j=1}^N A_{ij} x_j \] Das heißt, die $i$-te Komponente des Ergebnisses ist das Skalarprodukt der $i$-ten Zeile der Matrix mit dem Vektor $\vec{x}$. \textbf{Multiplikation von Matrizen.} Wenn $\mathbf{A} \in \C^{M \times K}$ und $\mathbf{B} \in \C^{K \times N}$, dann ist das Produkt $\mathbf{C} = \mathbf{A}\mathbf{B} \in \C^{M \times N}$: \[ C_{ij} = \sum_{k=1}^K A_{ik} B_{kj} \] Wichtig: Die Matrixmultiplikation ist im Allgemeinen \textbf{nicht kommutativ} --- $\mathbf{A}\mathbf{B} \neq \mathbf{B}\mathbf{A}$. Das ist einer der wichtigsten Unterschiede zur Multiplikation gewöhnlicher Zahlen! Alle diese Operationen --- Skalarprodukt, Matrix-Vektor-Multiplikation, Matrix-Matrix-Multiplikation --- haben ihre vollständigen Analoga in der Welt der Funktionen und Integraloperatoren. Nur überall dort, wo wir eine Summe $\sum$ hatten, erscheint ein Integral $\int$, und die Transposition wird durch einen komplizierteren Konjugationsoperator ersetzt. Aber das ist bereits eine Frage der Technik. \subsection{Die vier Hauptprobleme der linearen Algebra} Tatsächlich dreht sich, solange unsere Operatoren linear sind (das heißt, sie können in der Form (\ref{IntegralOperator}) oder $\mathbf{A}\vec{x}$ geschrieben werden), die ganze Welt der möglichen Probleme um einige typische Fragen. Listen wir sie auf. \textbf{Problem 1. Ein Skalarprodukt berechnen.} Gegeben zwei Vektoren --- finde die Zahl $\langle \vec{x}, \vec{y} \rangle$. Dies ist die grundlegende Operation, aus der wie aus Bausteinen alles andere aufgebaut wird. \textbf{Problem 2. Die Wirkung eines Operators berechnen.} Gegeben $\mathbf{A}$ und $\vec{x}$ --- finde $\mathbf{A}\vec{x}$. Oder gegeben $\mathbf{A}$ und $\mathbf{B}$ --- finde $\mathbf{A}\mathbf{B}$. Dies ist eine direkte Berechnung. \textbf{Problem 3. Ein lineares Gleichungssystem lösen.} Gegeben $\mathbf{A}$ und $\vec{b}$ --- finde $\vec{x}$ so, dass $\mathbf{A}\vec{x} = \vec{b}$. Dies ist das \textit{direkte} Problem: Der Operator und das Ergebnis sind bekannt, man muss finden, worauf er gewirkt hat. \textbf{Problem 4. Das Residuum minimieren (Methode der kleinsten Quadrate).} Aber was, wenn das System $\mathbf{A}\vec{x} = \vec{b}$ keine exakte Lösung hat? (Zum Beispiel gibt es mehr Gleichungen als Unbekannte --- ein überbestimmtes System, was für die Verarbeitung experimenteller Daten typisch ist.) Dann suchen wir ein $\vec{x}$, das das Residuum \textit{minimiert}: \[ \min_{\vec{x}} \|\mathbf{A}\vec{x} - \vec{b}\|_p \] Fast immer wird $p = 2$ als Norm verwendet --- dies ist die berühmte \textbf{Methode der kleinsten Quadrate (MKQ)}. Sie hat eine tiefe geometrische Interpretation: Wir projizieren den Vektor $\vec{b}$ auf den von den Spalten der Matrix $\mathbf{A}$ aufgespannten Unterraum. \subsection{Und was, wenn die rechte Seite null ist?} Und nun --- eine raffinierte Wendung. Stellen wir uns vor, dass $\vec{b} = \vec{0}$. Dann hat das Problem $\min_{\vec{x}} \|\mathbf{A}\vec{x}\|_2$ die triviale Lösung $\vec{x} = \vec{0}$. Aber das ist uninteressant! Stellen wir eine zusätzliche Bedingung: Wir suchen eine \textbf{von Null verschiedene} Lösung, zum Beispiel mit der Nebenbedingung $\|\vec{x}\|_2 = 1$. Das heißt, wir suchen die Richtung, in der die Matrix $\mathbf{A}$ den Raum am stärksten ``staucht'' (oder umgekehrt am stärksten streckt --- das ist bereits eine Frage des Vorzeichens). Und hier betreten \textbf{Eigenwerte und Eigenvektoren} die Bühne. Wir suchen solche speziellen Vektoren $\vec{v}$ und Zahlen $\lambda$, für die die Wirkung der Matrix sich einfach auf eine Streckung reduziert: \[ \mathbf{A}\vec{v} = \lambda \vec{v} \] Geometrisch: Ein Eigenvektor ist eine Richtung, die die Matrix \textit{nicht dreht}, sondern nur um den Faktor $\lambda$ streckt oder staucht. Die Aufgabe der Minimierung von $\|\mathbf{A}\vec{x}\|_2$ unter der Nebenbedingung $\|\vec{x}\|_2 = 1$ wird genau über die Eigenwerte von $\mathbf{A}^H \mathbf{A}$ oder die Singulärwerte von $\mathbf{A}$ gelöst: Das Minimum wird beim Eigenvektor der Matrix $\mathbf{A}^H\mathbf{A}$ erreicht, der dem \textit{betragsmäßig kleinsten} Eigenwert entspricht. Und das Maximum --- beim Vektor, der dem größten entspricht. \subsection{Und wenn der Operator selbst unbekannt ist?} Du wirst fragen: ``Und was, wenn uns der Operator unbekannt ist und wir nur die Eingabefunktionen und die Ergebnisse kennen?'' Hier ist die Antwort etwas mehrdeutig. Denk nach: In Matrixform haben wir nur $N$ Eingabeparameter (den Vektor $\vec{x}$), und wir wollen $N^2$ Parameter der Matrix $\mathbf{A}$ finden. Die Natur ist selten so eingerichtet, dass eine kleine Zahl von Eingaben eine größere Zahl von Parametern gut bestimmt. Dies ist ein klassisches \textbf{unterbestimmtes} Problem. Aber auch hier gibt es eine Lösung! Wenn wir sagen, dass dieselbe Matrix $\mathbf{A}$ auf mehrere verschiedene Vektoren $\vec{x}_1, \dots, \vec{x}_K$ mit bekannten Ergebnissen $\vec{b}_1, \dots, \vec{b}_K$ wirkt, dann können wir das Problem so formulieren: \begin{equation} \forall k = 1, \dots, K: \quad \mathbf{A}\vec{x}_k = \vec{b}_k, \quad \text{wobei } \mathbf{A} \text{ unbekannt ist.} \label{MatApprox} \end{equation} Wenn wir die Vektoren zu Matrizen $\mathbf{X} = (\vec{x}_1, \dots, \vec{x}_K)$ und $\mathbf{B} = (\vec{b}_1, \dots, \vec{b}_K)$ zusammenfassen, wird das Problem umgeschrieben als: \[ \mathbf{A}\mathbf{X} \approx \mathbf{B} \] Und dies ist, Achtung, \textbf{dasselbe Minimierungsproblem}, nur minimieren wir jetzt nicht über $\vec{x}$, sondern über $\mathbf{A}$: \[ \min_{\mathbf{A}} \|\mathbf{A}\mathbf{X} - \mathbf{B}\|_F \] (Hier ist $\|\cdot\|_F$ die Frobenius-Norm für Matrizen, die wir im vorigen Kapitel besprochen haben.) Transponieren wir zur Bequemlichkeit: $\mathbf{X}^T \mathbf{A}^T \approx \mathbf{B}^T$. Und wir haben wieder ein Problem der Form ``minimiere $\|\mathbf{M}\vec{z} - \vec{c}\|_2$'' erhalten, nur für jede Spalte von $\mathbf{A}^T$ separat. Dies ist die klassische \textbf{lineare Regression} --- die Grundlage der Chemometrie, QSAR, Gerätekalibrierung und vielem mehr. \subsection{Singulärwertzerlegung: eine Brücke in große Welten} Und hier betritt eine der schönsten Konstruktionen der gesamten linearen Algebra die Bühne --- die \textbf{Singulärwertzerlegung (SVD, Singular Value Decomposition)}. Jede Matrix $\mathbf{A} \in \C^{M \times N}$ kann dargestellt werden als: \[ \mathbf{A} = \mathbf{U} \mathbf{\Sigma} \mathbf{V}^H \] wobei $\mathbf{U} \in \C^{M \times M}$ und $\mathbf{V} \in \C^{N \times N}$ unitäre Matrizen sind (ihre Spalten sind orthonormierte Vektoren) und $\mathbf{\Sigma}$ eine Diagonalmatrix mit nichtnegativen Zahlen $\sigma_1 \ge \sigma_2 \ge \dots \ge 0$ auf der Diagonalen ist. Diese Zahlen heißen \textbf{Singulärwerte}. Die geometrische Bedeutung der SVD ist atemberaubend: Jede Matrix ist einfach eine Drehung ($\mathbf{V}^H$), dann eine Streckung entlang der Achsen ($\mathbf{\Sigma}$), dann noch eine Drehung ($\mathbf{U}$). Das ist alles! Mehr können Matrizen nicht. SVD ist ein universelles Werkzeug. Durch sie werden gelöst: \begin{itemize} \item Probleme der kleinsten Quadrate, \item die Suche nach der Pseudoinversen, \item Datenkompression und Extraktion von Hauptkomponenten (PCA --- Principal Component Analysis, die Grundlage der Chemometrie!), \item Regularisierung schlecht konditionierter Probleme. \end{itemize} Hier gibt es einen wichtigen Punkt: Unitäre Matrizen haben eine sehr wichtige Eigenschaft, nämlich dass immer $\mathbf{V}^H \mathbf{V} = \mathbf{I}$ gilt, wobei $\mathbf{I}$ die Einheitsmatrix ist, also eine, die auf der Hauptdiagonalen Einsen hat und sonst mit Nullen gefüllt ist. \subsection{Eine Brücke zur Funktionalanalysis: Fredholm-Theorie} Und nun --- das Interessanteste. Erinnerst du dich an den Integraloperator (\ref{MatApprox})? Nun, auch für ihn gibt es ein Analogon der SVD! Es heißt \textbf{Singulärwertzerlegung eines kompakten Operators} oder, in klassischerer Formulierung, \textbf{Fredholm-Theorie}. Der Kern ist folgender: Der Kern des Integraloperators $K(x,y)$ kann in eine Reihe nach ``Singulärfunktionen'' $u_k(x)$ und $v_k(y)$ entwickelt werden: \[ K(x,y) = \sum_{k=1}^\infty \sigma_k \, u_k(x) \overline{v_k(y)} \] wobei $\sigma_k$ die Singulärwerte sind (gegen Null abnehmend) und $u_k, v_k$ orthonormierte Funktionensysteme sind. Dies ist genau das Analogon der SVD für Matrizen, nur im unendlich-dimensionalen Raum! Und alle Ideen, die wir an Matrizen verstanden haben, übertragen sich hierher fast wörtlich. Aber überladen wir unser Bewusstsein jetzt nicht mit Fredholm und Funktionalanalysis. Die gute Nachricht ist, dass \textbf{der größte Teil dessen, was wir in Chemie und Signalverarbeitung brauchen, auf Matrixebene erledigt werden kann}. Und wo es nicht mehr geht (zum Beispiel in der Quantenmechanik oder in der strengen Theorie der Integralgleichungen), erinnern wir uns einfach, dass ``irgendwo da draußen Fredholm ist'', und fragen bei Bedarf die KI nach Details --- sie wird uns gerne über die Fredholm-Theorie, über Hilbert--Schmidt-Kerne und über die Spektren kompakter Operatoren erzählen. Die Hauptsache ist, die \textit{Struktur} zu verstehen. Und die Struktur ist überall dieselbe: Operator, Raum, Wirkung, Entwicklung nach ``Basisrichtungen''. %%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Kondition: wenn eine Matrix uns ``täuschen will''} \subsection{Wenn eine Lösung existiert, aber es scheint, als gäbe es sie nicht} Also, wir haben eine quadratische Matrix $\mathbf{A}$ und ein lineares System $\mathbf{A}\vec{x} = \vec{b}$. Was können wir darüber sagen? Wir erinnern uns aus dem Schulkurs, dass ein Gleichungssystem manchmal so ist, dass, wenn wir mit dem Einsetzen beginnen, entweder keine Lösungen existieren oder wir am Ende zwei oder mehr identische Gleichungen erhalten --- und dann gibt es unendlich viele Lösungen. In der Chemie, und überhaupt in jeder industriellen Aufgabe, kommen die Ausgangsdaten für Matrixgleichungen oft aus Messungen mit Geräten. Und hier kann uns gleich auf drei Arten Pech widerfahren: \begin{enumerate} \item Unser System wird entartet sein (keine eindeutige Lösung). \item Unser System wird viele Lösungen haben. \item \textbf{Der heimtückischste Fall:} Wir befinden uns beim Lösen in Maschinenarithmetik \textit{sehr nahe} am schlechten Fall, bemerken es aber nicht. Und wir erhalten eine Antwort, die vernünftig aussieht, aber tatsächlich völliger Unsinn ist. \end{enumerate} Wann kann das passieren? Klären wir das an einem lebendigen Beispiel. \subsection{Ein Beispiel aus der Chromatographie: die Falle des Großen und des Kleinen} Angenommen, wir haben zwei Chromatogramme, in denen zwei Substanzen vor dem Hintergrund eines Lösungsmittels aufgenommen wurden. Der Peak des Lösungsmittels kam sehr breit heraus und ``kroch'' auf die Peaks der gesuchten Substanzen, während die Antworten der gesuchten Substanzen sehr schwach ausfielen. Tatsächlich haben wir im Peak der gesuchten Substanzen: \begin{itemize} \item Erstes Chromatogramm: $(a + b_1)$, wobei $a$ eine riesige Antwort des Lösungsmittels ist und $b_1$ eine sehr kleine Zahl von der gesuchten Substanz. \item Zweites Chromatogramm: $(a + b_2)$, aber da die Temperatur etwas höher wurde, ist diese Messung leicht ungenau, sagen wir, sie erhöhte sich um einen gewissen Wert $\varepsilon$ \textit{relativ zum Peak des Lösungsmittels}, also $(1+\varepsilon)(a+b_2)$. \end{itemize} Wenn wir das zweite vom ersten subtrahieren wollen, um $b_1 - b_2$ zu berechnen, erhalten wir stattdessen: \[ (a+b_1) - (1+\varepsilon)(a+b_2) = b_1 - b_2 - \varepsilon(a+b_2) \approx b_1 - b_2 - \varepsilon a \] Das heißt, wenn $\varepsilon a$ größenordnungsmäßig mit $b_1 - b_2$ vergleichbar ist, können wir nicht nur einen großen Fehler erhalten, sondern auch das \textbf{entgegengesetzte Vorzeichen}! Dies ist das Wesen der schlechten Kondition: wenn in einer Zahl gleichzeitig eine sehr große und eine sehr kleine Komponente ``leben'' und wir versuchen, die kleine durch Subtraktion zu extrahieren. \begin{warningbox}[Regel für die Subtraktion naher Zahlen] Subtrahiere niemals zwei betragsmäßig nahe Zahlen, wenn du relative Genauigkeit brauchst. Der Verlust signifikanter Stellen ist kein Fehler des Computers, sondern eine fundamentale Eigenschaft der Arithmetik. \end{warningbox} \subsection{Wie die SVD den Mechanismus der Katastrophe aufdeckt} Dasselbe geschieht beim Lösen von Gleichungssystemen. Und wir haben eine schöne Möglichkeit, \textit{im Voraus} abzuschätzen, was zu tun ist. Angenommen, wir haben ein System $\mathbf{A}\vec{x} = \vec{b}$, wobei wir $\vec{x}$ suchen und alles andere gegeben ist. Erinnern wir uns an die Singulärwertzerlegung $\mathbf{A} = \mathbf{U} \mathbf{\Sigma} \mathbf{V}^H$. Dann kann die Lösung in mehreren Schritten umgeschrieben werden: \begin{align*} \mathbf{U} \mathbf{\Sigma} \mathbf{V}^H \vec{x} &= \vec{b} \\ \mathbf{\Sigma} \mathbf{V}^H \vec{x} &= \mathbf{U}^H \vec{b} \equiv \vec{t}_1 \\ \mathbf{V}^H \vec{x} &= \mathbf{\Sigma}^{-1} \vec{t}_1 \equiv \vec{t}_2 \\ \vec{x} &= \mathbf{V} \vec{t}_2 \end{align*} Hier nutzen wir die Eigenschaft, dass das Lösen eines Systems mit unitären Matrizen $\mathbf{U}, \mathbf{V}$ sehr einfach ist --- es genügt, mit $\mathbf{U}^H$ oder $\mathbf{V}^H$ zu multiplizieren, da $\mathbf{U}^H \mathbf{U} = \mathbf{I}$. Und ein System mit einer Diagonalmatrix zu lösen ist ohnehin trivial. Zeichnen wir, wie das für eine $3 \times 3$-Matrix aussieht: \[ \begin{pmatrix} \sigma_1 & 0 & 0 \\ 0 & \sigma_2 & 0 \\ 0 & 0 & \sigma_3 \end{pmatrix} \begin{pmatrix} t_{2,1} \\ t_{2,2} \\ t_{2,3} \end{pmatrix} = \begin{pmatrix} t_{1,1} \\ t_{1,2} \\ t_{1,3} \end{pmatrix} \quad \Rightarrow \quad \begin{cases} \sigma_1 t_{2,1} = t_{1,1} \\ \sigma_2 t_{2,2} = t_{1,2} \\ \sigma_3 t_{2,3} = t_{1,3} \end{cases} \quad \Rightarrow \quad t_{2,k} = \frac{t_{1,k}}{\sigma_k} \] Siehst du? Jede Gleichung ist einfach eine Division einer Komponente durch den entsprechenden Singulärwert. Und hier --- das Wichtigste. Stellen wir uns vor, wir haben die Singulärwerte der Ausgangsmatrix berechnet und bemerkt, dass der größte Singulärwert $\sigma_1$ sehr, sehr viel größer ist als der kleinste $\sigma_N$. Dann geschieht beim Schritt $\vec{t}_2 = \mathbf{\Sigma}^{-1} \vec{t}_1$ Folgendes: \begin{itemize} \item Im Vektor $\vec{t}_1$ ist der Messfehler mehr oder weniger gleichmäßig über alle Komponenten verteilt. \item Nach der Multiplikation mit $\mathbf{\Sigma}^{-1}$ wächst die Komponente, die dem \textit{kleinen} $\sigma_N$ entspricht, um den Faktor $\sigma_1/\sigma_N$! \item Und wenn $\sigma_N$ sogar null ist? Dann teilen wir durch null --- und die Lösung existiert nicht. \end{itemize} Genau deshalb hat man beschlossen, das Verhältnis $\sigma_1 / \sigma_N$ als \textbf{Konditionszahl} der Matrix zu bezeichnen: \[ \cond(\mathbf{A}) = \frac{\sigma_{\max}}{\sigma_{\min}} \] Tatsächlich ist dies eine Charakteristik dafür, \textit{wie sehr wir die Lösung mit einer solchen Matrix verderben können}. Wenn $\cond(\mathbf{A}) \approx 10^k$, dann verlieren wir beim Lösen des Systems etwa $k$ signifikante Stellen an Genauigkeit. Für double precision (16 signifikante Stellen) bedeutet dies, dass bei $\cond(\mathbf{A}) \approx 10^{16}$ die Antwort vollständig aus Rauschen besteht. \begin{tipbox}[Wichtige Bemerkung] Die Konditionszahl ist eine Charakteristik der \textit{Matrix}, nicht des Algorithmus. Kein noch so genialer Algorithmus kann ein schlecht konditioniertes System genauer lösen, als $\cond(\mathbf{A})$ erlaubt. Dies ist eine fundamentale Beschränkung des Problems, nicht unserer Berechnungsmethoden. Gleichzeitig beeinflusst die Kondition \textit{nicht} die Berechnung der Eigenvektoren symmetrischer/hermitescher Matrizen oder die Suche nach Singulärvektoren --- diese Probleme sind in der Regel viel stabiler. \end{tipbox} \subsection{Wo kommt das alles in der Chemie und darüber hinaus vor?} Wir hatten viel trockene Theorie, und es wäre gut, sich daran zu erinnern, wo das alles in der Chemie angewendet wird. \textbf{Quantenchemie: die Schrödinger-Gleichung.} Ja, die Schrödinger-Wellengleichung und alle ihre Näherungen --- die Hartree--Fock-Gleichungen, Density Functional Theory (DFT) --- all dies sind Eigenwertprobleme: die Suche nach Vektoren $\vec{x}$ und Werten $\lambda$, die die Gleichung $\mathbf{A}\vec{x} = \lambda \vec{x}$ erfüllen. Die Matrizen hier --- Fock-Matrizen oder Hamilton-Matrizen --- haben oft Dimensionen von Tausenden und Zehntausenden, und ihre Kondition bestimmt direkt, ob wir überhaupt eine physikalisch sinnvolle Antwort erhalten können. \textbf{Massenspektrometrie und Spektroskopie.} Die Suche nach einem aufgenommenen Spektrum in einer Datenbank von Antworten in einem Massenspektrometer --- das sind Algorithmen, die auf der Lösung linearer Systeme basieren. Wenn du eine Mischung von Substanzen hast und verstehen willst, woraus sie besteht, löst du das System $\mathbf{A}\vec{x} = \vec{b}$, wobei $\mathbf{A}$ die Matrix der Referenzspektren ist, $\vec{b}$ das gemessene Spektrum der Mischung und $\vec{x}$ die Konzentrationen der Komponenten. Und wenn die Spektren der Komponenten ähnlich sind --- ist die Matrix schlecht konditioniert, und die Konzentrationen werden ``springen''. \textbf{Kryo-Elektronenmikroskopie.} Die Rekonstruktion der 3D-Struktur von Proteinen aus Tausenden von 2D-Projektionen ist ein gigantisches dünnbesetztes Gleichungssystem. Die Kondition ist dort ein Schlüsselfaktor, der die Auflösung der resultierenden Struktur bestimmt. \textbf{NMR-Spektroskopie.} Die Entfaltung (Dekonvolution) ist ein klassisches schlecht konditioniertes Problem. Genau deshalb verwendet man in der NMR die Tikhonov-Regularisierung und die Fourier-Transformation --- sie verwandeln ein schlecht konditioniertes Problem in ein diagonales. \textbf{Kalibrierung analytischer Geräte.} Die PLS-Regression (Partial Least Squares) ist im Wesentlichen eine SVD mit Regularisierung. Die Kondition der Matrix der Spektren beeinflusst direkt die Genauigkeit der Konzentrationsvorhersage. \textbf{Internetsuche (PageRank).} Sogar die allererste Google-Suche war so eingerichtet, dass alle zuvor gefundenen Dokumente in eine riesige Matrix eingeordnet wurden und für sie die Singulärwertzerlegung berechnet wurde (oder, was äquivalent ist, der Haupt-Eigenvektor der Übergangsmatrix gesucht wurde), was half, das gesuchte Dokument sofort anhand einer passenden Wortkombination zu finden. In der modernen KI wird SVD sehr, sehr oft verwendet --- von der Kompression der Gewichte neuronaler Netze bis zur Dimensionsreduktion in Empfehlungssystemen. \textbf{Signal- und Bildverarbeitung.} Die Dekonvolution verschwommener Bilder in der Mikroskopie, die Rauschunterdrückung in EKG/EEG, die Kompression von Audio (MP3) und Video --- all dies sind Probleme, bei denen die Kondition eine Schlüsselrolle spielt. Es gab eine Zeit, da waren diese drei Probleme --- das Lösen linearer Systeme, die Suche nach Eigenwerten und die SVD --- ein Stolperstein bei der Lösung großer industrieller Aufgaben. Es gab sogar Firmen, die \textit{nur} solche Solver entwickelten --- und sonst nichts. Allerdings war das noch in den 1990er Jahren. \subsection{Was kostet das? Rechenkomplexität} Grob gesagt, wenn wir eine quadratische $N \times N$-Matrix haben und entweder ein lineares System mit ihr lösen oder ihre Eigen- oder Singulärwerte und -vektoren suchen, dann müssen wir etwa $\mathcal{O}(N^3)$ arithmetische Operationen aufwenden, wenn wir die Struktur oder Eigenschaften dieser Matrix nicht speziell nutzen. Dies ist der ``Preis'' des vollständigen Problems. Aber wenn die Matrix spezielle Eigenschaften hat, kann der Preis stark gesenkt werden --- manchmal bis auf $\mathcal{O}(N)$ oder $\mathcal{O}(N \log N)$. Und darüber --- weiter unten. % === VERBINDUNG SVD UND EIGENWERTE === \subsection{Der Zusammenhang zwischen Singulärwertzerlegung und Eigenwerten} Und nun --- eine der schönsten Tatsachen der gesamten linearen Algebra. Angenommen, wir haben die Singulärwertzerlegung: \[ \mathbf{A} = \mathbf{U} \mathbf{\Sigma} \mathbf{V}^H \] Betrachten wir, was passiert, wenn wir $\mathbf{A}$ mit ihrer hermitesch konjugierten $\mathbf{A}^H$ multiplizieren: \begin{align*} \mathbf{A} \mathbf{A}^H &= (\mathbf{U} \mathbf{\Sigma} \mathbf{V}^H)(\mathbf{V} \mathbf{\Sigma} \mathbf{U}^H) \\ &= \mathbf{U} \mathbf{\Sigma} (\mathbf{V}^H \mathbf{V}) \mathbf{\Sigma} \mathbf{U}^H \\ &= \mathbf{U} \mathbf{\Sigma}^2 \mathbf{U}^H \end{align*} Hier haben wir verwendet, dass $\mathbf{V}^H \mathbf{V} = \mathbf{I}$ (Unitärität von $\mathbf{V}$) und $\mathbf{\Sigma}^T = \mathbf{\Sigma}$ (Diagonalmatrix). Analog: \begin{align*} \mathbf{A}^H \mathbf{A} &= (\mathbf{V} \mathbf{\Sigma} \mathbf{U}^H)(\mathbf{U} \mathbf{\Sigma} \mathbf{V}^H) \\ &= \mathbf{V} \mathbf{\Sigma} (\mathbf{U}^H \mathbf{U}) \mathbf{\Sigma} \mathbf{V}^H \\ &= \mathbf{V} \mathbf{\Sigma}^2 \mathbf{V}^H \end{align*} \textbf{Was bedeutet das?} \begin{itemize} \item Die Singulärvektoren der Matrix $\mathbf{A}$ sind \textit{genau} die Eigenvektoren der Matrizen $\mathbf{A}\mathbf{A}^H$ und $\mathbf{A}^H\mathbf{A}$. \item Die Singulärwerte $\sigma_k$ der Matrix $\mathbf{A}$ sind die Quadratwurzeln aus den Eigenwerten $\lambda_k$ der Matrizen $\mathbf{A}\mathbf{A}^H$ oder $\mathbf{A}^H\mathbf{A}$: \[ \sigma_k = \sqrt{\lambda_k} \] \end{itemize} Dies ist ein tiefer Zusammenhang zwischen zwei scheinbar verschiedenen Problemen: der Singulärwertzerlegung und der Eigenwertsuche. \begin{warningbox}[Vorsicht: das Quadrat der Kondition!] Aber es gibt ein heimtückisches Detail. Die Konditionszahl der Matrizen $\mathbf{A}\mathbf{A}^H$ und $\mathbf{A}^H\mathbf{A}$ ist das \textbf{Quadrat} der Konditionszahl der Ausgangsmatrix: \[ \cond(\mathbf{A}\mathbf{A}^H) = \cond(\mathbf{A}^H\mathbf{A}) = \cond(\mathbf{A})^2 \] Wenn $\cond(\mathbf{A}) = 10^8$, dann ist $\cond(\mathbf{A}^H\mathbf{A}) = 10^{16}$ --- und wir verlieren alle 16 signifikanten Stellen an Genauigkeit! Deshalb muss der Übergang vom Problem der Singulärwertzerlegung zum Eigenwertproblem für $\mathbf{A}^H\mathbf{A}$ \textbf{äußerst vorsichtig} erfolgen. In modernen Bibliotheken (LAPACK) wird die SVD direkt berechnet, ohne explizite Bildung von $\mathbf{A}^H\mathbf{A}$ --- genau um diese Katastrophe zu vermeiden. \end{warningbox} %%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Klassifikation von Matrizen: Wer ist wer} Nun, da klar ist, dass praktisch alles auf diesen ``einfachen'' mathematischen Problemen basiert, systematisieren wir ein wenig unser Wissen über solche Lösungen. Denn wir haben immer eine Matrix, und von den Eigenschaften dieser Matrix hängt sehr viel ab. \subsection{Nach der Form} \begin{tabularx}{\textwidth}{l l X} \toprule \textbf{Typ} & \textbf{Größe} & \textbf{Beschreibung} \\ \midrule Quadratisch & $N \times N$ & Gleiche Anzahl von Zeilen und Spalten. Nur für solche Matrizen sind Determinante, Eigenwerte und inverse Matrix definiert. \\ ``Stehend'' (hoch) & $M \times N, \; M > N$ & Mehr Gleichungen als Unbekannte. Hat meist keine exakte Lösung --- wir lösen über MKQ. \\ ``Liegend'' (breit) & $M \times N, \; M < N$ & Mehr Unbekannte als Gleichungen. Es gibt unendlich viele Lösungen --- wir suchen die mit minimaler Norm. \\ \bottomrule \end{tabularx} \subsection{Nach Symmetrie und Struktur der Elemente} \begin{tabularx}{\textwidth}{l X l} \toprule \textbf{Typ} & \textbf{Definition und Eigenschaften} & \textbf{Komplexität} \\ \midrule Symmetrisch reell & $\mathbf{A} = \mathbf{A}^T$, Elemente $\in \R$. Eigenwerte \textbf{reell}, Eigenvektoren --- orthogonal. & $\mathcal{O}(N^3/3)$ \\ Hermitesch komplex & $\mathbf{A} = \mathbf{A}^H$, Elemente $\in \C$. Eigenwerte \textbf{reell}, Eigenvektoren --- orthonormiert. & $\mathcal{O}(4N^3/3)$ \\ Schiefsymmetrisch & $\mathbf{A} = -\mathbf{A}^T$. Eigenwerte --- rein imaginär oder null. & $\mathcal{O}(N^3/3)$ \\ Orthogonal / Unitär & $\mathbf{A}^T\mathbf{A} = \mathbf{I}$ / $\mathbf{A}^H\mathbf{A} = \mathbf{I}$. Erhält Längen und Winkel. Eigenwerte betragsmäßig gleich 1. & $\mathcal{O}(N^2)$ für Multiplikation \\ Einheitsmatrix & $\mathbf{A} = \mathbf{I}$. Diagonale aus Einsen, Rest --- Nullen. & $\mathcal{O}(N)$ \\ Diagonal & $A_{ij} = 0$ für $i \neq j$. Wird als Vektor aus $N$ Zahlen gespeichert. & $\mathcal{O}(N)$ \\ \bottomrule \end{tabularx} \subsection{Nach positiver Definitheit} Dies ist eine Unterklasse symmetrischer/hermitescher Matrizen, und sie ist für Optimierung und Statistik von entscheidender Bedeutung. \begin{tabularx}{\textwidth}{l X} \toprule \textbf{Typ} & \textbf{Definition und Eigenschaften} \\ \midrule Positiv definit ($\mathbf{A} \succ 0$) & $\forall \vec{x} \neq 0: \; \vec{x}^H \mathbf{A} \vec{x} > 0$. Alle Eigenwerte streng positiv. Kondition --- Verhältnis des größten zum kleinsten Eigenwert. Wird mit dem Cholesky-Verfahren in $\mathcal{O}(N^3/3)$ gelöst. \\ Positiv semidefinit ($\mathbf{A} \succeq 0$) & $\forall \vec{x}: \; \vec{x}^H \mathbf{A} \vec{x} \ge 0$. Alle Eigenwerte $\ge 0$. Kann entartet sein. \\ \bottomrule \end{tabularx} \subsection{Nach der Struktur der Anordnung der Nullelemente} \begin{tabularx}{\textwidth}{l X l} \toprule \textbf{Typ} & \textbf{Definition und Eigenschaften} & \textbf{Komplexität} \\ \midrule Bandmatrix & Nullelemente nur in der Nähe der Hauptdiagonalen: $A_{ij} = 0$ für $|i-j| > k$. Wird als $N \times (2k+1)$ gespeichert. & $\mathcal{O}(N k^2)$ \\ Toeplitz-Matrix & $A_{ij}$ hängt nur von $i-j$ ab. Jede Diagonale ist eine Konstante. Tritt bei Problemen mit stationären Prozessen auf. & $\mathcal{O}(N^2)$, über FFT --- $\mathcal{O}(N \log N)$ \\ Zirkulante Matrix & Toeplitz + Periodizität: Jede Zeile ist eine zyklische Verschiebung der vorherigen. Wird durch die Fourier-Transformation diagonalisiert. & $\mathcal{O}(N \log N)$ über FFT \\ Blockmatrix & Besteht aus Blöcken, von denen jeder eine Matrix ist. Erlaubt rekursive Algorithmen. & Hängt von der Blockstruktur ab \\ Dünnbesetzte Matrix (sparse) & Die meisten Elemente sind Nullen. Wird in speziellen Formaten gespeichert (CSR, CSC). Grundlage großer Berechnungen. & $\mathcal{O}(\text{nnz})$, wobei nnz die Anzahl der Nichtnullen ist \\ \bottomrule \end{tabularx} \subsection{Entartung} Eine \textbf{entartete (singuläre) Matrix} ist eine Matrix, deren Determinante null ist, oder, was äquivalent ist, die mindestens einen Singulärwert gleich null hat. Für eine solche Matrix: \begin{itemize} \item existiert keine inverse Matrix, \item hat das System $\mathbf{A}\vec{x} = \vec{b}$ entweder keine Lösungen oder unendlich viele, \item ist $\rank(\mathbf{A}) < N$. \end{itemize} In der Praxis sind Matrizen fast nie \textit{exakt} entartet --- aber sie können \textit{fast} entartet sein, also mit einem sehr kleinen $\sigma_{\min}$. Und dies ist ein viel heimtückischerer Fall, weil formal eine Lösung existiert, sie aber vollständig durch das Rauschen in den Daten bestimmt wird. % === ERGÄNZUNG ZUR KLASSIFIKATION DER MATRIZEN === \subsection{Zusätzliche wichtige Matrizentypen} \begin{tabularx}{\textwidth}{l X l} \toprule \textbf{Typ} & \textbf{Definition und Eigenschaften} & \textbf{Komplexität} \\ \midrule Obere Dreiecksmatrix & $A_{ij} = 0$ für $i > j$. Alle Nichtnullen oberhalb oder auf der Hauptdiagonalen. Lösung des Systems --- Rückwärtseinsetzen. & $\mathcal{O}(N^2)$ \\ Untere Dreiecksmatrix & $A_{ij} = 0$ für $i < j$. Alle Nichtnullen unterhalb oder auf der Hauptdiagonalen. Lösung des Systems --- Vorwärtseinsetzen. & $\mathcal{O}(N^2)$ \\ Permutationsmatrix & In jeder Zeile und jeder Spalte genau eine Eins, der Rest --- Nullen. Multiplikation mit einer solchen Matrix ist eine Permutation von Zeilen/Spalten. $\mathbf{P}^T = \mathbf{P}^{-1}$. & $\mathcal{O}(N)$ \\ \bottomrule \end{tabularx} Permutationsmatrizen sind nicht nur eine Abstraktion. Sie sind entscheidend für die numerische Stabilität von Algorithmen (zum Beispiel bei der LU-Zerlegung mit Pivotierung), und wir werden ihnen noch begegnen. \subsection{Übersichtstabelle: Was und wie lösen} \begin{longtable}{l c c c} \toprule \textbf{Problem} & \textbf{Allgemeiner Fall} & \textbf{Symmetrisch} & \textbf{Dünnbesetzt} \\ \midrule Lineares System $\mathbf{A}\vec{x} = \vec{b}$ & LU, $\mathcal{O}(N^3)$ & Cholesky, $\mathcal{O}(N^3/3)$ & Iterativ, $\mathcal{O}(\text{nnz} \cdot k)$ \\ MKQ $\min \|\mathbf{A}\vec{x} - \vec{b}\|_2$ & QR, $\mathcal{O}(MN^2)$ & --- & LSQR, iterativ \\ Eigenwerte & $\mathcal{O}(N^3)$ & $\mathcal{O}(N^3/3)$, alle reell & Lanczos, $\mathcal{O}(\text{nnz} \cdot k)$ \\ SVD & $\mathcal{O}(MN^2)$ & Über Eigenwerte von $\mathbf{A}^T\mathbf{A}$ & Truncated SVD, iterativ \\ \bottomrule \end{longtable} Hier ist $k$ die Anzahl der Iterationen, die von der gewünschten Genauigkeit und der Kondition abhängt. \begin{tipbox}[Hauptschlussfolgerung] Bevor du ein Problem löst, \textit{schau dir die Matrix an}. Ihre Eigenschaften --- Symmetrie, Dünnbesetztheit, Bandstruktur --- können die Komplexität von kubisch auf linear senken. Und ihre Kondition wird dir sagen, ob es überhaupt lohnt, der erhaltenen Antwort zu vertrauen. \end{tipbox} %%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Wie lineare Algebra-Probleme gelöst werden: von der Theorie zur Praxis} \subsection{Es gibt keine einzige Lösung für alle Fälle} Wie also werden diese linearen Algebra-Probleme gelöst? Sicher gibt es eine oder vielleicht ein paar der einfachsten und zuverlässigsten Lösungen? Leider gibt es auch jetzt keine einzige Lösung für alle Lebenslagen, und es ist keine in Sicht. Der Autor dieses Skripts hat einst an der Entwicklung von mehr als hundert verschiedenen solcher Algorithmen mitgewirkt. Ja, es gibt sehr viele solcher Algorithmen, und man muss zumindest kurz, auch ohne es jemals implementiert zu haben, verstehen, wann und was man verwenden kann. \subsection{Klassifikation der Lösungsmethoden} Die Lösungsmethoden lassen sich etwa so klassifizieren: \subsubsection{Direkte Zerlegungsmethoden} \textbf{LU-Zerlegung:} $\mathbf{A} = \mathbf{L}\mathbf{U}$, wobei $\mathbf{L}$ eine untere und $\mathbf{U}$ eine obere Dreiecksmatrix ist. Mit einer solchen Faktorisierung reduziert sich die Lösung des Systems $\mathbf{A}\vec{x} = \vec{b}$ auf zwei einfache Einsetzungen: \begin{enumerate} \item Wir lösen $\mathbf{L}\vec{y} = \vec{b}$ (Vorwärtseinsetzen, $\mathcal{O}(N^2)$) \item Wir lösen $\mathbf{U}\vec{x} = \vec{y}$ (Rückwärtseinsetzen, $\mathcal{O}(N^2)$) \end{enumerate} Aber im Allgemeinen gilt $\cond(\mathbf{L}) \cdot \cond(\mathbf{U}) \ge \cond(\mathbf{A})$, und wenn man keine Permutationen durchführt, kann das Produkt der Konditionen wesentlich größer werden als $\cond(\mathbf{A})$. Das heißt, wir können unser Ausgangsproblem verderben! Deshalb verwendet man in der Praxis fast immer die \textbf{LUP-Zerlegung}: $\mathbf{P}\mathbf{A} = \mathbf{L}\mathbf{U}$, wobei $\mathbf{P}$ eine Permutationsmatrix ist. Die Permutationen werden so gewählt, dass auf der Diagonalen von $\mathbf{U}$ die größtmöglichen Elemente stehen (Pivotierung). Dies garantiert numerische Stabilität. \textbf{Für spezielle Matrizen:} \begin{itemize} \item Wenn $\mathbf{A}$ eine Bandmatrix ist, dann verlassen $\mathbf{L}$ und $\mathbf{U}$ nicht das Band. Komplexität $\mathcal{O}(N k^2)$, wobei $k$ die Bandbreite ist. \item Wenn $\mathbf{A}$ dünnbesetzt ist, können $\mathbf{L}$ und $\mathbf{U}$ wesentlich mehr Nullelemente enthalten (fill-in), und dies kann ein großes Problem sein. \item Für eine symmetrische Matrix --- $\mathbf{A} = \mathbf{L}\mathbf{D}\mathbf{L}^H$ oder $\mathbf{P}\mathbf{A}\mathbf{P}^H = \mathbf{L}\mathbf{D}\mathbf{L}^H$. \item Für eine positiv definite --- $\mathbf{A} = \mathbf{L}\mathbf{L}^H$, das sogenannte \textbf{Cholesky-Verfahren}. Aber das Problem des Wachstums der Konditionszahl existiert auch für positiv definite Matrizen. \end{itemize} \textbf{QR-Zerlegung:} $\mathbf{A} = \mathbf{Q}\mathbf{R}$, wobei $\mathbf{Q}$ unitär ($\mathbf{Q}^H\mathbf{Q} = \mathbf{I}$) und $\mathbf{R}$ eine obere Dreiecksmatrix ist. Diese Methode \textbf{garantiert keine Erhöhung der Konditionszahl} der Lösung, da unitäre Transformationen die Längen von Vektoren erhalten. Sie ist gut für dichte Matrizen ohne Struktur. Aber wenn die Ausgangsmatrix dünnbesetzt ist, dann wird $\mathbf{Q}$ fast immer dicht sein (nur für sehr spezielle Fälle und Methoden), was die Anwendbarkeit dieser Methode auf riesige Probleme stark einschränkt, wenn die Matrix selbst in den Speicher passt, ihre dichte Version $N^2$ aber nicht. \textbf{Schur-Zerlegung:} $\mathbf{A} = \mathbf{Q}\mathbf{G}\mathbf{Q}^H$, wobei $\mathbf{G}$ eine obere Dreiecksmatrix ist (für reelle Matrizen --- eine obere quasi-Dreiecksmatrix mit $2 \times 2$-Blöcken auf der Diagonalen). Sie wird zur Suche nach Eigenwerten und Eigenvektoren verwendet, die man durch Lösen des Eigenwertproblems für die bereits dreieckige Matrix $\mathbf{G}$ findet (und für eine Dreiecksmatrix sind die Eigenwerte einfach die Diagonalelemente!). \subsubsection{Iterative Methoden} Dies sind Methoden zur Lösung eines linearen Systems oder eines Eigenwertproblems, wenn wir aufgrund seiner Struktur effizient mit der Matrix multiplizieren können und die Lösung iterativ aufbauen. Es gibt hier viele Methoden, aber wichtig ist: \textbf{die Konvergenz iterativer Methoden hängt von der Konditionszahl ab}. Je größer sie ist, desto langsamer konvergieren sie. (Es gibt Ausnahmen: Wenn die Matrix nur wenige sehr große und sehr kleine Singulärwerte hat, kann die Konvergenz wesentlich schneller sein.) Alle diese Methoden lassen sich zudem in Unterklassen einteilen: \begin{enumerate} \item \textbf{Spezielle Methoden für positiv definite Matrizen} mit geringem zusätzlichem Speicheraufwand und vernünftiger Iterationszahl. Die bekannteste ist das \textbf{Verfahren der konjugierten Gradienten (Conjugate Gradients, CG)}. \item \textbf{Rückführung eines nichtsymmetrischen Problems auf ein symmetrisches:} Statt $\mathbf{A}\vec{x} = \vec{b}$ löst man $\mathbf{A}^H\mathbf{A}\vec{x} = \mathbf{A}^H\vec{b}$. Wenn die Konditionszahl $\cond(\mathbf{A}^H\mathbf{A})$ nicht so enorm im Vergleich zu $\cond(\mathbf{A})$ ist, dann ist das vernünftig. Aber denk daran: $\cond(\mathbf{A}^H\mathbf{A}) = \cond(\mathbf{A})^2$, das ist also ein zweischneidiges Schwert. \item \textbf{Vorkonditionierer (preconditioners):} spezielle Matrizen $\mathbf{P} \approx \mathbf{A}^{-1}$, sodass wir statt $\mathbf{A}\vec{x} = \vec{b}$ das System $\mathbf{P}\mathbf{A}\vec{x} = \mathbf{P}\vec{b}$ lösen, in der Annahme, dass $\cond(\mathbf{P}\mathbf{A}) \ll \cond(\mathbf{A})$, was die Lösung durch solche iterativen Methoden stark beschleunigt. Übrigens werden für große dünnbesetzte Matrizen Vorkonditionierer oft als \textbf{unvollständige LU-Zerlegung} konstruiert: Man konstruiert für die Ausgangsmatrix $\mathbf{A} \approx \mathbf{L}\mathbf{U}$, versucht aber, neue Nullelemente bei der Konstruktion von $\mathbf{L}\mathbf{U}$ ``wegzuwerfen''. In diesem Fall kann man diese Matrix als etwas Ähnliches wie die Ausgangsmatrix noch anwenden, und indem man mit ihr löst, vorkonditioniert man das Ausgangsproblem und verbessert die Konvergenz. \end{enumerate} \subsection{Zwei Welten: dichte und dünnbesetzte Probleme} Tatsächlich lassen sich die meisten Lösungsalgorithmen in zwei Klassen einteilen: \begin{enumerate} \item Wenn wir bereit sind, $\mathcal{O}(N^3)$ arithmetische Operationen aufzuwenden, und $\mathcal{O}(N^2)$ Speicher haben. \item Wenn wir uns einschränken müssen und die Matrix so riesig ist und eine spezielle Struktur hat, dass wir mit ihr multiplizieren können, sie aber in voller Form zu nehmen sehr schwierig ist. \end{enumerate} Im ersten Fall --- das sind \textbf{Zerlegungsmethoden} (LU, QR, Cholesky). Im zweiten --- das sind \textbf{iterative Methoden}. Da moderne Prozessoren und Speicher eine spezifische Struktur haben, haben iterative Methoden eine etwas schlechtere Leistung als vollständige Zerlegungsmethoden. Deshalb ist die Anwendung iterativer Methoden auf sehr kleine Matrizen praktisch nie gerechtfertigt. \subsection{Beispiele aus der Quantenchemie: DFT} \textbf{Beispiel 1: Ein kleines System.} Ein kleines System mit einigen Elektronen nach der DFT-Methode, und die Größe des Hamilton-Operators ist etwa $1000 \times 1000$. Wir passen in den Speicher (das sind nur $\sim 8$ MB für double precision). Wir können alle Eigenwerte und die Menge der uns interessierenden Eigenvektoren finden, und dafür ist die offensichtliche Wahl eine direkte Methode (zum Beispiel die Schur-Zerlegung oder der QR-Algorithmus). \textbf{Beispiel 2: Eine große chemische Struktur.} Eine ziemlich große chemische Struktur nach der DFT-Methode, in der es etwa tausend Elektronenpaare in den äußeren Orbitalen gibt. Dafür haben wir den Hamilton-Operator der Schrödinger-Gleichung aufgebaut, und wir müssen für jedes Elektronenpaar einen eigenen Eigenvektor finden. Und der Hamilton-Operator ist so gegeben, dass wir etwa hundert Basisfunktionen pro Orbital genommen haben, das heißt, die Dimension dieses Hamilton-Operators ist etwa $100\,000 \times 100\,000$. Das scheint nicht sehr viel, aber allein die Matrix eines solchen Hamilton-Operators belegt bereits etwa \textbf{80 GB} im Arbeitsspeicher, und die Lösung selbst wird noch fast ein Gigabyte belegen. Hier ist eine iterative Lösung viel naheliegender (zum Beispiel das Lanczos-Verfahren oder Davidson), die nur die wenigen benötigten Eigenvektoren findet, ohne mit der vollen Matrix zu arbeiten. \subsection{Erfinde das Rad nicht neu: Bibliotheken} Heutzutage gibt es in vielen Programmiersprachen gut und über Jahre ausgefeilte Bibliotheken zur Lösung solcher Gleichungssysteme, zur Suche nach Singulärwerten und -vektoren sowie nach Eigenvektoren und -werten. Es genügt, in der Dokumentation von \textbf{NumPy} für Python und in \textbf{LAPACK/BLAS} für C/C++ nachzuschauen --- und alles wird sofort klar. Darüber hinaus ist der Code sehr gut für moderne Computer optimiert und ziemlich komplex. Zum Beispiel enthält die moderne Version der SVD in LAPACK etwa \textbf{eine halbe Million Zeilen Code}, und sie zu wiederholen oder gar besser zu machen ist wirklich eine fast unmögliche Aufgabe. \begin{tipbox}[Praktischer Rat] Schreibe niemals eigene Solver für lineare Algebra, es sei denn, es ist eine Übungsaufgabe. Verwende: \begin{itemize} \item \textbf{Python:} NumPy (\texttt{numpy.linalg}), SciPy (\texttt{scipy.linalg}) \item \textbf{C/C++:} LAPACK, BLAS, Eigen \item \textbf{Fortran:} LAPACK (das ist sein natürliches Element) \item \textbf{MATLAB:} eingebaute Funktionen (\texttt{eig}, \texttt{svd}, \texttt{lu}, \texttt{qr}) \end{itemize} Diese Bibliotheken haben Jahrzehnte der Optimierung und des Testens hinter sich. Sie wissen über den Prozessor-Cache, über SIMD-Instruktionen, über Multithreading --- all das, was du nicht von Hand wiederholen kannst. \end{tipbox} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Nichtlineare Operatoren: wenn die Welt aufhört, gerade zu sein} \subsection{Was ist ein nichtlinearer Operator?} Bisher haben wir über lineare Operatoren gesprochen --- solche, die als $\mathbf{A}\vec{x}$ oder $\int K(x,y)g(y)\diff y$ geschrieben werden können. Aber die reale Welt ist nichtlinear. Ein \textbf{nichtlinearer Operator} ist eine Abbildung $\mathcal{F}$, die auf eine Funktion oder einen Vektor $\vec{x}$ wirkt, aber die Eigenschaften der Additivität und Homogenität \textit{nicht} besitzt: \[ \mathcal{F}(\alpha \vec{x} + \beta \vec{y}) \neq \alpha \mathcal{F}(\vec{x}) + \beta \mathcal{F}(\vec{y}) \] Grob gesagt, definieren wir nun, dass wir eine Funktion $f(\vec{x})$ haben, die von einem oder mehreren Parametern abhängt. Und wir wollen entweder: \begin{itemize} \item Die Gleichung $f(\vec{x}) = 0$ lösen (ein System nichtlinearer Gleichungen), \item Das Minimum $\min_{\vec{x}} f(\vec{x})$ finden (ein Optimierungsproblem). \end{itemize} \subsection{Ein Beispiel aus der Chromatographie: das Teller-Modell} Ein klassisches Beispiel ist das System von Bilanzgleichungen und Dissoziationskonstanten auf jedem chromatographischen Teller. Wir wollen einen Chromatographen modellieren, genauer gesagt seine Säule. In ihr findet eine Trennung statt, und die mobile Phase erzeugt beim Wandern entlang der Säule bei jedem Schritt faktisch ein solches Gleichgewicht. Du wirst sagen --- oh, hier gibt es doch nur etwa 4--5 Gleichungen und ebenso viele Unbekannte! Ja, ich stimme zu, nicht viele. Aber: \begin{enumerate} \item Es gibt viele Teller (sagen wir tausend) --- das ist eins. \item Wir zerlegen die Durchgangszeit der Substanz durch die Säule in viele Zeitschritte --- das ist zwei. \end{enumerate} Das heißt, bei jedem Zeitschritt haben wir tausendmal dasselbe Gleichungssystem mit unterschiedlichen Parametern der Eingangskonzentrationen. Und es wird wesentlich mehr solcher Zeitschritte geben als Teller, denn die Substanz wird ja von der Säule zurückgehalten. Das heißt, wir müssen so etwa \textbf{zehn Millionen Mal} ein System aus nur 5 Gleichungen mit 5 Unbekannten lösen. Aber wenn wir es auch nur einmal nicht lösen --- können wir keine weitere Lösung erhalten. Leider ist ein solches System nicht linear. Die Bilanzgleichungen sind linear, aber die Gleichungen, die die Dissoziationskonstanten verknüpfen, enthalten Produkte und Verhältnisse der Unbekannten zueinander. Und, überraschenderweise, kann ein solches System manchmal wirklich sehr schlecht lösbar sein. \subsection{Klassifikation nichtlinearer Probleme} Bevor wir über Methoden sprechen, klassifizieren wir die Probleme selbst: \begin{tabularx}{\textwidth}{l X} \toprule \textbf{Eigenschaft} & \textbf{Beschreibung} \\ \midrule \textbf{Glattheit:} & \\ \quad Glatt & Die Funktion hat stetige Ableitungen (mindestens die erste, besser auch die zweite). Beispiel: $f(x) = x^2 + \sin x$. \\ \quad Fast glatt & Die Funktion ist fast überall glatt, aber es gibt Knickstellen oder Unstetigkeiten der Ableitungen. Beispiel: $f(x) = |x|$. \\ \quad Unstetig & Die Funktion hat Unstetigkeiten oder ist sehr verrauscht. Beispiel: experimentelle Daten mit Artefakten. \\ \midrule \textbf{Anzahl der Minima:} & \\ \quad Unimodal & Es gibt nur ein globales Minimum (oder Maximum). Beispiel: konvexe Funktionen. \\ \quad Multimodal & Es gibt mehrere lokale Minima, und die Aufgabe ist, das globale zu finden. Beispiel: Potentiale komplexer Moleküle. \\ \bottomrule \end{tabularx} \subsection{Methoden zur Lösung nichtlinearer Probleme} Nun systematisieren wir die Lösungsmethoden. Jede Methode erfordert etwas, hängt von etwas ab, und jede hat ihre Stärken und Schwächen. \subsubsection{Monte-Carlo-Methoden} \textbf{Idee:} Zufällige Suche. Wir erzeugen zufällige Punkte im Parameterraum, werten die Funktion dort aus und wählen die beste. \textbf{Anforderungen:} Keine. Sie kann sogar mit unstetigen Funktionen arbeiten. \textbf{Vorteile:} \begin{itemize} \item Bleibt nicht in lokalen Minima hängen (bei ausreichender Iterationszahl), \item Einfach zu implementieren, \item Parallelisiert perfekt. \end{itemize} \textbf{Nachteile:} \begin{itemize} \item Sehr langsame Konvergenz: Der Fehler nimmt ab wie $\mathcal{O}(1/\sqrt{N})$, wobei $N$ die Anzahl der Iterationen ist, \item Nutzt keine Informationen über die Struktur der Funktion. \end{itemize} \textbf{Wann verwenden:} Wenn die Funktion sehr verrauscht, unstetig ist oder wenn man einfach ``irgendeine'' Lösung zur Initialisierung anderer Methoden finden muss. \subsubsection{Gradientenverfahren (Gradient Descent)} \textbf{Idee:} Sich in die dem Gradienten entgegengesetzte Richtung bewegen: \[ \vec{x}_{k+1} = \vec{x}_k - \alpha_k \nabla f(\vec{x}_k) \] wobei $\alpha_k$ die Lernrate ist. \textbf{Anforderungen:} Die Funktion muss differenzierbar (glatt) sein. \textbf{Vorteile:} \begin{itemize} \item Einfach zu implementieren, \item Garantierte Konvergenz zu einem lokalen Minimum (bei richtiger Wahl von $\alpha$), \item Funktioniert gut in hohen Dimensionen. \end{itemize} \textbf{Nachteile:} \begin{itemize} \item Bleibt in lokalen Minima hängen, \item Kann in ``Schluchten'' sehr langsam konvergieren (wenn die Eigenwerte der Hesse-Matrix stark differieren), \item Erfordert die Wahl der Schrittweite $\alpha$ (zu groß --- divergiert, zu klein --- langsam). \end{itemize} \subsubsection{Verfahren der konjugierten Richtungen (Conjugate Gradient für Optimierung)} \textbf{Idee:} Suchrichtungen so konstruieren, dass sie bezüglich der Hesse-Matrix der Funktion ``konjugiert'' sind. Dies erlaubt, das Minimum einer quadratischen Funktion in $N$ Schritten zu finden (wobei $N$ die Dimension ist). \textbf{Wichtige Bemerkung:} Das Verfahren der konjugierten Gradienten zur Lösung linearer Systeme iterativ \textit{heißt nur} genauso wie das Verfahren der konjugierten Gradienten zur nichtlinearen Minimierung. Tatsächlich ist es derselbe Algorithmus, nur auf verschiedene Probleme angewendet: \begin{itemize} \item Für lineare Systeme $\mathbf{A}\vec{x} = \vec{b}$ --- das ist die Minimierung der quadratischen Funktion $f(\vec{x}) = \frac{1}{2}\vec{x}^T\mathbf{A}\vec{x} - \vec{b}^T\vec{x}$, \item Für nichtlineare Optimierung --- das ist die Verallgemeinerung derselben Idee auf beliebige Funktionen. \end{itemize} \textbf{Anforderungen:} Eine glatte Funktion, möglichst mit stetigem Gradienten. \textbf{Vorteile:} \begin{itemize} \item Konvergiert schneller als gewöhnlicher Gradientenabstieg, \item Erfordert keine Speicherung der Hesse-Matrix (im Gegensatz zum Newton-Verfahren), \item Speicher --- $\mathcal{O}(N)$. \end{itemize} \textbf{Nachteile:} \begin{itemize} \item Bleibt in lokalen Minima hängen, \item Für nichtquadratische Funktionen ist ein ``Neustart'' (reset) der Richtungen erforderlich. \end{itemize} \subsubsection{Newton-Verfahren} \textbf{Idee:} Nicht nur den Gradienten, sondern auch die zweiten Ableitungen (die Hesse-Matrix) verwenden. Wir entwickeln die Funktion in eine Taylor-Reihe bis zur zweiten Ordnung: \[ f(\vec{x} + \vec{h}) \approx f(\vec{x}) + \nabla f(\vec{x})^T \vec{h} + \frac{1}{2} \vec{h}^T \mathbf{H}(\vec{x}) \vec{h} \] und finden das Minimum dieser quadratischen Approximation. Wir erhalten die Iteration: \[ \vec{x}_{k+1} = \vec{x}_k - \mathbf{H}^{-1}(\vec{x}_k) \nabla f(\vec{x}_k) \] wobei $\mathbf{H}$ die Matrix der zweiten Ableitungen (die Hesse-Matrix) ist. \textbf{Anforderungen:} Die Funktion muss zweimal differenzierbar sein. Die Hesse-Matrix muss berechenbar sein. \textbf{Vorteile:} \begin{itemize} \item \textbf{Quadratische Konvergenz:} Wenn man nahe an der Lösung startet, verdoppelt sich die Anzahl der korrekten Ziffern bei jeder Iteration, \item Unempfindlich gegenüber ``Schluchten'' (im Gegensatz zum Gradientenabstieg). \end{itemize} \textbf{Nachteile:} \begin{itemize} \item Erfordert die Berechnung der Hesse-Matrix --- $\mathcal{O}(N^2)$ Elemente, \item Erfordert die Lösung eines linearen Systems mit der Hesse-Matrix --- $\mathcal{O}(N^3)$ Operationen, \item Kann divergieren, wenn man weit von der Lösung startet, \item Die Hesse-Matrix kann entartet oder nicht positiv definit sein. \end{itemize} \subsubsection{Das Newton--Raphson-Verfahren und Block-Newton für die Quantenmechanik} Das \textbf{Newton--Raphson-Verfahren} ist eine Variante des Newton-Verfahrens zur Lösung von Systemen nichtlinearer Gleichungen $\vec{F}(\vec{x}) = \vec{0}$: \[ \vec{x}_{k+1} = \vec{x}_k - \mathbf{J}^{-1}(\vec{x}_k) \vec{F}(\vec{x}_k) \] wobei $\mathbf{J}$ die Jacobi-Matrix ist (die Matrix der ersten Ableitungen $\partial F_i / \partial x_j$). \textbf{Block-Newton für die Quantenmechanik:} In quantenchemischen Problemen (zum Beispiel im Hartree--Fock-Verfahren oder DFT) entsteht oft ein Gleichungssystem, das in Blöcke aufgeteilt werden kann: \begin{itemize} \item Gleichungen für die Orbitale (Entwicklungskoeffizienten), \item Gleichungen für die Dichte, \item Gleichungen für die Energie. \end{itemize} Das Block-Newton-Verfahren berücksichtigt diese Struktur: Die Hesse-Matrix wird in Blockform dargestellt, und jeder Block wird separat behandelt. Dies erlaubt: \begin{itemize} \item Effiziente Nutzung der Struktur des Problems, \item Parallelisierung der Berechnungen, \item Anwendung verschiedener Methoden auf verschiedene Blöcke (zum Beispiel Newton für einen Block, Gradientenabstieg für einen anderen). \end{itemize} \subsubsection{BFGS- und L-BFGS-Verfahren} \textbf{BFGS (Broyden--Fletcher--Goldfarb--Shanno)} ist ein quasi-Newton-Verfahren. Die Idee: Die Hesse-Matrix nicht explizit zu berechnen, sondern sie bei jeder Iteration unter Verwendung von Informationen über die Änderung des Gradienten anzunähern. \textbf{Anforderungen:} Eine glatte Funktion mit stetigem Gradienten. \textbf{Vorteile:} \begin{itemize} \item Konvergenz fast wie bei Newton, aber ohne Berechnung der Hesse-Matrix, \item Garantierte positive Definitheit der Hesse-Approximation, \item Funktioniert in der Praxis gut --- eine der beliebtesten Methoden. \end{itemize} \textbf{Nachteile:} \begin{itemize} \item Erfordert die Speicherung einer $N \times N$-Matrix (Hesse-Approximation) --- $\mathcal{O}(N^2)$ Speicher. \end{itemize} \textbf{L-BFGS (Limited-memory BFGS)} ist eine Modifikation von BFGS für große Probleme. Die Idee: Nicht die volle Matrix zu speichern, sondern nur die letzten $m$ Paare von Vektoren $(\vec{s}_k, \vec{y}_k)$, wobei: \[ \vec{s}_k = \vec{x}_{k+1} - \vec{x}_k, \quad \vec{y}_k = \nabla f_{k+1} - \nabla f_k \] \textbf{Vorteile:} \begin{itemize} \item Speicher --- $\mathcal{O}(mN)$, wobei $m \sim 5\text{--}20$ (hängt nicht von $N^2$ ab!), \item Ideal für große Probleme ($N > 1000$), \item Sehr beliebt im maschinellen Lernen und in der Quantenchemie. \end{itemize} \textbf{Nachteile:} \begin{itemize} \item Konvergiert etwas langsamer als vollständiges BFGS, \item Erfordert die Einstellung des Parameters $m$. \end{itemize} \subsubsection{Simplex-Verfahren (Nelder--Mead)} \textbf{Idee:} Wir konstruieren ein Simplex (im $N$-dimensionalen Raum sind das $N+1$ Punkte), werten die Funktion an den Ecken aus und spiegeln, ziehen zusammen oder erweitern das Simplex iterativ, um uns dem Minimum zu nähern. \textbf{Anforderungen:} Keine! Die Funktion kann nichtglatt, verrauscht, sogar unstetig sein. \textbf{Vorteile:} \begin{itemize} \item Erfordert keine Gradienten, \item Einfach zu implementieren, \item Funktioniert gut für kleine Dimensionen ($N < 10$). \end{itemize} \textbf{Nachteile:} \begin{itemize} \item Sehr langsam für große $N$, \item Kann hängen bleiben, \item Keine theoretischen Konvergenzgarantien. \end{itemize} \textbf{Wann verwenden:} Wenn die Funktion sehr verrauscht ist oder wenn $N$ klein ist und man zu faul ist, Gradienten abzuleiten. \subsubsection{Simulated Annealing (Simulierte Abkühlung)} \textbf{Idee:} Inspiriert vom physikalischen Prozess des Ausglühens von Metallen. Wir beginnen mit einer hohen ``Temperatur'' $T$, die es dem Algorithmus erlaubt, über lokale Minima ``hinwegzuspringen''. Wir senken $T$ allmählich, und der Algorithmus ``friert'' im globalen Minimum ein. Bei jedem Schritt: \begin{enumerate} \item Wir schlagen eine zufällige Änderung $\vec{x} \to \vec{x}'$ vor, \item Wir berechnen $\Delta f = f(\vec{x}') - f(\vec{x})$, \item Wenn $\Delta f < 0$ --- akzeptieren wir die Änderung, \item Wenn $\Delta f > 0$ --- akzeptieren wir mit Wahrscheinlichkeit $P = \exp(-\Delta f / T)$. \end{enumerate} \textbf{Anforderungen:} Keine. \textbf{Vorteile:} \begin{itemize} \item Kann das globale Minimum finden (bei richtigem ``Abkühlungsplan''), \item Bleibt in frühen Phasen nicht in lokalen Minima hängen, \item Einfach zu implementieren. \end{itemize} \textbf{Nachteile:} \begin{itemize} \item Sehr langsam, \item Erfordert die Einstellung des Plans $T(t)$, \item Keine Konvergenzgarantien in vernünftiger Zeit. \end{itemize} \textbf{Wann verwenden:} Wenn das Problem multimodal ist (viele lokale Minima) und das globale gefunden werden muss. \subsection{Berechnung von Gradienten} Alle Gradientenverfahren erfordern die Berechnung von $\nabla f(\vec{x})$. Wie macht man das? \subsubsection{Finite Differenzen} \textbf{Idee:} Wir approximieren die Ableitung: \[ \frac{\partial f}{\partial x_j} \approx \frac{f(\vec{x} + h \vec{e}_j) - f(\vec{x})}{h} \] wobei $\vec{e}_j$ der $j$-te Basisvektor ist. \textbf{Probleme:} \begin{enumerate} \item \textbf{Teuer:} Zur Berechnung des Gradienten im $N$-dimensionalen Raum braucht man $N+1$ Funktionsauswertungen (oder $2N$ für die zentrale Differenz). \item \textbf{Instabil:} Wenn $h$ zu groß ist --- großer Approximationsfehler. Wenn $h$ zu klein ist --- Genauigkeitsverlust durch Subtraktion naher Zahlen. Optimal ist $h \sim \sqrt{\epsilon_{\text{mach}}}$, wobei $\epsilon_{\text{mach}}$ die Maschinengenauigkeit ist. \item \textbf{Rauschen:} Wenn die Funktion mit Rauschen berechnet wird (zum Beispiel experimentelle Daten), dann verstärken finite Differenzen das Rauschen. \end{enumerate} \subsubsection{Automatisches Differenzieren (AD): analytische Berechnung des Gradienten} Und nun --- Magie! Es stellt sich heraus, dass man den Gradienten einer Funktion \textit{analytisch} berechnen kann, indem man automatisches Differenzieren verwendet, und das in nur \textbf{einigen wenigen Malen} so teuer wie die Berechnung der Funktion selbst! \textbf{Idee:} Jede Funktion $f(\vec{x})$ ist eine Komposition elementarer Operationen ($+$, $-$, $\times$, $\div$, $\sin$, $\cos$, $\exp$, $\log$ usw.). Wir können: \begin{enumerate} \item Die Berechnung der Funktion als \textbf{Berechnungsgraphen} darstellen, \item Die \textbf{Kettenregel} in umgekehrter Reihenfolge anwenden (reverse mode). \end{enumerate} \textbf{Ergebnis:} Der Gradient wird zu $\mathcal{O}(1) \times$ den Kosten der Funktion berechnet (in der Praxis --- 3--5 Mal teurer, aber nicht $N$ Mal!). \textbf{Wie es im Kern funktioniert:} \begin{enumerate} \item Vorwärtspass: Wir berechnen $f(\vec{x})$ und speichern alle Zwischenergebnisse. \item Rückwärtspass: Wir gehen den Graphen in umgekehrter Richtung durch und wenden die Kettenregel an: \[ \frac{\partial f}{\partial x_j} = \sum_{\text{Pfade}} \frac{\partial f}{\partial u_1} \frac{\partial u_1}{\partial u_2} \cdots \frac{\partial u_k}{\partial x_j} \] \end{enumerate} \textbf{Wo es verwendet wird:} \begin{itemize} \item \textbf{PyTorch, TensorFlow:} Die Grundlage aller modernen neuronalen Netze ist genau das reverse-mode automatic differentiation. \item \textbf{Quantenchemie:} Programme wie PySCF, Psi4 verwenden AD zur Berechnung von Energiegradienten bezüglich der Kernkoordinaten. \item \textbf{Optimierung:} Alle modernen Bibliotheken (SciPy, JAX) verwenden AD. \end{itemize} \begin{successbox}[Hauptschlussfolgerung] Berechne niemals Gradienten mit finiten Differenzen, wenn automatisches Differenzieren verwendet werden kann! Es ist schneller, genauer und stabiler. \end{successbox} \subsection{Ein lebendiges Beispiel: Kraftfelder und Konformere} Ein weiteres lebendiges Beispiel aus der Chemie ist die Modellierung von Molekülen. Zur Modellierung kleiner Moleküle verwendet man oft eine Darstellung, in der die Gesamtenergie als Summe einfacher Wechselwirkungen geschrieben wird: \begin{enumerate} \item \textbf{Van-der-Waals-Wechselwirkungen} zwischen nicht gebundenen Atomen (Lennard-Jones-Potential): \[ E_{\text{vdW}} = \sum_{i < j} 4\varepsilon_{ij} \left[ \left(\frac{\sigma_{ij}}{r_{ij}}\right)^{12} - \left(\frac{\sigma_{ij}}{r_{ij}}\right)^6 \right] \] \item \textbf{Harmonische Bindungen} zwischen gebundenen Atomen: \[ E_{\text{bond}} = \sum_{\text{Bindungen}} \frac{1}{2} k_b (r - r_0)^2 \] \item \textbf{Winkelwechselwirkungen} zwischen drei aufeinanderfolgenden Atomen: \[ E_{\text{angle}} = \sum_{\text{Winkel}} \frac{1}{2} k_\theta (\theta - \theta_0)^2 \] \item \textbf{Torsions- (Diederv-)Wechselwirkungen} zwischen vier aufeinanderfolgenden Atomen: \[ E_{\text{torsion}} = \sum_{\text{Dieder}} \frac{V_n}{2} [1 + \cos(n\phi - \gamma)] \] \end{enumerate} Es wird nicht viele geben, aber etwa das Quadrat der Anzahl der Atome (für Van-der-Waals). Jede solche Funktion ist eine einfache Formel: manchmal mit Sinus, manchmal mit Exponentialfunktionen, manchmal mit Potenzen. Und insgesamt hat ein Molekül $3N$ Freiheitsgrade, da jedes Atom 3 Raumkoordinaten hat. \subsubsection{Das Konformerproblem} Es scheint, die Aufgabe ist klar: Man muss das Minimum dieser Funktion $E(\vec{x})$ finden, wobei $\vec{x} \in \R^{3N}$ die Koordinaten aller Atome sind. Aber hier beginnt das Interessanteste. Die Energiefunktion eines Moleküls ist eine \textbf{multimodale Landschaft} mit einer enormen Anzahl lokaler Minima. Jedes lokale Minimum entspricht einer eigenen \textbf{Konformation} (Konformer) des Moleküls --- einer eigenen Art, das Molekül im Raum zu verdrehen. Zum Beispiel kann für ein Protein aus 100 Aminosäuren die Anzahl möglicher Konformere $10^{100}$ erreichen --- das ist das berühmte \textbf{Levinthal-Paradoxon}. Und das globale Minimum (die native Struktur) ist nur eine von $10^{100}$ Möglichkeiten. \subsubsection{Warum einfache Minimierung nicht funktioniert} Wenn wir eine zufällige Anfangskonformation nehmen und gewöhnlichen Gradientenabstieg oder BFGS starten, konvergieren wir sehr schnell (innerhalb einiger hundert Iterationen) zum \textit{nächstgelegenen} lokalen Minimum. Aber das wird fast sicher \textit{nicht} das globale Minimum sein --- also nicht die Konformation, die das Molekül in der Realität annimmt. \subsubsection{Lösung: eine Kombination von Methoden} Deshalb verwendet man in der Praxis eine Kombination von Methoden: \begin{enumerate} \item \textbf{Simulated Annealing (simulierte Abkühlung):} Wir beginnen mit einer hohen ``Temperatur'', die es dem Algorithmus erlaubt, über Energiebarrieren zu springen und verschiedene Bereiche des Konformationsraums zu erkunden. Wir senken die Temperatur allmählich und ``frieren'' das System in einem energiearmen Bereich ein. \item \textbf{Lokale Minimierung (BFGS):} Nachdem Simulated Annealing einen ``guten'' Bereich gefunden hat, starten wir eine schnelle lokale Methode (BFGS oder das Verfahren der konjugierten Gradienten) für die präzise Konvergenz zu einem lokalen Minimum. \item \textbf{Wiederholung:} Wir starten diese Kombination viele Male mit verschiedenen Anfangsbedingungen und wählen das Konformer mit der niedrigsten Energie. \end{enumerate} \begin{tipbox}[Beliebte Programme] Dieser Ansatz ist in beliebten Programmen der Molekulardynamik implementiert: \begin{itemize} \item \textbf{AMBER, GROMACS, NAMD:} Verwenden Kraftfelder (AMBER, CHARMM, OPLS) und Methoden wie Simulated Annealing + lokale Minimierung. \item \textbf{Rosetta:} Verwendet für die Vorhersage der Proteinstruktur einen Fragmentansatz + Monte Carlo + lokale Minimierung. \end{itemize} \end{tipbox} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Die schnelle Fourier-Transformation: wenn $O(N^2)$ zu $O(N \log N)$ wird} \subsection{Das Problem der Mustersuche} Manchmal müssen wir ein bestimmtes Fragment in einer riesigen Sequenz finden. Wenn es ein sehr kurzes Fragment ist, dann ist eine einfache Durchmusterung der ganzen Sequenz und der Vergleich mit diesem kurzen Fragment eine recht gute Strategie. So verfährt man oft bei der Suche nach kurzen Fragmenten in Proteinstrukturen. Aber manchmal sind die Dimensionen dessen, was verglichen werden soll, fast gleich. Zum Beispiel haben wir eine periodische Kristallstruktur, und wir wissen, dass es in ihr eine gewisse Abweichung gibt, und wir verstehen sogar, dass die Abweichung nur teilweise dem ähnelt, womit wir vergleichen wollen. Formal: Die Ausgangsdaten sind $v_i$, $i = 1, \dots, N$, und das, womit verglichen werden soll, ist $w_j$, $j = 1, \dots, M$, wobei $M < N$. Wir können explizit in einer Schleife berechnen: \[ \forall k = 0, \dots, N-M: \quad s_k = \sum_{j=1}^M w_j v_{j+k} \] und das Maximum über $s$ finden. Aber die Rechenkomplexität einer solchen Suche ist quadratisch in $N$: Wenn $M \simeq N/2$, dann müssen $N^2/2$ arithmetische Operationen ausgeführt werden. Für $N = 10^6$ sind das $5 \times 10^{11}$ Operationen --- zu viele! \subsection{Matrixformulierung} Hier haben sich die Leute ausgedacht, wie man solche $s_k$ viel schneller finden kann. Stellen wir die Berechnung dieser $s$ in Matrixform dar. Angenommen, wir haben die Matrix: \[ \mathbf{W} = \begin{pmatrix} w_1 & w_2 & \cdots & w_M & 0 & \cdots & 0 \\ 0 & w_1 & w_2 & \cdots & w_M & \ddots & \vdots \\ \vdots & \ddots & \ddots & \ddots & & \ddots & 0 \\ 0 & \cdots & 0 & w_1 & w_2 & \cdots & w_M \end{pmatrix} \] der Größe $(N-M+1) \times N$. Wenn wir sie mit unserem Vektor $\vec{v}$ multiplizieren, erhalten wir genau unsere ersehnten $s_k$: \[ \vec{s} = \mathbf{W} \vec{v} \] Aber wie machen wir das schnell? \subsection{Zirkulante Matrix} Wenn wir diese Matrix von unten um weitere $M-1$ Zeilen der Form ergänzen: \begin{align*} & w_M, 0, \dots, 0, w_1, \dots, w_{M-1} \\ & w_{M-1}, w_M, 0, \dots, 0, w_1, \dots, w_{M-2} \\ & \vdots \\ & w_2, \dots, w_M, 0, \dots, 0, w_1 \end{align*} dann erhalten wir nach der Multiplikation denselben Vektor $\vec{s}$, aber am Ende werden ihm noch $M-1$ weitere Zahlen angehängt, die wir verwerfen können. Aber diese erweiterte Matrix ist eine \textbf{zirkulante Matrix} $\mathbf{C}$! Sie hat eine bemerkenswerte Eigenschaft: Jede nachfolgende Zeile entsteht durch zyklische Verschiebung der vorherigen. \subsection{Diagonalisierung der zirkulanten Matrix} Eine zirkulante Matrix hat eine interessante Darstellung über die Fourier-Matrix: \[ \mathbf{C} = \frac{1}{N} \mathbf{F}^H \diag(\mathbf{F} \vec{c}) \mathbf{F} \] wobei: \begin{itemize} \item $\vec{c}$ die erste Spalte der Ausgangsmatrix $\mathbf{C}$ ist, \item $\mathbf{F}$ die Matrix der diskreten Fourier-Transformation (DFT) ist, \item $\mathbf{F}^H$ die hermitesch konjugierte (komplex konjugierte und transponierte) Matrix ist. \end{itemize} Die Fourier-Matrix ist definiert als: \[ \mathbf{F} = \{f_{jk}\}_{j,k=0}^{N-1}, \quad f_{jk} = e^{-2\pi i j k / N} \] \subsection{Schnelle Multiplikation über FFT} Nun das Interessanteste. Wenn wir die Fourier-Matrix schnell mit einem Vektor multiplizieren können, dann können wir dieses $\vec{s}$ schneller berechnen: \[ \vec{s} = \mathbf{C} \vec{v} = \frac{1}{N} \mathbf{F}^H \diag(\mathbf{F} \vec{c}) \mathbf{F} \vec{v} \] Diese Berechnung besteht aus drei Schritten: \begin{enumerate} \item $\vec{a} = \mathbf{F} \vec{v}$ --- direkte DFT des Vektors $\vec{v}$, \item $\vec{b} = \diag(\mathbf{F} \vec{c}) \vec{a}$ --- elementweise Multiplikation (das ist $\mathcal{O}(N)$), \item $\vec{s} = \frac{1}{N} \mathbf{F}^H \vec{b}$ --- inverse DFT. \end{enumerate} \subsection{Zerlegung der Fourier-Matrix} Wie multipliziert man $\mathbf{F}$ schnell mit einem Vektor? Es stellt sich heraus, dass $\mathbf{F}$ als Produkt spezieller Matrizen dargestellt werden kann: \begin{enumerate} \item Einer Permutationsmatrix (bit-reversal permutation), \item Mehrere Paare von Matrizen: \begin{itemize} \item Blockmatrizen der Form $\begin{pmatrix} \mathbf{I} & \mathbf{I} \\ \mathbf{I} & -\mathbf{I} \end{pmatrix}$, \item Diagonalmatrizen mit Elementen $e^{-2\pi i k / N}$ (die sogenannten ``twiddle factors''). \end{itemize} \end{enumerate} Schreiben wir diese Matrizen explizit für $N = 8$ auf: \textbf{Permutationsmatrix (bit-reversal):} \[ \mathbf{P} = \begin{pmatrix} 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \end{pmatrix} \] (Zeilen gemäß Bit-Umkehrung permutiert: $0, 4, 2, 6, 1, 5, 3, 7$) \textbf{Blockmatrix (erste Stufe):} \[ \mathbf{B}_1 = \begin{pmatrix} 1 & 0 & 0 & 0 & 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 & 0 & 1 \\ 1 & 0 & 0 & 0 & -1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & -1 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 & -1 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 & 0 & -1 \end{pmatrix} \] \textbf{Diagonalmatrix (twiddle factors, erste Stufe):} \[ \mathbf{D}_1 = \diag\left(1, 1, 1, 1, 1, e^{-2\pi i / 8}, e^{-4\pi i / 8}, e^{-6\pi i / 8}\right) \] Und so weiter für $\log_2 N$ Stufen. Dann erfordert jede solche Multiplikation nur $N$ und $2N$ arithmetische Operationen, und die Gesamtzahl solcher Schritte ist $\log_2 N$. \textbf{Insgesamt:} Die Multiplikation der Fourier-Matrix mit einem Vektor kann in $\mathcal{O}(N \log_2 N)$ arithmetischen Operationen durchgeführt werden! Und damit können auch alle $s_k$ in denselben $\mathcal{O}(N \log_2 N)$ berechnet werden, und nicht in $\mathcal{O}(N^2)$. Für $N = 10^6$: \begin{itemize} \item Naiver Algorithmus: $10^{12}$ Operationen, \item FFT: $2 \times 10^7$ Operationen (50\,000 Mal schneller!). \end{itemize} \subsection{Die schnelle Fourier-Transformation (FFT)} Dies ist die berühmte \textbf{schnelle Fourier-Transformation} (Fast Fourier Transform, FFT), entdeckt von Cooley und Tukey im Jahr 1965 (obwohl Gauss sie bereits 1805 kannte!). Sie ist einer der gefragtesten Algorithmen in der modernen Chemie. \subsubsection{Anwendung in der NMR} Mit ihrer Hilfe wird zum Beispiel das ursprüngliche FID (Free Induction Decay) aus der NMR in den Spektralbereich transformiert. Wenn man sich jede Zeile der Fourier-Matrix ansieht, kann man in ihr eine oszillierende Funktion erkennen: \begin{itemize} \item Je näher die Zeile an der Mitte der Matrix liegt, desto größer die Oszillation, \item Die erste Zeile hat überhaupt keine --- sie ist einfach eine Konstante, \item Die unterste (oder zweite) Zeile hat eine Oszillationsperiode gleich der Länge der Zeile. \end{itemize} Genau deshalb erhalten wir, nachdem wir die Fourier-Matrix mit dem FID multipliziert haben, Maxima dort, wo die Resonanzfrequenzen sind: Es ist einfach so, dass eine solche Zeile in der Oszillationsfrequenz mit der Resonanzfrequenz übereinstimmte. \begin{tipbox}[Warum ist FFT für die NMR so wichtig] Moderne FIDs sind sehr lang --- sie werden oft in Hunderttausende und sogar Millionen von Zahlen digitalisiert. Wenn es die schnelle Fourier-Transformation nicht gäbe, würde die Multiplikation mit einer solchen Matrix auf modernen Computern sogar länger dauern als die NMR-Aufnahme selbst --- also wirklich Minuten und Zehner von Minuten selbst auf modernen Workstations. Dank der FFT dauert die Transformation \textbf{Bruchteile einer Sekunde}. \end{tipbox} \subsubsection{Weitere Anwendungen der FFT in der Chemie} \begin{itemize} \item \textbf{Kryo-Elektronenmikroskopie:} Rekonstruktion von 3D-Strukturen aus 2D-Projektionen. \item \textbf{Röntgenstrukturanalyse:} Umwandlung des Beugungsbildes in Elektronendichte. \item \textbf{Molekulardynamik:} Berechnung elektrostatischer Wechselwirkungen über die PPPM-Methode (Particle-Particle Particle-Mesh). \item \textbf{Signalverarbeitung:} Rauschfilterung, Datenkompression. \item \textbf{Chemometrie:} Schnelle Faltung und Korrelation von Spektren. \end{itemize} \begin{successbox}[Fazit] FFT ist einer jener Algorithmen, die die Welt verändert haben. Ohne ihn gäbe es keine moderne Spektroskopie, Signalverarbeitung, Audio- und Videokompression und vieles mehr. Es ist ein Muss für jeden Chemiker, der mit Daten arbeitet. \end{successbox} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Wenn die FFT das Offensichtliche nicht sieht: Spektrales Leck und die Prony-Methode} \subsection{Das Paradoxon: Das Auge sieht eine Sinuskurve, die FFT nicht} Obwohl die FFT ein Muss-Algorithmus fast überall in der experimentellen Chemie ist, hat auch sie ihre Besonderheiten, auf die man unbedingt achten muss. Betrachten wir ein Beispiel. Wir haben ein Spektrum aufgenommen, das Signal oszilliert, und das ist mit bloßem Auge im Diagramm sichtbar. Aber wir konnten nur etwas mehr als eine Periode aufzeichnen. Als wir die FFT anwendeten, erhielten wir etwas sehr Merkwürdiges: Der Peak scheint da zu sein, aber auch wieder nicht --- er ist irgendwie sehr verschmiert. Und wenn wir viel anderes Rauschen hatten, erhielten wir in einem anderen Teil des Spektrums etwas sehr Ähnliches, und ja, zufällig erwies sich das Rauschen als ``heller'' als das Signal, und unser Programm lieferte ein völlig falsches Ergebnis. Aber wir sehen diese Sinuskurve doch mit unseren Augen! Warum sieht die FFT sie nicht? \subsection{Spektrales Leck (Spectral Leakage)} Die FFT ist dann gut, und \textit{nur dann}, wenn wir sie auf ein Signal angewendet haben, in das \textbf{genau eine ganze Zahl von Perioden} dieses Signals passt. Warum? Weil die FFT implizit annimmt, dass unser endliches Signal \textit{periodisch fortgesetzt} wird über das Messfenster hinaus. Wenn genau $k$ volle Perioden in das Fenster passen, dann ist die periodische Fortsetzung glatt, und die FFT liefert einen klaren Peak. Aber wenn das Signal zum Beispiel $1.4$ Perioden enthält, dann entsteht bei der periodischen Fortsetzung ein \textbf{Sprung} an der Grenze des Fensters. Die FFT ist gezwungen, diesen Sprung durch eine Vielzahl von Harmonischen zu approximieren --- und Energie ``leckt'' aus dem Hauptpeak in alle anderen Frequenzen. Dieses Phänomen heißt \textbf{spektrales Leck} (spectral leakage). Mathematisch: Wenn die wahre Frequenz des Signals \textit{zwischen} zwei benachbarte Frequenzbins der FFT fällt, dann erhalten wir statt einer Delta-Funktion einen verschmierten ``Buckel'' (eine Funktion der Form $\sin(x)/x$, die sogenannte \textit{sinc}). \begin{warningbox}[FFT-Regel] Die FFT sieht nur die Frequenzen, die ganzzahlig oft in das Messfenster passen. Alles andere ist Leck. \end{warningbox} \subsection{Zero-Padding: eine einfache, aber nicht ideale Lösung} Was sollen wir dann tun? Die einfachste Variante ist, viele Nullen an das Ende des Signals anzuhängen (zero-padding) und zu hoffen, dass wir insgesamt wesentlich mehr Glück haben. Denn wenn wir zum Beispiel Nullen so anhängen, dass die Ausgangsdaten um den Faktor 5 gestreckt werden, dann erhalten wir sicher 7 Perioden (statt 1.4). Und wenn wir sie zum Beispiel um den Faktor 7 vergrößern, dann sind es 9.8 Perioden --- auch sehr nahe an einer ganzen Zahl. Oft hängt man eine unterschiedliche Anzahl von Nullen an und versucht dann, die erhaltenen Spektren zu kombinieren, indem man errät, welche Peaks in den entsprechenden erweiterten Spektren genauer waren. Das hilft wirklich! \begin{tipbox}[Wichtige Nuance des Zero-Padding] Zero-Padding \textit{erhöht nicht die tatsächliche Frequenzauflösung} --- es interpoliert nur das Spektrum und macht es glatter. Die Auflösung wird durch die \textit{Länge des Ausgangssignals} bestimmt, nicht durch die Länge des mit Nullen aufgefüllten. Aber in der Praxis hilft Zero-Padding, eine ganze Zahl von Perioden zu ``treffen'' und das Leck zu verringern. \end{tipbox} \subsection{Frequenzdrift: wenn das Signal ``davon driftet''} Aber es gibt noch ein Problem, mit dem Zero-Padding nicht fertig wird. Manchmal messen wir lange und mühsam ein Spektrum, und es ``ist leicht'' davongedriftet. Das heißt, am Anfang hatten wir die Frequenz $\omega$, und am Ende $\omega + \varepsilon$, und diese kleine Korrektur genügt, damit wir mit der FFT völlig danebenliegen und statt eines klaren Einzelpeaks etwas sehr Verschwommenes erhalten. Aber mit dem Auge sieht man es doch --- da ist die Sinuskurve, hier und dort, am Anfang und am Ende! Nur die FFT sieht sie wieder nicht. Der Grund ist, dass die FFT \textbf{Stationarität} des Signals annimmt --- also dass sich die Frequenzen mit der Zeit nicht ändern. Wenn die Frequenz driftet, kann keine Harmonische der FFT das gesamte Signal gut approximieren. Was sollen wir tun? \subsection{Die Prony-Methode: die Idee der Verschiebungen} Auf diese Frage hat uns vor einigen Jahrhunderten Gaspard de Prony (1795) geantwortet. Genauer gesagt, er hat sich ausgedacht, wie man sein eigenes Problem löst (die Entwicklung eines Signals in eine Summe von Exponentialfunktionen), aber sie löst gerade auch unser Problem --- wenn das Spektrum in der Zeit leicht davondriftet. \textbf{Idee:} Wir nehmen das Signal und speichern es in einem Vektor $\vec{s} = (s_0, s_1, \dots, s_{N-1})^T$. Dann verschieben wir diesen Vektor um einen Abtastwert nach unten, dann noch einmal, und stellen ihn wieder daneben. Wir machen das $L$ Mal. Was haben wir hier? Wenn wir ein periodisches Signal mit einer unbekannten Sinuskurve $\sin(\omega t)$ hatten, dann erhalten wir nach einer Verschiebung um $p$ Abtastwerte $\sin(\omega t + \omega p)$. Unter Erinnerung an die Additionsformeln aus Kapitel 1: \[ \sin(\omega t + \omega p) = \sin(\omega t)\cos(\omega p) + \cos(\omega t)\sin(\omega p) \] Da $\cos(\omega p)$ und $\sin(\omega p)$ \textit{Konstanten} für den gesamten Vektor sind (sie hängen nicht von $t$ ab), wird die Menge solcher verschobenen Vektoren nur \textbf{so viele von Null verschiedene Singulärwerte haben, wie wir das Doppelte der Anzahl verschiedener oszillierender Harmonischer haben}. Darüber hinaus, wenn eine Frequenz mit der Zeit ein wenig ``davongegangen'' ist, \textit{hat das keinerlei Einfluss} auf eine solche Approximation, zerstört aber das Ergebnis der FFT völlig, wie wir oben bemerkt haben. \subsection{Prony-Methode, Variante 1: lineare Prädiktion} Unser Signal werde als Summe von $K$ komplexen Exponentialfunktionen modelliert: \[ s_n = \sum_{m=1}^{K} c_m z_m^n, \quad n = 0, 1, \dots, N-1 \] wobei $z_m = e^{(\alpha_m + \i \omega_m)\Delta t}$ komplexe ``Frequenzen'' sind (die sowohl die Dämpfung $\alpha_m$ als auch die Oszillation $\omega_m$ enthalten) und $c_m$ komplexe Amplituden. \textbf{Schlüsseltatsache:} Wenn das Signal aus $K$ Exponentialfunktionen besteht und kein Rauschen vorhanden ist, dann erfüllt es eine \textbf{lineare Rekursionsbeziehung} der Ordnung $K$: \[ s_{n} + a_1 s_{n-1} + a_2 s_{n-2} + \cdots + a_K s_{n-K} = 0 \quad \forall n \ge K \] Das bedeutet, dass der $(n+1)$-te Abtastwert \textit{exakt vorhergesagt} werden kann aus $K$ vorhergehenden! Die Koeffizienten $a_1, \dots, a_K$ sind die Koeffizienten des charakteristischen Polynoms, dessen Wurzeln die $z_m$ sind: \[ P(z) = z^K + a_1 z^{K-1} + a_2 z^{K-2} + \cdots + a_K = \prod_{m=1}^{K}(z - z_m) \] \textbf{Wie findet man die Koeffizienten?} Schreiben wir die Rekursionsbeziehung für alle verfügbaren Abtastwerte in Matrixform: \[ \underbrace{ \begin{pmatrix} s_{K-1} & s_{K-2} & \cdots & s_0 \\ s_{K} & s_{K-1} & \cdots & s_1 \\ \vdots & \vdots & \ddots & \vdots \\ s_{N-2} & s_{N-3} & \cdots & s_{N-K-1} \end{pmatrix} }_{\mathbf{S} \;\; (N-K) \times K} \underbrace{ \begin{pmatrix} a_1 \\ a_2 \\ \vdots \\ a_K \end{pmatrix} }_{\vec{a}} = - \underbrace{ \begin{pmatrix} s_K \\ s_{K+1} \\ \vdots \\ s_{N-1} \end{pmatrix} }_{\vec{s}_{\text{future}}} \] Dies ist ein überbestimmtes System (wenn $N > 2K$), und wir lösen es mit der Methode der kleinsten Quadrate: \[ \vec{a} = -(\mathbf{S}^H \mathbf{S})^{-1} \mathbf{S}^H \vec{s}_{\text{future}} \] Nachdem wir $\vec{a}$ gefunden haben, suchen wir die Wurzeln des Polynoms $P(z)$ --- das sind unsere $z_m$, aus denen die Frequenzen $\omega_m$ und Dämpfungen $\alpha_m$ extrahiert werden. \subsection{Prony-Methode, Variante 2: SVD der Verschiebungsmatrix} Die zweite Variante ist rauschresistenter und eleganter. \textbf{Schritt 1: Wir konstruieren die Hankel-Matrix aus den Verschiebungen des Signals.} \[ \mathbf{H} = \begin{pmatrix} s_0 & s_1 & \cdots & s_{L-1} \\ s_1 & s_2 & \cdots & s_L \\ \vdots & \vdots & \ddots & \vdots \\ s_{N-L} & s_{N-L+1} & \cdots & s_{N-1} \end{pmatrix} \] der Größe $(N-L+1) \times L$, wobei $L > K$ (wir wählen $L$ sicher größer als die erwartete Anzahl von Harmonischen). \textbf{Schritt 2: Wir machen die SVD.} \[ \mathbf{H} = \mathbf{U} \mathbf{\Sigma} \mathbf{V}^H \] Wenn das Signal aus $K$ Exponentialfunktionen besteht und kein Rauschen vorhanden ist, dann ist $\rank(\mathbf{H}) = K$, und die letzten $L - K$ Singulärwerte sind null: \[ \sigma_1 \ge \sigma_2 \ge \cdots \ge \sigma_K > 0, \quad \sigma_{K+1} = \cdots = \sigma_L = 0 \] Bei Rauschen sind die kleinen Singulärwerte nicht null, aber sie sind \textit{wesentlich kleiner} als die ersten $K$. Wir verwerfen sie (dies nennt man \textbf{truncated SVD} oder \textbf{Regularisierung}). \textbf{Schritt 3: Wir extrahieren das Polynom aus dem Nullraum.} Der Singulärvektor $\vec{v}_{\min}$, der dem \textit{kleinsten} Singulärwert entspricht (die letzte Spalte von $\mathbf{V}$), liegt im (approximativen) Nullraum der Matrix $\mathbf{H}$. Das bedeutet: \[ \mathbf{H} \vec{v}_{\min} \approx \vec{0} \] Wenn wir die Komponenten $\vec{v}_{\min} = (v_0, v_1, \dots, v_{L-1})^T$ schreiben, dann ist diese Beziehung äquivalent zu: \[ \sum_{j=0}^{L-1} v_j s_{n+j} \approx 0 \quad \forall n \] Das heißt, $\vec{v}_{\min}$ enthält die Koeffizienten des Polynoms: \[ Q(z) = v_0 + v_1 z + v_2 z^2 + \cdots + v_{L-1} z^{L-1} \] \textbf{Schritt 4: Wir suchen die Wurzeln des Polynoms.} Die Wurzeln von $Q(z)$ enthalten $K$ ``Signal''-Wurzeln $z_m = e^{(\alpha_m + \i \omega_m)\Delta t}$ (die innerhalb oder auf dem Einheitskreis liegen) und $L - K$ ``Rausch''-Wurzeln (zufällig verstreut). Aus den Signalwurzeln extrahieren wir: \[ \omega_m = \frac{\arg(z_m)}{\Delta t}, \quad \alpha_m = \frac{\ln|z_m|}{\Delta t} \] \begin{tipbox}[Warum ist die SVD-Variante besser?] Die SVD-Variante der Prony-Methode ist rauschresistenter, weil: \begin{enumerate} \item Die Abschneidung kleiner Singulärwerte automatisch das Rauschen filtert. \item Wir lösen das überbestimmte System nicht direkt (was schlecht konditioniert sein kann), sondern verwenden eine orthogonale Zerlegung. \item Die Anzahl der Harmonischen $K$ wird automatisch aus dem ``Sprung'' im Spektrum der Singulärwerte bestimmt. \end{enumerate} \end{tipbox} \subsection{Grenzen der Prony-Methode} Aber Prony ist kein Allheilmittel. Sobald im Ausgangssignal eine Komponente auftaucht, die sowohl durch die Funktion selbst als auch durch ihre Verschiebungen gut approximiert wird (zum Beispiel weißes Rauschen oder ein sehr breitbandiges Signal), erhalten wir sofort eine ``Wurzel'' dafür, und dann stellt sich heraus, dass sie \textbf{falsch} ist. Weitere Einschränkungen: \begin{itemize} \item Die Methode setzt voraus, dass das Signal eine Summe von \textit{Exponentialfunktionen} ist (einschließlich Sinuskurven als Spezialfall). Wenn das Signal eine andere Struktur hat (zum Beispiel Rechteckimpulse), wird die Methode schlecht funktionieren. \item Die Anzahl der Harmonischen $K$ muss im Voraus bekannt sein oder erraten werden (obwohl die SVD dabei hilft). \item Die Suche nach den Wurzeln eines Polynoms hohen Grades ($L > 50$) kann selbst numerisch instabil sein. \item Die Methode ist empfindlich gegenüber starkem Rauschen: Bei einem niedrigen Signal-Rausch-Verhältnis können falsche Wurzeln sich als echte ``maskieren''. \end{itemize} Gleichzeitig wird diese Methode in der NMR-Spektroskopie sehr aktiv und seit langem verwendet und ergänzt die FFT gut. In der Praxis verwendet man oft einen \textbf{Hybridansatz}: \begin{enumerate} \item FFT für einen schnellen Überblick über das Spektrum und die Bestimmung der Anzahl der Peaks, \item Die Prony-Methode (oder ihre modernen Varianten: MUSIC, ESPRIT, Matrix Pencil) zur genauen Bestimmung von Frequenzen, Dämpfungen und Amplituden einzelner Peaks. \end{enumerate} \begin{successbox}[Fazit] Die FFT ist ein mächtiges, aber ``blindes'' Werkzeug: Sie sieht nur, was in ihr Frequenzraster passt. Die Prony-Methode ``sieht'' Frequenzen zwischen den Bins und ist robust gegenüber Drift, erfordert aber mehr Berechnungen und Vorsicht. In der realen NMR-Spektroskopie arbeiten sie im Tandem. \end{successbox} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Kleinste Quadrate, Residuen und Tikhonov-Regularisierung} \subsection{Problemstellung} Also, wir lösen oft Probleme, die wir Probleme der kleinsten Quadrate nennen. Betrachten wir einige von ihnen, um genauer zu verstehen, wie man das lösen kann. \subsection{Ein Beispiel aus der Medizin: Computertomographie (CT)} Nehmen wir ein Beispiel aus einem benachbarten Bereich --- aus der Medizin. Angenommen, wir haben einen Röntgen-CT-Scanner. Wir haben den Patienten (oder das, was wir untersuchen) unbeweglich auf eine Plattform gelegt und stellen um ihn herum auf gegenüberliegenden Seiten einen Sender und Empfänger von Röntgenstrahlen auf. Wir führen Messungen durch, wie stark die Strahlung vom untersuchten Körper absorbiert wurde, und positionieren diese Empfänger und Sender in verschiedenen Winkeln. Normalerweise liegt der Patient auf einer Liege, und Empfänger und Sender befinden sich auf der Oberfläche eines Rings. Dieser Ring ist so positioniert, dass seine Achse etwa entlang der Liege verläuft, und der Ring selbst rotiert um seine Achse und bewegt sich entlang der Liege. Was geschieht hier? Weiches Gewebe absorbiert Röntgenstrahlen kaum, während Knochen ziemlich stark absorbieren. Wir können den gesamten Raum, in dem sich der Patient befindet, gedanklich in kleine (normalerweise gleiche) Würfel unterteilen --- \textbf{Voxel} (volume elements). Wenn wir die Position von Sender und Empfänger genau kennen, können wir eine Linie (den Röntgenstrahl) ziehen, die unsere Würfel schneidet. Jeder solche Würfel habe eine Unbekannte --- die Intensität der Absorption von Röntgenstrahlen. Wir nummerieren alle diese Würfel und schreiben die Intensität ihrer Absorption als unbekannten Vektor $\vec{x} = (x_1, \dots, x_N)^T$. Dann liefert uns eine gegenseitige Anordnung von Sender und Empfänger eine ziemlich dünnbesetzte Menge von Koeffizienten --- die Längen des Durchgangs des Strahls durch den entsprechenden Würfel. All dies sei in einer Zeile der Matrix $\mathbf{A}$ geschrieben, und die Anzahl der Spalten in ihr sei $N$ (die Anzahl der Voxel), und die Anzahl der Zeilen sei gerade gleich der Anzahl der durchgeführten Messungen (sie sei $M$). Die gemessenen Werte selbst aber (wir können zum Beispiel eine Röntgenquelle und viele Empfänger haben; dann ist jede gegenseitige Anordnung von Empfänger und Quelle eine Zeile) werden im Vektor $\vec{b}$ geschrieben, ebenfalls der Länge $M$. Dann können wir das Problem der kleinsten Quadrate formulieren als: \[ \min_{\vec{x}} \|\mathbf{A}\vec{x} - \vec{b}\|_2 \] das heißt, wir wollen, dass jede Zeile der Matrix $\mathbf{A}$, multipliziert mit $\vec{x}$ (die Gesamtabsorptionsintensität), gleich dem ist, was wir messen. \subsection{Normalgleichungen} Wir erinnern uns, dass man ein solches Problem lösen kann, aber betrachten wir es sorgfältig. Die Residuenfunktion: \[ E = \|\mathbf{A}\vec{x} - \vec{b}\|_2^2 = (\mathbf{A}\vec{x} - \vec{b})^H (\mathbf{A}\vec{x} - \vec{b}) = \vec{x}^H \mathbf{A}^H \mathbf{A} \vec{x} - 2 \vec{x}^H \mathbf{A}^H \vec{b} + \|\vec{b}\|_2^2 \] Wenn wir die Ableitung von $E$ nach allen $x_j$ finden und sie null setzen (im Minimum), erhalten wir das Gleichungssystem: \[ 2 \mathbf{A}^H \mathbf{A} \vec{x} - 2 \mathbf{A}^H \vec{b} = 0 \quad \Rightarrow \quad (\mathbf{A}^H \mathbf{A}) \vec{x} = \mathbf{A}^H \vec{b} \] Dies sind die sogenannten \textbf{Normalgleichungen}. Hier scheint alles gut: Wir erhalten eine positiv semidefinite Matrix $\mathbf{A}^H \mathbf{A}$ und eine rechte Seite, mit der wir dieses Problem lösen können. \subsection{Das Entartungsproblem} Aber was, wenn wir solche Würfel konstruiert haben, aber keine einzige Messung hatten, bei der der Röntgenstrahl durch zum Beispiel den $i$-ten Würfel ging? Offensichtlich wird dann in der Ausgangsmatrix $\mathbf{A}$ die entsprechende $i$-te Spalte nur Nullen enthalten, und die Matrix $\mathbf{A}^H \mathbf{A}$ wird Nullen sowohl in der $i$-ten Spalte als auch in der $i$-ten Zeile haben, also garantiert einen Singulärwert null haben. Angenommen, dieser Würfel war irgendwo im Raum, und in ihm ist kein Teil des untersuchten Körpers, das heißt --- zu unserem Glück --- brauchen wir diesen Würfel nicht. Aber eine solche Formulierung \textit{verdirbt} die Lösung: Wir können dieses Problem einfach nicht lösen, da die Matrix $\mathbf{A}^H \mathbf{A}$ entartet sein wird. Darüber hinaus können wir nicht garantieren, dass selbst wenn $\mathbf{A}^H \mathbf{A}$ keine solchen Nullspalten und -zeilen hat, ihre Kondition gut sein wird. Das heißt, wir können statt einer Lösung völlig falsche Zahlen erhalten. \subsection{Lösung über SVD und Pseudoinverse} Wie soll man in diesem Fall vorgehen? Kehren wir zur Singulärwertzerlegung der Matrix $\mathbf{A} = \mathbf{U} \mathbf{\Sigma} \mathbf{V}^H$ zurück. Wir können bemerken, dass: \begin{itemize} \item Wenn in $\mathbf{\Sigma}$ Nullen auf der Diagonalen sind, beziehen sie sich auf jene Würfel, durch die der Röntgenstrahl nicht ging. \item Wenn etwas sehr Kleines vorhanden ist, bezieht es sich auf solche Experimente, bei denen der Röntgenstrahl nur einige Male und möglicherweise nur streifend in einen Würfel oder eine Menge von Würfeln gelangte. \end{itemize} Das heißt, wir sollten solche Daten nicht berücksichtigen, aber sie sind sehr schwer bereits bei der Bildung der Matrix $\mathbf{A}$ herauszufiltern. Aber da diese Daten wenig informativ sind, lass sie uns einfach ``wegwerfen''! Das heißt, führen wir diese SVD aus und nullen die kleinen Diagonalwerte in $\mathbf{\Sigma}$ und rekonstruieren $\mathbf{A}$. Das ist natürlich gut, aber wir wissen immer noch nicht, wie wir ein solches Problem lösen sollen, da $\mathbf{A}^H \mathbf{A}$ ebenfalls Null-Singulärwerte enthalten wird. Aber wenn man für ein quadratisches nicht entartetes System $\mathbf{A}\vec{x} = \vec{b}$ die Lösung darstellen kann als: \[ \mathbf{A} = \mathbf{U} \mathbf{\Sigma} \mathbf{V}^H \quad \Rightarrow \quad \vec{x} = \mathbf{V} \mathbf{\Sigma}^{-1} \mathbf{U}^H \vec{b} \] dann könnten wir unser Problem der kleinsten Quadrate wohl analog lösen. Zwar werden wir dort, wo $\mathbf{\Sigma}$ Nullen hat, in $\mathbf{\Sigma}^{-1}$ sie nicht invertieren, sondern dort ebenfalls eine Null eintragen. Dann wäre es falsch, sie invers zu nennen, und wir nennen sie \textbf{Pseudoinverse}, bezeichnet als $\mathbf{\Sigma}^+$ (oder $\mathbf{A}^+$ für die gesamte Matrix). \subsection{Das Größenproblem: wenn die SVD nicht in den Speicher passt} Alles wäre gut, nur ist in realen Problemen die Matrix $\mathbf{A}$ oft riesig. Häufig werden ihre Elemente analytisch berechnet, und es besteht keine Notwendigkeit, sie zu speichern. Und die dichte Matrix $\mathbf{U}$, deren Dimension gleich der von $\mathbf{A}$ ist, passt normalerweise nicht in den Speicher. Zum Beispiel verwendet man für einen gewöhnlichen CT-Scan Würfel von etwa 1 mm Größe. Pro Sekunde erfolgen etwa 100--500 Messungen an 100--500 Empfängern, und die gesamte Messung dauert etwa eine Minute. Das heißt, $N$ und $M$ können leicht in der Größenordnung von 10--100 Millionen liegen, und die vollständige Singulärmatrix würde 100 Petabyte Speicher erfordern, was unvernünftig viel ist. \subsection{Tikhonov-Regularisierung} Man kann bemerken, dass, wenn man zur entarteten Matrix $\mathbf{A}^H \mathbf{A}$ die Einheitsmatrix $\lambda \mathbf{I}$ addiert, sodass $\lambda$ im Bereich jener sehr kleinen Singulärwerte liegt, die wir genullt haben, dann diese Addition geschrieben werden kann als: \begin{align*} \mathbf{A}^H \mathbf{A} + \lambda \mathbf{I} &= \mathbf{V} \mathbf{\Sigma} \mathbf{U}^H \mathbf{U} \mathbf{\Sigma} \mathbf{V}^H + \lambda \mathbf{V} \mathbf{V}^H \\ &= \mathbf{V} (\mathbf{\Sigma}^2 + \lambda \mathbf{I}) \mathbf{V}^H \end{align*} Und wir sehen, dass wir dieses Lambda zu jedem Singulärwert addieren und so eine entartete oder schlecht konditionierte Matrix in eine gut konditionierte ``verwandeln''. Darüber hinaus ändern wir bei kleinem $\lambda$ (in der Größenordnung der kleinen Singulärwerte) die großen Singulärwerte fast nicht, ersetzen aber die kleinen einfach durch $\lambda$. Wenn wir dann das System mit einer solchen Matrix lösen, ``verderben'' wir die großen Singulärwerte einfach nicht und reduzieren die Antwort der kleinen Singulärwerte wesentlich. Also ist die Kondition der Matrix wesentlich kleiner geworden, und wir ``verderben'' die Lösung nicht durch die Matrix. Da wir wissen, dass die Matrix dünnbesetzt ist, können wir einige iterative Methoden anwenden (sogar das Verfahren der konjugierten Gradienten) und sehr schnell konvergieren. Tatsächlich ist ein solches $\lambda$ der sogenannte \textbf{Tikhonov-Regularisierer}, denn wir können statt des Ausgangsproblems das folgende lösen: \[ \min_{\vec{x}} \|\mathbf{A}\vec{x} - \vec{b}\|_2^2 + \lambda \|\vec{x}\|_2^2 \] wobei wir faktisch verlangen, dass sowohl das Residuum als auch die Norm der Lösung gleichzeitig minimiert werden und ihr gegenseitiges Verhältnis gerade durch den Wert dieses $\lambda$ geregelt wird. Hier gibt es viel schöne Theorie, die einst in den 1970er Jahren von Tikhonov und seinen Nachfolgern hergeleitet wurde, aber der Kern bleibt einfach: Ein solcher zusätzlicher ``Regularisierer'' hilft, \textbf{schlecht gestellte} (ill-posed) Probleme zu lösen. \subsection{Iterative Verkleinerung von $\lambda$} Darüber hinaus setzt man oft zunächst $\lambda$ ziemlich groß, dann konvergiert der konjugierte Gradient sehr schnell (die Kondition der Matrix wird sehr klein sein), und dann verkleinert man $\lambda$, wobei man jedes Mal die vorherige Lösung als Startwert verwendet. Dies erlaubt: \begin{itemize} \item Schnell eine ``grobe'' Lösung mit großem $\lambda$ zu erhalten, \item Sie allmählich zu verfeinern, indem man $\lambda$ verkleinert, \item Probleme mit der Kondition in frühen Phasen zu vermeiden. \end{itemize} \begin{successbox}[Fazit] Die Tikhonov-Regularisierung ist ein mächtiges Werkzeug zur Lösung schlecht konditionierter Probleme der kleinsten Quadrate. Sie fügt eine Strafe für die Norm der Lösung hinzu, was die Lösung stabilisiert und die Verwendung iterativer Methoden auch für entartete Systeme erlaubt. \end{successbox} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Basisfunktionen: von finiten Differenzen zu Splines} \subsection{Die Idee: das Unbekannte durch Bekanntes approximieren} Bisher haben wir mit diskreten Vektoren und Matrizen gearbeitet. Aber reale Probleme --- Differentialgleichungen, Integrale, Wellenfunktionen --- leben im kontinuierlichen Raum. Wie reduzieren wir ein unendlich-dimensionales Problem auf ein endlich-dimensionales? Antwort: \textbf{Basisfunktionen}. Wir stellen die unbekannte Funktion $f(x)$ als Linearkombination bekannter Basisfunktionen $\phi_k(x)$ dar: \[ f(x) \approx \sum_{k=1}^N c_k \phi_k(x) \] und suchen die Koeffizienten $c_k$. Dies verwandelt das Problem der Suche nach einer Funktion in das Problem der Suche nach einem Vektor $\vec{c} \in \R^N$. \subsection{Das Ritz-Verfahren (Variationsverfahren)} Einer der allgemeinsten Ansätze ist das \textbf{Ritz-Verfahren}. Angenommen, wir haben ein Funktional $E[f]$, das wir minimieren wollen (zum Beispiel die Energie in der Quantenmechanik). Wir setzen die Entwicklung ein: \[ f(x) \approx \sum_{k=1}^N c_k \phi_k(x) \] und erhalten eine Funktion der Koeffizienten: \[ E(\vec{c}) = E\left[\sum_{k=1}^N c_k \phi_k\right] \] Nun minimieren wir $E(\vec{c})$ über $\vec{c}$ --- das ist bereits ein gewöhnliches Optimierungsproblem in $\R^N$. \subsection{Finite Differenzen (Finite Difference Method, FDM)} Die einfachste Wahl von Basisfunktionen sind \textbf{lokale Konstanten} oder \textbf{lineare Funktionen} auf einem äquidistanten Gitter. \subsubsection{Ein Beispiel aus der Chemie: Diffusion} Angenommen, wir haben die Diffusionsgleichung für die Konzentration eines Stoffes $C(x,t)$: \[ \frac{\partial C}{\partial t} = D \frac{\partial^2 C}{\partial x^2} \] Wir diskretisieren den Raum mit Schrittweite $h$: $x_j = jh$, und die Zeit mit Schrittweite $\Delta t$: $t_n = n\Delta t$. Die zweite Ableitung approximieren wir durch die zentrale Differenz: \[ \frac{\partial^2 C}{\partial x^2}\bigg|_{x_j} \approx \frac{C_{j+1} - 2C_j + C_{j-1}}{h^2} \] wobei $C_j^n \approx C(x_j, t_n)$. Dann das Euler-Schema in der Zeit: \[ \frac{C_j^{n+1} - C_j^n}{\Delta t} = D \frac{C_{j+1}^n - 2C_j^n + C_{j-1}^n}{h^2} \] Dies ergibt eine explizite Rekursionsformel: \[ C_j^{n+1} = C_j^n + \frac{D\Delta t}{h^2} (C_{j+1}^n - 2C_j^n + C_{j-1}^n) \] \begin{warningbox}[Stabilität] Das explizite Schema ist nur stabil für $\frac{D\Delta t}{h^2} \le \frac{1}{2}$. Wenn der Zeitschritt zu groß ist --- ``zerfällt'' die Lösung. \end{warningbox} \subsection{Finite Elemente (Finite Element Method, FEM)} Ein komplizierterer, aber flexiblerer Ansatz sind \textbf{finite Elemente}. Wir unterteilen das Gebiet in kleine Elemente (Dreiecke, Tetraeder) und approximieren auf jedem Element die Funktion durch \textbf{lokale Polynome} (üblicherweise linear oder quadratisch). Die Basisfunktionen sind \textbf{stückweise lineare ``Hüte''} (hat functions), von denen jede an einem Knoten gleich 1 und an allen anderen gleich 0 ist. Vorteile der FEM: \begin{itemize} \item Kann mit komplexen Geometrien arbeiten, \item Kann das Gitter anpassen (verdichten, wo sich die Lösung schnell ändert), \item Gut theoretisch fundiert. \end{itemize} \subsection{Basisfunktionen in der Quantenchemie: Gauß-Orbitale} Und nun --- das Interessanteste für Chemiker. In der Quantenchemie (Hartree--Fock-Verfahren, DFT) lösen wir die Schrödinger-Gleichung: \[ \hat{H} \Psi = E \Psi \] wobei $\hat{H}$ der Hamilton-Operator und $\Psi$ die Wellenfunktion ist. Wir stellen die Molekülorbitale $\psi_i(\vec{r})$ als Linearkombination von \textbf{Atomorbitalen} dar (LCAO --- Linear Combination of Atomic Orbitals): \[ \psi_i(\vec{r}) = \sum_{\mu=1}^K c_{\mu i} \phi_\mu(\vec{r}) \] wobei $\phi_\mu(\vec{r})$ die an den Atomen zentrierten Basisfunktionen sind. \subsubsection{Gauß-Funktionen (Gaussian Type Orbitals, GTO)} In den meisten modernen Programmen (Gaussian, ORCA, Q-Chem) verwendet man als Basisfunktionen \textbf{Gauß-Orbitale}: \[ \phi_\mu(\vec{r}) = (x - X_A)^l (y - Y_A)^m (z - Z_A)^n \exp(-\alpha |\vec{r} - \vec{R}_A|^2) \] wobei: \begin{itemize} \item $\vec{R}_A = (X_A, Y_A, Z_A)$ die Koordinaten des Atoms $A$ sind, \item $l, m, n$ die Drehimpulsquantenzahlen ($s, p, d, f$-Orbitale) sind, \item $\alpha$ der Exponent ist (kontrolliert die ``Größe'' des Orbitals). \end{itemize} Warum gerade Gauß-Funktionen? \begin{itemize} \item \textbf{Das Produkt von Gauß-Funktionen ist wieder eine Gauß-Funktion:} Dies ist entscheidend für die Berechnung von Mehrzentrenintegralen (Elektron-Elektron-Abstoßung), \item \textbf{Analytische Integrale:} Alle Integrale werden analytisch berechnet, \item \textbf{Kontraktionen:} Reale Atomorbitale werden durch eine Summe mehrerer Gauß-Funktionen (Kontraktionen) approximiert, was die Genauigkeit erhöht. \end{itemize} \subsubsection{Das Ritz-Verfahren in der Quantenchemie} Wir setzen die LCAO-Entwicklung in das Hartree--Fock- oder DFT-Energiefunktional ein und minimieren über die Koeffizienten $c_{\mu i}$. Dies führt auf das Eigenwertproblem: \[ \mathbf{F} \vec{c}_i = \varepsilon_i \mathbf{S} \vec{c}_i \] wobei: \begin{itemize} \item $\mathbf{F}$ die Fock-Matrix (oder die Kohn--Sham-Matrix in DFT) ist, \item $\mathbf{S}$ die Überlappungsmatrix ist ($S_{\mu\nu} = \langle \phi_\mu | \phi_\nu \rangle$), \item $\varepsilon_i$ die Orbitalenergien sind, \item $\vec{c}_i$ die Entwicklungskoeffizienten des $i$-ten Molekülorbitals sind. \end{itemize} Dies ist ein verallgemeinertes Eigenwertproblem, und es wird iterativ gelöst (mit der Methode des selbstkonsistenten Feldes, SCF). \begin{tipbox}[Basissätze] Beliebte Basissätze in der Quantenchemie: \begin{itemize} \item \textbf{STO-3G:} Minimalbasis (3 Gauß-Funktionen pro Atomorbital), \item \textbf{6-31G*:} Double-Zeta-Qualität mit Polarisationsfunktionen, \item \textbf{cc-pVTZ:} Korrelationskonsistente Triple-Zeta (hohe Genauigkeit). \end{itemize} Je größer der Basissatz --- desto genauer das Ergebnis, aber desto teurer die Berechnungen ($\mathcal{O}(N^4)$ für Hartree--Fock, wobei $N$ die Anzahl der Basisfunktionen ist). \end{tipbox} \subsection{Splines: glatte stückweise polynomielle Funktionen} Und nun --- über Splines. Dies sind Basisfunktionen, die in der Signalverarbeitung, Computergrafik und numerischen Methoden weit verbreitet sind. \subsubsection{Was ist ein Spline?} Ein \textbf{Spline} ist eine stückweise polynomielle Funktion, die an den Verbindungsstellen \textbf{glatt} ist (stetige Ableitungen bis zu einer bestimmten Ordnung hat). Der beliebteste ist der \textbf{kubische Spline} (Grad 3, Glattheit $C^2$ --- die Funktion, die erste und die zweite Ableitung sind stetig). \subsubsection{B-Splines (Basis Splines)} \textbf{B-Splines} sind ein spezieller Satz von Basissplines mit \textbf{kompaktem Träger} (compact support). Jeder B-Spline ist nur auf einem kleinen Intervall von null verschieden. Vorteile: \begin{itemize} \item \textbf{Lokalität:} Die Änderung eines Koeffizienten beeinflusst nur einen kleinen Bereich, \item \textbf{Numerische Stabilität:} Die Matrizen sind gut konditioniert, \item \textbf{Rekursive Definition:} B-Splines werden rekursiv konstruiert (Formel von de Boor). \end{itemize} \subsubsection{Zweidimensionale Splines} B-Splines lassen sich leicht auf den zweidimensionalen Fall über das \textbf{Tensorprodukt} verallgemeinern: \[ B_{ij}(x,y) = B_i(x) \cdot B_j(y) \] Dies wird verwendet in: \begin{itemize} \item Bildverarbeitung (Glättung, Interpolation), \item Finiten Elementen (höhere Ordnungen), \item Computergrafik (Oberflächen). \end{itemize} \subsubsection{Bézier-Splines} \textbf{Bézier-Splines} sind parametrische Kurven, die durch \textbf{Kontrollpunkte} definiert werden. Die Kurve verläuft nicht durch die Kontrollpunkte, sondern wird von ihnen ``angezogen''. Die Formel einer Bézier-Kurve vom Grad $n$: \[ \vec{B}(t) = \sum_{i=0}^n \binom{n}{i} (1-t)^{n-i} t^i \vec{P}_i, \quad t \in [0,1] \] wobei $\vec{P}_i$ die Kontrollpunkte sind und $\binom{n}{i} (1-t)^{n-i} t^i$ die \textbf{Bernstein-Polynome}. Anwendungen: \begin{itemize} \item Computergrafik (Vektorgrafik, TrueType-Schriften), \item CAD-Systeme (AutoCAD, SolidWorks), \item Animation und Trajektorien. \end{itemize} \subsubsection{Zusammenhang mit finiten Elementen} Interessanterweise gehen B-Splines mit kompaktem Träger \textbf{fließend} in Basisfinite Elemente über. Nimmt man B-Splines vom Grad 0 --- das sind stückweise konstante Funktionen (wie in den einfachsten finiten Elementen). Grad 1 --- stückweise lineare ``Hüte''. Grad 2 und höher --- glattere Basen. Dies erlaubt die Konstruktion \textbf{isoparametrischer finiter Elemente} höherer Ordnung, die bei weniger Knoten eine genauere Approximation liefern. \begin{successbox}[Fazit] Basisfunktionen sind die Brücke zwischen der kontinuierlichen Welt der Differentialgleichungen und der diskreten Welt der Computer. Finite Differenzen sind die einfachsten, finite Elemente die flexibelsten, Gauß-Orbitale sind auf die Quantenchemie spezialisiert, und Splines sind für glatte Interpolation und Grafik. Das Verständnis ihrer Eigenschaften ist der Schlüssel zur Wahl der richtigen Methode für dein Problem. \end{successbox} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Integration: von der Analytik zu Monte Carlo} \subsection{Wozu brauchen wir Integration?} Das Integral ist eine der fundamentalsten Operationen in Mathematik und Physik. Wir berechnen ständig: \begin{itemize} \item Flächen und Volumina, \item Erwartungswerte und Varianzen, \item Normierungskonstanten in der Quantenmechanik, \item Mehrdimensionale Integrale in der statistischen Physik, \item Faltungen in der Signalverarbeitung. \end{itemize} Und wenn wir in der Schule gelernt haben, Integrale analytisch zu berechnen, dann endet die Analytik im wirklichen Leben --- besonders in Chemie und Physik --- oft sehr schnell. \subsection{Analytische Integration: wann sie funktioniert} Erinnern wir uns an die Grundlagen. Wenn wir eine analytisch gegebene Funktion $f(x)$ haben und ihre Stammfunktion $F(x)$ kennen, dann gilt: \[ \int_a^b f(x)\,dx = F(b) - F(a) \] Zum Beispiel: \[ \int_0^1 x^2\,dx = \frac{x^3}{3}\bigg|_0^1 = \frac{1}{3} \] Schön, exakt, schnell. Aber: \begin{itemize} \item Nicht alle Funktionen haben eine elementare Stammfunktion (zum Beispiel $e^{-x^2}$ --- das Fehlerintegral), \item Nicht alle Funktionen sind analytisch gegeben --- oft haben wir nur eine Menge von Punkten (experimentelle Daten), \item In mehrdimensionalen Fällen ist die Analytik fast immer machtlos. \end{itemize} \subsection{Numerische Integration in 1D: das Simpson-Verfahren} Was tun, wenn das Integral analytisch nicht berechnet werden kann? Antwort: die Funktion durch etwas Einfaches approximieren und dieses Einfache integrieren. Im vorigen Kapitel haben wir gelernt, kubische Splines zu konstruieren. Und hier konstruieren wir ebenfalls Polynome und integrieren mit solchen Stücken. \subsubsection{Rechteckverfahren} Der einfachste Ansatz: Wir teilen das Intervall $[a,b]$ in $N$ gleiche Teile der Breite $h = (b-a)/N$ und approximieren die Funktion auf jedem Teilintervall durch eine Konstante (den Wert in der Mitte): \[ \int_a^b f(x)\,dx \approx h \sum_{k=0}^{N-1} f\left(a + \left(k + \frac{1}{2}\right)h\right) \] Der Fehler ist $\mathcal{O}(h^2)$. \subsubsection{Trapezverfahren} Etwas besser: Wir approximieren die Funktion auf jedem Intervall durch ein lineares Polynom: \[ \int_a^b f(x)\,dx \approx \frac{h}{2} \left[ f(a) + 2\sum_{k=1}^{N-1} f(a+kh) + f(b) \right] \] Der Fehler ist $\mathcal{O}(h^2)$. \subsubsection{Simpson-Verfahren} Noch besser: Wir approximieren die Funktion auf jedem Paar von Intervallen durch ein \textbf{quadratisches Polynom}: \[ \int_a^b f(x)\,dx \approx \frac{h}{3} \left[ f(a) + 4\sum_{k=1,3,5,\dots}^{N-1} f(a+kh) + 2\sum_{k=2,4,6,\dots}^{N-2} f(a+kh) + f(b) \right] \] wobei $N$ eine gerade Zahl ist. Der Fehler ist $\mathcal{O}(h^4)$! Das ist bereits sehr gut für glatte Funktionen. \begin{tipbox}[Warum ist Simpson so gut?] Das Simpson-Verfahren ist im Wesentlichen die Integration stückweise quadratischer Splines. Wenn die Funktion glatt ist (stetige Ableitungen bis zur 4. Ordnung), dann nimmt der Fehler wie $h^4$ ab, was viel schneller ist als bei Trapezen. Für sehr glatte Funktionen gibt es noch genauere Methoden --- Gauß-Quadraturen, die optimal gewählte Knoten und Gewichte verwenden. \end{tipbox} \subsection{Der Fluch der Dimension} Aber was, wenn wir eine 2-, 3-, 4-, 10-dimensionale Funktion integrieren? Dann brauchen wir $N^{10}$ Integrationspunkte, wobei $N$ mindestens 100 ist... Das ist irgendwie viel! Rechnen wir. Für ein 10-dimensionales Integral mit $N = 100$ Punkten entlang jeder Achse: \[ 100^{10} = 10^{20} \text{ Punkte} \] Selbst wenn jeder Punkt in 1 Nanosekunde berechnet wird, dauert das: \[ 10^{20} \text{ ns} = 10^{11} \text{ s} \approx 3000 \text{ Jahre} \] Dies ist der \textbf{Fluch der Dimension} (curse of dimensionality). \subsubsection{Wo kommt das in der Chemie vor?} \begin{itemize} \item \textbf{Quantenchemie:} Berechnung von Mehrelektronenintegralen (Elektron-Elektron-Abstoßung) --- das sind 6-dimensionale Integrale (3 Koordinaten pro Elektron). \item \textbf{Statistische Mechanik:} Die Zustandssumme ist ein Integral über den Phasenraum aller Teilchen. Für $N$ Teilchen --- das ist ein $6N$-dimensionales Integral (3 Koordinaten + 3 Impulse pro Teilchen). \item \textbf{Molekulardynamik:} Mittelung über Konfigurationen --- mehrdimensionale Integration. \item \textbf{Maschinelles Lernen:} Bayessche Inferenz --- Integration über den Raum der Modellparameter. \end{itemize} \subsection{Die Monte-Carlo-Methode zur Integration} Und hier betritt die Monte-Carlo-Methode die Bühne --- einer der elegantesten und überraschendsten Algorithmen der numerischen Mathematik. \subsubsection{Die Idee im Kern} Stellen wir uns vor, wir wollen das Integral berechnen: \[ I = \int_{[a,b]^d} f(\vec{x})\,d\vec{x} \] wobei $d$ die Dimension ist (kann sehr groß sein). Die Monte-Carlo-Methode sagt: ``Lass uns einfach $N$ zufällige Punkte $\vec{x}_1, \vec{x}_2, \dots, \vec{x}_N$ gleichmäßig in das Integrationsgebiet werfen, die Funktion an ihnen auswerten und mitteln''. Die Formel: \[ I \approx \frac{V}{N} \sum_{k=1}^N f(\vec{x}_k) \] wobei $V = (b-a)^d$ das Volumen des Integrationsgebietes ist und $\vec{x}_k$ zufällige, in $[a,b]^d$ gleichmäßig verteilte Punkte sind. \subsubsection{Warum funktioniert das?} Dies ist einfach das Gesetz der großen Zahlen. Der Erwartungswert der Zufallsvariablen $f(\vec{X})$, wobei $\vec{X}$ in $[a,b]^d$ gleichmäßig verteilt ist, ist: \[ \E[f(\vec{X})] = \frac{1}{V} \int_{[a,b]^d} f(\vec{x})\,d\vec{x} = \frac{I}{V} \] Nach dem Gesetz der großen Zahlen konvergiert der Mittelwert $\frac{1}{N}\sum f(\vec{x}_k)$ gegen $\E[f(\vec{X})]$ für $N \to \infty$. \subsubsection{Konvergenzgeschwindigkeit} Und hier --- Magie! Der Fehler der Monte-Carlo-Methode nimmt ab wie: \[ \text{Fehler} \sim \frac{\sigma}{\sqrt{N}} \] wobei $\sigma$ die Standardabweichung der Funktion $f(\vec{x})$ im Integrationsgebiet ist. \textbf{Wichtig:} Diese Konvergenzgeschwindigkeit \textit{hängt nicht von der Dimension} $d$ ab! Zum Vergleich: \begin{itemize} \item Simpson-Verfahren in 1D: Fehler $\sim h^4 \sim N^{-4}$, \item Simpson-Verfahren im $d$-dimensionalen Fall: Fehler $\sim N^{-4/d}$ (katastrophal langsam für große $d$), \item Monte-Carlo-Methode in beliebiger Dimension: Fehler $\sim N^{-1/2}$. \end{itemize} Ja, Monte Carlo konvergiert langsamer als Simpson in 1D. Aber im 10-dimensionalen Fall: \begin{itemize} \item Simpson: $N^{-4/10} = N^{-0.4}$, \item Monte Carlo: $N^{-0.5}$. \end{itemize} Monte Carlo ist \textit{schneller}! \begin{warningbox}[Der Preis von Monte Carlo] Die Monte-Carlo-Methode konvergiert wie $N^{-1/2}$ --- das bedeutet, um die Genauigkeit um den Faktor 10 zu erhöhen, braucht man 100-mal mehr Punkte. Das ist nach absoluten Maßstäben langsam, aber es ist die \textit{einzige} Methode, die in hohen Dimensionen funktioniert. \end{warningbox} \subsection{Ein Beispiel aus der Quantenchemie: Berechnung von Orbitalen} Betrachten wir, wie die Monte-Carlo-Methode in der realen Quantenchemie angewendet wird. \subsubsection{Aufgabe: Normierung der Wellenfunktion} In der Quantenmechanik beschreibt die Wellenfunktion $\Psi(\vec{r}_1, \vec{r}_2, \dots, \vec{r}_N)$ den Zustand von $N$ Elektronen. Sie muss normiert sein: \[ \int |\Psi(\vec{r}_1, \vec{r}_2, \dots, \vec{r}_N)|^2 \,d\vec{r}_1\,d\vec{r}_2\cdots d\vec{r}_N = 1 \] Das ist ein $3N$-dimensionales Integral! Für ein Wassermolekül ($N = 10$ Elektronen) --- das ist ein 30-dimensionales Integral. \subsubsection{Anwendung von Monte Carlo} Wir erzeugen $M$ zufällige Elektronenkonfigurationen: \[ \{\vec{r}_1^{(k)}, \vec{r}_2^{(k)}, \dots, \vec{r}_{10}^{(k)}\}, \quad k = 1, \dots, M \] wobei jede Koordinate $\vec{r}_i^{(k)} = (x_i^{(k)}, y_i^{(k)}, z_i^{(k)})$ aus einer Verteilung (zum Beispiel einer Gauß-Verteilung) gewählt wird. Wir berechnen $|\Psi|^2$ in jeder Konfiguration und mitteln: \[ \int |\Psi|^2 \approx \frac{V^{30}}{M} \sum_{k=1}^M |\Psi(\vec{r}_1^{(k)}, \dots, \vec{r}_{10}^{(k)})|^2 \] \subsubsection{Variational Monte Carlo (VMC)} Aber das ist noch nicht alles! In der Quantenchemie gibt es die \textbf{Variational Monte Carlo (VMC)}-Methode, die Monte Carlo zur \textit{Minimierung} der Energie verwendet. Wir wählen eine Testwellenfunktion $\Psi_T(\vec{r}_1, \dots, \vec{r}_N; \vec{\alpha})$ mit Parametern $\vec{\alpha}$ und berechnen die Energie: \[ E(\vec{\alpha}) = \frac{\int \Psi_T^* \hat{H} \Psi_T \,d\vec{r}_1\cdots d\vec{r}_N}{\int |\Psi_T|^2 \,d\vec{r}_1\cdots d\vec{r}_N} \] Beide Integrale sind $3N$-dimensional, und wir berechnen sie mit Monte Carlo. Dann minimieren wir $E(\vec{\alpha})$ über die Parameter $\vec{\alpha}$ --- und erhalten eine approximative Wellenfunktion. \subsubsection{Quantum Monte Carlo (QMC)} Es gibt noch fortgeschrittenere Methoden --- \textbf{Diffusion Monte Carlo (DMC)}, die die Schrödinger-Gleichung direkt lösen, indem sie die ``Diffusion'' von Elektronen in imaginärer Zeit modellieren. Diese Methoden liefern \textit{fast exakte} Lösungen für kleine Moleküle und werden als Referenz zur Prüfung anderer Methoden verwendet. \begin{tipbox}[Wo wird Monte Carlo in der Chemie verwendet?] \begin{itemize} \item \textbf{Quantum Monte Carlo (QMC):} Exakte Lösung der Schrödinger-Gleichung für kleine Systeme, \item \textbf{Monte-Carlo-Molekulardynamik:} Sampling von Konfigurationen in der statistischen Mechanik, \item \textbf{Integration in DFT:} Numerische Integration des Austausch-Korrelations-Funktionals, \item \textbf{Bayessche Optimierung:} Suche nach dem globalen Minimum der Energie eines Moleküls. \end{itemize} \end{tipbox} \subsection{Verbesserungen der Monte-Carlo-Methode} Die grundlegende Monte-Carlo-Methode ist einfach eine gleichmäßige Zufallswanderung. Aber es gibt viele Verbesserungen: \subsubsection{Importance Sampling} Statt einer gleichmäßigen Verteilung verwenden wir eine Verteilung proportional zu $|f(\vec{x})|$. Dann fallen die Punkte häufiger in Bereiche, in denen die Funktion groß ist, und der Fehler verringert sich. \subsubsection{Metropolis-Verfahren (Metropolis-Hastings)} Für komplizierte mehrdimensionale Integrale verwenden wir Markov-Ketten: Jeder nächste Punkt hängt vom vorherigen ab. Dies erlaubt eine effiziente Erkundung des Raums, auch wenn er sehr groß ist. \subsubsection{Quasi-Monte-Carlo} Statt zufälliger Punkte verwenden wir \textbf{deterministische} Folgen mit geringer Diskrepanz (Sobol-, Halton-Folgen). Sie bedecken den Raum ``gleichmäßiger'' als zufällige Punkte und konvergieren schneller: Fehler $\sim N^{-1} (\log N)^d$. \begin{successbox}[Fazit] Die Monte-Carlo-Methode ist die einzige praktikable Möglichkeit zur Berechnung mehrdimensionaler Integrale. Ihre Konvergenzgeschwindigkeit hängt nicht von der Dimension ab, was sie in der Quantenchemie, der statistischen Physik und im maschinellen Lernen unverzichtbar macht. Obwohl sie in niedrigen Dimensionen langsamer konvergiert als deterministische Methoden, hat sie in hohen Dimensionen einfach keine Konkurrenz. \end{successbox} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Statistik: wie man Information vom Rauschen trennt} Bisher haben wir hauptsächlich mathematische Aufgaben betrachtet, in denen die Zahlen als bekannt gelten. Aber die experimentelle Chemie ist anders eingerichtet. Wenn wir die Konzentration eines Stoffes, die Lage einer Spektrallinie oder die Intensität eines Signals messen, gibt das Gerät uns niemals den ``wahren'' Wert mit unendlicher Genauigkeit. Jede Messung unterscheidet sich ein wenig von der anderen. Deshalb haben wir es statt mit einer einzigen Zahl mit einer \textbf{Zufallsvariablen} zu tun. Angenommen, wir messen dieselbe Größe mehrere Male. Wir erhalten: $(x_1,x_2,\ldots,x_N)=\vec{x}$, wobei $\vec{x}=\mu+\vec{\varepsilon},$ und \begin{itemize} \item $\mu$ der unbekannte wahre oder mittlere Wert ist; \item $\vec{\varepsilon}$ der zufällige Messfehler ist. \end{itemize} Diese einfache Gleichung ist eine der grundlegenden Ideen der Statistik. Statistik beginnt dort, wo wir verstehen: \textit{Wir beobachten Daten, interessieren uns aber für verborgene Größen.} \subsection{Der Mittelwert} Die einfachste Schätzung von $\mu$ ist das arithmetische Mittel: % $$ \mu = \frac{1}{N} \sum_{i=1}^N x_i. $$ % Warum gerade dieses? Wenn die Fehler einen Mittelwert nahe null haben, dann minimiert das arithmetische Mittel die Summe der quadrierten Abweichungen: $\min_{\mu} \sum_{i=1}^N |x_i - \mu|^2$, und die Lösung ist genau durch die obige Formel gegeben. Übrigens erklärt dies bereits, warum die Wiederholung eines Experiments gewöhnlich die Genauigkeit des Ergebnisses erhöht. \subsection{Varianz und Standardabweichung} Der Mittelwert sagt uns, wo das Zentrum der Daten liegt, sagt aber nichts darüber, wie stark die Ergebnisse streuen. Dafür führt man die Varianz ein --- ein Maß dafür, wie stark unsere Messungen um den wahren Wert (zu dem der Mittelwert strebt) ``streuen'': $$ \sigma^2 = \frac{1}{N}\sum_{i=1}^N |x_i - \mu|^2 $$ Ja, wir haben sie doch gerade minimiert, oder? Richtig, aber der Wert des Minimums wird genau dann null sein, wenn alle Werte $x_i$ gleich waren; andernfalls --- ist dieser Wert gerade die Varianz. Und um sie in denselben Einheiten zu messen, erinnern wir uns einfach, dass die zweite Norm eine Wurzel enthält, und schreiben $$ \sigma = \sqrt{\frac{1}{N}\sum_{i=1}^N |x_i - \mu|^2} $$ Das heißt, die Varianz ist im Wesentlichen das \textbf{Quadrat der Länge des Vektors der Abweichungen vom Mittelwert}. Die Statistik verwandelt sich wieder in lineare Algebra. \subsection{Die Normalverteilung} In vielen physikalischen und chemischen Messungen werden zufällige Fehler näherungsweise durch die Normalverteilung beschrieben: $$ p(x) = \frac{1}{\sqrt{2\pi\sigma^2}} \exp\left( -\frac{(x-\mu)^2}{2\sigma^2} \right). $$ Es ist nicht notwendig, sich diese Formel zu merken. Es ist viel wichtiger, die Bedeutung der beiden Parameter zu verstehen: $\mu$ und $\sigma$. Der erste bestimmt die Lage des Zentrums der Verteilung, der zweite --- ihre Breite, das heißt, je größer $\sigma$, desto stärker die Streuung der Messungen. \begin{tipbox}[Intuition] \begin{itemize} \item Der Mittelwert beantwortet die Frage: \textit{``Wo befindet sich das Ergebnis?''} \item Die Standardabweichung beantwortet die Frage: \textit{``Wie stark streuen die Ergebnisse?''} \end{itemize} \end{tipbox} \subsection{Warum wird der Mittelwert bei Wiederholung des Experiments genauer?} Angenommen, wir haben $N$ unabhängige Messungen gemacht. Für unabhängige Fehler ist die Varianz des Mittelwerts $\displaystyle \frac{\sigma^2}{N}$, folglich \textbf{verringert sich der Fehler des Mittelwerts wie} $\frac{\sigma}{\sqrt N}$. Deshalb erhöht die Vergrößerung der Anzahl der Messungen die Genauigkeit, aber nicht linear (wie bei Monte Carlo!). \begin{tipbox} Um den zufälligen Fehler um den Faktor $10$ zu verringern, braucht man im idealisierten Fall etwa $100$ Mal mehr unabhängige Messungen. \end{tipbox} \subsection{Zufälliger Fehler und systematischer Fehler} Hier muss eine sehr wichtige Unterscheidung getroffen werden. Wenn das Gerät $x=\mu+\varepsilon$ liefert, wobei $\varepsilon$ zufällig um null schwankt, können wiederholte Messungen den Einfluss dieses Fehlers verringern. Aber stellen wir uns vor, das Gerät überschätzt das Ergebnis systematisch: $x=\mu+b+\varepsilon$, wobei $b$ ein konstanter Versatz ist. Dann ergibt die Mittelung $\bar{x}\approx\mu+b$. Und so oft wir das Experiment auch wiederholen, der Versatz $b$ wird nicht verschwinden. \begin{warningbox}[Wichtig] Die Wiederholung von Messungen verringert den zufälligen Fehler, beseitigt aber nicht den systematischen Fehler. Die Statistik kann ein falsch kalibriertes Gerät nicht korrigieren. \end{warningbox} Dies ist einer der Gründe, warum in der Chemie Kalibrierung, Kontrollproben und unabhängige Messmethoden so wichtig sind. \subsection{Zwei Größen können zusammenhängen} Angenommen, zwei Größen werden gleichzeitig gemessen: $x_i$ und $y_i$. Zum Beispiel: % \begin{itemize} \item die Konzentration eines Stoffes und die Intensität eines Signals; \item Temperatur und Reaktionsgeschwindigkeit; \item Druck und Volumen; \item zwei spektrale Charakteristiken. \end{itemize} Uns kann die Frage interessieren: \textit{Ändern sich $x$ und $y$ übereinstimmend?} Dafür verwendet man die Kovarianz: $\displaystyle \operatorname{Cov}(x,y) = \frac{1}{N}(x-\mu_x)^T (y-\mu_y)$ Wenn große Werte von $x$ gewöhnlich großen Werten von $y$ entsprechen, ist die Kovarianz positiv. Wenn große Werte von $x$ kleinen $y$ entsprechen, ist sie negativ. Wenn es keinen linearen Zusammenhang gibt, kann die Kovarianz nahe null sein. \subsection{Korrelation} Die Kovarianz hängt von den Maßeinheiten ab. Deshalb verwendet man oft die normierte Größe: $\displaystyle \rho_{xy} = \frac{\operatorname{Cov}(x,y)} {\sigma_x\sigma_y}.$ Für sie gilt $\displaystyle -1\leq\rho_{xy}\leq1$. Ein Wert nahe $1$ bedeutet einen starken positiven linearen Zusammenhang, ein Wert nahe $-1$ --- einen starken negativen, und ein Wert nahe $0$ bedeutet das Fehlen einer ausgeprägten linearen Korrelation. Aber hier muss man sich eine sehr wichtige Sache merken: \begin{warningbox}[Korrelation bedeutet nicht Kausalität] Wenn zwei Größen korrelieren, bedeutet das noch nicht, dass die eine die andere verursacht. Die Korrelation spricht von einem statistischen Zusammenhang zwischen den Daten. Um einen kausalen Zusammenhang festzustellen, sind zusätzliche physikalische, chemische oder experimentelle Argumente nötig. \end{warningbox} \subsection{Die Kovarianzmatrix} Wenn wir nicht zwei, sondern $p$ gemessene Größen haben, $x_1,x_2,\ldots,x_p$, dann können die Kovarianzen aller Paare in einer einzigen Matrix gesammelt werden: $$ \Sigma= \begin{pmatrix} \operatorname{Cov}(x_1,x_1) & \operatorname{Cov}(x_1,x_2) & \cdots\\ \operatorname{Cov}(x_2,x_1) & \operatorname{Cov}(x_2,x_2) & \cdots\\ \vdots & \vdots & \ddots \end{pmatrix}. $$ Diese Matrix ist symmetrisch: $\Sigma=\Sigma^T$, und sie speichert nicht nur statistische Information. Sie ist ein mathematisches Objekt der linearen Algebra. Genau deshalb tauchen Eigenwerte und Eigenvektoren, denen wir bereits früher begegnet sind, wieder in der Statistik auf. \subsection{PCA: Statistik trifft lineare Algebra} Stellen wir uns vor, wir haben $N$ chemische Proben, für die jeweils $p$ Merkmale gemessen wurden. Wir erhalten eine Datenmatrix: $X\in\mathbb{R}^{N\times p}$. Einige Merkmale können stark miteinander zusammenhängen. Wenn zum Beispiel zwei gemessene Größen sich fast immer gemeinsam ändern, ist die Information über sie teilweise redundant. Es wäre wünschenswert, neue Koordinaten zu finden, in denen: \begin{itemize} \item die erste Koordinate die maximal mögliche Variation der Daten enthält; \item die zweite --- die maximal mögliche verbleibende Variation; \item die dritte --- die nächste, und so weiter. \end{itemize} Dies ist die Grundidee der \textbf{Hauptkomponentenanalyse} (Principal Component Analysis, PCA). Mathematisch ist PCA eng mit den Eigenvektoren der Kovarianzmatrix verbunden: $\displaystyle \Sigma v_i=\lambda_i v_i$. Die Eigenvektoren $v_i$ geben neue Richtungen an, und die Eigenwerte $\lambda_i$ zeigen, wie viel Variation der Daten auf die entsprechende Richtung entfällt. Wenn die ersten wenigen Eigenwerte wesentlich größer sind als die übrigen, $$ \lambda_1,\lambda_2,\ldots,\lambda_k \gg \lambda_{k+1},\ldots,\lambda_p, $$ dann kann der größte Teil der Information durch nur wenige Koordinaten dargestellt werden. Zum Beispiel, $$ 5000\text{ gemessene Parameter} \quad\longrightarrow\quad 20\text{ Hauptkomponenten}. $$ Das bedeutet nicht, dass wir chemische Information ``weggeworfen'' haben. Wir haben eine kompaktere mathematische Darstellung der Daten gefunden. Und hier erscheint wieder die SVD, der wir bereits begegnet sind: $$ X=U\Sigma V^T. $$ PCA und SVD erweisen sich als zwei Seiten derselben linear-algebraischen Idee. \subsection{Regression: wenn wir ein Modell bauen wollen} Angenommen, wir haben die Konzentration $c_i$ und das entsprechende Signal $y_i$ gemessen. Wir nehmen ein lineares Modell an: $$ y_i=ac_i+b+\varepsilon_i. $$ Da die Messungen Fehler enthalten, liegen die Punkte gewöhnlich nicht genau auf einer einzigen Geraden. Deshalb suchen wir solche $a$ und $b$, die die Summe der Fehlerquadrate minimieren: $$ \min_{a,b} \sum_{i=1}^N \left(y_i-ac_i-b\right)^2. $$ Dies ist die \textbf{Methode der kleinsten Quadrate}. Somit ist die Regression nichts völlig Neues. Wir kennen diese mathematische Konstruktion bereits: $$ \boxed{ \text{Daten} \rightarrow \text{Modell} \rightarrow \text{Residuum} \rightarrow \text{Minimierung} } $$ Genau deshalb sind lineare Algebra und Optimierung für die Statistik so wichtig. \subsection{Überanpassung} Nun entsteht ein komplizierteres Problem. Wenn es wenige Daten gibt, können wir eine sehr komplizierte Funktion wählen, die praktisch durch jeden experimentellen Punkt verläuft. Zum Beispiel kann man statt einer Geraden ein Polynom hohen Grades verwenden: $$ y=a_0+a_1x+a_2x^2+\cdots+a_kx^k. $$ Bei genügend großem $k$ kann man die vorhandenen Messungen fast perfekt beschreiben. Aber das bedeutet nicht unbedingt, dass wir das richtige physikalische Gesetz gefunden haben. Wir haben uns möglicherweise einfach an das Rauschen angepasst. Dies nennt man \textbf{Überanpassung} (overfitting). Genau deshalb interessiert uns in Statistik und maschinellem Lernen nicht nur die Frage: \textit{``Wie gut beschreibt das Modell die bekannten Daten?''} sondern auch: \textit{``Wie gut funktioniert es auf neuen Daten?''} \subsection{Training, Validierung und Test} Wenn wir ein Modell aus den vorhandenen Daten bauen, ist es nützlich, die Daten in Teile aufzuteilen: $$ \boxed{ \text{Training} \quad+\quad \text{Validierung} \quad+\quad \text{Test} } $$ Auf dem Trainingssatz wird das Modell trainiert. Der Validierungssatz wird zur Wahl der Parameter des Modells verwendet. Der Testsatz muss unabhängig bleiben und wird zur abschließenden Bewertung der Fähigkeit des Modells verwendet, auf neuen Daten zu arbeiten. Diese Idee wird besonders wichtig, wenn wir zum maschinellen Lernen kommen. \subsection{Was die Statistik wirklich zu tun versucht} Man kann sagen, dass sich die Statistik nicht so sehr mit dem ``Zählen von Mittelwerten'' befasst als mit einer allgemeineren Aufgabe: \textbf{zuverlässige Information aus unvollkommenen Daten zu extrahieren.} Wir beobachten: $$ \text{Daten} = \text{Struktur} + \text{Rauschen}. $$ Unsere Aufgabe besteht darin, aus den Daten die uns interessierende Struktur zu rekonstruieren und gleichzeitig abzuschätzen, wie sehr wir ihr vertrauen können. Im einfachsten Fall sieht das so aus: $$ x_1,x_2,\ldots,x_N \longrightarrow \bar{x}\pm\text{Unsicherheit}. $$ Im komplizierteren Fall: $$ X \longrightarrow \text{PCA} \longrightarrow \text{niedrigdimensionale Darstellung}. $$ Und noch komplizierter: $$ X \longrightarrow \text{Modell} \longrightarrow \text{Vorhersage} \longrightarrow \text{Validierung auf neuen Daten}. $$ Genau von hier aus wächst auf natürliche Weise das moderne maschinelle Lernen. \begin{successbox}[Was wichtig zu merken ist] \begin{enumerate} \item Eine reale Messung enthält einen zufälligen und möglicherweise einen systematischen Fehler. \item Der Mittelwert beschreibt den zentralen Wert der Daten. \item Die Standardabweichung beschreibt ihre Streuung. \item Der zufällige Fehler des Mittelwerts verringert sich wie $1/\sqrt N$. \item Der systematische Fehler wird durch Mittelung nicht beseitigt. \item Die Kovarianz beschreibt die gemeinsame Änderung von Größen. \item PCA sucht die informativsten Richtungen in mehrdimensionalen Daten. \item Die Regression baut aus Messungen ein mathematisches Modell. \item Ein gutes Modell muss nicht nur auf bekannten Daten, sondern auch auf neuen Daten funktionieren. \end{enumerate} \end{successbox} \subsection{Und wieder lineare Algebra} Und vielleicht die angenehmste Beobachtung für uns ist, dass sich die Statistik als kein so fremdes Fach erwiesen hat. Wir begannen mit Messungen: $\displaystyle x_1,x_2,\ldots,x_N, $ gingen zu Vektoren über: $\displaystyle \vec x, $ dann zu Normen: $\displaystyle \|\vec x\|_2, $ zu Kovarianzmatrizen: $\displaystyle \Sigma, $ zu Eigenwerten: $\displaystyle \Sigma v=\lambda v, $ zur SVD: $\displaystyle X=U\Sigma V^T, $ und schließlich zur Optimierung: $\displaystyle \min \|Ax-b\|_2^2. $ Das heißt, die Statistik zerstört unser mathematisches Bild nicht. Sie zeigt, wie dieselbe Mathematik beginnt, mit \textbf{realen, unvollkommenen und verrauschten Daten} zu arbeiten. Und das ist genau die Situation, mit der ein moderner Chemiker fast jeden Tag konfrontiert ist. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Von SVD zu Compressed Sensing: wenn weniger Daten vorhanden sind als nötig} \subsection{Eine praktische Aufgabe: Trennung von Spektren} Betrachten wir eine praktische Aufgabe. In unser Labor wurde eine Masse bunter organischer Substanzen gebracht, und zur Hand war nur ein VIS/IR-Spektrometer für Fluoreszenzspektren. Man sagte uns, dass die Proben mehrere Reaktionsprodukte mit vermutlich unterschiedlichen Spektren enthalten, und bat uns, sowohl die Konzentrationen als auch die Spektren der reinen Substanzen zu bestimmen. Wir wählten eine nicht sehr starke UV-Quelle, damit die Fluoreszenz sehr hell war und unterschiedliche Spektren für verschiedene Substanzen lieferte. Das scheint eine Aufgabe für SVD zu sein! Denn wenn wir jedes solche Spektrum in seine eigene Spalte der Matrix $\mathbf{A}$ setzen, annehmen, dass wir hypothetisch eine Matrix reiner Spektren dieser Substanzen $\mathbf{P}$ haben, und für jede Probe Konzentrationskoeffizienten der Substanzen selbst $\mathbf{Q}$, dann: \[ \mathbf{A} = \mathbf{P}\mathbf{Q} \] und wir können $\mathbf{P}$ als die linken Singulärvektoren der Matrix $\mathbf{A}$ erhalten und $\mathbf{Q}$ als die rechten Singulärvektoren, von links mit der Diagonalen multipliziert. Wir können die Singulärwertzerlegung auch immer schreiben als: \[ \forall i, j: \quad a_{ij} = \sum_{r=1}^R d_r u_{ir} v_{rj} \] Darüber hinaus, wenn wir eine einzige Substanz hätten ($R=1$), wäre es genau so. \subsection{Aber etwas ging schief...} Statt schöner reiner Spektren mit positiven Werten erhielten wir in $\mathbf{P}$ irgendwelche seltsamen Diagramme --- mit negativen Werten, Oszillationen, ``Geister''-Peaks. Tatsächlich --- alles war genau so, wie es sein sollte, die Spektren \textit{können} so nicht erhalten werden. Die Spektren selbst sind positiv definite Diagramme, also ist das Skalarprodukt eines beliebigen Paares, obwohl es sehr nahe bei null liegen kann, nicht notwendigerweise null. In der Singulärwertzerlegung verlangen wir, dass alle Spalten von $\mathbf{U}$ zueinander \textbf{orthogonal} sind. Und reale Spektren von Substanzen sind nicht orthogonal! Sie sind einfach ``ähnlich'' oder ``nicht ähnlich''. Wir brauchen etwas anderes, um diese Spektren so auseinanderzuziehen, dass unsere Mathematik sie getrennt sieht. \subsection{Die Magie der dritten Dimension: das Kruskal-Theorem} Wenn wir das Experiment so aufsetzen, dass wir 3D-Daten statt 2D wie im obigen Beispiel erhalten, dann kann mathematisch bewiesen werden, dass auch eine \textbf{nicht orthogonale} Zerlegung existiert. Wenn wir zum Beispiel Fluoreszenzspektren bei mehreren verschiedenen Anregungswellenlängen aufnehmen können. Dann haben wir bereits dreidimensionale Daten: \[ \forall i, j, k: \quad a_{ijk} = \sum_{r=1}^R \alpha_r \, b_{ir} \, c_{jr} \, d_{kr} \] und nach dem \textbf{Kruskal-Theorem} (Kruskal, 1977) haben wir keine Anforderungen an die Orthogonalität jeder dieser Matrizen! \subsubsection{Was das Kruskal-Theorem besagt} Das Kruskal-Theorem ist einer der Eckpfeiler der mehrdimensionalen Datenanalyse. Es gibt Bedingungen für die \textbf{Eindeutigkeit} der Tensorzerlegung (der sogenannten CANDECOMP/PARAFAC-Zerlegung). Für einen Tensor 3. Ordnung $\mathcal{A} = \sum_{r=1}^R \vec{a}_r \circ \vec{b}_r \circ \vec{c}_r$ (wobei $\circ$ das äußere Produkt von Vektoren ist), ist die Zerlegung \textbf{eindeutig} bis auf Permutation und Skalierung der Komponenten, wenn: \[ k_A + k_B + k_C \geq 2R + 2 \] wobei $k_A$ der \textbf{k-Rang} (Kruskal-Rang) der Matrix $\mathbf{A}$ ist, also die maximale Zahl $k$, sodass jede Teilmenge von $k$ Spalten der Matrix $\mathbf{A}$ linear unabhängig ist. \begin{tipbox}[Warum ist das so wichtig?] In der gewöhnlichen SVD haben wir \textit{unendlich viele} Zerlegungen $\mathbf{A} = \mathbf{U}\mathbf{\Sigma}\mathbf{V}^H$ --- jede orthogonale Transformation im Inneren ergibt eine neue gültige Zerlegung. Deshalb kann die SVD nicht die ``physikalischen'' Komponenten isolieren --- nur mathematisch bequeme (orthogonale). Das Kruskal-Theorem sagt: Für Tensoren 3. und höherer Ordnung ist bei Erfüllung der Bedingung an die k-Ränge die Zerlegung \textit{eindeutig}! Das bedeutet, wir können genau die physikalischen Komponenten wiederherstellen, die in den Daten waren --- ohne die Forderung nach Orthogonalität. \end{tipbox} \subsubsection{Lösungsmethoden: PARAFAC und ALS} In der Praxis sucht man die Tensorzerlegung mit iterativen Methoden. Die beliebteste ist \textbf{ALS} (Alternating Least Squares): \begin{enumerate} \item Wir fixieren $\mathbf{B}$ und $\mathbf{C}$ und lösen das Problem der kleinsten Quadrate für $\mathbf{A}$, \item Wir fixieren $\mathbf{A}$ und $\mathbf{C}$ und lösen für $\mathbf{B}$, \item Wir fixieren $\mathbf{A}$ und $\mathbf{B}$ und lösen für $\mathbf{C}$, \item Wir wiederholen bis zur Konvergenz. \end{enumerate} Dies ist ein klassisches nichtlineares Minimierungsproblem (wir haben es im Kapitel über nichtlineare Operatoren besprochen), und es kann in lokalen Minima hängen bleiben. Deshalb startet man ALS in der Praxis viele Male mit verschiedenen Anfangsapproximationen. \begin{orangebox}[Anwendbarkeit in der Chemie] PARAFAC/CANDECOMP wird weit verbreitet verwendet in: \begin{itemize} \item Fluoreszenzspektroskopie (EEM --- Excitation-Emission Matrix), \item Chromatographie-Massenspektrometrie, \item NMR-Spektroskopie, \item Analyse metabolischer Daten (Metabolomik). \end{itemize} \end{orangebox} \subsection{Mehrdimensionale NMR-Spektren} Und nun --- das Interessanteste für Strukturchemiker. In der NMR-Spektroskopie von Proteinen kann man sehr mehrdimensionale Spektren konstruieren: 2D, 3D, 4D und sogar 5D. Jede Dimension ist ihre eigene Frequenz (die chemische Verschiebung eines bestimmten Kerns: $^1$H, $^{13}$C, $^{15}$N). Kreuzpeaks in solchen Spektren zeigen, welche Kerne im Raum nahe beieinander liegen, was die Rekonstruktion der 3D-Struktur des Proteins erlaubt. Aber hier ist das Problem: Wenn wir ein 4D-Spektrum mit Auflösung $1024 \times 256 \times 256 \times 64$ Punkten aufnehmen wollen, dann ist die Gesamtzahl der Punkte $\sim 4 \times 10^9$. Und jede Messung braucht Zeit (um die Relaxation abzuwarten), und insgesamt sind das \textbf{Wochen} Spektrometerbetrieb! Es gab Zeiten, in denen man, um alle Kreuzpeaks von Ubiquitin (ein Peptid aus 76 Aminosäuren) zu bestimmen, etwa \textbf{drei Wochen} ein Spektrum aufnehmen musste. Das Protein konnte in dieser Zeit degradieren! \subsection{Sparse Sampling: weniger bedeutet mehr} Aber man kann nur zufällig dünnbesetzte eindimensionale Spektren aufnehmen --- und zwar nicht 10\%, nicht 1\%, sondern wirklich sehr, sehr wenige, und dieselbe Zerlegung anwenden, nur für dünnbesetzte Daten (sparse data), und ebenso zuverlässige Ergebnisse erhalten! Indem der Autor vor 20 Jahren zusammen mit seinen Kollegen solche Methoden anwendete, beschleunigte er die Aufnahme eines solchen mehrdimensionalen Spektrums ohne Informationsverlust von drei Wochen auf \textbf{15 Minuten}. Die Idee ist einfach: Wenn wir wissen, dass das Spektrum in einer gewissen Basis \textit{dünnbesetzt} ist (das heißt, aus einer kleinen Zahl von Peaks besteht), dann müssen wir nicht alle Punkte messen. Es genügt, eine zufällige Teilmenge zu messen und dann die fehlenden Punkte durch Optimierung \textit{wiederherzustellen}. \subsection{Compressed Sensing: eine Revolution in der Messtechnik} Dies bringt uns zu einer der schönsten Ideen in der modernen Mathematik und Signalverarbeitung --- \textbf{Compressed Sensing} (komprimierte Messung, oder compressive sampling). \subsubsection{Klassische Theorie: Nyquist--Shannon} Die traditionelle Theorie der Diskretisierung (Nyquist--Shannon) besagt: Um ein Signal mit maximaler Frequenz $f_{\max}$ wiederherzustellen, muss man es mit einer Frequenz von mindestens $2f_{\max}$ messen. Für ein 4D-NMR-Spektrum mit Auflösung $1024 \times 256 \times 256 \times 64$ bedeutet dies, dass wir \textit{verpflichtet} sind, alle $4 \times 10^9$ Punkte zu messen. Weniger --- unmöglich, sonst gibt es Aliasing (Frequenzüberlappung). \subsubsection{Die Revolution: Compressed Sensing} Aber in den Jahren 2004--2006 zeigten Donoho, Candès, Tao und andere Mathematiker: Wenn ein Signal in einer gewissen Basis \textbf{dünnbesetzt} (sparse) ist, dann kann es aus einer \textbf{wesentlich geringeren} Anzahl von Messungen wiederhergestellt werden! Formal: Das Signal $\vec{x} \in \R^N$ habe eine dünnbesetzte Darstellung $\vec{x} = \mathbf{\Psi}\vec{s}$, wobei $\vec{s}$ nur $K \ll N$ von Null verschiedene Elemente hat. Dann können wir $\vec{y} = \mathbf{\Phi}\vec{x}$ messen, wobei $\mathbf{\Phi} \in \R^{M \times N}$ die Messmatrix mit $M \sim K \log(N/K) \ll N$ ist, und $\vec{x}$ durch Lösung des Problems rekonstruieren: \[ \min_{\vec{s}} \|\vec{s}\|_1 \quad \text{unter der Bedingung} \quad \vec{y} = \mathbf{\Phi}\mathbf{\Psi}\vec{s} \] \begin{tipbox} \textbf{Warum die L1-Norm und nicht L0?} Die Idee ist, dass wir die \textit{am dünnsten besetzte} Lösung finden wollen (minimiere $\|\vec{s}\|_0$ --- die Anzahl der von Null verschiedenen Elemente). Aber die Minimierung der L0-Norm ist ein NP-schweres Problem (Durchmusterung aller Kombinationen). Die Magie von Compressed Sensing ist, dass unter bestimmten Bedingungen an die Matrix $\mathbf{\Phi}$ (die sogenannte Restricted Isometry Property, RIP) die Minimierung der L1-Norm \textit{genau dasselbe Ergebnis} liefert wie die Minimierung von L0! Und die L1-Minimierung ist ein konvexes Problem, das effizient gelöst wird (es ist lineare Programmierung). \end{tipbox} \subsubsection{Anwendbarkeitsbedingungen} Für eine erfolgreiche Rekonstruktion sind zwei Bedingungen nötig: \begin{enumerate} \item \textbf{Dünnbesetztheit:} Das Signal muss in einer gewissen Basis dünnbesetzt sein (zum Beispiel ist ein NMR-Spektrum eine Menge von Delta-Funktionen, also im Frequenzbereich sehr dünnbesetzt). \item \textbf{Inkohärenz:} Die Messmatrix $\mathbf{\Phi}$ muss ``inkohärent'' mit der Basis $\mathbf{\Psi}$ sein. In der Praxis bedeutet dies, dass die Messpunkte \textbf{zufällig} (oder pseudozufällig) gewählt werden müssen. \end{enumerate} \subsection{Beispiele für Compressed Sensing in der Chemie und darüber hinaus} \subsubsection{Schnelle MRT (Magnetic Resonance Imaging)} Eine der bekanntesten Anwendungen ist die Beschleunigung der MRT in der Medizin. Die klassische MRT erfordert eine langwierige Abtastung (der k-Raum wird zeilenweise gefüllt). Mit Compressed Sensing kann man nur \textbf{zufällige} Zeilen des k-Raums messen und das vollständige Bild durch L1-Minimierung rekonstruieren. Dies verkürzt die Scanzeit um das 5--10-fache --- entscheidend für Patienten, die nicht lange still liegen können. \subsubsection{Dünnbesetzte NMR-Spektroskopie} In der mehrdimensionalen NMR von Proteinen erlaubt Compressed Sensing: \begin{itemize} \item Nur eine zufällige Teilmenge der Punkte im mehrdimensionalen k-Raum zu messen, \item Das vollständige Spektrum durch L1-Minimierung zu rekonstruieren, \item Die Experimentzeit von Wochen auf Stunden oder Minuten zu verkürzen. \end{itemize} Das ist genau das, was der Autor des Skripts vor 20 Jahren tat und was heute zum Standard in modernen NMR-Spektrometern geworden ist (Methoden wie sparse sampling, Poisson-gap sampling usw.). \subsubsection{Single-Pixel-Kamera (Rice-Kamera)} Ein erstaunliches Beispiel: eine Kamera, die \textbf{keine Pixelmatrix} hat, sondern nur einen einzigen Detektor. Das Bild wird durch komprimierte Messungen mit einem DMD (Digital Micromirror Device) rekonstruiert. Dies funktioniert, weil reale Bilder in der Wavelet-Basis dünnbesetzt sind. \subsubsection{Massenspektrometrie} In der Massenspektrometrie wird Compressed Sensing zur Beschleunigung von FT-ICR (Fourier Transform Ion Cyclotron Resonance) und Orbitrap verwendet --- man kann das FID kürzer messen und dennoch eine hohe Auflösung erhalten. \subsubsection{Datenkompression} JPEG verwendet die diskrete Kosinustransformation (ein naher Verwandter der Fourier-Transformation), und Bilder sind in dieser Basis dünnbesetzt --- deshalb komprimiert JPEG so gut. Compressed Sensing ist die nächste Stufe: Wir komprimieren nicht nur bereits gemessene Daten, sondern messen gleich weniger. \subsection{Zusammenhang mit früheren Kapiteln} Betrachten wir, wie Compressed Sensing alles verbindet, was wir durchgegangen sind: \begin{itemize} \item \textbf{Lineare Algebra:} Wir lösen ein überbestimmtes System $\mathbf{\Phi}\mathbf{\Psi}\vec{s} = \vec{y}$ durch L1-Norm-Minimierung. \item \textbf{Kondition:} Die Matrix $\mathbf{\Phi}$ muss gut konditioniert sein und die Restricted Isometry Property (RIP) erfüllen. \item \textbf{SVD und Tensoren:} Für mehrdimensionale Daten verwenden wir Tensorzerlegungen (PARAFAC), die Eindeutigkeit ohne Orthogonalität liefern. \item \textbf{Nichtlineare Optimierung:} L1-Minimierung ist ein konvexes, aber nicht glattes Problem (im Nullpunkt gibt es keine Ableitung). Wir verwenden Methoden wie ISTA (Iterative Shrinkage-Thresholding Algorithm) oder ADMM. \item \textbf{FFT:} Der Operator $\mathbf{\Psi}$ ist oft die Fourier-Transformation, und wir verwenden FFT zur schnellen Multiplikation. \end{itemize} \begin{successbox}[Hauptschlussfolgerung] Compressed Sensing ist nicht nur ein weiterer Algorithmus. Es ist eine \textit{Philosophie des Messens}: Wenn wir wissen, dass das Signal einfach (dünnbesetzt) ist, dann können wir es wesentlich weniger messen, als die klassische Theorie verlangt. Dies hat MRT, NMR-Spektroskopie, Massenspektrometrie und viele andere Bereiche revolutioniert. Und all dies basiert auf schöner Mathematik: konvexe Optimierung, Wahrscheinlichkeitstheorie (Zufallsmatrizen) und lineare Algebra. \end{successbox} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Informationskompression: von Archiven bis JPEG} \subsection{Zwei Welten der Kompression} Informationskompression ist eine der praktischsten Bereiche der angewandten Mathematik. Wir komprimieren und dekomprimieren ständig Daten: Archive, Bilder, Videos, Musik. Aber nicht alle Kompression ist gleich. Es gibt zwei grundlegend verschiedene Ansätze: \begin{itemize} \item \textbf{Verlustfreie Kompression (lossless):} Wir können die Ausgangsdaten \textit{Bit für Bit} wiederherstellen. Beispiele: ZIP, PNG, FLAC. \item \textbf{Verlustbehaftete Kompression (lossy):} Wir opfern einen Teil der Information für eine stärkere Kompression. Beispiele: JPEG, MP3, MP4. \end{itemize} \subsection{Shannons Informationsentropie} Bevor wir über Methoden sprechen, verstehen wir, \textit{wie viel} Daten überhaupt komprimiert werden können. Claude Shannon führte 1948 den Begriff der \textbf{Informationsentropie} ein. Wenn wir ein Alphabet aus $n$ Symbolen haben und Symbol $i$ mit Wahrscheinlichkeit $p_i$ vorkommt, dann ist die Entropie (die durchschnittliche Informationsmenge pro Symbol): \[ H = -\sum_{i=1}^n p_i \log_2 p_i \quad \text{[Bit]} \] Dies ist die \textbf{theoretische Grenze} der verlustfreien Kompression. Wir können Daten nicht stärker komprimieren als $H$ Bit pro Symbol. \begin{tipbox}[Beispiel] In englischem Text kommt der Buchstabe ``e'' am häufigsten vor ($p \approx 0.13$), ``z'' --- sehr selten ($p \approx 0.001$). Wenn wir alle Buchstaben gleich kodieren würden (8 Bit pro Symbol), würden wir zu viele Bits für seltene Buchstaben verschwenden. Die Entropie englischen Textes beträgt etwa 4.2 Bit pro Symbol, also können wir ihn theoretisch etwa um den Faktor 2 komprimieren. \end{tipbox} \subsection{Verlustfreie Kompression: Huffman-Kodierung} Die Idee ist einfach: Wir bemerken, dass einige Symbole und Sequenzen von 2--3 Symbolen häufiger vorkommen, legen eine Tabelle solcher Symbole an, und je häufiger ein Symbol oder eine Sequenz vorkommt, desto weniger Bits verwenden wir zu seiner Kodierung. \subsubsection{Der Huffman-Algorithmus} \begin{enumerate} \item Wir zählen die Häufigkeiten aller Symbole im Text. \item Wir bauen einen Binärbaum: Wir vereinigen die zwei seltensten Symbole zu einem Knoten mit der Summenhäufigkeit und wiederholen, bis wir einen Baum erhalten. \item Der Kode eines Symbols ist der Pfad von der Wurzel zum Blatt (0 --- links, 1 --- rechts). \end{enumerate} Häufige Symbole landen näher an der Wurzel --- ihre Kodes sind kürzer. Seltene --- weiter weg, Kodes länger. \begin{tipbox}[Beispiel] Wir haben die Symbole A (Häufigkeit 0.5), B (0.25), C (0.125), D (0.125). Der Huffman-Baum ergibt die Kodes: \begin{itemize} \item A: 0 (1 Bit) \item B: 10 (2 Bit) \item C: 110 (3 Bit) \item D: 111 (3 Bit) \end{itemize} Mittlere Länge: $0.5 \times 1 + 0.25 \times 2 + 0.125 \times 3 + 0.125 \times 3 = 1.75$ Bit pro Symbol --- nahe an der Entropie! \end{tipbox} \subsubsection{Wörterbuchmethoden: LZ77, LZ78, LZW} Algorithmen der LZ-Familie (Lempel--Ziv) gehen weiter: Sie suchen nicht nur häufige Symbole, sondern auch \textbf{sich wiederholende Sequenzen}. Idee: Wenn wir die Zeichenkette ``abrakadabra'' bereits gesehen haben, dann schreiben wir sie beim erneuten Auftreten nicht neu, sondern verweisen auf das vorherige Vorkommen: ``(zurück 11, Länge 11)''. Dies ist die Grundlage der Formate ZIP, GZIP, PNG. \subsection{Ein chemisches Beispiel: Suche nach Proteinsequenzen} Und nun --- ein lebendiges Beispiel aus der Biologie. Wir haben eine Datenbank von Proteinsequenzen (Millionen von Aminosäuren) und müssen herausfinden, ob sie Homologe (ähnliche Proteine) zu unserem neuen Protein enthält. Der direkte Vergleich ``unser Protein vs. jedes Protein in der Datenbank'' ist zu langsam. Wir brauchen einen cleveren Algorithmus für Kompression und Suche. \subsubsection{BLAST (Basic Local Alignment Search Tool)} BLAST ist einer der meistzitierten Algorithmen in der Biologie. Die Idee: \begin{enumerate} \item Wir zerlegen unser Protein in kurze ``Wörter'' (k-Mere, üblicherweise $k = 3$ für Proteine). \item Für jedes Wort suchen wir ähnliche Wörter in der Datenbank (mit kleinen Mutationen). \item Wir erweitern die Übereinstimmungen in beide Richtungen, bis die Ähnlichkeit unter eine Schwelle fällt. \end{enumerate} Dies ist im Wesentlichen eine Wörterbuchkompressionsmethode + schnelle Suche über eine Hash-Tabelle. BLAST erlaubt es, Homologe in Sekunden zu finden, obwohl ein vollständiges Alignment Stunden dauern würde. \begin{orangebox}[Warum ist das wichtig?] BLAST und seine Varianten (PSI-BLAST, BLASTP, BLASTN) sind die Grundlage der modernen Bioinformatik. Ohne sie gäbe es weder Genomentschlüsselung noch Wirkstoffentwicklung noch Evolutionsforschung. \end{orangebox} \subsection{Verlustbehaftete Kompression: JPEG und Niedrigrang-Approximation} Nun --- verlustbehaftete Kompression. Die Idee: Wir opfern einen Teil der Information, den das menschliche Auge (oder Ohr) sowieso nicht bemerkt. \subsubsection{JPEG über die diskrete Kosinustransformation (DCT)} JPEG funktioniert so: \begin{enumerate} \item Wir zerlegen das Bild in Blöcke von $8 \times 8$ Pixeln. \item Auf jeden Block wenden wir die \textbf{diskrete Kosinustransformation} (DCT) an --- ein naher Verwandter der Fourier-Transformation. \item Die DCT zerlegt den Block in 64 Frequenzkomponenten (von niedrigen Frequenzen --- der allgemeine Hintergrund, bis zu hohen Frequenzen --- feine Details). \item Wir quantisieren die Koeffizienten: Wir teilen durch eine Quantisierungsmatrix und runden. Hohe Frequenzen (feine Details) werden gröber quantisiert --- wir ``verlieren'' sie. \item Die verbleibenden von Null verschiedenen Koeffizienten werden verlustfrei komprimiert (Huffman). \end{enumerate} \subsubsection{Niedrigrang-Approximation über SVD} Eine alternative Sicht auf die Bildkompression ist über SVD. Wir haben die Bildmatrix $\mathbf{A} \in \R^{M \times N}$. SVD: \[ \mathbf{A} = \mathbf{U} \mathbf{\Sigma} \mathbf{V}^H = \sum_{k=1}^{\min(M,N)} \sigma_k \vec{u}_k \vec{v}_k^H \] Wir können $\mathbf{A}$ durch die ersten $r$ Singulärwerte approximieren: \[ \mathbf{A}_r = \sum_{k=1}^r \sigma_k \vec{u}_k \vec{v}_k^H \] Dies ist die \textbf{Niedrigrang-Approximation} vom Rang $r$. Nach dem Eckart--Young-Theorem ist dies die \textit{beste} Approximation vom Rang $r$ im Sinne der L2-Norm. \begin{tipbox}[Kompressionsrate] Ausgangsmatrix: $M \times N$ Zahlen. Approximation vom Rang $r$: $r(M + N + 1)$ Zahlen. Wenn $r \ll \min(M,N)$, ist die Kompression enorm! Zum Beispiel für ein Bild $1000 \times 1000$ und $r = 50$: Kompression um $1000 \times 1000 / (50 \times 2001) \approx 10$ Mal. \end{tipbox} \subsection{Wavelets: lokale Frequenzen} Sowohl DCT als auch SVD sind globale Transformationen: Sie zerlegen das ganze Bild in Frequenzen. Aber Bilder haben \textit{lokale} Besonderheiten: Kanten von Objekten, Texturen, Punkte. \textbf{Wavelets} sind Basisfunktionen, die sowohl im Raum als auch in der Frequenz lokalisiert sind. Sie erlauben die Analyse eines Signals auf verschiedenen Skalen. \subsubsection{Die Wavelet-Transformation} Statt Sinuskurven (wie bei Fourier) verwenden wir ``kleine Wellen'' (wavelets) --- Funktionen, die schnell abklingen. Die beliebteste ist das \textbf{Haar-Wavelet}: \[ \psi(t) = \begin{cases} 1, & 0 \leq t < 0.5 \\ -1, & 0.5 \leq t < 1 \\ 0, & \text{sonst} \end{cases} \] Die Wavelet-Transformation zerlegt ein Signal nach Skalen (Frequenzen) und Positionen. Dies erlaubt: \begin{itemize} \item Bilder besser zu komprimieren als JPEG (das Format JPEG2000 verwendet Wavelets), \item Kanten zu detektieren (Kanten von Objekten), \item Rauschen zu entfernen und wichtige Details zu bewahren. \end{itemize} \begin{successbox}[Fazit] Informationskompression ist ein Gleichgewicht zwischen mathematischer Theorie (Entropie, SVD, Wavelets) und Psychophysik (was das Auge sieht, was das Ohr hört). Von Archiven bis JPEG --- überall eine Idee: Struktur in den Daten finden und sie für eine kompakte Darstellung nutzen. \end{successbox} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Bildvergleich: von der Faltung zu neuronalen Netzen} \subsection{Die Aufgabe: ein Objekt in einem Bild finden} Wie vergleicht man zwei Bilder und versteht, dass sie dasselbe Objekt enthalten? Wir haben bereits zwei Sequenzen verglichen, indem wir eine zirkulante Matrix konstruiert haben (Kapitel über FFT). Wahrscheinlich kann man Bilder genauso vergleichen? Sie seien 2D. Aber nicht alles ist hier so gut, wie es scheint. Wenn ein Bild relativ zum anderen einfach \textbf{verschoben} ist entlang einer oder beider Achsen --- ja, alles wird funktionieren. Wir können die 2D-Faltung (über FFT) verwenden und das Maximum finden --- das ist die Position des Objekts. Aber wenn es eine \textbf{Drehung} gibt? Oder eine \textbf{Streckung}? Oder eine \textbf{Änderung der Beleuchtung}? Einfache Faltung hilft nicht mehr. \subsection{Edge Detection: Kantenerkennung} Der erste Schritt zum Verständnis des Bildinhalts ist die Erkennung von \textbf{Kanten} (edges). Kanten sind Stellen, an denen sich die Pixelintensität stark ändert. \subsubsection{Der Sobel-Operator} Der einfachste Weg ist die Berechnung des Intensitätsgradienten. Der Sobel-Operator verwendet zwei $3 \times 3$-Masken zur Berechnung der Ableitungen nach $x$ und $y$: \[ G_x = \begin{pmatrix} -1 & 0 & 1 \\ -2 & 0 & 2 \\ -1 & 0 & 1 \end{pmatrix} * I, \quad G_y = \begin{pmatrix} -1 & -2 & -1 \\ 0 & 0 & 0 \\ 1 & 2 & 1 \end{pmatrix} * I \] wobei $*$ die Faltung und $I$ das Bild ist. Gradientenbetrag: $G = \sqrt{G_x^2 + G_y^2}$. Wo $G$ groß ist --- dort ist eine Kante. \subsubsection{Canny-Kantendetektor} Ein fortgeschrittenerer Algorithmus: \begin{enumerate} \item Wir glätten das Bild mit einem Gauß-Filter (Rauschen entfernen). \item Wir berechnen den Gradienten (wie bei Sobel). \item Wir unterdrücken Nichtmaxima (nur die stärksten Kanten behalten). \item Wir wenden Hysterese an: Wenn eine Kante über der oberen Schwelle liegt --- behalten, unter der unteren --- entfernen, dazwischen --- behalten, wenn sie mit einer starken Kante verbunden ist. \end{enumerate} Das Ergebnis --- dünne, zusammenhängende Linien von Kanten. \subsection{Keypoints: Schlüsselpunkte eines Bildes} Kanten sind gut, aber wir brauchen etwas ``punktartigeres'', um Bilder zu vergleichen. Dafür sucht man \textbf{Keypoints} (interest points) --- Stellen, die leicht zu finden und stabil gegenüber Transformationen sind. \subsubsection{Harris-Eckendetektor} Idee: Wir suchen Ecken --- Stellen, an denen eine Kante die Richtung ändert. Für jedes Pixel berechnen wir die Matrix der zweiten Momente des Gradienten: \[ \mathbf{M} = \sum_{x,y} w(x,y) \begin{pmatrix} G_x^2 & G_x G_y \\ G_x G_y & G_y^2 \end{pmatrix} \] wobei $w(x,y)$ ein Fenster ist (zum Beispiel Gauß). Wenn beide Eigenwerte von $\mathbf{M}$ groß sind --- das ist eine Ecke. Wenn einer groß und der andere klein ist --- das ist eine Kante. Wenn beide klein sind --- das ist ein flacher Bereich. \subsubsection{SIFT (Scale-Invariant Feature Transform)} SIFT ist einer der beliebtesten Algorithmen (Lowe, 1999). Er findet Keypoints, die stabil sind gegenüber: \begin{itemize} \item Skalierung (das Objekt kann näher oder weiter sein), \item Drehung, \item Änderung der Beleuchtung, \item kleinen Änderungen des Blickwinkels. \end{itemize} Algorithmus: \begin{enumerate} \item Wir bauen eine \textbf{Skalenpyramide}: Wir verwischen das Bild mit einem Gauß-Filter mit unterschiedlichem $\sigma$ und verkleinern es. \item Wir suchen \textbf{Extrema} in der Differenz von Gauß-Funktionen (DoG --- Difference of Gaussians) --- das ist eine Approximation des Laplace-Operators. \item Für jeden Keypoint berechnen wir einen \textbf{Deskriptor}: ein Histogramm der Gradienten in einer $16 \times 16$-Pixel-Nachbarschaft. Das ist ein Vektor aus 128 Zahlen. \item Wir vergleichen die Deskriptoren zwischen Bildern (euklidischer Abstand). \end{enumerate} \subsubsection{SURF, ORB, AKAZE} Es gibt viele Verbesserungen von SIFT: \begin{itemize} \item \textbf{SURF} (Speeded-Up Robust Features) --- schneller, verwendet Integralbilder. \item \textbf{ORB} (Oriented FAST and Rotated BRIEF) --- sehr schnell, patentfrei. \item \textbf{AKAZE} --- verwendet nichtlineare Diffusionsskalen. \end{itemize} \subsection{Finden der Transformation: RANSAC} Wir haben Keypoints in zwei Bildern gefunden und sie einander zugeordnet (feature matching). Aber nicht alle Zuordnungen sind korrekt --- es gibt Ausreißer (outliers). Wie findet man die Transformation (Drehung, Skalierung, Verschiebung), die ein Bild in das andere überführt, wenn Ausreißer vorhanden sind? \subsubsection{RANSAC (Random Sample Consensus)} RANSAC ist ein eleganter Algorithmus zur robusten Parameterschätzung: \begin{enumerate} \item Wir wählen zufällig die minimale Anzahl von Punkten (für eine affine Transformation --- 3 Punktpaare). \item Aus diesen Punkten berechnen wir die Parameter der Transformation. \item Wir zählen, wie viele \textit{aller} Punkte mit dieser Transformation übereinstimmen (inliers) --- Abstand kleiner als eine Schwelle. \item Wir wiederholen $N$ Mal, wählen die Transformation mit der maximalen Anzahl von Inliers. \item Schließlich verfeinern wir die Parameter über alle Inliers (Methode der kleinsten Quadrate). \end{enumerate} \begin{tipbox}[Warum funktioniert RANSAC?] Wenn der Anteil der Ausreißer nicht zu groß ist (sagen wir, $< 50\%$), dann ist die Wahrscheinlichkeit, zufällig nur Inliers zu wählen, von Null verschieden. Wenn wir es oft genug wiederholen, finden wir fast sicher die ``richtige'' Transformation. \end{tipbox} \subsection{Partikelschwarm-Methoden und Optimierung} Manchmal haben wir keine expliziten Keypoints, aber wir wissen, dass ein Bild eine transformierte Version des anderen ist. Wie finden wir die Transformationsparameter? Dies ist ein Optimierungsproblem: Wir maximieren eine \textbf{Ähnlichkeitsmetrik} (zum Beispiel die Mutual Information) über die Transformationsparameter. \subsubsection{Particle Swarm Optimization (PSO)} Die Partikelschwarm-Methode ist eine heuristische Optimierungsmethode, inspiriert vom Verhalten eines Vogelschwarms oder Fischschwarms: \begin{enumerate} \item Wir initialisieren einen ``Schwarm'' aus $N$ Partikeln --- jedes Partikel ist eine Menge von Transformationsparametern (zum Beispiel Drehwinkel, Skalierung, Verschiebungen entlang $x$ und $y$). \item Jedes Partikel hat eine ``Geschwindigkeit'' und eine ``Position''. \item Bei jedem Schritt ``merkt'' sich das Partikel seine beste Position (personal best) und kennt die beste Position des gesamten Schwarms (global best). \item Die Geschwindigkeit des Partikels wird als gewichtete Summe aktualisiert: Trägheit (weiter in dieselbe Richtung bewegen), kognitive Komponente (Streben zum personal best) und soziale Komponente (Streben zum global best). \item Wir wiederholen bis zur Konvergenz. \end{enumerate} PSO ist besonders nützlich, wenn: \begin{itemize} \item Die Ähnlichkeitsfunktion nichtglatt ist oder viele lokale Maxima hat, \item Der Parameterraum groß ist (zum Beispiel eine 3D-Transformation mit 12 Parametern), \item Wir den Gradienten nicht berechnen können (zum Beispiel ist die Mutual Information nicht differenzierbar). \end{itemize} \subsection{Neuronale Netze: von handgefertigten Merkmalen zu gelernten} Alle oben besprochenen Methoden (Sobel, Harris, SIFT) sind \textbf{handgefertigte Merkmale} (hand-crafted features). Wir selbst denken uns aus, was eine ``Kante'', eine ``Ecke'', ein ``interessanter Punkt'' ist und wie man daraus einen Deskriptor baut. Aber was, wenn wir den Computer \textit{selbst} lernen lassen, Merkmale zu finden? \subsubsection{Faltungsneuronale Netze (CNN)} \textbf{Convolutional Neural Networks} sind ein spezieller Typ neuronaler Netze, der für die Arbeit mit Bildern geschaffen wurde. Die Idee: \begin{enumerate} \item \textbf{Faltungsschichten:} Wir wenden eine Menge trainierbarer Filter an (wie Sobel-Masken, aber die Parameter werden aus den Daten gelernt). Jeder Filter extrahiert sein eigenes Merkmal --- von einfachen Kanten in den ersten Schichten bis zu komplexen Texturen und Objektteilen in den tiefen Schichten. \item \textbf{Pooling-Schichten:} Wir verkleinern die Feature-Map (max pooling --- das Maximum in einem $2 \times 2$-Fenster nehmen). Dies verleiht Invarianz gegenüber kleinen Verschiebungen. \item \textbf{Vollständig verbundene Schichten:} Am Ende --- ein gewöhnliches neuronales Netz, das klassifiziert oder regressiert. \end{enumerate} \subsubsection{Hierarchie der Merkmale} Erstaunlicherweise lernen CNNs selbst dieselbe Hierarchie, die wir von Hand konstruiert haben: \begin{itemize} \item \textbf{Erste Schichten:} Einfache Kanten, Gradienten (wie Sobel), \item \textbf{Mittlere Schichten:} Texturen, Ecken, einfache Formen (wie SIFT-Deskriptoren), \item \textbf{Tiefe Schichten:} Objektteile (Augen, Räder, Blätter), \item \textbf{Letzte Schichten:} Ganze Objekte (Gesicht, Auto, Baum). \end{itemize} Dies ist \textbf{Feature Learning} statt handgefertigter Merkmale. \subsection{Vortrainierte Modelle und Transfer Learning} Das Training eines CNN von Grund auf erfordert Millionen annotierter Bilder und Tage der Berechnung auf einer GPU. Aber es gibt einen Trick --- \textbf{Transfer Learning}. Idee: Wir nehmen ein Netz, das auf einer riesigen Datenbank vortrainiert wurde (zum Beispiel ImageNet --- 1.4 Millionen Bilder, 1000 Klassen), und verwenden es als \textbf{Feature-Extraktor}. \subsubsection{Beliebte Architekturen} \begin{itemize} \item \textbf{VGG (2014):} Einfach, tief (16--19 Schichten), gute Merkmale. \item \textbf{ResNet (2015):} Sehr tief (bis zu 152 Schichten) mit ``skip connections'' --- löst das Problem der verschwindenden Gradienten. \item \textbf{EfficientNet (2019):} In allen Parametern optimiert (Genauigkeit, Geschwindigkeit, Größe). \item \textbf{Vision Transformer (ViT, 2020):} Verwendet die Transformer-Architektur (wie in NLP) für Bilder. \end{itemize} \subsubsection{Wie man ein vortrainiertes Modell verwendet} \begin{enumerate} \item \textbf{Feature Extraction:} Wir nehmen ein vortrainiertes Netz, schneiden die letzte Klassifikationsschicht ab und verwenden die vorletzte Schicht als Feature-Extraktor. Wir erhalten einen Vektor von, sagen wir, 2048 Zahlen für jedes Bild. Wir vergleichen die Vektoren (Kosinus-Abstand). \item \textbf{Fine-Tuning:} Wir nehmen ein vortrainiertes Netz und trainieren es auf unseren Daten weiter (zum Beispiel auf mikroskopischen Bildern von Zellen). Die ersten Schichten werden ``eingefroren'' (sie funktionieren bereits gut), die letzten --- trainiert. \item \textbf{Siamese Networks:} Zwei identische Netze mit gemeinsamen Gewichten, trainiert so, dass ähnliche Bilder nahe Feature-Vektoren haben und unterschiedliche --- weit entfernte. \end{enumerate} \subsection{Anwendungen in Chemie und Biologie} \subsubsection{Kryo-Elektronenmikroskopie (cryo-EM)} In der cryo-EM erhalten wir Tausende von 2D-Projektionen von Proteinmolekülen in zufälligen Orientierungen. Die Aufgabe: \begin{enumerate} \item Alle Moleküle auf den Mikrographien finden (particle picking) --- das ist eine Objektdetektionsaufgabe, \item Die Orientierung jeder Projektion bestimmen --- das ist eine Bildvergleichsaufgabe, \item Die 3D-Struktur rekonstruieren --- das ist eine inverse Tomographieaufgabe. \end{enumerate} Moderne Programme (RELION, cryoSPARC) verwenden CNNs für particle picking und Klassifikation. Dies hat eine ``Resolution Revolution'' ermöglicht --- Proteinstrukturen mit atomarer Auflösung zu bestimmen. \subsubsection{Mikroskopie und Zellanalyse} In der biologischen Forschung muss man: \begin{itemize} \item Zellen in mikroskopischen Bildern zählen, \item Zelltypen klassifizieren, \item Zellbewegungen über die Zeit verfolgen, \item Anomalien finden (Krebszellen). \end{itemize} CNN (insbesondere U-Net für die Segmentierung) ist der Standard in diesem Bereich. \subsubsection{Chemometrie und Wirkstoffforschung} In der Wirkstoffforschung muss man molekulare Strukturen vergleichen. Aber Moleküle sind keine Bilder! Man kann sie jedoch als 2D-Bilder darstellen (molekulare Graphen, auf ein Gitter gezeichnet) und CNN zur Aktivitätsvorhersage verwenden. \subsection{Zusammenhang mit früheren Kapiteln} Betrachten wir, wie der Bildvergleich alles verbindet, was wir durchgegangen sind: \begin{itemize} \item \textbf{Lineare Algebra:} Faltung ist die Multiplikation einer Matrix mit einem Vektor (oder Tensor). CNN ist eine Folge linearer Operationen mit Nichtlinearitäten. \item \textbf{SVD und Niedrigrang-Approximation:} CNN können über SVD komprimiert werden (Niedrigrang-Faktorisierung von Faltungskernen). \item \textbf{FFT:} Faltung über FFT ist $\mathcal{O}(N \log N)$ statt $\mathcal{O}(N^2)$. In CNN ist dies entscheidend für die Geschwindigkeit. \item \textbf{Nichtlineare Optimierung:} Das Training eines CNN ist die Minimierung einer Verlustfunktion (cross-entropy, MSE) über Backpropagation + SGD/Adam (Varianten des Gradientenabstiegs). \item \textbf{Compressed Sensing:} In MRT und cryo-EM verwenden wir komprimierte Messungen + CNN zur Rekonstruktion. \end{itemize} \begin{successbox}[Hauptschlussfolgerung] Bildvergleich ist von einfach zu komplex: \begin{enumerate} \item Faltung (für Verschiebungen) --- über FFT, \item Keypoints (SIFT) + RANSAC (für Drehungen und Skalierung), \item Optimierung (PSO) für komplexe Transformationen, \item Neuronale Netze (CNN) --- für semantische Ähnlichkeit. \end{enumerate} Jede Methode ist ein Kompromiss zwischen Universalität, Geschwindigkeit und Genauigkeit. Und alle basieren auf der Mathematik, die wir durchgegangen sind: lineare Algebra, Optimierung, Fourier-Analyse, Wahrscheinlichkeitstheorie. \end{successbox} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Maschinelles Lernen: vom Perzeptron zu AlphaFold} \subsection{Was ist maschinelles Lernen?} Maschinelles Lernen (Machine Learning, ML) ist ein Zweig der künstlichen Intelligenz, in dem der Computer \textit{selbst lernt}, Aufgaben auf der Grundlage von Daten zu lösen, statt nach fest vorgeschriebenen Regeln. Statt ein Programm ``wenn Temperatur > 100, dann siedet Wasser'' zu schreiben, geben wir dem Computer Tausende von Beispielen ``Temperatur --- Zustand des Wassers'' und bitten ihn, das Muster zu finden. \subsubsection{Drei Arten des Lernens} \begin{itemize} \item \textbf{Überwachtes Lernen (supervised):} Es gibt annotierte Daten (Eingabe $\to$ Ausgabe). Die Aufgabe ist, die Ausgabe für neue Eingaben vorherzusagen. Beispiele: Bildklassifikation, Vorhersage von Moleküleigenschaften. \item \textbf{Unüberwachtes Lernen (unsupervised):} Daten ohne Annotation. Die Aufgabe ist, Struktur, Cluster, verborgene Muster zu finden. Beispiele: Clustering von Spektren, PCA. \item \textbf{Verstärkungslernen (reinforcement):} Ein Agent interagiert mit der Umgebung und erhält Belohnungen/Strafen. Die Aufgabe ist, eine Strategie zu lernen, die die Belohnung maximiert. Beispiele: Schachspielen, Robotersteuerung. \end{itemize} \subsection{Vom Perzeptron zu neuronalen Netzen} \subsubsection{Das Perzeptron (1958, Rosenblatt)} Das einfachste ``neuronale Netz'' ist das Perzeptron. Es berechnet: \[ y = \sigma\left(\sum_{i=1}^n w_i x_i + b\right) \] wobei $x_i$ die Eingaben, $w_i$ die Gewichte, $b$ der Bias und $\sigma$ die Aktivierungsfunktion ist (zum Beispiel eine Stufe oder eine Sigmoid-Funktion). Das Perzeptron kann nur \textbf{linear trennbare} Probleme lösen (zum Beispiel logisches ODER). Für XOR (exklusives ODER) genügt ein Perzeptron nicht --- man braucht ein mehrschichtiges Netz. \subsubsection{Mehrschichtiges Perzeptron (MLP)} Ein Netz aus mehreren Schichten: Eingabeschicht, eine oder mehrere verborgene Schichten, Ausgabeschicht. Jedes Neuron ist mit allen Neuronen der nächsten Schicht verbunden. Formel für die $l$-te Schicht: \[ \vec{h}^{(l)} = \sigma\left(\mathbf{W}^{(l)} \vec{h}^{(l-1)} + \vec{b}^{(l)}\right) \] wobei $\mathbf{W}^{(l)}$ die Gewichtsmatrix, $\vec{b}^{(l)}$ der Bias-Vektor und $\sigma$ eine nichtlineare Aktivierungsfunktion ist (ReLU, tanh, sigmoid). \subsubsection{Training: Backpropagation} Wie trainiert man ein neuronales Netz? Man minimiert die Verlustfunktion $L$ (zum Beispiel MSE für Regression oder cross-entropy für Klassifikation) über Gradientenabstieg. \textbf{Backpropagation} ist ein effizienter Algorithmus zur Berechnung von Gradienten über die Kettenregel. Wir gehen von der Ausgabe zur Eingabe und berechnen $\partial L / \partial w_{ij}$ für jedes Gewicht und aktualisieren die Gewichte: \[ w_{ij} \leftarrow w_{ij} - \eta \frac{\partial L}{\partial w_{ij}} \] wobei $\eta$ die Lernrate ist. \subsection{Faltungsneuronale Netze (CNN)} Für Bilder ist MLP ineffizient: zu viele Parameter (jedes Pixel ist mit allen Neuronen verbunden). CNN verwenden \textbf{Faltungsschichten}: Ein kleiner Filter (zum Beispiel $3 \times 3$) gleitet über das Bild und berechnet eine Faltung. Dies ergibt: \begin{itemize} \item \textbf{Lokalität der Verbindungen:} Jedes Neuron schaut nur auf einen kleinen Bereich, \item \textbf{Gewichtsteilung:} Derselbe Filter wird auf das gesamte Bild angewendet, \item \textbf{Verschiebungsinvarianz:} Das Objekt wird überall gefunden. \end{itemize} Beliebte Architekturen: LeNet (1998), AlexNet (2012), VGG (2014), ResNet (2015), EfficientNet (2019). \subsection{Rekurrente neuronale Netze (RNN)} Für Sequenzen (Text, Zeitreihen) braucht man Netze mit \textbf{Gedächtnis}. RNN übertragen den verborgenen Zustand von einem Schritt zum nächsten: \[ \vec{h}_t = \sigma\left(\mathbf{W}_h \vec{h}_{t-1} + \mathbf{W}_x \vec{x}_t + \vec{b}\right) \] Das Problem gewöhnlicher RNN sind \textbf{verschwindende Gradienten} (vanishing gradients): Das Netz lernt bei langen Sequenzen schlecht. Lösungen: \begin{itemize} \item \textbf{LSTM} (Long Short-Term Memory, 1997): Fügt ``Tore'' (gates) hinzu, die den Informationsfluss steuern. \item \textbf{GRU} (Gated Recurrent Unit, 2014): Eine vereinfachte Version von LSTM. \end{itemize} \subsection{Die Revolution: Transformer und Attention} Im Jahr 2017 erschien die Arbeit ``Attention Is All You Need'' (Vaswani et al.), die NLP und vieles mehr auf den Kopf stellte. \subsubsection{Der Attention-Mechanismus} Idee: Statt die ganze Sequenz in einen Vektor zu komprimieren (wie bei RNN), erlauben wir jedem Element, auf alle anderen Elemente zu ``schauen'' und sie nach Wichtigkeit zu gewichten. Für eine Eingabesequenz $\vec{x}_1, \dots, \vec{x}_n$ berechnen wir: \[ \text{Attention}(\mathbf{Q}, \mathbf{K}, \mathbf{V}) = \text{softmax}\left(\frac{\mathbf{Q}\mathbf{K}^T}{\sqrt{d_k}}\right)\mathbf{V} \] wobei $\mathbf{Q}$ (query), $\mathbf{K}$ (key), $\mathbf{V}$ (value) lineare Projektionen der Eingaben sind. \subsubsection{Transformer} Die Transformer-Architektur verzichtet vollständig auf Rekurrenz und Faltungen. Stattdessen wird nur der Attention-Mechanismus (self-attention) und vollständig verbundene Schichten verwendet. Schlüsselkomponenten: \begin{itemize} \item \textbf{Multi-Head Attention:} Mehrere parallele Attention-Mechanismen, von denen jeder lernt, sich auf verschiedene Aspekte der Daten zu konzentrieren. \item \textbf{Positional Encoding:} Da der Transformer keine Rekurrenz hat, kennt er die Reihenfolge der Elemente nicht. Wir fügen Positionskodierungen hinzu (Sinus/Kosinus oder trainierbare Vektoren). \item \textbf{Layer Normalization:} Normalisierung der Aktivierungen für Trainingsstabilität. \item \textbf{Residual Connections:} ``Skip connections'' zur Bekämpfung verschwindender Gradienten (wie in ResNet). \end{itemize} Transformer wurde zur Grundlage für: \begin{itemize} \item \textbf{BERT} (2018): Bidirectional Encoder Representations from Transformers --- Vortraining auf maskierten Wörtern. \item \textbf{GPT} (2018--2023): Generative Pre-trained Transformer --- autoregressive Textgenerierung. GPT-3 (175 Milliarden Parameter), GPT-4 (multimodal). \item \textbf{T5, BART}: Encoder-Decoder-Architekturen für Übersetzung, Zusammenfassung. \end{itemize} \subsection{Maschinelles Lernen in der Chemie: Evolution} \subsubsection{Die QSAR-Ära (1990er -- 2010er)} QSAR (Quantitative Structure-Activity Relationship) ist ein klassischer Ansatz: Wir berechnen \textbf{molekulare Deskriptoren} (numerische Eigenschaften eines Moleküls) und bauen eine Regression/Klassifikation. Deskriptoren: \begin{itemize} \item Physikochemische: Molmasse, logP (Lipophilie), polare Oberfläche, \item Topologische: Konnektivitätsindizes, Formen des Molekülgraphen, \item Elektronische: Atomladungen, Orbitalenergien, \item 3D-Deskriptoren: Trägheitsmomente, Gyrationsradien. \end{itemize} Methoden: PLS (Partial Least Squares), Random Forest, SVM (Support Vector Machines). \begin{tipbox}[Beispiel] Vorhersage der Toxizität eines Moleküls: Wir berechnen 200 Deskriptoren, sammeln eine Datenbank von 10\,000 Molekülen mit bekannter Toxizität, trainieren einen Random Forest. Genauigkeit --- 80--85\%. \end{tipbox} \subsubsection{Graphneuronale Netze (2015 -- heute)} Ein Molekül ist ein \textbf{Graph}: Atome sind Knoten, Bindungen sind Kanten. Graphneuronale Netze (Graph Neural Networks, GNN) arbeiten direkt mit Graphen. \textbf{Message Passing Neural Networks (MPNN):} \begin{enumerate} \item Jedes Atom hat eine Anfangsdarstellung (Merkmalsvektor: Atomtyp, Ladung, Hybridisierung). \item Bei jedem Schritt ``tauschen Atome Nachrichten'' mit Nachbarn über Bindungen aus. \item Nach $K$ Schritten ``kennt'' jedes Atom seine Nachbarschaft vom Radius $K$ Bindungen. \item Die finale Darstellung des Moleküls ist eine Aggregation aller atomaren Darstellungen. \end{enumerate} Beliebte Architekturen: \begin{itemize} \item \textbf{GCN} (Graph Convolutional Network): Ein Analogon der Faltung für Graphen. \item \textbf{GAT} (Graph Attention Network): Attention zwischen Atomen. \item \textbf{SchNet, DimeNet, SphereNet}: Berücksichtigen 3D-Koordinaten der Atome und Winkel. \end{itemize} \begin{successbox}[Warum sind GNN so gut für die Chemie?] \begin{itemize} \item \textbf{Permutationsinvarianz:} Die Reihenfolge der Atome ist egal, \item \textbf{Strukturbewusstsein:} GNN sehen Bindungen und Geometrie, \item \textbf{Interpretierbarkeit:} Man kann schauen, welche Atome/Bindungen für die Vorhersage wichtig sind. \end{itemize} \end{successbox} \subsubsection{Vorhersage von Moleküleigenschaften} Moderne Muss-Modelle für Chemiker: \begin{tabularx}{\textwidth}{l X l} \toprule \textbf{Aufgabe} & \textbf{Modell} & \textbf{Genauigkeit} \\ \midrule Molekülenergie & SchNet, DimeNet++ & MAE $\sim$1 kcal/mol \\ Solvatation & SolTranX & MAE $\sim$0.5 kcal/mol \\ pKa & ChemProp, GNN & MAE $\sim$0.3 \\ LogP & MolCLR, GNN & MAE $\sim$0.2 \\ Toxizität & GraphCL, GNN & AUC $\sim$0.9 \\ Wirkstoffähnlichkeit & MolBERT & AUC $\sim$0.85 \\ \bottomrule \end{tabularx} \subsection{Die AlphaFold-Revolution: Vorhersage der Proteinstruktur} \subsubsection{Das Problem} Die Vorhersage der 3D-Struktur eines Proteins aus seiner Aminosäuresequenz ist eine der großen Herausforderungen der Biologie. Experimentelle Methoden (X-ray, cryo-EM, NMR) sind teuer und langsam. \subsubsection{AlphaFold (2020, DeepMind)} AlphaFold 2 ist ein Durchbruch, der das Problem der Proteinstrukturvorhersage mit atomarer Genauigkeit gelöst hat. Architektur: \begin{enumerate} \item \textbf{Evoformer:} Ein Transformer, der verarbeitet: \begin{itemize} \item Multiple Sequenzalignment (MSA) --- evolutionäre Information, \item Pair representation --- paarweise Abstände zwischen Resten. \end{itemize} \item \textbf{Structure Module:} Baut iterativ 3D-Koordinaten der Atome auf und minimiert die Verlustfunktion. \item \textbf{Confidence Estimation:} Sagt pLDDT (per-residue confidence) voraus --- wie sicher das Modell bei jedem Rest ist. \end{enumerate} Ergebnisse auf CASP14 (Critical Assessment of Structure Prediction): \begin{itemize} \item Mittleres GDT\_TS (Global Distance Test) --- 92.4 (experimentelle Genauigkeit --- $\sim$90), \item Für 2/3 der Proteine --- Genauigkeit $< 1$ Å (atomare Genauigkeit). \end{itemize} \subsubsection{AlphaFold 3 (2024)} Erweiterung auf: \begin{itemize} \item Protein--Ligand-Komplexe, \item DNA/RNA-Strukturen, \item posttranslationale Modifikationen, \item ionische Wechselwirkungen. \end{itemize} \subsection{Ein Fundamentmodell für Konformere} \subsubsection{Das Problem des Konformationsraums} Ein Molekül ist keine statische Struktur. Es ``atmet'' ständig, dreht sich um Bindungen, ändert die Konformation. Für das Wirkstoffdesign ist es entscheidend, nicht eine Struktur zu kennen, sondern ein \textbf{Ensemble} niedrigenergetischer Konformere. Klassische Methoden: \begin{itemize} \item \textbf{Molecular Dynamics (MD):} Wir modellieren die Bewegung der Atome in der Zeit, aber das ist langsam (Nanosekunden --- Mikrosekunden). \item \textbf{Monte Carlo:} Zufällige Wanderung durch den Konformationsraum, aber ineffizient für große Moleküle. \item \textbf{Systematic/Rotamer search:} Durchmusterung aller Kombinationen von Torsionswinkeln, aber exponentielles Wachstum. \end{itemize} \subsubsection{ML-Ansatz: Konformer-Vorhersage} Moderne Modelle lernen, die Verteilung der Konformere direkt aus Daten vorherzusagen. \textbf{GeoLDM (Geometric Latent Diffusion Models, 2023):} \begin{enumerate} \item \textbf{Encoder:} Kodiert eine 3D-Konformation in einen latenten Raum. \item \textbf{Diffusion Process:} Fügt Rauschen zu den latenten Vektoren hinzu (wie in DALL-E, Stable Diffusion). \item \textbf{Decoder:} Erzeugt neue Konformationen aus Rauschen über den umgekehrten Diffusionsprozess. \item \textbf{Energy Model:} Filtert die erzeugten Konformationen nach Energie. \end{enumerate} \textbf{ConfGF (Conformation Graph Factor, 2022):} \begin{itemize} \item Verwendet GNN zur Vorhersage von 3D-Koordinaten, \item Erzeugt ein Ensemble von Konformeren durch Sampling, \item Berücksichtigt die Boltzmann-Verteilung (niedrigenergetische Konformere sind wahrscheinlicher). \end{itemize} \subsubsection{Fundamentmodell: Uni-Mol (2023)} Uni-Mol ist ein ``Foundation Model'' für Moleküle, vortrainiert auf Millionen von Konformationen. Architektur: \begin{enumerate} \item \textbf{3D Transformer:} Arbeitet mit Atomkoordinaten und Merkmalen. \item \textbf{Pre-training tasks:} \begin{itemize} \item Vorhersage von Abständen zwischen Atomen, \item Vorhersage von Winkeln und Diederwinkeln, \item Maskierte Atomvorhersage (wie BERT), \item Kontrastives Lernen (ähnliche Konformere --- nah). \end{itemize} \item \textbf{Fine-tuning:} Weiteres Training auf konkreten Aufgaben (Energie, Eigenschaften, Konformere). \end{enumerate} Ergebnisse: \begin{itemize} \item Energievorhersage: MAE $\sim$0.5 kcal/mol (besser als DFT für einige Klassen), \item Konformer-Generierung: RMSD $< 0.5$ Å für 90\% der Moleküle, \item Eigenschaftsvorhersage: State-of-the-Art auf den meisten Benchmarks. \end{itemize} \subsection{Muss-Werkzeuge für den modernen Chemiker} \subsubsection{Software} \begin{tabularx}{\textwidth}{l X} \toprule \textbf{Werkzeug} & \textbf{Zweck} \\ \midrule RDKit & Chemoinformatik: Deskriptoren, Fingerprints, Ähnlichkeit \\ DeepChem & ML für Wirkstoffforschung, GNN \\ PyTorch Geometric & Graphneuronale Netze \\ OpenMM & Molekulardynamik auf GPU \\ ASE (Atomic Simulation Environment) & Quantenchemie + ML \\ SchNetPack & SchNet und andere GNN für Chemie \\ AlphaFold (ColabFold) & Proteinstrukturvorhersage \\ Uni-Mol & Foundation Model für Moleküle \\ \bottomrule \end{tabularx} \subsubsection{Typischer Workflow} \begin{enumerate} \item \textbf{Datensammlung:} PubChem, ChEMBL, PDB (Protein Data Bank). \item \textbf{Vorverarbeitung:} RDKit für Deskriptoren, Generierung von Konformeren. \item \textbf{Modellierung:} GNN zur Eigenschaftsvorhersage, AlphaFold für Struktur. \item \textbf{Validierung:} Cross-Validation, externer Testsatz, experimentelle Verifikation. \item \textbf{Interpretation:} Attention-Gewichte, SHAP-Werte zum Verständnis der Vorhersagen. \end{enumerate} \subsection{Zusammenhang mit früheren Kapiteln} Betrachten wir, wie ML alles verbindet, was wir durchgegangen sind: \begin{itemize} \item \textbf{Lineare Algebra:} Neuronale Netze sind eine Folge von Matrixmultiplikationen $\mathbf{W}\vec{x} + \vec{b}$. Backpropagation ist die Kettenregel für Matrizen. \item \textbf{SVD und Tensoren:} Modellkompression durch Niedrigrang-Faktorisierung. Tensorzerlegungen für effiziente Berechnungen. \item \textbf{FFT:} Faltungsschichten verwenden FFT zur Beschleunigung. Attention über FFT (Linear Attention). \item \textbf{Nichtlineare Optimierung:} Das Training neuronaler Netze ist die Minimierung eines Verlusts über SGD, Adam (adaptive Methoden), L-BFGS (für Fine-Tuning). \item \textbf{Compressed Sensing:} Sparse Training, Pruning neuronaler Netze. \item \textbf{Graphen und Tensoren:} GNN arbeiten mit Molekülgraphen. 3D-Modelle verwenden Tensoren von Koordinaten. \item \textbf{Monte Carlo:} Diffusionsmodelle sind stochastische Prozesse. Variational Inference --- ein Bayesscher Ansatz. \end{itemize} \begin{successbox}[Hauptschlussfolgerung] Maschinelles Lernen ist kein Ersatz für klassische Chemie, sondern ein \textbf{mächtiges Werkzeug}, das: \begin{itemize} \item Berechnungen um das Tausendfache beschleunigt (Eigenschaftsvorhersage vs. DFT), \item neue Möglichkeiten eröffnet (Proteinstrukturvorhersage, Molekülgenerierung), \item Verständnis der Mathematik erfordert (lineare Algebra, Optimierung, Wahrscheinlichkeitstheorie). \end{itemize} Der moderne Chemiker muss wissen: \begin{itemize} \item Grundlagen des ML (neuronale Netze, GNN, Transformer), \item Werkzeuge (RDKit, PyTorch, DeepChem), \item Wann ML funktioniert und wann nicht (physikalische Grenzen, Interpretierbarkeit). \end{itemize} \end{successbox} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Moderne Rechentechnik: vom Transistor zum Supercomputer} \subsection{Zahlen, die man kennen sollte} Ein moderner Computer --- eine Workstation --- ist: \begin{itemize} \item 128--1024 GB Arbeitsspeicher, \item $\sim$500 Gflop/s Spitzenleistung (Milliarden Gleitkommaoperationen pro Sekunde), \item mehrere Prozessoren (von 4 bis 128 Kernen). \end{itemize} Aber hier ist ein Paradoxon: Die Geschwindigkeit des Speicherzugriffs ist \textbf{Hunderte Male} langsamer als die Geschwindigkeit arithmetischer Operationen. Und wenn der Zugriff zufällig ist (nicht sequentielles Lesen-Schreiben), dann noch einmal \textbf{hundert Mal} langsamer. Warum ist das so? Klären wir das. \subsection{Wie ein Prozessor aufgebaut ist} Alle sagen, dass ein Prozessor ``viele Transistoren'' hat. Aber was tun sie? Tatsächlich ist ein Prozessor eine Entität, die auf \textbf{Instruktionen} zugreift und auf \textbf{Daten} zugreift. Beides belegt Platz im Speicher. \subsubsection{Prozessorbefehle} Der Prozessor liest Befehle, die üblicherweise erlauben: \begin{itemize} \item Auf Daten im Speicher zuzugreifen und sie in \textbf{Register} zu legen --- einen temporären, sehr, sehr schnellen Speicher, \item Eine Operation auszuführen: Addition, Subtraktion, Multiplikation, Vergleich, \item Einen \textbf{Sprung} auszuführen: Standardmäßig führt der Prozessor Befehl für Befehl aus, aber wenn er auf eine Sprunginstruktion trifft, kann er vorwärts oder rückwärts springen. \end{itemize} Sprünge sind: \begin{itemize} \item \textbf{Unbedingt:} Der Sprung erfolgt immer, \item \textbf{Bedingt:} Wenn wir zuvor Zahlen verglichen haben und die Antwort ``ja'' oder ``nein'' haben, dann erfolgt der Sprung nur bei Erfüllung der Bedingung. \end{itemize} \subsubsection{SIMD: ein Befehl --- viele Daten} Ein moderner Prozessorbefehl kann mit mehreren Zahlen gleichzeitig arbeiten. Zum Beispiel kann eine AVX-512-Instruktion 16 Paare von Additionen gleichzeitig ausführen! Dies nennt man \textbf{SIMD} (Single Instruction, Multiple Data) --- eine Instruktion, viele Daten. Wenn die Daten bereits auf dem Prozessor selbst sind, ist das ziemlich einfach. \begin{tipbox}[Woher kommt SIMD?] Die Leute bemerkten, dass man oft identische Operationen ausführen muss: zum Beispiel Pixel identisch transformieren, um ein Bild zu zeichnen, Array-Elemente identisch verarbeiten. Warum dieselbe Instruktion 16 Mal ausführen, wenn man es einmal tun kann --- aber für 16 Zahlen gleichzeitig? \end{tipbox} Aber um mit jedem Befehl solche Operationen mit \textit{verschiedenen} Daten auszuführen, bräuchten wir ein sehr leistungsfähiges System zum Pumpen von Daten aus dem Speicher in den Prozessor und zurück. \subsubsection{Das Memory-Wall-Problem} Leider würde ein solches System Hunderte und Tausende Male mehr Transistoren belegen, als für einen arithmetischen Block nötig sind. Und die Leute bemerkten, dass wir in Programmen oft mit einer kleinen Menge von Zahlen arbeiten und dann zu einer anderen Menge übergehen, und so weiter. Einen kleinen Speicherblock mit schnellem Zugriff zu bauen ist nicht so kostspielig an Transistoren. Nur den Programmierer zu zwingen, Daten in diesen Block aus dem Arbeitsspeicher und zurück zu schleppen, wäre schwierig. Deshalb hat man das automatisiert --- und der \textbf{Prozessor-Cache} entstand. Es ist schneller Speicher, es gibt nicht viel davon, aber der Prozessor ``merkt'' sich selbst, was du gerade verwendest, und hält diese Daten einige Zeit in seinem Cache, in der Erwartung, dass sie später vielleicht wieder gebraucht werden. \subsection{Speicherhierarchie} Ein moderner Computer ist eine mehrschichtige Speicherstruktur: \begin{tabularx}{\textwidth}{l l l X} \toprule \textbf{Ebene} & \textbf{Größe} & \textbf{Geschwindigkeit} & \textbf{Beschreibung} \\ \midrule Register & 32--128 $\times$ 64 Bit & $\sim$0.3 ns & Zellen im Prozessor, mit denen er direkt arbeitet \\ L1-Cache & $\sim$32--64 KB & $\sim$1 ns & Getrennt für Instruktionen und Daten, eigener pro Kern \\ L2-Cache & $\sim$256 KB -- 1 MB & $\sim$3--5 ns & Eigener pro Kern oder gemeinsam für ein Kernpaar \\ L3-Cache & $\sim$8--64 MB & $\sim$10--20 ns & Gemeinsam für alle Kerne des Prozessors \\ RAM (DRAM) & 128--1024 GB & $\sim$50--100 ns & Hauptspeicher \\ Festplatte (SSD) & 1--10 TB & $\sim$10--100 $\mu$s & Langzeitspeicherung \\ \bottomrule \end{tabularx} Vergleichen wir, wie viel eine moderne Workstation in einer Nanosekunde durchschnittlich leisten kann: \begin{itemize} \item \textbf{Arithmetik:} Alle ihre Prozessoren können unter Berücksichtigung langer SIMD-Instruktionen etwa \textbf{1000 Gleitkommaoperationen} ausführen. \item \textbf{L1-Cache:} Man kann 2--3 Zahlen pro Prozessor lesen und schreiben, insgesamt $\sim$200 Zahlen. \item \textbf{L2/L3-Cache:} Diese Größe fällt um Zehner. \item \textbf{RAM:} Nur wenige Lese- und Schreibvorgänge pro Nanosekunde, und nur, wenn dies in großen Blöcken von etwa 512 Bytes etwa sequenziell geschieht. \item \textbf{Zufälliger Zugriff:} Wenn wir jedes Mal auf verschiedene Blöcke zugreifen, fällt die Geschwindigkeit um weitere Zehner --- mehrere Nanosekunden können pro Zahl nötig sein. \end{itemize} \begin{warningbox}[Memory Wall] Der zufällige Speicherzugriff ist \textbf{Zehntausende Male} langsamer als die reale Rechenleistung eines modernen Prozessors! Dies ist der Haupt-Bottleneck (``Flaschenhals'') moderner Berechnungen. Genau deshalb: \begin{itemize} \item Matrixalgorithmen schneller arbeiten, wenn die Daten ``dicht'' im Speicher liegen, \item dünnbesetzte Matrizen langsam sind (viele zufällige Zugriffe), \item cache-orientierte Algorithmen eine eigene Wissenschaft sind. \end{itemize} \end{warningbox} \subsection{Grafikprozessoren (GPU)} Aber das war den Leuten nicht genug. Als sie 3D-Objekte der Computergrafik auf dem Bildschirm zeichneten, bemerkten sie, dass all dies fast vollständig \textbf{parallel} erledigt werden kann. Grob: Wenn wir 1024 Prozessoren hätten, würden wir den ganzen Bildschirm in ein Gitter von $32 \times 32$ Quadraten aufteilen, und jeder Prozessor würde in seinem Quadrat zeichnen. So entstanden Grafikkarten --- sie haben eine Menge kleiner spezialisierter Prozessoren, die vollständig parallel arbeiten. \subsubsection{Von der Grafik zu Berechnungen} Ende der 1990er Jahre wurden solche Grafikkarten zur Beschleunigung der Grafikdarstellung verwendet. Um 2005 begann man, solche Karten so herzustellen, dass kleine Code-Stücke auf Tausenden solcher kleiner Prozessoren ausgeführt werden konnten --- so entstanden \textbf{CUDA} (NVIDIA) und \textbf{OpenCL} (ein offener Standard). Bis 2010 begannen alle numerischen Methoden auf Grafikkarten überzugehen. Allmählich fügte man zur Single- und Double-Gleitkomma-Arithmetik die sogenannten \textbf{Halb-} (FP16) und \textbf{Viertel-} (FP8, INT8) Gleitkommazahlen hinzu --- weil mit ihnen Berechnungen schneller waren und sie weniger Speicher belegten. In modernen NVIDIA-Grafikkarten (H100, B200) ist es möglich, etwa \textbf{eine Million Operationen} mit solchen kleinen Objekten in derselben einen Nanosekunde auszuführen --- und dies hat es ermöglicht, solche Grafikkarten für moderne Systeme der künstlichen Intelligenz zu verwenden. \subsubsection{GPU-Architektur} Eine moderne GPU (zum Beispiel NVIDIA H100) ist: \begin{itemize} \item \textbf{Streaming Multiprocessors (SM):} $\sim$132 Multiprocessoren, \item \textbf{CUDA-Kerne:} $\sim$16\,896 Kerne für FP32-Arithmetik, \item \textbf{Tensor-Kerne:} $\sim$528 spezialisierte Kerne für Matrixoperationen (FP16, FP8, INT8), \item \textbf{Speicher:} 80 GB HBM3 (High Bandwidth Memory) mit einer Bandbreite von $\sim$3 TB/s. \end{itemize} \subsubsection{Threads und Warps} Eine GPU arbeitet mit \textbf{Threads}. Tausende von Threads laufen parallel, aber sie sind in Gruppen organisiert: \begin{itemize} \item \textbf{Warp:} Eine Gruppe von 32 Threads, die \textbf{synchron} ausgeführt werden --- alle 32 Threads führen dieselbe Instruktion zum selben Zeitpunkt aus. \item \textbf{Block:} Eine Gruppe von Warps (bis zu 1024 Threads), die Daten über einen gemeinsamen \textbf{shared memory} austauschen können. \item \textbf{Grid:} Alle auf der GPU gestarteten Blöcke. \end{itemize} \begin{orangebox}[Regel der GPU-Effizienz] Damit eine GPU effizient arbeitet, dürfen die Programme auf jedem Multiprocessor für jeden Thread \textbf{nicht auseinanderlaufen}. Das bedeutet: \begin{itemize} \item Alle 32 Threads in einem Warp müssen \textit{dieselbe} Instruktion ausführen (ohne bedingte Sprünge, die die Threads trennen --- ``warp divergence''), \item Alle Threads müssen auf \textit{aufeinanderfolgende} Speicheradressen zugreifen (coalesced memory access), \item Es müssen genügend Threads vorhanden sein, um die Speicherlatenzen zu ``verstecken''. \end{itemize} Wenn diese Bedingungen nicht erfüllt sind --- arbeitet die GPU um ein Vielfaches langsamer. \end{orangebox} \subsection{Paralleles Rechnen: Flynns Taxonomie} Im Jahr 1966 schlug Michael Flynn eine Klassifikation paralleler Architekturen nach Instruktions- und Datenströmen vor: \begin{tabularx}{\textwidth}{l X l} \toprule \textbf{Typ} & \textbf{Beschreibung} & \textbf{Beispiel} \\ \midrule \textbf{SISD} & Single Instruction, Single Data. Ein Instruktionsstrom, ein Datenstrom. Klassischer sequentieller Prozessor. & Alte CPUs \\ \textbf{SIMD} & Single Instruction, Multiple Data. Eine Instruktion auf viele Daten angewendet. & Vektorinstruktionen (AVX), GPU \\ \textbf{MISD} & Multiple Instruction, Single Data. Verschiedene Instruktionen auf dieselben Daten. Selten. & Einige spezialisierte Architekturen \\ \textbf{MIMD} & Multiple Instruction, Multiple Data. Viele Prozessoren, jeder führt seine eigenen Instruktionen auf seinen eigenen Daten aus. & Moderne Multiprozessoren, Cluster \\ \bottomrule \end{tabularx} \subsubsection{Gemeinsamer und verteilter Speicher} In MIMD-Systemen gibt es zwei Ansätze: \begin{itemize} \item \textbf{Gemeinsamer Speicher (shared memory):} Alle Prozessoren haben Zugriff auf einen einzigen Speicher. Beispiel: Multi-Core-CPU. Leicht zu programmieren (alle sehen dieselben Daten), aber schwer zu skalieren (Speicherzugriffskonflikte). \item \textbf{Verteilter Speicher (distributed memory):} Jeder Prozessor hat seinen eigenen lokalen Speicher, Datenaustausch --- über ein Netzwerk. Beispiel: Cluster, Supercomputer. Schwieriger zu programmieren (man muss Daten explizit übertragen --- MPI), aber skaliert bis zu Millionen von Kernen. \end{itemize} \subsection{TOP500: die leistungsstärksten Supercomputer der Welt} Die TOP500-Liste ist eine Rangliste der 500 leistungsstärksten Supercomputer der Welt, die zweimal jährlich aktualisiert wird (Juni und November). Die Leistung wird im LINPACK-Benchmark gemessen --- der Lösung eines großen linearen Gleichungssystems. \subsubsection{Moderne Spitzenreiter (2024--2025)} \begin{tabularx}{\textwidth}{l l l X} \toprule \textbf{Name} & \textbf{Rang} & \textbf{Leistung} & \textbf{Architektur} \\ \midrule Frontier & \#1 & $\sim$1.2 Eflop/s & AMD EPYC + AMD Instinct MI250X \\ Aurora & \#2 & $\sim$1 Eflop/s & Intel Xeon + Intel Max PVC \\ Eagle & \#3 & $\sim$561 Pflop/s & AMD EPYC + NVIDIA H100 \\ Fugaku & \#4 & $\sim$442 Pflop/s & Fujitsu A64FX (ARM) \\ LUMI & \#5 & $\sim$231 Pflop/s & AMD EPYC + AMD MI250X \\ \bottomrule \end{tabularx} \begin{tipbox}[Maßstab] 1 Eflop/s = $10^{18}$ Operationen pro Sekunde. Frontier ist $\sim$1.2 Millionen Milliarden Operationen pro Sekunde. Zum Vergleich: Dein Laptop ist $\sim$0.1 Tflop/s = $10^{11}$ Operationen pro Sekunde. Der Unterschied ist 10 Millionen Mal! \end{tipbox} \subsubsection{Wie ein moderner Supercomputer aufgebaut ist} Typische Architektur: \begin{enumerate} \item \textbf{Knoten (nodes):} $\sim$10\,000--100\,000 Server, jeder mit 2--4 CPUs + 4--8 GPUs, \item \textbf{Netzwerk:} Hochgeschwindigkeits-Interconnect (InfiniBand, Slingshot) mit einer Bandbreite von 200--400 Gbit/s pro Knoten, \item \textbf{Speicher:} Paralleles Dateisystem (Lustre, GPFS) mit einer Bandbreite von $\sim$1--10 TB/s, \item \textbf{Kühlung:} Flüssigkeitskühlung (Wasser oder sogar zweiphasig), da die Leistung 20--50 MW beträgt. \end{enumerate} \subsection{Digitalisierer: eine Brücke zwischen der analogen und der digitalen Welt} Aber der Prozessor --- er muss die Daten für die Berechnung doch irgendwoher nehmen. Und hier betreten \textbf{Analog-Digital-Umsetzer} (ADC --- Analog-to-Digital Converters) die Bühne. Sie werden durch zwei Parameter bestimmt: \begin{itemize} \item \textbf{Genauigkeit (Bit-Tiefe):} Wie viele Bits zur Kodierung eines Abtastwerts verwendet werden, \item \textbf{Geschwindigkeit (Abtastrate):} Wie viele Abtastwerte pro Sekunde. \end{itemize} \subsubsection{Der Kompromiss: Genauigkeit vs. Geschwindigkeit} Es gibt einen fundamentalen Kompromiss: Je höher die Genauigkeit, desto niedriger die Geschwindigkeit. \begin{tabularx}{\textwidth}{l l l X} \toprule \textbf{Bit-Tiefe} & \textbf{Geschwindigkeit} & \textbf{Minimales Bit} & \textbf{Anwendung} \\ \midrule 12 Bit & $\sim$10 GS/s & $\sim$500 $\mu$V & Oszilloskope, schnelle Bildgebung \\ 16 Bit & $\sim$1 GS/s & $\sim$10--50 $\mu$V & Audio, NMR, präzise Messungen \\ 24 Bit & $\sim$4 MS/s & $\sim$100--500 nV & Hochauflösendes Audio, Seismographie \\ \bottomrule \end{tabularx} \begin{warningbox}[Warum werden Nanovolt nicht direkt gemessen?] 24 Bit ergeben einen Bereich von Zehnern und Hunderten von Nanovolt. Aber Nanovolt werden von niemandem direkt gemessen --- Radarstörungen von allem um uns herum betragen etwa Zehner von Mikrovolt! Das heißt, man misst dort üblicherweise \textit{Ströme}, nicht Spannungen --- oder verwendet spezielle abgeschirmte Kammern. \end{warningbox} \subsection{Physik des Schaltens: vom Transistor zum Megavolt} In einem modernen Prozessor funktioniert alles, indem Transistoren bei etwa 0.7 V ein- und ausgeschaltet werden. Die Geschwindigkeit eines solchen Schaltens beträgt Zehner von Gigahertz. Aber die Ströme dort sind sehr klein (Mikroampere pro Transistor). \subsubsection{Der Kompromiss: Spannung vs. Geschwindigkeit} Wenn man die Spannung erhöht, beginnt die Schaltgeschwindigkeit zu fallen: \begin{tabularx}{\textwidth}{l l X} \toprule \textbf{Spannung} & \textbf{Schaltzeit} & \textbf{Technologie} \\ \midrule 0.7 V & $\sim$0.05 ns (20 GHz) & Moderne CMOS-Transistoren \\ 40 V & $\sim$3--5 ns & Leistungs-MOSFETs (gewöhnliche Nennwerte) \\ 1000 V & $\sim$3--5 ns & Leistungs-MOSFETs (Spitzennennwerte) \\ 4.5--5 kV & $\sim$20--30 ns & IGBT, Hochspannungs-MOSFETs \\ $>$10 kV & $\sim$100 ns -- 1 $\mu$s & Ein einzelner Halbleiter schafft es nicht mehr \\ 1 MV & $\sim$100--500 $\mu$s & Röhren, Halbleiterbaugruppen, Vervielfacher \\ \bottomrule \end{tabularx} \textit{Obwohl moderne CMOS-Transistoren bei Frequenzen um 20 GHz schalten, muss man für die Ausführung eines Prozessortakts einen Spielraum von mehreren solchen Schaltvorgängen haben, deshalb liegt die Taktfrequenz von Prozessoren im Bereich von etwa 3-4 GHz.} \subsubsection{Die Grenze des Spannungsanstiegs} Über $\sim$5 kV schafft es ein einzelner Halbleiter nicht mehr --- man muss übergehen auf: \begin{itemize} \item \textbf{Gute alte Röhren} (!) --- Vakuumtrioden, Thyratrons, \item \textbf{Halbleiterbaugruppen} --- Reihenschaltung mehrerer Transistoren, \item \textbf{Spannungsvervielfacher} --- Kaskadenschaltungen mit Röhren oder Funkenstrecken. \end{itemize} Selbst 1 Megavolt zu erzeugen ist nicht sehr schwierig (Marx-Generator, Cockcroft--Walton-Kaskaden). Eine andere Sache ist, dieses Megavolt in einer Nanosekunde ein- und auszuschalten. Dafür braucht man wirklich Hunderte von Mikrosekunden. \begin{successbox}[Fundamentale Grenze] Tatsächlich ist $\sim 10^{12}$ V/s die moderne Grenze des Spannungsanstiegs praktisch im gesamten Bereich möglicher Spannungen. Dies ist eine physikalische Beschränkung, verbunden mit: \begin{itemize} \item der Geschwindigkeit der Bewegung von Ladungsträgern in Halbleitern (Grenze $\sim 10^7$ cm/s), \item parasitären Kapazitäten und Induktivitäten, \item der Ausbreitungsgeschwindigkeit elektromagnetischer Wellen im Medium. \end{itemize} \end{successbox} \subsection{Zusammenhang mit früheren Kapiteln} Betrachten wir, wie alles, was wir durchgegangen sind, mit der ``Hardware'' zusammenhängt: \begin{itemize} \item \textbf{Lineare Algebra:} Matrixoperationen sind die Grundlage aller Berechnungen. GPUs sind genau dafür optimiert (Tensor Cores). \item \textbf{FFT:} Die schnelle Fourier-Transformation ist $\mathcal{O}(N \log N)$ Operationen, und sie ist entscheidend für NMR, Signalverarbeitung, Datenkompression. GPUs beschleunigen FFT um Zehner. \item \textbf{Iterative Methoden:} Das Verfahren der konjugierten Gradienten, GMRES --- das ist die Grundlage zur Lösung großer dünnbesetzter Systeme. Sie arbeiten auf Supercomputern mit verteiltem Speicher (MPI + OpenMP + CUDA). \item \textbf{Maschinelles Lernen:} Das Training neuronaler Netze sind Milliarden von Matrixmultiplikationen. Ohne GPUs und Tensor Cores wären moderne Modelle (GPT-4, AlphaFold) unmöglich. \item \textbf{Compressed Sensing:} Die L1-Minimierung ist eine iterative Methode, die Tausende von Iterationen erfordert. GPUs beschleunigen sie um das 100--1000-fache. \item \textbf{Digitalisierer:} ADC ist die Brücke zwischen der analogen Welt (Spektren, Signale) und der digitalen (Prozessor). Die Genauigkeit des ADC bestimmt, welche Mathematik wir anwenden können. \end{itemize} \begin{successbox}[Hauptschlussfolgerung] Moderne Rechentechnik ist: \begin{itemize} \item \textbf{Speicherhierarchie:} Register $\to$ L1 $\to$ L2 $\to$ L3 $\to$ RAM $\to$ Festplatte. Jede Stufe ist 10--100 Mal langsamer, aber 10--1000 Mal größer. \item \textbf{Parallelismus:} SIMD (Vektorinstruktionen), Mehrkernigkeit (CPU), massiver Parallelismus (GPU), Cluster (Supercomputer). \item \textbf{Spezialisierung:} Tensor Cores für KI, FPGA für spezifische Aufgaben, ASIC für Mining. \item \textbf{Physikalische Grenzen:} Memory Wall, Power Wall, Lichtgeschwindigkeit --- all dies begrenzt das Wachstum der Leistung. \end{itemize} Für einen Chemiker bedeutet dies: \begin{itemize} \item Zu verstehen, wo der ``Flaschenhals'' ist (Speicher oder Berechnungen), \item Code schreiben zu können, der Cache und GPU effizient nutzt, \item Zu wissen, wann eine Aufgabe auf einer Workstation gelöst wird und wann ein Supercomputer nötig ist. \end{itemize} \end{successbox} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Nachwort: Wie man dieses Wissen nutzt} Früher, nach dem Lesen eines ähnlichen Buches, entstand, wenn man auf eine reale Aufgabe stieß, die Notwendigkeit, eine Lösung zu implementieren. Dafür musste man eine oder mehrere Programmiersprachen recht gut beherrschen und verstehen, wie man fremde Softwaresysteme und -pakete verwendet. Im modernen Zeitalter der Entwicklung künstlicher Intelligenz kann man all dies mit Unterstützung von KI-Systemen tun --- und es wird deutlich effektiver sein. Aber dafür muss man ein wenig die wichtigsten Anwendungen von Programmiersprachen und den sogenannten \textbf{Workflow} verstehen --- wie Programmentwicklung abläuft. \subsection{Zwei Welten der Programmiersprachen} Obwohl es eine enorme Anzahl von Programmiersprachen gibt, lassen sie sich in zwei grundlegend verschiedene Klassen einteilen: \begin{tabularx}{\textwidth}{l X} \toprule \textbf{Kompiliert} & \textbf{Interpretiert / Bytecode} \\ \midrule C, C++, CUDA, Fortran, Rust & Python, JavaScript, Java, C\#, R \\ \midrule Code wird in Maschineninstruktionen für einen bestimmten Prozessor kompiliert & Code wird zur Laufzeit interpretiert oder in Bytecode kompiliert \\ \midrule Maximale Leistung, direkter Zugriff auf die ``Hardware'' & Flexibilität, Plattformunabhängigkeit, schnelle Entwicklung \\ \midrule \textbf{Wann verwenden:} & \textbf{Wann verwenden:} \\ \quad --- Schwere Berechnungen (DFT, MD, ML) & \quad --- Prototyping, Datenanalyse \\ \quad --- GPU-Beschleunigung (CUDA) & \quad --- Visualisierung, Berichte \\ \quad --- Eingebettete Systeme & \quad --- Web-Schnittstellen, Skripte \\ \bottomrule \end{tabularx} \subsubsection{Warum interpretierte Sprachen bequem sind} Interpretierte Sprachen erlauben es, größere Flexibilität beim Wechsel von einer Plattform zur anderen zu erreichen. Zum Beispiel öffnen wir dieselbe Webseite auf einem Desktop, einem Smartphone und einem Tablet --- und sehen denselben Inhalt. Dieser Inhalt wird üblicherweise mit der interpretierten Sprache JavaScript erzeugt, und wenn wir unseren Code für jeden der Prozessoren jedes Mal senden müssten, müsste der Server alle Varianten solchen Codes im Voraus haben. Dies ist manchmal sogar unmöglich vorherzusehen --- da sogar verschiedene Android-Smartphones unterschiedliche Prozessorarchitekturen haben können. \subsubsection{Typischer Workflow eines Chemikers} Um einige Ergebnisse schnell anzuzeigen, ist es oft einfacher, interpretierte Programmiersprachen zu verwenden (Python --- der absolute Standard in der Wissenschaft). Aber wenn die Aufgabe komplex ist und viele Rechenressourcen erfordert, wird es notwendig, kompilierte Programmiersprachen zu verwenden. Ein typischer Workflow sieht so aus: \begin{enumerate} \item \textbf{Prototyp in Python:} Schnell die Idee prüfen, die Daten visualisieren, verstehen, ob der Algorithmus funktioniert. \item \textbf{Optimierung:} Wenn der Code funktioniert, aber langsam ist --- die ``Engpässe'' in C++/CUDA umschreiben oder fertige Bibliotheken verwenden (NumPy, SciPy, PyTorch). \item \textbf{Produktion:} Wenn die Aufgabe regelmäßig gelöst wird --- in ein Paket packen, Tests, Dokumentation hinzufügen. \end{enumerate} \subsection{Arbeit mit KI-Assistenten} Heutzutage kann praktisch jeder Algorithmus implementiert werden, indem man einen korrekten Prompt mit Hilfe von Chats schreibt. Aber es ist zu bemerken: Jeder Chat ist fast wie ein Mensch. Und wenn man ihm eine nicht vollständig durchdachte Aufgabe gibt, kann er auch Unsinn machen oder, ohne zu verstehen, etwas nicht ganz Richtiges programmieren. \subsubsection{Goldene Regeln für die Arbeit mit KI} \begin{enumerate} \item \textbf{Zerlege die Aufgabe in kleine Blöcke.} Bitte nicht ``schreibe ein Programm zur Lösung der Schrödinger-Gleichung''. Bitte: ``schreibe eine Funktion, die die Fock-Matrix für eine gegebene Basis berechnet''. \item \textbf{Bitte um Tests.} Für jeden Block bitte den Chat, Testprüfprogramme zu schreiben, solche Tests auszuführen und sicherzustellen, dass alles in Ordnung ist. \item \textbf{Kombiniere schrittweise.} Erst nachdem jeder Block funktioniert, kombiniere sie, um das Endergebnis zu erhalten. \item \textbf{Verlange Erklärungen.} Wenn der Chat etwas geschrieben hat, das du nicht verstehst --- bitte ihn, es zu erklären. Wenn die Erklärung unklar ist --- bitte um eine einfachere. Das ist dein Skript, du musst jede Zeile verstehen. \item \textbf{Prüfe an bekannten Beispielen.} Wenn du ein Gleichungssystem löst --- prüfe an einem System, für das du die Antwort kennst. Wenn du eine Fourier-Transformation durchführst --- prüfe an einer Sinuskurve mit bekannter Frequenz. \end{enumerate} \begin{warningbox}[Typische Fehler] \begin{itemize} \item \textbf{Blindes Vertrauen:} Die KI kann Code schreiben, der ``richtig aussieht'', aber subtile Fehler enthält (zum Beispiel falsche Normierung, falsche Array-Grenzen, Speicherlecks). \item \textbf{Ignorieren des Kontexts:} Die KI kennt deine konkrete Aufgabe nicht --- du musst sie ausführlich erklären, mit Beispielen, mit Einschränkungen. \item \textbf{Fehlen von Tests:} Code ohne Tests ist Code, der im unpassendsten Moment brechen wird. \end{itemize} \end{warningbox} \subsection{Muss-Werkzeuge für einen Chemiker} Hier ist die minimale Menge von Werkzeugen, die du kennen solltest: \begin{tabularx}{\textwidth}{l X} \toprule \textbf{Sprache} & \textbf{Anwendung} \\ \midrule \textbf{Python} & Der absolute Standard: Datenanalyse, ML, Visualisierung, Skripte \\ \textbf{R} & Statistische Analyse, Bioinformatik \\ \textbf{MATLAB} & Ingenieurberechnungen, Signalverarbeitung (Alternative zu Python) \\ \textbf{C++/CUDA} & Hochleistungsrechnen, GPU \\ \textbf{Bash/Shell} & Automatisierung, Arbeit mit Clustern \\ \bottomrule \end{tabularx} \begin{tabularx}{\textwidth}{l X} \toprule \textbf{Bibliothek} & \textbf{Zweck} \\ \midrule \textbf{NumPy, SciPy} & Lineare Algebra, Optimierung, Integration \\ \textbf{Pandas} & Arbeit mit tabellarischen Daten \\ \textbf{Matplotlib, Seaborn} & Visualisierung \\ \textbf{PyTorch, TensorFlow} & Maschinelles Lernen, neuronale Netze \\ \textbf{RDKit} & Chemoinformatik, Moleküle \\ \textbf{ASE, PySCF} & Quantenchemie \\ \textbf{OpenMM, GROMACS} & Molekulardynamik \\ \bottomrule \end{tabularx} \subsection{Letzter Rat} Mathematik ist keine Sammlung von Formeln, die man auswendig lernen muss. Sie ist eine \textbf{Denkweise}. Sie ist die Fähigkeit, Struktur im Chaos zu sehen, Muster zu finden, Modelle zu bauen, Hypothesen zu prüfen. Wenn du auf eine reale Aufgabe stößt --- egal, ob es die Vorhersage einer Proteinstruktur, die Analyse eines NMR-Spektrums, die Optimierung einer Synthese oder etwas völlig Neues ist --- erinnere dich an dieses Skript. Erinnere dich daran, dass: \begin{itemize} \item jede Aufgabe entweder ein Gleichungssystem oder ein Optimierungsproblem oder ein Approximationsproblem ist, \item es für jede Aufgabe bewährte mathematische Methoden gibt, \item moderne Werkzeuge (Python, GPU, KI) es erlauben, Aufgaben zu lösen, die vor 20 Jahren unmöglich waren, \item die Hauptsache ist, das Wesen der Aufgabe zu verstehen, und der Rest ist Technik. \end{itemize} Wir glauben an euch. Ihr werdet es schaffen. Und wenn dieses Skript euch eines Tages hilft, eine Aufgabe zu lösen, die die Welt verändert --- werden wir unendlich stolz sein. \begin{successbox}[Letztes Wort] Hab keine Angst, Fehler zu machen. Hab keine Angst zu fragen. Hab keine Angst zu experimentieren. Mathematik ist keine Prüfung, sie ist ein Abenteuer. Und ihr steht am Anfang des interessantesten Weges. Viel Erfolg! \end{successbox} \vfill \begin{center} \textit{Mit Liebe, \\ Papa und Mama} \\[0.5cm] \textit{September 2026} \end{center} \end{document}