Physics-Informed Broad Learning System

• PIBLS is a physics-informed PDE solver that replaces deep backpropagation with a broad, shallow learning system and a direct least-squares approach.
• It reformulates PDE enforcement into a single-step linear algebra problem, offering rapid convergence and machine-precision accuracy.
• PIBLS maintains mesh-free modeling while ensuring theoretical soundness through a universal approximation guarantee in Sobolev spaces.

Searching arXiv for the specified paper and closely related work to ground the article. arxiv_search(query="Physics-Informed Broad Learning System partial differential equations", max_results=5) arxiv_search(query="(Yu et al., 18 Jun 2026)", max_results=10) Physics-Informed Broad Learning System (PIBLS) is a physics-informed, backpropagation-free solver for partial differential equations (PDEs) built on the Broad Learning System (BLS) architecture. It preserves the physics-informed idea of Physics-Informed Neural Networks (PINNs)—enforcing the PDE, boundary conditions, and initial conditions through residuals—but replaces the usual deep network and iterative gradient descent with a single-step least-squares solve for the output weights. In the reported formulation, PIBLS uses a shallow-but-wide BLS, keeps the internal random basis fixed, trains only the final linear output weights, and solves those weights by least squares rather than gradient-based optimization. It is presented as mesh-free like PINNs but not iterative, with the stated aim of combining the physics enforcement of PINNs, the speed of closed-form/randomized solvers, and stronger expressive power through the broad basis of BLS (Yu et al., 18 Jun 2026).

1. Conceptual position among PDE solvers

PIBLS is introduced against two established baselines for PDE computation: traditional mesh-based numerical solvers and PINNs. Classical solvers such as FEM, FVM, and FDM require mesh generation or discretization, incur high cost for fine resolution, and are typically PDE-specific and mesh-dependent. PINNs offer a mesh-free alternative by approximating the PDE solution $u(\boldsymbol{x})$ with a deep neural network $\hat{u}(\boldsymbol{x};\theta)$ and training $\theta$ by minimizing a physics-informed loss using backpropagation, often with Adam and L-BFGS. PIBLS retains the residual-based enforcement of governing equations and constraints, but changes the optimization mechanism fundamentally (Yu et al., 18 Jun 2026).

Relative to PINNs, PIBLS uses a shallow-but-wide BLS instead of a deep network, keeps the internal random basis fixed, trains only the final linear output weights, and solves those weights by least squares rather than gradient-based optimization. The reported motivation is to eliminate PINN issues identified in the formulation: slow convergence, sensitivity to optimization hyperparameters, pathological loss landscapes, and gradient vanishing or imbalance. The conceptual novelty is summarized as follows: PIBLS turns PDE solving into a direct linear algebra problem over a rich random basis, with physics enforced through residual equations.

This positioning suggests that PIBLS occupies an intermediate methodological space. It remains physics-informed and mesh-free, as in PINNs, but its computational core is closer to randomized linear-algebraic solvers than to deep-network training.

2. Physics-informed direct least-squares formulation

The reported framework considers a general PDE with differential, boundary, and initial operators: $$\begin{aligned} \mathcal{N}[u(\boldsymbol{x})] &= f(\boldsymbol{x}), \quad \boldsymbol{x} \in \Omega,\ \mathcal{B}[u(\boldsymbol{x})] &= g(\boldsymbol{x}), \quad \boldsymbol{x} \in \partial \Omega_B,\ \mathcal{I}[u(\boldsymbol{x})] &= s(\boldsymbol{x}), \quad \boldsymbol{x} \in \partial \Omega_I. \end{aligned}$$

PIBLS predicts

$$\hat{u}(\boldsymbol{x}) = \mathbf{A}(\boldsymbol{x}) \mathbf{W}_{\text{out}},$$

where $\mathbf{A}$ is the BLS basis matrix and $\mathbf{W}_{\text{out}}$ are the only trainable parameters. At collocation points, the residuals are defined by

