Relaxed Greedy Algorithmus (RGA) – Theorie & Experimente

Worum geht es beim RGA?

Der Relaxed Greedy Algorithmus (RGA) minimiert eine konvexe Energie $$\mathcal R:H\to\mathbb R$$ über der konvexen Hülle eines Wörterbuchs $ \mathbb D \subset H $. Die Approximation $$u_n \in \Sigma_{n,M}(\mathbb D) = \left\{\sum_{i=1}^n a_i d_i : d_i \in \mathbb D,\ \sum |a_i|\le M\right\}$$ wird schrittweise aufgebaut: in jedem Schritt wird das “beste” Atom $g_n\in\mathbb D$ bzgl. der Richtungsableitung der Energie gewählt und mit einer Relaxation gemischt.

Im Unterschied zu vielen Deep-Ritz-Implementierungen gibt es für den RGA unter geeigneten Annahmen eine explizite $\mathcal O(1/n)$-Konvergenztheorie in der Energielücke. Für quadratische PDE-Energien übersetzt sich das über eine Poincaré-Ungleichung in eine $\mathcal O(n^{-1/2})$-Rate im $L^2$-Fehler.

Mathematisches Setup (Poisson & Wörterbuch)

Poisson mit Dirichlet-Randbedingungen

Referenzproblem: eindimensionales Poisson-Problem auf $[0,\pi]$ mit homogener Dirichlet-Randbedingung und Quelle $f\equiv 1$:

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

Das Dirichlet-Energiefunktional lautet $$ \mathcal R(u) = \tfrac12\!\int_0^\pi |u'(x)|^2\,dx - \int_0^\pi f(x)\,u(x)\,dx, $$ mit $H = H_0^1(0,\pi)$ und Skalarprodukt $\langle u,v\rangle_H = \int_0^\pi u'(x)v'(x)\,dx$.

Für die numerische Behandlung von Dirichlet-Randbedingungen wird oft eine Strafvariante verwendet:

$$ \mathcal R_\delta(u) = \tfrac12\!\int_0^\pi |u'|^2\,dx - \int_0^\pi f u\,dx + \frac{1}{2\delta}\big(u(0)^2 + u(\pi)^2\big), $$

mit $\delta>0$ und der zugehörigen Norm $$ \|u\|_\delta^2 = \int_0^\pi |u'|^2\,dx + \delta^{-1}\big(u(0)^2+u(\pi)^2\big). $$

Dictionary & Barron-Raum

Das Wörterbuch besteht typischerweise aus flachen Netzwerkatomen:

ReLU$^k$: $g_{w,b}(x) = \sigma_k(wx+b) = (wx+b)_+^k$, z. B. ReLU$^2$.
tanh-Atome: $g_{w,b}(x) = \tanh(wx+b)$.

Die (abgeschlossene) konvexe Hülle $B_1(\mathbb D)$ und das zugehörige Minkowski-Funktional $\|\cdot\|_{\mathcal K_1(\mathbb D)}$ definieren den Barron-Raum $$ \mathcal K_1(\mathbb D) = \{ f : \|f\|_{\mathcal K_1(\mathbb D)} <\infty\}, $$ der als Regularitätsklasse in der Analyse des RGA auftaucht. Funktionen mit $\|u^\ast\|_{\mathcal K_1(\mathbb D)}\le M$ sind für das Greedy-Verfahren grundsätzlich gut approximierbar.

Konvergenztheorie des RGA

K-Glattheit & Konvexität

Für die Konvergenz wird angenommen:

$\mathcal R$ ist konvex und Gâteaux-differenzierbar.
Der Gradient ist $K$-Lipschitz in $\|\cdot\|_H$:
$ \|\nabla\mathcal R(u)-\nabla\mathcal R(v)\|_H \le K\|u-v\|_H.$
Das Wörterbuch $\mathbb D$ ist symmetrisch ($g\in\mathbb D\Rightarrow -g\in\mathbb D$) und gleichmäßig beschränkt: $\sup_{g\in\mathbb D}\|g\|_H \le C$.

Im Poisson-Setting (auch mit Penalty) gilt:

Gradientenformel: $\nabla\mathcal R_\delta(u) = u - w_\delta$ für ein fixes $w_\delta$.
1-Glattheit in der $\delta$-Norm: $ \|\nabla\mathcal R_\delta(u)-\nabla\mathcal R_\delta(v)\|_\delta = \|u-v\|_\delta.$
$\mathcal R_\delta$ ist stark konvex in $\|\cdot\|_\delta$.

RGA-Algorithmus

Mit Budget $M>0$ und Schrittweitenplan $$\alpha_n = \min\!\bigl(1,\,\tfrac{2}{n}\bigr)$$ lautet ein Iterationsschritt:

Greedy-Auswahl: $$g_n \in \arg\max_{g\in\mathbb D}\ \langle g,\nabla\mathcal R(u_{n-1})\rangle_H,$$ mit Vorzeichenkorrektur über die Symmetrie.
Relaxiertes Update: $$u_n = (1-\alpha_n)u_{n-1} - M\alpha_n g_n.$$

Konvergenzrate

Unter den genannten Annahmen gilt (Siegel et al.):

$$\mathcal R(u_n) - \inf_{\|v\|_{\mathcal K_1(\mathbb D)}\le M} \mathcal R(v) \;\le\; \frac{32\,(CM)^2 K}{n}.$$

Für PDE-Energien mit $K=1$ und Poisson-Geometrie liefert eine Poincaré-Ungleichung:

$$\|u_n-u^\ast\|_{L^2(\Omega)}^2 \;\le\; 2C_P^2\bigl(\mathcal R(u_n)-\mathcal R(u^\ast)\bigr) \;\lesssim\; n^{-1}.$$

Für den $L^2$-Fehler bedeutet das: $\|u_n-u^\ast\|_{L^2} = \mathcal{O}(n^{-1/2})$.

Fehlerzerlegung: Modellierung, Diskretisierung, Optimierung

Praktisch wird die kontinuierliche Energie $\mathcal R$ durch ein diskretes Funktional $\mathcal R_N$ (Quadratur/Monte-Carlo) ersetzt und nur näherungsweise minimiert. Für eine numerische Lösung $\bar u_{\Theta,N}$ ergibt sich:

$$\mathcal R(\bar u_{\Theta,N})-\mathcal R(u) \;\le\; \underbrace{\bigl[\mathcal R(u_\Theta)-\mathcal R(u)\bigr]}_{\text{Modellierungsfehler}} + \underbrace{2\sup_{v\in\mathcal F_\Theta}|\mathcal R(v)-\mathcal R_N(v)|}_{\text{Diskretisierung}} + \underbrace{\bigl[\mathcal R_N(\bar u_{\Theta,N})-\mathcal R_N(u_{\Theta,N})\bigr]}_{\text{Optimierung}}.$$

Modellierungsfehler: hängt von der Approximationseigenschaft des Dictionaries $\mathbb D$ (Barron-Regularität) und dem Budget $M$ ab. Für flache ReLU$^k$-Netze gibt es Raten $\mathcal O(n^{-1/2-\gamma})$ im $H^m$.
Diskretisierungsfehler: Monte-Carlo: $\mathcal O(N^{-1/2})$ via Rademacherkomplexität; hochgradige Quadratur: $\mathcal O(N^{-(k+1)/d})$ für hinreichend glatte Integranden.
Optimierungsfehler: RGA liefert $\mathcal O(1/n)$-Abnahme in der Energie, sofern das $\arg\max$ ausreichend gut gelöst wird.

Für eine balancierte Wahl von $N$ und $n$ (z. B. $N\sim n^2$ im 1D-Fall) können alle drei Beiträge auf die Skala $\mathcal O(1/n)$ gebracht werden; dann bestimmt RGA asymptotisch den Gesamterror.

Große Konstanten: ReLU, Randstrafe & Sichtbarkeit

In der Schranke $$\mathcal R(u_n)-\inf_{\|v\|_{\mathcal K_1}\le M}\mathcal R(v) \le \dfrac{32(CM)^2K}{n}$$ steckt eine große Vorfaktorkonstante:

$C = \sup_{g\in\mathbb D}\|g\|_\delta$ hängt vom Wörterbuch ab.
Mit Dirichlet-Strafe wächst $\|g\|_\delta^2 \sim \delta^{-1}$ wegen der Randwerte.

Beispiel für ReLU$^2$-Atome $g_{w,b}(x)=(wx+b)_+^2$ auf $[0,\pi]$:

$$\|g_{w,b}\|_\delta^2 \le \frac{28}{3}\pi^3 + \frac{17}{\delta}\pi^4 =: C(\delta)^2.$$

Für $\delta = 10^{-3}$ liegt $C(\delta)$ schon im Bereich $10^3$. Bei Budget $M=10$ führt das zu $$32(CM)^2 \approx 5\cdot 10^9,$$ und die $\mathcal O(1/n)$-Kurve wird erst für astronomische $n$ sichtbar. Numerisch äußert sich das als langes Plateau, bevor irgendwann wieder ein Abfall einsetzt (wenn man sehr viele Iterationen zulässt).

Glatte Dictionaries (tanh)

Für $\tanh$-Atome $g_{w,b}(x)=\tanh(wx+b)$ mit $|w|\le W_{\max}$ ergibt sich dagegen eine viel mildere Schranke:

$$\|g_{w,b}\|_\delta^2 \le \tfrac{4}{3}|w| + \tfrac{2}{\delta} \le \tfrac{4}{3}W_{\max} + \tfrac{2}{\delta}.$$

Die abhängigen Konstanten sind deutlich kleiner, die $\mathcal O(1/n)$-Rate wird bei moderaten Iterationszahlen sichtbar. Der Preis: $\tanh$-Dictionaries haben andere (teilweise schlechtere) Approximationseigenschaften als polynomiale ReLU$^2$-Atome.

Normierung & andere Wörterbücher

Eine natürliche Idee: alle Atome in $\|\cdot\|_\delta$ auf Norm $1$ skalieren, also $ \widetilde g = g / \|g\|_\delta $. Dann ist $\sup_{\widetilde g}\|\widetilde g\|_\delta = 1$ und die Schranke wird formal “schön”: $C=1$.

Aber: die Variationsnorm $\|\cdot\|_{\mathcal K_1}$ hängt vom Wörterbuch ab. Normierung ändert die Geometrie des Budgetballs:

$$m\,\|v\|_{\mathcal K_1(\mathbb D)} \le \|v\|_{\mathcal K_1(\widetilde{\mathbb D})} \le C\,\|v\|_{\mathcal K_1(\mathbb D)},$$

und die Budgetmengen erfüllen

$$\mathcal B_{\mathbb D}(M) \subset \mathcal B_{\widetilde{\mathbb D}}(C M),\qquad \mathcal B_{\widetilde{\mathbb D}}(M) \subset \mathcal B_{\mathbb D}(M/m).$$

Um denselben modellierbaren Funktionsraum zu erhalten, muss das Budget nach Normierung um den Faktor $C$ wachsen. Normierung alleine “heilt” die Konstanten also nicht – sie verschiebt das Problem ins Budget $M$.

Numerisch: zu kleines Budget im normierten Setting ⇒ schneller Modellierungs-Stopp, Plateau. Erst ab hinreichend großem $M$ (oft deutlich höher als im unnormierten Fall) setzt wieder ein klarer Abfall des Fehlers ein.

Hüllfunktionen: Randbedingungen ohne Penalty

Alternative zu Randstrafen: jedes Atom incorporates die Randbedingung bereits über eine Formfunktion $S$, z. B. $S(x)=\sin x$ oder $S(x)=x(\pi-x)$. Statt $g(x)=\phi(wx+b)$ wird

$$g_{w,b}^S(x) = S(x)\,\phi(wx+b)$$

verwendet. Dann gilt:

$g_{w,b}^S \in H_0^1$, d. h. Dirichlet-Randwerte werden exakt erfüllt.
Kein $\delta$-Term und kein $\delta^{-1/2}$ mehr in den Konstanten.
Atome bleiben gleichmäßig beschränkt: $\|g_{w,b}^S\|_H^2 \lesssim \|S'\|_{L^2}^2 + W_{\max}^2 \|S\|_{L^2}^2$.

Die Wahl von $S$ hat aber großen Einfluss:

Günstig: $S(x)=\sin x$ oder $S(x)=x(\pi-x)$ – qualitativ ähnlich zur Lösung $\Rightarrow$ kleiner Modellierungsbias, schneller Fehlerabfall.
Ungünstig: $S(x)=\sin(2x)$ (zusätzliche innere Nullstelle) oder schräge Polynome $\Rightarrow$ der Funktionsraum wird unphysikalisch eingeschränkt, Plateau auf höherem Niveau, größeres Budget nötig.

Greedy Line-Search (nach Jaggi 2013)

Da $\mathcal R_\delta$ quadratisch ist, kann entlang jeder Richtung $g$ eine exakte eindimensionale Line-Search berechnet werden:

$$\mathcal R_\delta(u-\lambda g) = \mathcal R_\delta(u) - \lambda \langle \nabla\mathcal R_\delta(u), g\rangle_\delta + \frac{\lambda^2}{2}\|g\|_\delta^2,$$ $$\lambda^\star = \frac{\langle \nabla\mathcal R_\delta(u), g\rangle_\delta}{\|g\|_\delta^2}.$$

Greedy-Variante:

Wähle $g_k\in\pm\mathbb D$ mit maximaler Korrelation wie im RGA.
Setze $\lambda_k^\star$ wie oben.
Clippe $\lambda_k$, so dass die Variationsnorm $k_1(u_k)$ das Budget $M$ nicht überschreitet.
Update: $u_{k+1}=u_k - \lambda_k g_k$.

Diese Variante ist eng mit Frank-Wolfe-Methoden mit exakter Line-Search verwandt (vgl. Jaggi, 2013). In den Experimenten:

Sehr schneller anfänglicher Fehlerabfall.
Plateaus durch Penalty-Bias und Budgetgrenzen.
Für kleinere $\delta$ sinkt das Plateau-Niveau, teilweise mit zweiter “Abstiegsphase”.

Matrix- & GPU-Implementierung des RGA

Der $\arg\max$-Schritt dominiert die Kosten: für jedes Atom $g_j$ wird $$r_j = \mathcal R_\delta'(u)[g_j]$$ benötigt. Diskret am Gitter $(x_i)_{i=1}^N$ mit Gewichten $(w_i)$ gilt (schematisch)

$$r_j \approx \sum_i w_i\,u'_i\,g'_{j,i} - \sum_i w_i\,f_i\,g_{j,i} + \frac{1}{\delta}\bigl(u(0)g_j(0)+u(\pi)g_j(\pi)\bigr).$$

In Matrixform mit

$G = [g_1,\dots,g_P]\in\mathbb R^{N\times P}$,
$G'=[g_1',\dots,g_P']$,
Diagonalmatrix $W=\mathrm{diag}(w_1,\dots,w_N)$,
Randprojektor $E = e_0e_0^\top+e_\pi e_\pi^\top$,

wird daraus

$$r = (G')^\top W u' - G^\top W f + \frac{1}{\delta} G^\top E u.$$

Alle Scores $r_j$ entstehen also durch zwei Matrix-Vektor-Produkte plus einen dünn besetzten Randterm – perfekt für GPU-Beschleunigung.

Die Atome selbst werden vektorisiert erzeugt:

Voraktivierungen $Z_{i,j}=w_j x_i + b_j$.
Aktivierungen $G_{i,j}=\phi(Z_{i,j})$, Ableitungen $G'_{i,j}=\phi'(Z_{i,j})\,w_j$.

Dadurch lassen sich Millionen Iterationen fahren, sodass die asymptotische $\mathcal O(1/n)$-Phase trotz großer Konstanten numerisch erreichbar und sichtbar wird.

Numerische Experimente – Überblick

Die Experimente prüfen systematisch:

Verhältnis Quadraturfeinheit $N$ vs. Iterationszahl $n$ (kein Quadratur-Floor bis $n\ll N^2$).
Effekt exakter vs. approximativer $\arg\max$-Suche.
Konstanten $C(\delta)$ für ReLU$^2$ vs. $\tanh$-Dictionary und deren Einfluss auf Plateaus.
Normierung der Atome und notwendige Budget-Reskalierung.
Formfunktionen $S(x)=\sin x$, $\sin(2x)$, $x^2(x-\pi)$
Greedy Line-Search vs. klassischer RGA.

Typische Plots (3×2-Grid):

$|\langle g,\nabla\mathcal R(u_k)\rangle|$ vs. $k$.
Variationsnorm $k_1(u_k)$ relativ zu $M$.
$L^2$-Fehler mit Referenzneigungen ($n^{-1},n^{-1/2}$).
Profilvergleich $u_k$ vs. $u^\ast$.
Argmax-Spuren der Parameter (z. B. $(\omega,b)$).

Auf dieser Basis lässt sich sauber unterscheiden, ob ein Plateau von Optimierung (argmax, Konstanten, Budget) oder von Diskretisierung/Modellierung (Quadratur, Penalty-Bias, Wahl von $S$) stammt.

Code & Repositories

Die in dieser Seite beschriebenen RGA-Varianten sowie GPU-vektorisierte 1D-Implementierungen sind in folgenden Git-Repositories verfügbar:

RGA 1D, sin-Hülle, GPU-basiertes b-Grid: rga1d_sin_bgrid_gpu
Komplettes PDE-RGA-Projekt (RGA, OGA, Line-Search, 1D/2D, Debug-Plots): pde-relaxed-greedy

Clone-Befehle

git clone https://git.numexp.org/ferdinandkruppa/rga1d_sin_bgrid_gpu.git

git clone https://git.numexp.org/ferdinandkruppa/pde-relaxed-greedy.git

Quellen

Siegel, J. W.; Hong, Q.; Jin, X.; Hao, W.; Xu, J.: Greedy Training Algorithms for Neural Networks and Applications to PDEs, arXiv:2107.04466 (rev. 2023).
Siegel, J. W.; Xu, J.: Approximation Rates for Neural Networks in Barron Spaces, Variation Spaces Preprint, 2023.
Jaggi, M.: Revisiting Frank-Wolfe: Projection-Free Sparse Convex Optimization, ICML 2013.
E, W.; Yu, B.: The Deep Ritz Method: A Deep Learning-Based Numerical Algorithm for Solving Variational Problems, 2018.