Collocation-Point Resampling Methods

• Collocation-point resampling is a family of techniques that transform the collocation set using explicit barycentric interpolation and adaptive repositioning to enhance numerical accuracy.
• The method adapts both spectral and PINN-based approaches by transferring values via resampling matrices or concentrating collocation points in regions with high residuals or curvature.
• Practical implementations leverage fixed low-discrepancy sets, Hessian- and NTK-based metrics, and oversampling strategies to balance computational cost with improved convergence and error reduction.

Searching arXiv for papers on collocation-point resampling and closely related methods. Collocation-point resampling denotes a family of techniques that modify where a collocation method enforces its governing equations. In one line of work, resampling is an interpolation operation: values known on a source collocation grid are transferred to delayed, shifted, or otherwise transformed target points through a barycentric resampling matrix, yielding global spectral discretizations for functional and delay differential equations (Hale, 2024). In another line of work, especially for physics-informed neural networks (PINNs), resampling means replacing or moving the collocation set itself during training so that residual evaluations are concentrated in regions judged most informative by residual magnitude, curvature, Hessian structure, Neural Tangent Kernel criteria, reduced-order residual dynamics, or direct gradient ascent on the squared residual (Florido et al., 2024, Caradot et al., 1 Apr 2025, Visser et al., 2024, Lau et al., 2024, Celaya et al., 13 Jan 2025). Related uses also appear in rectangular collocation for the Kohn–Sham equation, where machine learning reduces the collocation set size while retaining mHartree accuracy, and in boundary element collocation, where oversampling produces rectangular systems that improve robustness and sometimes convergence (Ku et al., 2019, Maierhofer et al., 2021).

1. Definitions and major variants

The unifying idea is to treat the collocation set as an object that can be transformed rather than fixed. For classical spectral collocation, this transformation is algebraic. Let $\{x_j\}_{j=0}^N$ be the Chebyshev–Gauss–Lobatto nodes on $[-1,1]$,

$$x_j=-\cos\Bigl(\frac{j\pi}{N}\Bigr),\qquad j=0,\dots,N.$$

Given nodal data $y_j=y(x_j)$, the degree-$N$ interpolant can be written in the second barycentric form,

$$p(x)=\frac{\displaystyle\sum_{j=0}^N \frac{w_j}{x-x_j}\,y_j}{\displaystyle\sum_{j=0}^N \frac{w_j}{x-x_j}},$$

with barycentric weights

$$w_0=w_N=\frac12(-1)^j,\qquad w_j=(-1)^j,\quad j=1,\dots,N-1.$$

Resampling consists of evaluating this interpolant at target points $\xi_i$ through a matrix $R$ whose entries are $R_{i,j}=\ell_j(\xi_i)$, where $[-1,1]$0 are the Lagrange basis functions (Hale, 2024).

For PINNs, the term usually refers to adaptive replacement of the set $[-1,1]$1 of collocation points used in the residual loss

$[-1,1]$2

The objective is to relocate a fixed number of points toward regions where the current network estimate is least accurate, thereby improving the objective used in subsequent optimization (Florido et al., 2024). This can be implemented as probability-based resampling from a candidate pool, deterministic low-discrepancy initialization followed by adaptive replacement, or direct continuous movement of points by gradient-based optimizers (Florido et al., 2024, Visser et al., 2024).

A broader collocation literature uses the same principle in other ways. Rectangular collocation for the Kohn–Sham equation begins with a very large spatial grid and then filters and optimizes a much smaller collocation set via Gaussian process regression and a genetic algorithm (Ku et al., 2019). Oversampled boundary element collocation uses $[-1,1]$3 collocation points for an $[-1,1]$4-dimensional trial space and solves the resulting $[-1,1]$5 system in weighted least squares form (Maierhofer et al., 2021). This suggests that “collocation-point resampling” is best understood as a methodological umbrella rather than a single algorithm.

2. Spectral collocation resampling via barycentric interpolation