$$\begin{aligned} r_{\mathcal{N}}(\boldsymbol{x}_i) &= \mathcal{N}[\hat{u}(\boldsymbol{x}_i)] - f(\boldsymbol{x}_i),\ r_{\mathcal{B}}(\boldsymbol{x}_j) &= \mathcal{B}[\hat{u}(\boldsymbol{x}_j)] - g(\boldsymbol{x}_j),\ r_{\mathcal{I}}(\boldsymbol{x}_k) &= \mathcal{I}[\hat{u}(\boldsymbol{x}_k)] - s(\boldsymbol{x}_k). \end{aligned}$$

For linear PDEs, the key identity is

$$\mathcal{O}[\hat{u}(\boldsymbol{x})] = \mathcal{O}[\mathbf{A}\mathbf{W}_{\text{out}}] = \mathcal{O}[\mathbf{A}]\,\mathbf{W}_{\text{out}}, \qquad \mathcal{O}\in\{\mathcal{N},\mathcal{B},\mathcal{I}\}.$$

All conditions are then stacked into a single linear system,

$$\mathbf{A}_{\text{sys}}\mathbf{W}_{\text{out}} = \mathbf{T},$$

with

$\hat{u}(\boldsymbol{x};\theta)$0

The PDE solve therefore becomes the least-squares problem

$\hat{u}(\boldsymbol{x};\theta)$1

This formulation is the core “physics-informed direct solve” idea. Instead of iteratively minimizing a composite loss over all network parameters, PIBLS reduces the physics-constrained approximation problem to direct least-squares optimization in the output layer.

3. Broad Learning System architecture

The BLS architecture used in PIBLS consists of feature nodes, enhancement nodes, and a final linear readout. For input $\hat{u}(\boldsymbol{x};\theta)$2, the feature layer maps the input into $\hat{u}(\boldsymbol{x};\theta)$3 groups of feature nodes, each with $\hat{u}(\boldsymbol{x};\theta)$4 nodes: $\hat{u}(\boldsymbol{x};\theta)$5 where $\hat{u}(\boldsymbol{x};\theta)$6, $\hat{u}(\boldsymbol{x};\theta)$7, $\hat{u}(\boldsymbol{x};\theta)$8, and the parameters are randomly initialized and fixed. The feature groups are concatenated as

$\hat{u}(\boldsymbol{x};\theta)$9

A second random nonlinear mapping produces the enhancement layer: $\theta$0 with $\theta$1, $\theta$2, and $\theta$3. The enhancement groups are concatenated as

$\theta$4

The full basis is

$\theta$5

and the surrogate solution is

$\theta$6

Because $\theta$7 is linear in $\theta$8, once $\theta$9 is constructed the solution reduces to solving a linear or nonlinear least-squares system; no gradient descent over hidden weights is required. The formulation also notes that feature weights and biases may optionally be fine-tuned with a sparse autoencoder, but the core PIBLS method still avoids backpropagation through the PDE loss (Yu et al., 18 Jun 2026).

This architecture is broad rather than deep: representational capacity is increased by width, grouped random mappings, and the composite basis $$\begin{aligned} \mathcal{N}[u(\boldsymbol{x})] &= f(\boldsymbol{x}), \quad \boldsymbol{x} \in \Omega,\ \mathcal{B}[u(\boldsymbol{x})] &= g(\boldsymbol{x}), \quad \boldsymbol{x} \in \partial \Omega_B,\ \mathcal{I}[u(\boldsymbol{x})] &= s(\boldsymbol{x}), \quad \boldsymbol{x} \in \partial \Omega_I. \end{aligned}$$0, rather than by repeated trainable depth.

4. Nonlinear PDE handling and the enhanced NLSQ-perturb algorithm

For nonlinear PDEs, the linear-system reduction no longer applies directly because $$\begin{aligned} \mathcal{N}[u(\boldsymbol{x})] &= f(\boldsymbol{x}), \quad \boldsymbol{x} \in \Omega,\ \mathcal{B}[u(\boldsymbol{x})] &= g(\boldsymbol{x}), \quad \boldsymbol{x} \in \partial \Omega_B,\ \mathcal{I}[u(\boldsymbol{x})] &= s(\boldsymbol{x}), \quad \boldsymbol{x} \in \partial \Omega_I. \end{aligned}$$1 is nonlinear in the solution field. The method therefore formulates a nonlinear least-squares objective,

