Orthogonal Greedy Algorithm (OGA)

Abstrakte Definition (Hilbertraum-Setting)

Der orthogonale Greedy-Algorithmus (OGA) ist ein iteratives Verfahren zur Approximation eines unbekannten Ziels $u$ in einem Hilbertraum

$$ (H,\langle\cdot,\cdot\rangle_H),\qquad \|v\|_H^2=\langle v,v\rangle_H. $$

Gegeben ist ein Dictionary $\mathbb D\subset H$, typischerweise symmetrisch ($g\in\mathbb D\Rightarrow -g\in\mathbb D$). Gesucht ist eine Folge $u_n\in H$, die sich aus endlich vielen Dictionary-Elementen zusammensetzt und $u$ schrittweise annähert.

$$ \begin{aligned} &u_0 = 0,\\[0.2em] &\text{für }n\ge1:\\ &\quad r_{n-1} = u - u_{n-1},\\ &\quad g_n = \arg\max_{g\in\mathbb D} \bigl|\langle g, r_{n-1}\rangle_H\bigr|,\\ &\quad V_n = \mathrm{span}\{g_1,\dots,g_n\},\\ &\quad u_n = P_{V_n}(u) = \arg\min_{v\in V_n}\|u-v\|_H. \end{aligned} $$

Der Verlust $L(v)=\tfrac12\|v-u\|_H^2$ hat Gradient $\nabla L(v)=v-u$. Damit ist das Residuum $r_{n-1}=u-u_{n-1}=-\nabla L(u_{n-1})$, und OGA wählt in jedem Schritt die Richtung $g_n$, die maximal mit dem Gradienten (bzw. Residuum) korreliert, und minimiert dann exakt im neuen Unterraum $V_n$.

Orthogonale Projektion im Unterraum

Die Projektion $u_n=P_{V_n}(u)$ lässt sich konkret formulieren, wenn $V_n=\mathrm{span}\{g_1,\dots,g_n\}$ ist. Schreibe

$$ u_n = \sum_{i=1}^n c_i g_i. $$

Die Orthogonalitätsbedingung $\langle u-u_n,g_j\rangle_H=0$ für alle $j=1,\dots,n$ führt auf das Gram-System

$$ G c = b,\qquad G_{ij} = \langle g_i,g_j\rangle_H,\quad b_i = \langle u,g_i\rangle_H. $$

Praktisch:

Bei jedem Schritt wird nur ein neuer Vektor $g_n$ hinzugefügt: $G$ wächst von $(n-1)\times(n-1)$ auf $n\times n$.
Man kann entweder
- das Gram-System bei jedem Schritt neu lösen (O(n³) in $n$), oder
- inkrementell arbeiten (z. B. Cholesky-Update) oder
- Gram–Schmidt-Orthonormalisierung durchführen und in einer ONB projizieren.

In Code ist es oft am einfachsten, die Matrix $G$ als dichte Matrix aufzubauen und pro Schritt ein kleines lineares Gleichungssystem zu lösen (z. B. numpy.linalg.solve). Die Dimension $n$ ist typischerweise viel kleiner als die Zahl der Gitter-/Quadraturpunkte.

Poisson-Fall: Diskretisierung der Energie

Für das Poisson-Problem

$$ -\Delta u = f \quad\text{in } \Omega,\qquad u = 0 \quad\text{auf } \partial\Omega $$

ist der natürliche Hilbertraum $H=H^1_0(\Omega)$ mit

$$ \langle v,w\rangle_H = \int_\Omega \nabla v\cdot\nabla w\,dx, $$

und das kontinuierliche Energie-Funktional (Dirichlet-Fall)

$$ \mathcal R(v) = \tfrac12 \int_\Omega |\nabla v|^2\,dx - \int_\Omega f\,v\,dx. $$

Zur diskreten OGA-Implementierung wählen wir Quadraturpunkte $x_i\in\Omega$ mit Gewichten $w_i>0$ (z. B. Trapezregel, Gauss–Legendre). Das diskrete Energie-Funktional ist dann

$$ \mathcal{R}_N(v) = \frac12\sum_{i=1}^{N} w_i\,|\nabla v(x_i)|^{2} - \sum_{i=1}^{N} w_i\,f(x_i)\,v(x_i). $$

Evaluationsoperator & diskretes Skalarprodukt

Für $v\in C^1(\Omega)$ speichern wir Funktionswerte und Gradienten in einem Vektor:

$$ I_{1,N}(v) = \bigl(v(x_1),\dots,v(x_N),\; \nabla v(x_1),\dots,\nabla v(x_N)\bigr)^\top \in \mathbb{R}^{N(1+d)}. $$

