Orthogonal Greedy Algorithm (OGA) for PDEs

This project implements the Orthogonal Greedy Algorithm (OGA) to solve 1D partial differential equations (PDEs), specifically the Poisson problem

-u'' = f \text{ on } [0,\pi].

The algorithm builds a sparse approximation of the solution using a dictionary of parameterized atoms (e.g. ReLU$^3$ activations). The implementation follows the theory in:

J. W. Siegel, Q. Hong, X. Jin, W. Hao, and J. Xu (2023), "Greedy training algorithms for neural networks and applications to PDEs," Journal of Computational Physics, 484, 112084.

Two strategies for homogeneous Dirichlet boundary conditions are supported and compared:

1. Shape functions (exact Dirichlet): Each atom is multiplied by a function that vanishes at the boundary.
2. Penalty method (Robin/penalized): Boundary values are discouraged via a Robin-type penalty term.

The main implementation is in OGA_PDE1d.py. Runtime options allow switching boundary treatment (use_form) and generating diagnostic plots (errors, convergence indicators, parameter traces).

---

1. Installation & Setup

git clone <repository-url>
cd OGA-PDE
python -m venv .venv
source .venv/bin/activate        # macOS/Linux
# or .venv\\Scripts\\activate   # Windows
pip install -r requirements.txt  # or: pip install numpy scipy matplotlib

Dependencies: numpy, scipy, matplotlib.

---

2. Problem Setting and Energy

The model problem is the 1D Poisson equation with homogeneous Dirichlet boundary conditions.

The energy functional is