In the Chebyshev-spectral setting, the resampling operator is explicit. If values are known on the source grid $[-1,1]$6 and interpolation is required at target points $[-1,1]$7, the barycentric resampling matrix is

$[-1,1]$8

The paper gives a short MATLAB routine, $[-1,1]$70 which forms the matrix entry-by-entry and handles the diagonal case $[-1,1]$9 by setting $$x_j=-\cos\Bigl(\frac{j\pi}{N}\Bigr),\qquad j=0,\dots,N.$$0 (Hale, 2024).

This resampling matrix is integrated with the pseudospectral differentiation matrix $$x_j=-\cos\Bigl(\frac{j\pi}{N}\Bigr),\qquad j=0,\dots,N.$$1, defined by differentiating the Lagrange basis at the collocation points,

$$x_j=-\cos\Bigl(\frac{j\pi}{N}\Bigr),\qquad j=0,\dots,N.$$2

For the delay differential equation

$$x_j=-\cos\Bigl(\frac{j\pi}{N}\Bigr),\qquad j=0,\dots,N.$$3

one maps Chebyshev points to $$x_j=-\cos\Bigl(\frac{j\pi}{N}\Bigr),\qquad j=0,\dots,N.$$4, forms $$x_j=-\cos\Bigl(\frac{j\pi}{N}\Bigr),\qquad j=0,\dots,N.$$5, the identity $$x_j=-\cos\Bigl(\frac{j\pi}{N}\Bigr),\qquad j=0,\dots,N.$$6, and the resampling matrix

$$x_j=-\cos\Bigl(\frac{j\pi}{N}\Bigr),\qquad j=0,\dots,N.$$7

The collocation-point discretization enforces

$$x_j=-\cos\Bigl(\frac{j\pi}{N}\Bigr),\qquad j=0,\dots,N.$$8

or in matrix form,

$$x_j=-\cos\Bigl(\frac{j\pi}{N}\Bigr),\qquad j=0,\dots,N.$$9

The boundary condition is imposed by row-replacement: the first row of $y_j=y(x_j)$0 is replaced by $y_j=y(x_j)$1 and $y_j=y(x_j)$2 by $y_j=y(x_j)$3, after which one solves the $y_j=y(x_j)$4 system $y_j=y(x_j)$5 (Hale, 2024).

The same construction extends to more general functional terms $y_j=y(x_j)$6 by building $y_j=y(x_j)$7, and it also extends to state-dependent delays by including $y_j=y(x_j)$8 and its Jacobian inside a Newton iteration (Hale, 2024). Because the method is a natural extension of standard spectral collocation based on polynomial interpolants, it can be incorporated into existing spectral discretizations such as Chebfun/Chebop (Hale, 2024).

The principal significance of this version of resampling is that it does not seek new collocation points; it transports an existing polynomial representation to functionally transformed arguments. A common misconception is therefore to equate all collocation-point resampling with adaptive point selection. In the spectral DDE/FDE literature, resampling is instead an interpolation operator embedded in the discretization itself.

3. Multidomain treatment, discontinuities, and convergence in spectral schemes

When the history function or the delay induces jumps in derivatives at known times $y_j=y(x_j)$9, the interval $N$0 is split into subintervals

$N$1

Each subinterval receives its own Chebyshev grid $N$2, barycentric weights $N$3, differentiation matrix $N$4, and identity $N$5 (Hale, 2024). The functional evaluation operator becomes a block matrix $N$6 whose $N$7 block is either zero or a local resampling matrix

$N$8

so that

$N$9

Continuity at an interface $$p(x)=\frac{\displaystyle\sum_{j=0}^N \frac{w_j}{x-x_j}\,y_j}{\displaystyle\sum_{j=0}^N \frac{w_j}{x-x_j}},$$0 is enforced by row replacement using

$$p(x)=\frac{\displaystyle\sum_{j=0}^N \frac{w_j}{x-x_j}\,y_j}{\displaystyle\sum_{j=0}^N \frac{w_j}{x-x_j}},$$1

