Core Abstractions (pypielm.core)

Core abstractions: BasePIELM, feature maps, and linear solvers.

Implemented across Steps 2 (core) and 5 (models).

BasePIELM

Abstract base class and shared configuration for all PIELM models.

All PIELM and PINN variants in PyPIELM inherit from BasePIELM, which is itself a torch.nn.Module. The ELM paradigm mandates that hidden-layer weights are randomly sampled and frozen; only the output layer weights beta are determined during BasePIELM.fit() (analytically, not via gradient descent).

class pypielm.core.base.PIELMConfig(hidden_dim=200, activation='tanh', ridge_lambda=1e-08, seed=42, device='cpu', dtype=<factory>)

Bases: object

Shared hyper-parameter configuration for PIELM models.

Parameters:

• hidden_dim (int) – Number of hidden neurons (random feature dimension H).

• activation (str) – Activation function name. Supported: 'tanh', 'sin', 'relu', 'sigmoid', 'softplus'.

• ridge_lambda (float) – Tikhonov regularisation coefficient λ used in the output-weight solve. Typically in [1e-12, 1e-4].

• seed (int) – Integer seed for reproducible hidden-layer weight sampling.

• device (str) – PyTorch device string (e.g. 'cpu', 'cuda', 'cuda:0').

• dtype (dtype) – Floating-point dtype. torch.float64 (double precision) is the default because PDE residuals require high numerical accuracy.

Example:

cfg = PIELMConfig(hidden_dim=500, activation="tanh", ridge_lambda=1e-10)

hidden_dim: int = 200

activation: str = 'tanh'

ridge_lambda: float = 1e-08

seed: int = 42

device: str = 'cpu'

dtype: dtype

class pypielm.core.base.BasePIELM(*args, **kwargs)

Bases: Module

Abstract base for all PIELM and PINN variants in PyPIELM.

Subclasses must implement fit(), predict(), score(), and get_feature_matrix(). Hidden-layer weights are frozen (ELM paradigm): they are initialised in __init__ and never updated by a gradient step. Output weights beta are solved analytically in fit().

The class inherits from torch.nn.Module to enable:

• model.to(device) — move all parameters to a device.

• torch.jit.script / torch.jit.trace — TorchScript export.

• torch.onnx.export — ONNX export.

• torch.compile — optional kernel fusion.

Scikit-learn compatibility:

• fit(dataset, ...) returns self (fluent chaining).

• predict(X) returns a torch.Tensor.

• score(X, y) returns a scalar float.

Initialize internal Module state, shared by both nn.Module and ScriptModule.

abstractmethod fit(dataset, *, pde_operator=None, bcs=None, ics=None, collocation_sampler=None)

Solve for output weights analytically.

Parameters:

• dataset (PIELMDataset) – A PIELMDataset containing X_train, y_train, and (optionally) validation splits.

• pde_operator (Any | None) – Optional differential operator applied at collocation points to build the physics residual block in the linear system. If None, the model falls back to pure data regression.

• bcs (list[Any] | None) – List of boundary condition objects (DirichletBC, etc.).

• ics (list[Any] | None) – List of initial condition objects.

• collocation_sampler (Any | None) – Overrides the default collocation sampler used to generate interior PDE points.

Returns:

model.fit(ds, pde_operator=op).predict(X_test).

Return type:

self — enables fluent chaining

abstractmethod predict(X)

Evaluate the surrogate solution at input coordinates X.

Parameters:

X (ndarray | Tensor) – Input coordinates of shape (N, d). Accepts both torch.Tensor and numpy.ndarray; the result is always a torch.Tensor.

Return type:

Tensor

Returns:

Predicted values of shape (N, 1) or (N, out_dim).

abstractmethod score(X, y, metric='relative_l2')

Compute a scalar accuracy metric on a held-out set.

Parameters:

• X (ndarray | Tensor) – Input coordinates of shape (N, d).

• y (ndarray | Tensor) – Reference values of shape (N,) or (N, out_dim).

• metric (str) – One of 'relative_l2', 'rmse', 'mae', 'r2', 'max_error'.