\mathcal{R}(v) = \tfrac12 \int_0^{\pi} (v'(x))^2,dx - \int_0^{\pi} f(x),v(x),dx.

The exact solution u_* is the unique minimizer of \mathcal{R} in H_0^1(0,\pi), i.e. functions with u(0)=u(\pi)=0.

Instead of solving a global discretization (FEM etc.), OGA builds a sparse expansion

u_n(x) = \sum_{k=1}^n c_k,g_k(x)

with atoms of the form

g_{\omega,b}(x) = \psi(x),\sigma(\omega x + b),

where \sigma is an activation (e.g. ReLU$^p$, \tanh, \sin), \omega, b are parameters chosen adaptively, and \psi controls boundary behaviour.

---

3. Orthogonal Greedy Algorithm (OGA)

Let H be a Hilbert space with inner product \langle \cdot,\cdot \rangle_H, and let \mathbb{D} \subset H be a symmetric dictionary. Denote the target solution by u_*. OGA constructs (u_n)_{n\ge0} as follows:

1. Initialization: u_0 = 0.
2. Greedy direction: choose

g_n = \operatorname*{argmax}*{g\in\mathbb{D}} \big|\langle g, r*{n-1} \rangle_H\big|,

where r_{n-1} = u_* - u_{n-1}. For quadratic energies, r_{n-1} is proportional to the negative gradient of \mathcal{R} at u_{n-1}, so this is equivalent to maximizing |\mathcal{R}'(u_{n-1})[g]|.

3. Span update: V_n = \operatorname{span}{g_1,\dots,g_n}.
4. Full correction (orthogonal projection):

u_n = P_{V_n}(u_*) = \operatorname*{argmin}*{v\in V_n} |u** - v|_H.

In practice this projection is obtained by solving a linear system A c = b, where the Gram matrix is

A_{ij} = \langle g_i, g_j \rangle_H,

and the right-hand side is

b_i = \langle u_*, g_i \rangle_H.

Key point: Unlike the Relaxed Greedy Algorithm (RGA), OGA re-optimizes all coefficients at every step ("full correction"), which accelerates convergence.

---

4. Specialization to Poisson in 1D

For the Poisson functional,

\mathcal{R}(v) = \tfrac12 a(v,v) - L(v)

with bilinear form and load functional

a(v,w) = \int_0^{\pi} v'(x) w'(x),dx,

L(w) = \int_0^{\pi} f(x) w(x),dx,

the Gâteaux-Differential ist

\mathcal{R}'(u)[\phi] = a(u,\phi) - L(\phi).

The greedy step in code is therefore

g_n = \operatorname*{argmax}*{g\in\mathbb{D}} |\mathcal{R}'(u*{n-1})[g]|.

After choosing g_n, the algorithm assembles

A_{ij} \approx \int_0^{\pi} g_i'(x) g_j'(x),dx ; + ; \text{(boundary terms if penalized)},

b_i \approx \int_0^{\pi} f(x) g_i(x),dx,

solves A c = b, and forms

u_n(x) = \sum_{i=1}^n c_i g_i(x).

Numerical integration is done by the trapezoidal rule on a uniform grid with N points in [0,\pi].

---

5. Boundary Handling

Boundary enforcement is controlled by the flag use_form.

5.1 Shape-function approach (use_form=True)

• Each atom is

g(x) = \psi(x),\sigma(\omega x + b), \quad \psi(x)=x(\pi-x).

• Because \psi(0)=\psi(\pi)=0, every atom satisfies homogeneous Dirichlet boundary conditions.
• Any linear combination of atoms lies in H_0^1(0,\pi).
• This exactly solves the correct Dirichlet problem.
• The penalty weight \delta is irrelevant in this mode.

Effect: no boundary bias, correct PDE model.

5.2 Penalty / Robin approach (use_form=False)

• Here \psi(x)\equiv 1, so boundary values are not forced to zero.
• The energy is modified to

$\mathcal{R}_\delta(u) = \tfrac12 \int_0^{\pi} (u'(x))^2,dx ; - \int_0^{\pi} f(x) u(x),dx ; + \tfrac{1}{2\delta}\big(u(0)^2 + u(\pi)^2\big).$

• The Euler--Lagrange system becomes interior Poisson -u''=f with Robin-type boundary conditions

u'(0) = \tfrac{1}{\delta} u(0), \quad u'(\pi) = -\tfrac{1}{\delta} u(\pi).

• As \delta \to 0, these approach homogeneous Dirichlet conditions.

Effect:

• Easy atoms (no shape multiplier).
• But for fixed \delta>0, the algorithm converges to the minimizer of \mathcal{R}*\delta, not to the true Dirichlet solution u**. This introduces an \mathcal{O}(\delta) bias in L^2 near the boundary.
• In practice this appears as an early error plateau whose level is proportional to \delta.

---

6. Convergence Behavior

The convergence rate of OGA is governed by the metric entropy \varepsilon_n(\mathbb{D})_H of the dictionary. For shallow neural-network-type dictionaries (affine parameters, ReLU$^p$ etc.) one typically has

\varepsilon_n(\mathbb{D})_H \sim n^{-1/2-\gamma} with \gamma>0.

This implies for the energy error

\mathcal{R}(u_n) - \mathcal{R}(u_*) = \mathcal{O}(n^{-1-2\gamma}),

and for the $L^2$-error

|u_n - u_*|_{L^2} = \mathcal{O}(n^{-1/2-\gamma}).

In computations:

• With use_form=True, the L^2 error decays extremely fast and eventually hits a floor dictated by quadrature accuracy (N large \Rightarrow floor very small).
• With use_form=False, and e.g. \delta=10^{-5}, the L^2 error quickly drops and then plateaus around 10^{-5}. This plateau is expected: it reflects the \mathcal{O}(\delta) boundary bias of the penalized Robin problem, not a failure of OGA.

---

7. Code Structure

• OGA_PDE1d.py

  • main class GreedyPDE1D_OGA_PenalizedTestPlots (OGA loop, boundary handling, projection, residual evaluation, logging, plotting).

• OGA_PDE_test.py

  • minimal driver / experiment configuration (problem setup, solver parameters, plotting).

• OGA_PDE_test.ipynb (optional)

  • interactive exploration and visualization.

• main.py (optional)

  • alternative scripted run.

---

8. Key Takeaways

• OGA builds a sparse expansion from a neural-network-like dictionary and re-optimizes all coefficients after every new atom (full correction), which accelerates convergence.
• use_form=True enforces u(0)=u(\pi)=0 exactly via the shape function \psi(x)=x(\pi-x) and solves the true Dirichlet Poisson problem in H_0^1(0,\pi).
• use_form=False uses a Robin-style penalty with weight \delta and therefore solves a modified boundary value problem; this produces an \mathcal{O}(\delta) error floor in L^2.
• Numerically, L^2 errors decay extremely fast, often outperforming classical greedy methods like RGA, which typically only achieve \mathcal{O}(n^{-1}) energy decay.

---

9. Citation

Please cite:

J. W. Siegel, Q. Hong, X. Jin, W. Hao, and J. Xu (2023). Greedy training algorithms for neural networks and applications to PDEs. Journal of Computational Physics, 484, 112084. doi: 10.1016/j.jcp.2023.112084