The resulting global blocked system is square and solves for all nodal values simultaneously (Hale, 2024).

The convergence theory stated in the paper is geometric under analyticity. If the exact solution is analytic in a Bernstein ellipse of radius $$p(x)=\frac{\displaystyle\sum_{j=0}^N \frac{w_j}{x-x_j}\,y_j}{\displaystyle\sum_{j=0}^N \frac{w_j}{x-x_j}},$$2 about $$p(x)=\frac{\displaystyle\sum_{j=0}^N \frac{w_j}{x-x_j}\,y_j}{\displaystyle\sum_{j=0}^N \frac{w_j}{x-x_j}},$$3, then the Chebyshev interpolant satisfies

$$p(x)=\frac{\displaystyle\sum_{j=0}^N \frac{w_j}{x-x_j}\,y_j}{\displaystyle\sum_{j=0}^N \frac{w_j}{x-x_j}},$$4

with analogous behavior for derivatives (Hale, 2024). This carries over to smooth spectral-collocation solutions of ODEs and DDEs built from such interpolants. If the solution is piecewise analytic with jumps in its $$p(x)=\frac{\displaystyle\sum_{j=0}^N \frac{w_j}{x-x_j}\,y_j}{\displaystyle\sum_{j=0}^N \frac{w_j}{x-x_j}},$$5-th derivative and subinterval boundaries are introduced at those jumps, the same $$p(x)=\frac{\displaystyle\sum_{j=0}^N \frac{w_j}{x-x_j}\,y_j}{\displaystyle\sum_{j=0}^N \frac{w_j}{x-x_j}},$$6 geometric convergence is recovered on each subinterval (Hale, 2024). By contrast, if no domain splitting is used in the presence of a derivative jump, the observed convergence is only algebraic, $$p(x)=\frac{\displaystyle\sum_{j=0}^N \frac{w_j}{x-x_j}\,y_j}{\displaystyle\sum_{j=0}^N \frac{w_j}{x-x_j}},$$7, because of Gibbs-type effects (Hale, 2024).

These statements place resampling within the standard spectral narrative: the resampling matrix itself is not a cure for nonsmoothness, but it preserves the global high-order structure once discontinuities are isolated. A plausible implication is that in this setting the decisive design variable is not merely the interpolation formula, but the compatibility between resampling and domain decomposition.

4. Adaptive resampling in PINNs: probability densities and guiding metrics

In PINNs, adaptive collocation-point resampling is usually framed as sampling from a candidate set $$p(x)=\frac{\displaystyle\sum_{j=0}^N \frac{w_j}{x-x_j}\,y_j}{\displaystyle\sum_{j=0}^N \frac{w_j}{x-x_j}},$$8 of size $$p(x)=\frac{\displaystyle\sum_{j=0}^N \frac{w_j}{x-x_j}\,y_j}{\displaystyle\sum_{j=0}^N \frac{w_j}{x-x_j}},$$9 according to a normalized probability density function derived from a guiding metric $$w_0=w_N=\frac12(-1)^j,\qquad w_j=(-1)^j,\quad j=1,\dots,N-1.$$0. A generic form is

$$w_0=w_N=\frac12(-1)^j,\qquad w_j=(-1)^j,\quad j=1,\dots,N-1.$$1

with hyper-parameters $$w_0=w_N=\frac12(-1)^j,\qquad w_j=(-1)^j,\quad j=1,\dots,N-1.$$2 (Florido et al., 2024). The strategies summarized in the literature include fixed uniform random sampling, fixed Hammersley sampling, residual-based sampling with either random or Hammersley initialization, solution geometric curvature with $$w_0=w_N=\frac12(-1)^j,\qquad w_j=(-1)^j,\quad j=1,\dots,N-1.$$3, and residual curvature with $$w_0=w_N=\frac12(-1)^j,\qquad w_j=(-1)^j,\quad j=1,\dots,N-1.$$4 (Florido et al., 2024).

For Burgers’ equation,