Return type:

float

Returns:

Scalar metric value. For error metrics (relative_l2, rmse, mae, max_error) lower is better; for R² higher is better.

abstractmethod get_feature_matrix(X)

Return the hidden-layer activation matrix Φ(X).

Parameters:

X (ndarray | Tensor) – Input coordinates of shape (N, d).

Return type:

Tensor

Returns:

Feature matrix H of shape (N, hidden_dim).

Feature Maps

Random and Fourier feature maps for ELM hidden layers (PyTorch).

All feature maps are torch.nn.Module subclasses with frozen parameters. They provide both a forward pass and analytic first/second partial derivatives (fast paths for tanh and sin/cos) as well as an autograd fallback for arbitrary activations.

class pypielm.core.feature_maps.RandomFeatureMap(input_dim, hidden_dim, activation='tanh', seed=42, w_scale=1.0, device='cpu', dtype=torch.float64)

Bases: Module

Random-weight ELM hidden layer with analytic partial derivatives.

The hidden representation is:

\[H_{ij} = \sigma\!\left(\sum_k x_{ik} W_{kj} + b_j\right)\]

where \(\mathbf{W} \in \mathbb{R}^{d \times H}\) and \(\mathbf{b} \in \mathbb{R}^H\) are sampled once at construction and never updated (frozen buffers).

Analytic first- and second-order partial derivatives are provided for 'tanh', 'sin', 'relu', 'sigmoid', and 'softplus'.

Parameters:

• input_dim (int) – Spatial dimension d of the input coordinates.

• hidden_dim (int) – Number of neurons H.

• activation (str) – Name of the activation function.

• seed (int) – Random seed for weight initialisation.

• w_scale (float) – Multiplicative scale applied to the random weights W.

• device (str | device) – Target device ('cpu', 'cuda', …).

• dtype (dtype) – Floating-point precision. torch.float64 recommended.

Example:

fm = RandomFeatureMap(input_dim=2, hidden_dim=200, seed=42)
H   = fm(X)          # shape (N, 200)
dH  = fm.d1(X, 0)    # ∂H/∂x₀,   shape (N, 200)
d2H = fm.d2(X, 0)    # ∂²H/∂x₀², shape (N, 200)

Initialize internal Module state, shared by both nn.Module and ScriptModule.

W: Tensor

b: Tensor

forward(X)

Compute hidden-layer activations H = σ(X @ W + b).

Parameters:

X (Tensor) – Input coordinates, shape (N, d).

Return type:

Tensor

Returns:

Feature matrix H, shape (N, H).

d1(X, axis)

Analytic first partial derivative ∂H/∂x_{axis}.