$$\begin{aligned} \mathcal{N}[u(\boldsymbol{x})] &= f(\boldsymbol{x}), \quad \boldsymbol{x} \in \Omega,\ \mathcal{B}[u(\boldsymbol{x})] &= g(\boldsymbol{x}), \quad \boldsymbol{x} \in \partial \Omega_B,\ \mathcal{I}[u(\boldsymbol{x})] &= s(\boldsymbol{x}), \quad \boldsymbol{x} \in \partial \Omega_I. \end{aligned}$$2

with the initial term omitted for steady-state problems. The reported solver combines Trust-Region Reflective (TRF) nonlinear least squares, exact Jacobians from analytical derivatives, and an Enhanced NLSQ-perturb algorithm (Yu et al., 18 Jun 2026).

Because the surrogate remains linear in the output weights, the residual Jacobian has the structure

$$\begin{aligned} \mathcal{N}[u(\boldsymbol{x})] &= f(\boldsymbol{x}), \quad \boldsymbol{x} \in \Omega,\ \mathcal{B}[u(\boldsymbol{x})] &= g(\boldsymbol{x}), \quad \boldsymbol{x} \in \partial \Omega_B,\ \mathcal{I}[u(\boldsymbol{x})] &= s(\boldsymbol{x}), \quad \boldsymbol{x} \in \partial \Omega_I. \end{aligned}$$3

This preserves analytical tractability even when the PDE residual is nonlinear in $$\begin{aligned} \mathcal{N}[u(\boldsymbol{x})] &= f(\boldsymbol{x}), \quad \boldsymbol{x} \in \Omega,\ \mathcal{B}[u(\boldsymbol{x})] &= g(\boldsymbol{x}), \quad \boldsymbol{x} \in \partial \Omega_B,\ \mathcal{I}[u(\boldsymbol{x})] &= s(\boldsymbol{x}), \quad \boldsymbol{x} \in \partial \Omega_I. \end{aligned}$$4.

The perturbation modification proceeds in five reported steps. First, a physics-informed initial guess $$\begin{aligned} \mathcal{N}[u(\boldsymbol{x})] &= f(\boldsymbol{x}), \quad \boldsymbol{x} \in \Omega,\ \mathcal{B}[u(\boldsymbol{x})] &= g(\boldsymbol{x}), \quad \boldsymbol{x} \in \partial \Omega_B,\ \mathcal{I}[u(\boldsymbol{x})] &= s(\boldsymbol{x}), \quad \boldsymbol{x} \in \partial \Omega_I. \end{aligned}$$5 is constructed by solving a linearized PDE. Second, the nonlinear least-squares problem is solved with TRF. Third, if the loss stagnates or local minima occur, the current weights are perturbed: $$\begin{aligned} \mathcal{N}[u(\boldsymbol{x})] &= f(\boldsymbol{x}), \quad \boldsymbol{x} \in \Omega,\ \mathcal{B}[u(\boldsymbol{x})] &= g(\boldsymbol{x}), \quad \boldsymbol{x} \in \partial \Omega_B,\ \mathcal{I}[u(\boldsymbol{x})] &= s(\boldsymbol{x}), \quad \boldsymbol{x} \in \partial \Omega_I. \end{aligned}$$6 Fourth, TRF is re-run from the perturbed guess. Fifth, the best solution among trials is retained.

The stated purpose of this perturbation loop is to reduce sensitivity to local minima and improve robustness on difficult nonlinear PDEs. A plausible implication is that the nonlinear PIBLS variant remains optimization-based, but in a lower-dimensional and more structured form than gradient-based training of all hidden and output parameters.

5. Universal approximation in Sobolev spaces