$$w_0=w_N=\frac12(-1)^j,\qquad w_j=(-1)^j,\quad j=1,\dots,N-1.$$5

the residual is

$$w_0=w_N=\frac12(-1)^j,\qquad w_j=(-1)^j,\quad j=1,\dots,N-1.$$6

The relevant guiding metrics are then the residual magnitude $$w_0=w_N=\frac12(-1)^j,\qquad w_j=(-1)^j,\quad j=1,\dots,N-1.$$7, the mixed derivative of the solution $$w_0=w_N=\frac12(-1)^j,\qquad w_j=(-1)^j,\quad j=1,\dots,N-1.$$8, and the mixed derivative of the residual $$w_0=w_N=\frac12(-1)^j,\qquad w_j=(-1)^j,\quad j=1,\dots,N-1.$$9, all computed by automatic differentiation (Florido et al., 2024). A typical cycle samples $\xi_i$0 candidate points uniformly in the domain, evaluates the guiding metric, constructs the normalized PDF, draws $\xi_i$1 new collocation points, replaces the previous set, and continues training for another $\xi_i$2 steps (Florido et al., 2024).

The empirical results in the comparative study are problem-dependent but structured. For Burgers’ with $\xi_i$3, $\xi_i$4, and 30 resamples, the averaged $\xi_i$5 errors over 20 runs were reported as $\xi_i$6 for fixed random, $\xi_i$7 for fixed Hammersley, $\xi_i$8 for PDE,R, $\xi_i$9 for PDE,H, $R$0 for $R$1, and $R$2 for $R$3 (Florido et al., 2024). For Allen–Cahn with $R$4, $R$5, and 50 resamples, the corresponding values were $R$6, $R$7, $R$8, $R$9, $R_{i,j}=\ell_j(\xi_i)$0, and $R_{i,j}=\ell_j(\xi_i)$1 (Florido et al., 2024). The same source states that all adaptive methods outperform fixed for small $R_{i,j}=\ell_j(\xi_i)$2, that $R_{i,j}=\ell_j(\xi_i)$3 often yields $R_{i,j}=\ell_j(\xi_i)$4–$R_{i,j}=\ell_j(\xi_i)$5 lower error than residual-only sampling at fixed $R_{i,j}=\ell_j(\xi_i)$6, and that there exists a problem-dependent threshold $R_{i,j}=\ell_j(\xi_i)$7 above which fixed Hammersley sampling suffices (Florido et al., 2024).

A central caution in this literature is explicit: no single approach is likely to be “optimal” (Florido et al., 2024). This is not merely a rhetorical qualification. Residual-based sampling directly targets the largest local contributions to the PINN loss, curvature-based metrics focus on shocks or interfaces, and fixed low-discrepancy sets preserve global coverage (Florido et al., 2024). The trade-off is therefore intrinsic to the choice of guiding information.

5. Hessian-, gradient-, NTK-, and DEIM-based variants

Several later methods refine the choice of guiding information. A Hessian-based method defines the interior residual $R_{i,j}=\ell_j(\xi_i)$8, the pointwise loss integrand $R_{i,j}=\ell_j(\xi_i)$9, and the Hessian matrix $[-1,1]$00. It then uses a local indicator such as $[-1,1]$01 or $[-1,1]$02 and samples from a density

$[-1,1]$03

optionally shifted by a small constant $[-1,1]$04 for exploration (Caradot et al., 1 Apr 2025). The motivation comes from a refined quadrature argument. In one dimension, if $[-1,1]$05 is partitioned into cells $[-1,1]$06 of length $[-1,1]$07, with

$[-1,1]$08

then Theorem 3.1 gives

$[-1,1]$09

and Theorem 3.2 states that the refined Hessian-based allocation yields a tighter worst-case error bound when $[-1,1]$10 is nonuniform (Caradot et al., 1 Apr 2025). Comparative experiments summarized in the paper report, for example, that Hessian-RAD converges in $[-1,1]$11 epochs to a ten-times smaller error on 1D Newton’s law of cooling, whereas the uniform PINN converges in $[-1,1]$12 epochs to $[-1,1]$13-error $[-1,1]$14 (Caradot et al., 1 Apr 2025).

