△ Mechanism confirmed, baseline not beaten 2026

Residual-Christoffel Collocation for Random-Feature PDE Networks

Implementation & benchmark of arXiv:2607.13382 — Residual-Christoffel Sampling for Random Feature Collocation of Linear PDEs

Usefulness8/10
Difficulty5/10
Novelty7/10

Source paper: Residual-Christoffel Sampling for Random Feature Collocation of Linear PDEs arXiv:2607.13382 ⓘ · analyzed Aug 30, 2026

AI-generated research hypothesis, automatically tested. Not peer-reviewed.

Idea description

Replace uniform collocation for a fixed random-feature neural PDE solver with sampling from the leverage-score density of the operator-applied features. Whiten the retained residual feature space before solving for output coefficients, so the sampled least-squares matrix has an identity-like expected Gram rather than inheriting severe anisotropy from the differential operator. The same construction can be used for a linearized neural network by treating Jacobian features as the trial functions.

Formulas

$$u_M(x)=\sum_{m=1}^{M}c_m\varphi_m(x),\qquad \psi(x)=\big(L\varphi_1(x),\ldots,L\varphi_M(x)\big)^\top.$$
$$G=\int_\Omega\psi(x)\psi(x)^\top\,d\mu(x)=U\Lambda U^\top,\qquad z(x)=\Lambda_r^{-1/2}U_r^\top\psi(x),\qquad q(x)=\frac{\|z(x)\|_2^2}{r}.$$
$$x_i\sim q(x)\,d\mu(x),\qquad A_i=\frac{z(x_i)^\top}{\sqrt{Nq(x_i)}},\qquad b_i=\frac{f(x_i)}{\sqrt{Nq(x_i)}},\qquad \mathbb{E}[A^\top A]=I_r.$$
$$G_\lambda=G+\lambda I,\qquad z_\lambda(x)=G_\lambda^{-1/2}\psi(x),\qquad d_\lambda=\operatorname{tr}\!\left(GG_\lambda^{-1}\right),\qquad q_\lambda(x)=\frac{\|z_\lambda(x)\|_2^2}{d_\lambda}.$$

Mathematical statement

For fixed trial functions $\{\varphi_m\}_{m=1}^M$, the paper defines the residual feature vector $\psi(x)=(L\varphi_1(x),\ldots,L\varphi_M(x))^\top$, where $L$ is the interior PDE operator. Let $\mu$ be a reference sampling measure on the interior, or a mixture of interior and boundary measures when boundary rows are included. The reference residual Gram is $G=\int_\Omega\psi(x)\psi(x)^\top d\mu(x)$. With eigendecomposition $G=U\Lambda U^\top$, retain the $r$ eigenvectors whose eigenvalues exceed threshold $\tau$; denote the retained factors by $U_r$ and $\Lambda_r$. Define whitened residual features $z(x)=\Lambda_r^{-1/2}U_r^\top\psi(x)$ and the residual-Christoffel density $q(x)=\|z(x)\|_2^2/r$ relative to $\mu$. Since $\int q\,d\mu=1$, sampling $x_i\sim q(x)d\mu(x)$ and multiplying each row by $(Nq(x_i))^{-1/2}$ gives $\mathbb{E}[A^\top A]=I_r$ for $A_i=(Nq(x_i))^{-1/2}z(x_i)^\top$. The transformed coefficients satisfy $c=U_r\Lambda_r^{-1/2}\theta$, so the PDE residual is $\psi(x)^\top c=z(x)^\top\theta$. A ridge-stable alternative uses $G_\lambda=G+\lambda I$, effective dimension $d_\lambda=\operatorname{tr}(GG_\lambda^{-1})$, and $q_\lambda(x)=\|G_\lambda^{-1/2}\psi(x)\|_2^2/d_\lambda$.

Implementation notes

1. Integration point: use a fixed random-feature MLP or ELM for a PDE, with trainable output coefficients only. For each candidate point, compute the feature vector $\varphi(x)$ and its differentiated vector $\psi(x)=L\varphi(x)$ using automatic differentiation. If enforcing boundary conditions strongly is impossible, concatenate interior rows $L\varphi(x)$ and boundary rows $B\varphi(x)$, and define $\mu$ as a mixture with an explicit interior/boundary sampling probability.

2. Pseudocode: draw a pilot set $\{\tilde x_j\}_{j=1}^P$ uniformly from $\mu$; form $\Psi_j=\psi(\tilde x_j)$ and estimate $G\leftarrow P^{-1}\sum_j\Psi_j\Psi_j^\top+\epsilon I$. Either eigendecompose $G$ and retain $\lambda_k\ge\tau$, then compute $z_j=\Lambda_r^{-1/2}U_r^\top\Psi_j$ and scores $s_j=\|z_j\|^2/r$, or use the ridge version $z_j=(G+\lambda I)^{-1/2}\Psi_j$ and $s_j=\|z_j\|^2/d_\lambda$. Normalize scores over the pilot set, sample collocation indices proportional to them, evaluate fresh rows, multiply row $i$ and target by $1/\sqrt{Nq_i}$, and solve the whitened least-squares problem for $\theta$. Recover output coefficients with $c=U_r\Lambda_r^{-1/2}\theta$.

3. Computed from the mathematics: the Gram estimate, eigenspace or ridge inverse square root, Christoffel scores, inverse-density weights, and whitened coefficient parameterization. Estimated empirically: the pilot approximation to $G$, the threshold $\tau$ or ridge $\lambda$, and the mixture weight between interior and boundary rows. Clip $q_i$ below by $q_{\min}$ to avoid exploding weights, and compare clipping levels in an ablation.