PIBLS is accompanied by a universal approximation theorem in Sobolev spaces. The assumptions are explicit: $$\begin{aligned} \mathcal{N}[u(\boldsymbol{x})] &= f(\boldsymbol{x}), \quad \boldsymbol{x} \in \Omega,\ \mathcal{B}[u(\boldsymbol{x})] &= g(\boldsymbol{x}), \quad \boldsymbol{x} \in \partial \Omega_B,\ \mathcal{I}[u(\boldsymbol{x})] &= s(\boldsymbol{x}), \quad \boldsymbol{x} \in \partial \Omega_I. \end{aligned}$$7 and $$\begin{aligned} \mathcal{N}[u(\boldsymbol{x})] &= f(\boldsymbol{x}), \quad \boldsymbol{x} \in \Omega,\ \mathcal{B}[u(\boldsymbol{x})] &= g(\boldsymbol{x}), \quad \boldsymbol{x} \in \partial \Omega_B,\ \mathcal{I}[u(\boldsymbol{x})] &= s(\boldsymbol{x}), \quad \boldsymbol{x} \in \partial \Omega_I. \end{aligned}$$8 are non-constant, bounded, and $$\begin{aligned} \mathcal{N}[u(\boldsymbol{x})] &= f(\boldsymbol{x}), \quad \boldsymbol{x} \in \Omega,\ \mathcal{B}[u(\boldsymbol{x})] &= g(\boldsymbol{x}), \quad \boldsymbol{x} \in \partial \Omega_B,\ \mathcal{I}[u(\boldsymbol{x})] &= s(\boldsymbol{x}), \quad \boldsymbol{x} \in \partial \Omega_I. \end{aligned}$$9; and all feature and enhancement parameters are i.i.d. sampled from continuous distributions with compact or bounded support. The target function satisfies

$$\hat{u}(\boldsymbol{x}) = \mathbf{A}(\boldsymbol{x}) \mathbf{W}_{\text{out}},$$0

where $$\hat{u}(\boldsymbol{x}) = \mathbf{A}(\boldsymbol{x}) \mathbf{W}_{\text{out}},$$1 is compact with Lipschitz boundary and

$$\hat{u}(\boldsymbol{x}) = \mathbf{A}(\boldsymbol{x}) \mathbf{W}_{\text{out}},$$2

The theorem states that for any compact domain with Lipschitz boundary $$\hat{u}(\boldsymbol{x}) = \mathbf{A}(\boldsymbol{x}) \mathbf{W}_{\text{out}},$$3 and any continuous function $$\hat{u}(\boldsymbol{x}) = \mathbf{A}(\boldsymbol{x}) \mathbf{W}_{\text{out}},$$4, $$\hat{u}(\boldsymbol{x}) = \mathbf{A}(\boldsymbol{x}) \mathbf{W}_{\text{out}},$$5, there exists a sequence of PIBLS approximants $$\hat{u}(\boldsymbol{x}) = \mathbf{A}(\boldsymbol{x}) \mathbf{W}_{\text{out}},$$6 such that

$$\hat{u}(\boldsymbol{x}) = \mathbf{A}(\boldsymbol{x}) \mathbf{W}_{\text{out}},$$7

Equivalently, for any $$\hat{u}(\boldsymbol{x}) = \mathbf{A}(\boldsymbol{x}) \mathbf{W}_{\text{out}},$$8,

$$\hat{u}(\boldsymbol{x}) = \mathbf{A}(\boldsymbol{x}) \mathbf{W}_{\text{out}},$$9

for sufficiently large $\mathbf{A}$0 and $\mathbf{A}$1 (Yu et al., 18 Jun 2026).

The proof proceeds in two stages. In the first, the random feature map $\mathbf{A}$2 is shown, for sufficiently large $\mathbf{A}$3, to be almost surely a smooth embedding, described as a diffeomorphism onto its image. Immersion is established through full-rank Jacobian arguments using

$\mathbf{A}$4

and injectivity is analyzed via the kernel

$\mathbf{A}$5

In the second stage, the residual

$\mathbf{A}$6

is approximated by a smooth function, pulled back through the feature-layer diffeomorphism, and then approximated using the enhancement layer. The argument uses equivalence of Sobolev norms under smooth coordinate changes: $\mathbf{A}$7

The stated interpretation is that PIBLS is not only numerically effective but also theoretically expressive enough to approximate sufficiently smooth PDE solutions arbitrarily well, provided the number of feature and enhancement nodes grows.

6. Benchmarks, baselines, and reported performance