Auf $\mathbb R^{N(1+d)}$ definieren wir ein Skalarprodukt, das nur die Gradient-Komponenten nutzt (um die $H^1_0$-Struktur nachzubilden):

$$ \langle z,y\rangle_{H,N} := \sum_{i=1}^N w_i\, z_{N+i:\,N+i+d-1}\cdot y_{N+i:\,N+i+d-1}. $$

Datenvektor $u_N$

Die schwache Form der PDE liefert $\langle u,g\rangle_H = \int_\Omega f\,g\,dx$. Nach Quadratur lassen sich diese Skalarprodukte vollständig aus $f(x_i)$ gewinnen. Man definiert den Datenvektor

$$ u_N = \bigl( 0,\dots,0,\; f(x_1),\dots,f(x_N) \bigr)^\top \in \mathbb R^{N(1+d)}. $$

Dann gilt formal

$$ \langle u,g\rangle_H \;\approx\; \langle u_N, I_{1,N}(g)\rangle_{H,N}. $$

Die Funktionswert-Slots sind Null, weil die Energie keinen $v^2$-Massenterm enthält; in den Gradient-Slots steckt direkt die Information über $f$.

Quadratische Darstellung der diskreten Energie

Mit obigen Definitionen lässt sich $\mathcal R_N$ als quadratische Form schreiben:

$$ \mathcal R_N(v) = \tfrac12\bigl\|I_{1,N}(v) - u_N\bigr\|_{H,N}^2 + C_N, \qquad C_N = \tfrac12\sum_{i=1}^N w_i f(x_i)^2. $$

Für OGA ist wichtig: Alle benötigten Skalarprodukte können mit $u_N$ und den ausgewerteten Atomen $I_{1,N}(g)$ berechnet werden, ohne die exakte Lösung $u$ zu kennen.

Diskreter OGA – Algorithmus & Implementierung

Identifiziere jedes Dictionary-Element $g$ mit seinem Evaluationsvektor

$$ v = I_{1,N}(g) \in \mathbb R^{N(1+d)}. $$

Dann läuft der OGA auf dem quadratischen Funktional $\mathcal R_N$ wie folgt:

$$ \boxed{ \begin{aligned} &U_{0,N} = 0 \in \mathbb R^{N(1+d)},\\ &\text{für }n \ge 1:\\ &\quad R_{n-1} = u_N - U_{n-1,N} \quad\text{(Residuum)},\\[0.2em] &\quad v_n = \arg\max_{g\in\mathbb D} \bigl|\langle I_{1,N}(g), R_{n-1}\rangle_{H,N}\bigr|,\\[0.2em] &\quad V_n = \mathrm{span}\{v_1,\dots,v_n\},\\ &\quad U_{n,N} = P_{V_n}^{(H,N)}(u_N) = \sum_{i=1}^n c_i v_i,\\ &\qquad G c = b,\quad G_{ij} = \langle v_i,v_j\rangle_{H,N},\quad b_i=\langle u_N,v_i\rangle_{H,N}. \end{aligned} } $$

Praktische Implementierung (Schritt für Schritt)

Gitter & Quadratur: Knoten $x_i$, Gewichte $w_i$ erzeugen.
Datenvektor $u_N$: $f(x_i)$ auswerten, Grad-Slots belegen.
Dictionary: Funktion $\theta \mapsto g_\theta$ + Ableitung programmieren.
Evaluationsvektor: $v = I_{1,N}(g_\theta)$ als concatenierte Werte+Gradienten.
Skalarprodukt: inner_HN, das nur Gradient-Slots nutzt.
Greedy-Schritt: $\theta$-Kandidaten durchlaufen, Score $|\langle v,R\rangle_{H,N}|$ maximieren.
Projektion: Gram-Matrix $G$ und Vektor $b$ aufbauen, Gleichung lösen.
Residuum aktualisieren: $R = u_N - U_n$ und weiter mit $n+1$.

OGA nutzt die Quadratik von $\mathcal R_N$ entscheidend aus: nur dann ist das Residuum tatsächlich (bis auf ein Vorzeichen) der Gradient. Für nichtquadratische PDE-Energien ist der Algorithmus in dieser Form nicht anwendbar.

1D-Poisson-Spezialfall auf $(0,\pi)$

Modellproblem:

$$ -u''(x)=f(x)\quad\text{auf }(0,\pi),\qquad u(0)=u(\pi)=0. $$

Hilbertraum:

$$ H = H^1_0(0,\pi),\qquad \langle v,w\rangle_H = \int_0^\pi v'(x)\,w'(x)\,dx. $$

Diskretisierung:

Knoten $x_i\in(0,\pi)$, Gewichte $w_i>0$ (z. B. Trapezregel),
$I_{1,N}(v) = (v(x_1),\dots,v(x_N), v'(x_1),\dots,v'(x_N))^\top\in\mathbb R^{2N}$,
$\langle z,y\rangle_{H,N} = \sum_{i=1}^N w_i z_{N+i}y_{N+i}$,
$u_N = (0,\dots,0,f(x_1),\dots,f(x_N))^\top$.

Dictionary-Beispiel:

$$ g_{\omega,b}(x) = \psi(x)\,\sigma(\omega x + b),\qquad \sigma(z)=\max\{z,0\}^k,\quad \omega\in\{\pm1\},\ b\in[-\pi,\pi], $$

mit Formfunktion $\psi(x)=x(\pi-x)$, so dass $g_{\omega,b}\in H_0^1(0,\pi)$ und die Dirichlet-Randbedingungen exakt erfüllt sind.

Konvergenz & Raten (theoretischer Überblick)

Die klassische OGA-Theorie (z. B. im Rahmen von Siegel et al.) arbeitet mit der Metric-Entropie des Wörterbuchs. Sei $B_1(\mathbb D)$ der konvexe Rumpf von $\mathbb D$ in $H$ und

$$ \varepsilon_n(B_1(\mathbb D))_H \lesssim n^{-\tfrac12-\gamma}, \qquad \gamma>0, $$

dann gilt stark vereinfacht:

$$ \|u_n - u\|_H^2 \;\le\; \|v-u\|_H^2 + K\,\|v\|_{\mathcal K_1(\mathbb D)}^2\,n^{-1-2\gamma}, $$

für alle Vergleichsfunktionen $v$ im Sparse-Konvexkörper $\mathcal K_1(\mathbb D)$, mit einer Konstanten $K>0$. Für Wörterbücher flacher ReLU$^k$-Modelle erhält man explizite $\gamma$-Werte und damit explizite Raten $\|u_n-u\|_H = \mathcal O(n^{-\tfrac12-\gamma})$.

Wichtiger Hinweis für das hier verwendete PDE-Setting

Die obige Konvergenztheorie setzt jedoch voraus, dass

der Konvexkörper $\mathcal K_1(\mathbb D)$ bzgl. $\|\cdot\|_H$ „kontrolliert“ ist (insbesondere: Vergleichsfunktionen $v$ mit endlicher $\|v\|_{\mathcal K_1(\mathbb D)}$), und
das Dictionary $\mathbb D$ bzw. sein konvexer Rumpf in geeigneter Weise normbeschränkt/kompakt ist (Entropiezahlen endlich).

In der hier betrachteten PDE-Implementierung des OGA sind diese Voraussetzungen im Allgemeinen nicht erfüllt:

Das verwendete Wörterbuch (z. B. parametrische Atome $g_{\omega,b}$) ist nicht a priori normbeschränkt,
es wird kein explizites Vergleichselement $v\in\mathcal K_1(\mathbb D)$ mit kontrollierter $\mathcal K_1$-Norm vorausgesetzt,
somit ist die notwendige Entropie-Abschätzung für $B_1(\mathbb D)$ nicht nachgewiesen.

Konsequenz: Für den konkret implementierten OGA in diesem PDE-Setting steht keine belastbare Konvergenztheorie zur Verfügung; die oben genannten Raten sind nur als theoretische Referenz für idealisierte, gut beschränkte Dictionaries zu verstehen und lassen sich nicht direkt auf die konkrete Implementierung übertragen.

Praktische Beobachtungen (1D-Tests, grob)

Mit Formfunktion (harte Randbedingung):
Dirichlet-Randwerte werden exakt erfüllt; OGA erreicht sehr kleine Endfehler (z. B. $L^2$-Fehler bis etwa $10^{-7}$ im 1D-Test).
Mit Strafterm (Penalty):
Dirichlet wird über einen $\delta$-Penalizer erzwungen; das Fehlerplateau liegt höher (z. B. Größenordnung $10^{-5}$), konsistent mit dem theoretischen $\mathcal O(\delta)$-Randfehler.

Diese Beobachtungen sind numerisch, aber – wie oben erläutert – nicht von einer vollständigen OGA-Konvergenztheorie im verwendeten Wörterbuch-Setting abgesichert.

Git & Quelle

Repository: git.numexp.org/ferdinandkruppa/pde-oga

Clone

git clone https://git.numexp.org/ferdinandkruppa/pde-oga.git

Literatur

Siegel, Hong, Jin, Hao, Xu: Greedy Training Algorithms for Neural Networks and Applications to PDEs, arXiv:2107.04466.