4. First experiment: solve 1D Poisson, 2D Poisson, and a variable-coefficient elliptic equation using the same 2-layer random-feature MLP and identical feature initialization. Compare uniform collocation, residual-magnitude adaptive sampling, residual-Christoffel sampling without whitening, and the full weighted-whitened method. Sweep the number of collocation points and record smallest singular value, condition number, QR or conjugate-gradient iterations, relative $L^2$ solution error, and PDE residual. The expected signal is a flatter condition-number curve and lower error at small sample counts; a compelling result would be reaching a fixed relative error with 2x fewer rows or reducing linear-solve iterations by at least 2x at equal feature count.

Verification

Mechanism confirmed, baseline not beaten

Stage 1 · Toy mechanism gate: Passed ✓

Stage 2 · Mechanism transferred to benchmark: Confirmed ✓

Stage 2 · Practical benchmark result: Significantly worse ✗

Methodology: 8 paired seeds, permutation test (p<0.05) against a learning-rate-tuned baseline. How verification works

Stage 1 — Mechanism check agent confidence 8/10

Built a fixed random-feature 1D Poisson solver with residual-Christoffel sampling, inverse-density weighting, eigenspace whitening, and a uniform-collocation baseline. The mechanism was confirmed: q integrated to 1.000, the population whitened Gram error was 3.97e-9, the direct q-sampling expectation error was 0.0286, and finite-sample Gram error followed the predicted N^-1/2 law with fitted slope -0.496. Conditioning improved dramatically, but PDE L2 solution error was worse than uniform sampling in this tiny setup, likely because boundary rows and aggressive residual weighting were not tuned.

Agent confidence
8/10
Baseline
Uniform raw residual condition median: 5.26e8 at N=14, 1.98e16 at N=28, and 2.02e12 at N=224; PDE L2 median: 2.48e-2, 1.58e-6, 4.27e-9, and 1.42e-9 for N=14,28,56,112.
Idea
Christoffel weighted-whitened condition median: 6.25e3, 18.3, 2.50, and 1.53 at N=14,28,56,112, with Gram-error slope -0.496; PDE L2 median: 3.58e-2, 1.70e-3, 9.15e-4, and 4.18e-4.

Limitations: Only a 1D Poisson problem and analytic tanh second derivatives were tested; no 2D or variable-coefficient PDE, autodiff Jacobian features, residual-magnitude adaptive baseline, q clipping ablation, timing/FLOP comparison, or independently estimated pilot Gram was evaluated. Boundary weighting was fixed rather than tuned, and the pilot distribution used a dense quadrature grid.

How to run: python3 experiment.py

Stage 2 — Benchmark vs. tuned baseline (arXiv:2607.13382)

Significantly WORSE than the tuned baseline (+51.1%, p=0.0081, wins 0 / 8 paired seeds).

Benchmark
Poisson_dirichlet_local
Model
mlp
Paired seeds
8
Baseline mean
0.2885
Idea mean
0.4359
Effect (Δ)
0.1474 (+51.1%; negative = idea better)
Wins
0 / 8 paired seeds
p-value
0.0081 (permutation test, 20 000 shuffles)
Smallest detectable effect
±16.5%
Mechanism
Confirmed ✓
Practical effect
Hurts
Baseline tuning
swept over 3 configs
Limitations:

The experiment used a small 1D Poisson proxy with 32 fixed random features, 96 collocation rows, 80 epochs, and three learning rates. It did not test 2D or variable-coefficient PDEs, larger networks, boundary-weight sweeps, q-clipping ablations, pilot-size sensitivity, or timing/FLOPs.

How to run: python3 bench_run.py

Verdict computed by deterministic test code from paired-seed statistics — not by the language model.

Stage 2 — Benchmark latest Failed ✗

Agent confidence: 10/10

Residual-Christoffel sampling with whitening was implemented and evaluated on the registered poisson_dirichlet_local PDE track using 8 paired seeds and a shared three-value learning-rate sweep. It reduced the trained-system median design condition from 3.34e10 to 3.26e4, confirming the mechanism, but significantly worsened independent test MSE; the benchmark verdict is idea worse (significant).

Baseline
Best lr=0.1, mean test MSE 0.2884851713 across 8 seeds.
Idea
Best lr=0.1, mean test MSE 0.4358671680 across 8 seeds; paired delta +0.1473819967, 0/8 wins, permutation p=0.0081.
Mechanism
Condition ratio idea/baseline=9.755321e-7; confirmed=true. Residual RMS was 5.61245 baseline versus 6.55377 idea.

Benchmark result

Significantly WORSE than the tuned baseline (+51.1%, p=0.0081, wins 0 / 8 paired seeds).

Benchmark
Poisson_dirichlet_local
Model
mlp
Paired seeds
8
Baseline mean
0.2885
Idea mean
0.4359
Effect (Δ)
0.1474 (+51.1%; negative = idea better)
Wins
0 / 8 paired seeds
p-value
0.0081 (permutation test, 20 000 shuffles)
Smallest detectable effect
±16.5%
Mechanism
Confirmed ✓
Practical effect
Hurts
Baseline tuning
swept over 3 configs
Limitations:

The experiment used a small 1D Poisson proxy with 32 fixed random features, 96 collocation rows, 80 epochs, and three learning rates. It did not test 2D or variable-coefficient PDEs, larger networks, boundary-weight sweeps, q-clipping ablations, pilot-size sensitivity, or timing/FLOPs.

How to run: python3 bench_run.py

Artifacts

Implementation overview ⬇ Download all as ZIP 7 files · code, reports and structured results