The reported evaluation covers 11 benchmark PDEs: 8 linear and 3 nonlinear. The tests include 1D, 2D, and time-dependent problems. The linear benchmarks TC-1 to TC-8 include advection, diffusion, advection-diffusion, 2D advection, two 2D Poisson problems, constant-coefficient wave/advection in space-time, and variable-coefficient wave/advection in space-time. The nonlinear benchmarks TC-9 to TC-11 are nonlinear Helmholtz, nonlinear spring oscillator, and viscous Burgers’ equation. Baselines are PINN #1, a parameter-matched shallow PINN trained with Adam for 5000 iterations; PINN #2, a deeper PINN with 10 layers and 50 nodes trained with Adam + L-BFGS; PIELM for linear cases; locELM for nonlinear cases; and FEM with high-resolution finite elements. Experiments ran in Python on an Intel Core i5-12600KF CPU with 32GB RAM in a WSL environment (Yu et al., 18 Jun 2026).

The reported metrics are max error, $\mathbf{A}$8 error, and training time. For the linear cases, PIBLS attains the following $\mathbf{A}$9 errors and times: TC-1, $\mathbf{W}_{\text{out}}$0 in $\mathbf{W}_{\text{out}}$1 s; TC-2, $\mathbf{W}_{\text{out}}$2 in $\mathbf{W}_{\text{out}}$3 s; TC-3, $\mathbf{W}_{\text{out}}$4 in $\mathbf{W}_{\text{out}}$5 s; TC-4, $\mathbf{W}_{\text{out}}$6 in $\mathbf{W}_{\text{out}}$7 s; TC-5, $\mathbf{W}_{\text{out}}$8 in $\mathbf{W}_{\text{out}}$9 s; TC-6, $$\begin{aligned} r_{\mathcal{N}}(\boldsymbol{x}_i) &= \mathcal{N}[\hat{u}(\boldsymbol{x}_i)] - f(\boldsymbol{x}_i),\ r_{\mathcal{B}}(\boldsymbol{x}_j) &= \mathcal{B}[\hat{u}(\boldsymbol{x}_j)] - g(\boldsymbol{x}_j),\ r_{\mathcal{I}}(\boldsymbol{x}_k) &= \mathcal{I}[\hat{u}(\boldsymbol{x}_k)] - s(\boldsymbol{x}_k). \end{aligned}$$0 in $$\begin{aligned} r_{\mathcal{N}}(\boldsymbol{x}_i) &= \mathcal{N}[\hat{u}(\boldsymbol{x}_i)] - f(\boldsymbol{x}_i),\ r_{\mathcal{B}}(\boldsymbol{x}_j) &= \mathcal{B}[\hat{u}(\boldsymbol{x}_j)] - g(\boldsymbol{x}_j),\ r_{\mathcal{I}}(\boldsymbol{x}_k) &= \mathcal{I}[\hat{u}(\boldsymbol{x}_k)] - s(\boldsymbol{x}_k). \end{aligned}$$1 s; TC-7, $$\begin{aligned} r_{\mathcal{N}}(\boldsymbol{x}_i) &= \mathcal{N}[\hat{u}(\boldsymbol{x}_i)] - f(\boldsymbol{x}_i),\ r_{\mathcal{B}}(\boldsymbol{x}_j) &= \mathcal{B}[\hat{u}(\boldsymbol{x}_j)] - g(\boldsymbol{x}_j),\ r_{\mathcal{I}}(\boldsymbol{x}_k) &= \mathcal{I}[\hat{u}(\boldsymbol{x}_k)] - s(\boldsymbol{x}_k). \end{aligned}$$2 in $$\begin{aligned} r_{\mathcal{N}}(\boldsymbol{x}_i) &= \mathcal{N}[\hat{u}(\boldsymbol{x}_i)] - f(\boldsymbol{x}_i),\ r_{\mathcal{B}}(\boldsymbol{x}_j) &= \mathcal{B}[\hat{u}(\boldsymbol{x}_j)] - g(\boldsymbol{x}_j),\ r_{\mathcal{I}}(\boldsymbol{x}_k) &= \mathcal{I}[\hat{u}(\boldsymbol{x}_k)] - s(\boldsymbol{x}_k). \end{aligned}$$3 s; and TC-8, $$\begin{aligned} r_{\mathcal{N}}(\boldsymbol{x}_i) &= \mathcal{N}[\hat{u}(\boldsymbol{x}_i)] - f(\boldsymbol{x}_i),\ r_{\mathcal{B}}(\boldsymbol{x}_j) &= \mathcal{B}[\hat{u}(\boldsymbol{x}_j)] - g(\boldsymbol{x}_j),\ r_{\mathcal{I}}(\boldsymbol{x}_k) &= \mathcal{I}[\hat{u}(\boldsymbol{x}_k)] - s(\boldsymbol{x}_k). \end{aligned}$$4 in $$\begin{aligned} r_{\mathcal{N}}(\boldsymbol{x}_i) &= \mathcal{N}[\hat{u}(\boldsymbol{x}_i)] - f(\boldsymbol{x}_i),\ r_{\mathcal{B}}(\boldsymbol{x}_j) &= \mathcal{B}[\hat{u}(\boldsymbol{x}_j)] - g(\boldsymbol{x}_j),\ r_{\mathcal{I}}(\boldsymbol{x}_k) &= \mathcal{I}[\hat{u}(\boldsymbol{x}_k)] - s(\boldsymbol{x}_k). \end{aligned}$$5 s. In the reported table, PIBLS is consistently many orders of magnitude more accurate than PINNs and usually better than PIELM and even FEM.