PACMANN formulates adaptive placement differently. It defines the pointwise residual $[-1,1]$15 and squared residual $[-1,1]$16, then moves the collocation points by gradient-based optimization guided by $[-1,1]$17 (Visser et al., 2024). Training is cast as a saddle-point problem in which the network parameters are minimized while the collocation set is chosen to maximize the residual term over the domain (Visser et al., 2024). The method updates points every $[-1,1]$18 training iterations, performs $[-1,1]$19 ascent steps using gradient ascent or Adam, and reinitializes points that leave the domain by uniform random draws (Visser et al., 2024). The paper states that the method matches the performance of state-of-the-art methods in low-dimensional problems and outperforms available approaches for high-dimensional problems, with the best performance observed for Adam (Visser et al., 2024). It also states a limitation: there is no formal guarantee of global optimality of the collocation-set placement, and the process may settle in local residual maxima (Visser et al., 2024).

PINNACLE replaces pointwise residual heuristics by an NTK-based score. For a candidate point $[-1,1]$20, the “convergence-degree” score is

$[-1,1]$21

where $[-1,1]$22 is the empirical NTK (Lau et al., 2024). The method forms a candidate pool over PDE and boundary point types, computes the eNTK via Nystrom approximation, selects collocation points either by sampling proportional to $[-1,1]$23 or by a K-means++ procedure in an NTK-derived embedding, and updates the fractions of PDE, boundary, and experimental points automatically (Lau et al., 2024). The theoretical result given in the summary relates the convergence degree of a training set to the PINN generalization error through the bound

$[-1,1]$24

where $[-1,1]$25 is defined from NTK eigenpairs and residual components (Lau et al., 2024).

QR-DEIM-based resampling uses residual dynamics over time rather than an instantaneous metric. Over a period $[-1,1]$26, residual snapshots on a fixed set $[-1,1]$27 are assembled into

$[-1,1]$28

with thin SVD

$[-1,1]$29

A reduced basis $[-1,1]$30 is selected by thresholding the singular values, and a column-pivoted QR factorization of $[-1,1]$31 identifies the DEIM interpolation points (Celaya et al., 13 Jan 2025). The resulting approximation

$[-1,1]$32

selects the most informative residual locations in the reduced basis sense (Celaya et al., 13 Jan 2025). The paper reports benchmark improvements on the wave, Allen–Cahn, and Burgers’ equations, including relative $[-1,1]$33 errors of $[-1,1]$34 for QR-DEIM on the wave equation, compared with $[-1,1]$35 for uniform random sampling, and $[-1,1]$36 for Burgers’, compared with $[-1,1]$37 for uniform random sampling (Celaya et al., 13 Jan 2025). The same source states that no PINN-specific convergence proof is given (Celaya et al., 13 Jan 2025).

6. Other collocation settings: rectangular collocation and oversampling

Outside the PINN literature, collocation-point resampling also appears as point-set compression or oversampling. In rectangular collocation for the Kohn–Sham equation,

$[-1,1]$38

one expands the wavefunction in a basis and enforces the equation at $[-1,1]$39 collocation points, typically with $[-1,1]$40 (Ku et al., 2019). The matrix form is

$[-1,1]$41

and after left-multiplication by $[-1,1]$42 one obtains the square generalized eigenproblem

$[-1,1]$43

The collocation set is first filtered from a large cube grid using a potential-based probability

$[-1,1]$44

then optimized by a Gaussian-process-regression fitness based on the RMSE of the predicted effective potential on a fixed test grid, with a genetic algorithm performing selection, uniform crossover, and pointwise mutation at probability $[-1,1]$45 (Ku et al., 2019). The reported outcome is a reduction of the collocation point set size by more than an order of magnitude, from about 51,000 points to 2,000 points, while retaining mHartree accuracy (Ku et al., 2019).