\[\frac{\partial H_{ij}}{\partial x_{i,\text{axis}}} = \sigma'(Z_{ij}) \cdot W_{\text{axis},j}\]

Parameters:

• X (Tensor) – Input coordinates, shape (N, d).

• axis (int) – Spatial dimension index (0-based).

Return type:

Tensor

Returns:

First-derivative matrix, shape (N, H).

d2(X, axis)

Analytic second partial derivative ∂²H/∂x_{axis}².

\[\frac{\partial^2 H_{ij}}{\partial x_{i,\text{axis}}^2} = \sigma''(Z_{ij}) \cdot W_{\text{axis},j}^2\]

Parameters:

• X (Tensor) – Input coordinates, shape (N, d).

• axis (int) – Spatial dimension index (0-based).

Return type:

Tensor

Returns:

Second-derivative matrix, shape (N, H).

laplacian(X, dims=None)

Compute the Laplacian feature matrix Σ_i ∂²H/∂x_i².

Parameters:

• X (Tensor) – Input coordinates, shape (N, d).

• dims (list[int] | None) – Dimensions to sum over. Defaults to all input dims.

Return type:

Tensor

Returns:

Laplacian feature matrix, shape (N, H).

class pypielm.core.feature_maps.FourierFeatureMap(input_dim, hidden_dim, freq_init='log_uniform', freq_min=1.0, freq_max=100.0, seed=42, w_scale=1.0, normalize_w=True, device='cpu', dtype=torch.float64)

Bases: Module

Generalised Fourier Feature ELM hidden layer.

Each hidden neuron computes:

\[\phi_j(\mathbf{x}) = \sqrt{2} \cos\!\left( \omega_j \, \mathbf{w}_j^\top \mathbf{x} + b_j \right)\]

where \(\mathbf{w}_j\) is a random direction, \(b_j \sim \mathcal{U}(0, 2\pi)\) is a random phase, and \(\omega_j\) is a per-neuron frequency coefficient.

Analytic derivatives:

\[ \begin{align}\begin{aligned}\frac{\partial \phi_j}{\partial x_i} = -\sqrt{2} \, \omega_j w_{ji} \sin(A_j)\\\frac{\partial^2 \phi_j}{\partial x_i^2} = -\sqrt{2} \, (\omega_j w_{ji})^2 \cos(A_j)\end{aligned}\end{align} \]

where \(A_j = \omega_j \mathbf{w}_j^\top \mathbf{x} + b_j\).

Parameters:

• input_dim (int) – Spatial dimension d.

• hidden_dim (int) – Number of neurons H.

• freq_init (Literal['uniform', 'log_uniform', 'auto']) – Frequency initialisation ('log_uniform', 'uniform').

• freq_min (float) – Minimum frequency.

• freq_max (float) – Maximum frequency.

• seed (int) – Random seed.

• w_scale (float) – Scale of Gaussian weight draw before normalisation.

• normalize_w (bool) – If True, unit-normalise each direction vector.

• device (str | device) – Target device.

• dtype (dtype) – Floating-point precision.

Initialize internal Module state, shared by both nn.Module and ScriptModule.

W: Tensor

b: Tensor

omega: Tensor

forward(X)

Compute φ(X) = √2 cos(ω (X @ W) + b).

Parameters:

X (Tensor) – Input coordinates, shape (N, d).

Return type:

Tensor

Returns:

Feature matrix, shape (N, H).

d1(X, axis)

Analytic ∂φ/∂x_{axis}.

Parameters:

• X (Tensor) – Input coordinates, shape (N, d).

• axis (int) – Spatial dimension index.

Return type:

Tensor

Returns:

First-derivative matrix, shape (N, H).

d2(X, axis)

Analytic ∂²φ/∂x_{axis}².

Parameters:

• X (Tensor) – Input coordinates, shape (N, d).

• axis (int) – Spatial dimension index.

Return type:

Tensor

Returns:

Second-derivative matrix, shape (N, H).

laplacian(X, dims=None)

Compute Σ_i ∂²φ/∂x_i².

Parameters:

• X (Tensor) – Input coordinates, shape (N, d).

• dims (list[int] | None) – Dimensions to include. Defaults to all input dims.

Return type:

Tensor

Returns:

Laplacian feature matrix, shape (N, H).

class pypielm.core.feature_maps.AutogradFeatureMap(input_dim, hidden_dim, activation_fn, seed=42, w_scale=1.0, device='cpu', dtype=torch.float64)

Bases: Module

Feature map with arbitrary activation, derivatives via autograd.

Wraps any user-supplied callable activation. Derivatives are computed via torch.autograd.grad() — slower than analytic but works for any smooth activation.

Parameters:

• input_dim (int) – Spatial dimension d.

• hidden_dim (int) – Number of neurons H.

• activation_fn (Callable[[Tensor], Tensor]) – Any differentiable callable f: Tensor → Tensor.

• seed (int) – Random seed.

• w_scale (float) – Weight scale.

• device (str | device) – Target device.

• dtype (dtype) – Floating-point precision.

Initialize internal Module state, shared by both nn.Module and ScriptModule.

W: Tensor

b: Tensor

forward(X)

H = activation(X @ W + b).

Return type:

Tensor

d1(X, axis)

Autograd ∂H/∂x_{axis}, shape (N, H).

Return type:

Tensor

d2(X, axis)

Autograd ∂²H/∂x_{axis}², shape (N, H).

Return type:

Tensor

Solvers

Linear solvers for ELM output-weight determination (PyTorch).

All solvers accept and return torch.Tensor objects and work on both CPU and CUDA devices. float64 (double precision) is strongly recommended for numerically stable PDE solves.

class pypielm.core.solver.BayesianSolveResult(beta_mean, beta_cov)

Bases: NamedTuple

Result of the Bayesian linear solve.

beta_mean: posterior mean of output weights, shape (H, out_dim). beta_cov: posterior precision matrix (inverse covariance), shape (H, H). Stored as precision (not covariance) for numerical stability; invert only when uncertainty quantification is needed.

Create new instance of BayesianSolveResult(beta_mean, beta_cov)

beta_mean: Tensor

Alias for field number 0

beta_cov: Tensor

Alias for field number 1

class pypielm.core.solver.WeightedLinearSystem(H, y, weight=1.0)

Bases: NamedTuple

One weighted observation block: y ≈ H @ beta + eps.

H: feature/design matrix, shape (N, H). y: target values, shape (N, out_dim) or (N,). weight: observation precision 1/σ²; rows are scaled by sqrt(weight) before the solve.

Create new instance of WeightedLinearSystem(H, y, weight)

H: Tensor

Alias for field number 0

y: Tensor

Alias for field number 1

weight: float

Alias for field number 2

pypielm.core.solver.ridge_solve(H, y, lam=1e-08)

Closed-form ridge regression output weights.

Solves:

\[\boldsymbol{\beta} = (\mathbf{H}^\top \mathbf{H} + \lambda \mathbf{I})^{-1} \mathbf{H}^\top \mathbf{y}\]

using torch.linalg.solve() for numerical stability.

Parameters:

• H (Tensor) – Feature matrix of shape (N, H).

• y (Tensor) – Target values of shape (N, out_dim) or (N,).

• lam (float) – Regularisation strength λ ≥ 0.

Return type:

Tensor

Returns:

Output weight matrix β of shape (H, out_dim).

pypielm.core.solver.rrqr_solve(H, y, tol=None)

Rank-revealing QR (divide-and-conquer SVD) least-squares solve.

Delegates to torch.linalg.lstsq() with driver='gelsd'.

Parameters:

• H (Tensor) – Feature matrix of shape (N, H).

• y (Tensor) – Target values of shape (N, out_dim) or (N,).

• tol (float | None) – Rank-truncation tolerance.

Return type:

Tensor

Returns:

Minimum-norm least-squares solution β of shape (H, out_dim).

pypielm.core.solver.bayesian_solve(blocks, prior_precision=0.0001)

Sequential Bayesian update over weighted observation blocks.

Computes the posterior:

\[ \begin{align}\begin{aligned}\boldsymbol{\Lambda}_{\text{post}} = \alpha \mathbf{I} + \sum_k \lambda_k \mathbf{H}_k^\top \mathbf{H}_k\\\boldsymbol{\beta}_{\text{post}} = \boldsymbol{\Lambda}_{\text{post}}^{-1} \sum_k \lambda_k \mathbf{H}_k^\top \mathbf{y}_k\end{aligned}\end{align} \]

Parameters:

• blocks (list[WeightedLinearSystem]) – List of WeightedLinearSystem objects.

• prior_precision (float) – Scalar α for the isotropic prior β ~ N(0, α⁻¹ I).

Return type:

BayesianSolveResult

Returns:

BayesianSolveResult with beta_mean and beta_cov (posterior precision matrix).

pypielm.core.solver.tikhonov_solve(H, y, L, lam=1e-08)

Generalised Tikhonov regularisation.

Solves:

\[\boldsymbol{\beta} = (\mathbf{H}^\top \mathbf{H} + \lambda \mathbf{L}^\top \mathbf{L})^{-1} \mathbf{H}^\top \mathbf{y}\]

Parameters:

• H (Tensor) – Feature matrix of shape (N, H).

• y (Tensor) – Target values of shape (N, out_dim) or (N,).

• L (Tensor) – Regularisation matrix of shape (r, H).

• lam (float) – Regularisation strength λ.

Return type:

Tensor

Returns:

Output weight matrix β of shape (H, out_dim).