For the nonlinear cases, PIBLS reports TC-9 with $$\begin{aligned} r_{\mathcal{N}}(\boldsymbol{x}_i) &= \mathcal{N}[\hat{u}(\boldsymbol{x}_i)] - f(\boldsymbol{x}_i),\ r_{\mathcal{B}}(\boldsymbol{x}_j) &= \mathcal{B}[\hat{u}(\boldsymbol{x}_j)] - g(\boldsymbol{x}_j),\ r_{\mathcal{I}}(\boldsymbol{x}_k) &= \mathcal{I}[\hat{u}(\boldsymbol{x}_k)] - s(\boldsymbol{x}_k). \end{aligned}$$6 in $$\begin{aligned} r_{\mathcal{N}}(\boldsymbol{x}_i) &= \mathcal{N}[\hat{u}(\boldsymbol{x}_i)] - f(\boldsymbol{x}_i),\ r_{\mathcal{B}}(\boldsymbol{x}_j) &= \mathcal{B}[\hat{u}(\boldsymbol{x}_j)] - g(\boldsymbol{x}_j),\ r_{\mathcal{I}}(\boldsymbol{x}_k) &= \mathcal{I}[\hat{u}(\boldsymbol{x}_k)] - s(\boldsymbol{x}_k). \end{aligned}$$7 s, TC-10 with $$\begin{aligned} r_{\mathcal{N}}(\boldsymbol{x}_i) &= \mathcal{N}[\hat{u}(\boldsymbol{x}_i)] - f(\boldsymbol{x}_i),\ r_{\mathcal{B}}(\boldsymbol{x}_j) &= \mathcal{B}[\hat{u}(\boldsymbol{x}_j)] - g(\boldsymbol{x}_j),\ r_{\mathcal{I}}(\boldsymbol{x}_k) &= \mathcal{I}[\hat{u}(\boldsymbol{x}_k)] - s(\boldsymbol{x}_k). \end{aligned}$$8 in $$\begin{aligned} r_{\mathcal{N}}(\boldsymbol{x}_i) &= \mathcal{N}[\hat{u}(\boldsymbol{x}_i)] - f(\boldsymbol{x}_i),\ r_{\mathcal{B}}(\boldsymbol{x}_j) &= \mathcal{B}[\hat{u}(\boldsymbol{x}_j)] - g(\boldsymbol{x}_j),\ r_{\mathcal{I}}(\boldsymbol{x}_k) &= \mathcal{I}[\hat{u}(\boldsymbol{x}_k)] - s(\boldsymbol{x}_k). \end{aligned}$$9 s, and TC-11 with $$\mathcal{O}[\hat{u}(\boldsymbol{x})] = \mathcal{O}[\mathbf{A}\mathbf{W}_{\text{out}}] = \mathcal{O}[\mathbf{A}]\,\mathbf{W}_{\text{out}}, \qquad \mathcal{O}\in\{\mathcal{N},\mathcal{B},\mathcal{I}\}.$$0 in $$\mathcal{O}[\hat{u}(\boldsymbol{x})] = \mathcal{O}[\mathbf{A}\mathbf{W}_{\text{out}}] = \mathcal{O}[\mathbf{A}]\,\mathbf{W}_{\text{out}}, \qquad \mathcal{O}\in\{\mathcal{N},\mathcal{B},\mathcal{I}\}.$$1 s. The comparative description states that PINNs often fail badly on nonlinear problems, especially TC-9 and TC-11; locELM performs well, but PIBLS is often more accurate and sometimes much faster; and for TC-11, PIBLS is dramatically better than PINNs and about $$\mathcal{O}[\hat{u}(\boldsymbol{x})] = \mathcal{O}[\mathbf{A}\mathbf{W}_{\text{out}}] = \mathcal{O}[\mathbf{A}]\,\mathbf{W}_{\text{out}}, \qquad \mathcal{O}\in\{\mathcal{N},\mathcal{B},\mathcal{I}\}.$$2 faster than locELM.