Oversampled boundary element collocation uses a different resampling logic. One chooses

$[-1,1]$46

forms the $[-1,1]$47 collocation matrix

$[-1,1]$48

and solves the overdetermined system through the normal equations

$[-1,1]$49

The summary states that oversampling at least lowers the error by a constant factor in most cases of practical interest, that oversampled collocation is much less affected by a poor choice of collocation points, and that in some cases convergence rates can be higher than the corresponding Galerkin method, although this requires oversampling at a faster than linear rate (Maierhofer et al., 2021). For equispaced points and Fourier-multiplier operators, Theorem 3.8 gives

$[-1,1]$50

from which the source concludes that for any fixed $[-1,1]$51 and $[-1,1]$52, the error constant decays like $[-1,1]$53 (Maierhofer et al., 2021).

These examples show that resampling is not restricted to online adaptation during optimization. It can also mean selecting a smaller but better collocation set before solving the discretized problem, or enforcing the same trial space at a larger set of collocation points to stabilize the method.

7. Practical issues, misconceptions, and open directions

Several practical recommendations recur across the literature. For adaptive PINN resampling, candidate-set sizes on the order of $[-1,1]$54–$[-1,1]$55 or $[-1,1]$56–$[-1,1]$57 per cycle are reported, with resampling every $[-1,1]$58–$[-1,1]$59 optimizer steps or every fixed number of epochs, and defaults such as $[-1,1]$60, $[-1,1]$61, or comparable small exploration constants (Florido et al., 2024, Caradot et al., 1 Apr 2025). Hammersley initialization is repeatedly used as a strong fixed baseline (Florido et al., 2024). PACMANN states that the extra cost remains a small constant factor when $[-1,1]$62, and its scaling is linear in both the number of residual points $[-1,1]$63 and the ambient dimension $[-1,1]$64 (Visser et al., 2024). QR-DEIM reports one SVD-plus-QR update every $[-1,1]$65 iterations and describes the overall overhead as approximately $[-1,1]$66–$[-1,1]$67 of training time if $[-1,1]$68 (Celaya et al., 13 Jan 2025).

Several misconceptions are explicitly contradicted by the cited work. First, adaptive resampling is not universally superior: fixed Hammersley sampling can suffice once the collocation budget exceeds a problem-dependent threshold $[-1,1]$69 (Florido et al., 2024). Second, residual magnitude is not the only meaningful guide. Curvature-based and Hessian-based indicators may better resolve shocks, interfaces, or steep regions, while NTK-based methods target training dynamics across point types rather than only instantaneous residuals (Florido et al., 2024, Caradot et al., 1 Apr 2025, Lau et al., 2024). Third, adaptive methods do not all scale the same way. The PACMANN abstract states that previous adaptive collocation point sampling methods have been reported to scale poorly to higher dimensions, which motivates its gradient-based point-movement formulation (Visser et al., 2024). Fourth, stronger empirical performance does not necessarily come with a formal PDE-specific guarantee: QR-DEIM explicitly lacks a PINN-specific convergence proof, and PACMANN explicitly lacks a formal guarantee of global optimality for the point placement (Celaya et al., 13 Jan 2025, Visser et al., 2024).

A plausible implication of the collected literature is that collocation-point resampling is increasingly being treated as a numerical design problem in its own right. In spectral DDE/FDE solvers, the design variable is the evaluation operator that maps global polynomial data to transformed arguments (Hale, 2024). In PINNs, it is the online distribution of residual evaluations under a fixed computational budget (Florido et al., 2024, Caradot et al., 1 Apr 2025, Visser et al., 2024, Lau et al., 2024, Celaya et al., 13 Jan 2025). In rectangular collocation and boundary element methods, it is the geometry and cardinality of the collocation set relative to the trial space (Ku et al., 2019, Maierhofer et al., 2021). What remains common across these settings is the premise that the accuracy, robustness, and efficiency of collocation methods depend not only on the approximation space and optimizer, but also on how the collocation points are transformed, selected, or oversampled.