The speed comparisons are also explicit: PINN #1 typically takes tens of seconds, PINN #2 takes hundreds to thousands of seconds, PIBLS typically takes fractions of a second, and the claimed speedup is 1 to 3 orders of magnitude faster than PINNs. The reported main empirical conclusion is that PIBLS attains machine-precision or near-machine-precision accuracy on many tests, sub-second to few-second solve times, and substantial advantages over gradient-trained PINNs.

7. Practical implications and stated limitations

The practical implications identified for scientific machine learning are fivefold. First, PIBLS offers fast PDE solving by replacing long training loops with direct least-squares computation, making it suitable for near-real-time simulation and design workflows. Second, it provides physics-informed mesh-free modeling, requiring no mesh generation and operating on scattered collocation points like PINNs. Third, it is presented as having better numerical stability than PINNs because it avoids backpropagation through a deep loss landscape and thereby avoids many optimization pathologies. Fourth, it has strong theoretical backing through the universal approximation proof in Sobolev spaces. Fifth, it shows excellent benchmark performance, especially on smooth linear problems and also on selected nonlinear PDEs with the enhanced solver (Yu et al., 18 Jun 2026).

The formulation also states that PIBLS is not a universal cure-all. Nonlinear PDEs are more sensitive, with performance depending on the initialization range $$\mathcal{O}[\hat{u}(\boldsymbol{x})] = \mathcal{O}[\mathbf{A}\mathbf{W}_{\text{out}}] = \mathcal{O}[\mathbf{A}]\,\mathbf{W}_{\text{out}}, \qquad \mathcal{O}\in\{\mathcal{N},\mathcal{B},\mathcal{I}\}.$$3, and nonlinear cases may require careful tuning and perturbation strategies. Random basis quality matters because the hidden parameters are fixed random draws, so representation quality depends on the random feature space, and too small or too large initialization can hurt nonlinear performance. Wide networks may increase cost: error saturates beyond some width while compute time keeps increasing, creating a trade-off between width and efficiency. PIBLS is also described as not fully adaptive like advanced PDE solvers because it still relies on collocation sampling and fixed random features; for complex geometries, discontinuities, or highly stiff problems, performance may be less predictable. Finally, the proof is asymptotic: universal approximation holds as $$\mathcal{O}[\hat{u}(\boldsymbol{x})] = \mathcal{O}[\mathbf{A}\mathbf{W}_{\text{out}}] = \mathcal{O}[\mathbf{A}]\,\mathbf{W}_{\text{out}}, \qquad \mathcal{O}\in\{\mathcal{N},\mathcal{B},\mathcal{I}\}.$$4, but practical accuracy depends on finite-width choices.

Taken together, these points place PIBLS within the broader development of scientific machine learning as a physics-informed, broad, random-feature-based solver that reformulates PDE enforcement into a least-squares problem. Its defining characteristics are the replacement of deep backpropagation with a BLS basis plus direct linear or nonlinear least-squares solve, the reported ability to achieve machine-level accuracy on many benchmark PDEs, and a Sobolev-space universal approximation guarantee.