accepted FastLSQ at the AI & PDE Workshop · ICLR 2026 read on OpenReview ↗
fastlsq v0.1.5 · differentiable digital twins, closed-form derivatives

A differentiable
digital twin, in closed form.

Born from solving PDEs: fastlsq fits a random-Fourier surrogate of your field in one least-squares call. You get the field and every derivative, ∂t, ∇, Δ, analytically, at any query point. A real digital twin you can drop into an optimiser, a Kalman filter, an inverse problem, or a control loop.

$ pip install fastlsq
view source read paper
§ 01 · live demo

Random sinusoids,
one solve, every derivative.

N features → least-squares → uₙ(x), u'ₙ(x), u''ₙ(x), …
N 4 ‖uₙ − u‖L²/‖u‖ · solve ·
target u(x) fit uₙ(x) features φⱼ
features N σ = 8.0
§ 02 · method

A differentiable surrogate
over a random Fourier basis.

u(x) ≈ Σⱼ βⱼ · sin(Wⱼ·x + bⱼ) · Dᵅu available in closed form

01 · the basis

The unknown solution is expanded in a fixed basis of random sinusoids: φⱼ(x) = sin(Wⱼ·x + bⱼ), with Wⱼ and bⱼ sampled once and never updated.

02 · why it works

Averaging enough random Fourier features recovers the Gaussian RBF kernel (Bochner). Smooth solutions live in that RKHS, so a few hundred sinusoids cover the space.

03 · the fit

Find β by a single least-squares call. Scalar or vector-valued targets: components stack into the same block-LSQ, one solve.

§ 03 · the differentiable twin

One β. Every derivative, exact.

solve once → query ∂ₓᵏ at any order, no re-fit

We fit β once against the target below. After that the twin is a closed-form function: differentiating it just multiplies each coefficient by Wⱼᵏ and shifts the phase by k·π⁄2. No re-solve, no finite differences, no autodiff graph. Step through the orders, the analytic curve (solid) lands exactly on the true derivative (dashed) every time.

∂ₓ⁰ uₙ(x) = N−½ Σ βⱼ sin(Wⱼx + bⱼ) max err · scale ·
true ∂ₓᵏu(x) analytic ∂ₓᵏuₙ(x)
derivative order k
β fixed · solved once
2D Poisson solve — same β, pick an operator:
true u(x, y)exact
true field
estimate uₙ(x, y)FastLSQ
FastLSQ estimate
error |uₙ − u|abs
absolute error
features N 1,500 operator applied analytically to β rel L² error 5.6 × 10⁻² re-solves 0 real fastlsq output
§ 04 · API

Compose any operator.
One formula handles it.

helmholtz_2d.py
helmholtz_2d.py · k = 10, [0, 1]² 
import numpy as np
import fastlsq as fl
from fastlsq.basis import SinusoidalBasis, Op

# Helmholtz 2D:  Δu + k² u = f   on  [0, 1]²
basis = SinusoidalBasis(d=2, N=1500, sigma=10.0)
L     = Op.laplacian(d=2) + (10.0**2) * Op.identity(d=2)

u, info = fl.solve(
    L      = L,
    f      = lambda x: 0.0,
    bc     = lambda x: np.sin(np.pi * x[..., 0]),
    domain = fl.Box(d=2),
    basis  = basis,
)

print(info.l2_error)   # 1.9e-6
print(info.wall_clock) # 0.08 s
status solved ✓ rel L² 1.9 × 10⁻⁶ wall-clock 0.08 s features 1500 scikit-fem P2 baseline 4.0 × 10⁻⁵
§ 05 · examples

Two physical demos.
One basis — solved once, or refit in the RL loop.

inverse heat-source localisation · radar stealth navigation
Animation: four hidden Gaussian heat sources localised from four sensors as L-BFGS-B drives the per-source position error down
recovered temperature field · sensor squares · source trajectories 1,200 iters · pos err 0.09

Inverse heat-source localisation.
Find the heaters from four sensors.

Four hidden heaters warm a square plate. Four sensors, one per quadrant, record a temperature time-series as the heat diffuses outward. From only those 4 × 60 = 240 numbers, recover where the heaters are — and their shape and strength: 24 unknown source parameters.

The forward model is the full space–time heat equation. Every step of the search needs a complete PDE solve, normally far too expensive to drop inside an optimiser that runs it a thousand times.

fastlsq fits the field in a sinusoidal random-Fourier basis whose ∂t and ∇² are closed-form. The (PDE + BC + IC) least-squares operator is assembled and Cholesky-factored once; each forward solve is then a single back-substitution, and the loss gradient over all 24 parameters rides back through the same factor in one adjoint solve.

Animation: a drone surfs the low-intensity corridors of a radar interference field from start to goal while a FastLSQ digital twin, refit every few steps, reconstructs the field along its path
radar field · drone comet-trail · twin refit every 3 steps 362 steps · undetected

Stealth navigation.
A world model refit in the loop.

A drone must cross a radar interference field corner to corner without the intensity it senses ever tripping a detector. It has no map — only a tiny five-point cross of field samples around its current position, streaming in as it moves.

From that trickle it refits a FastLSQ surrogate of the field every few steps — one Tikhonov least-squares solve, about 8 ms — then reads the surrogate's analytic gradient to surf the low-intensity corridors between the bright fringes toward the goal.

This is FastLSQ as a differentiable digital twin inside a control loop: cheap enough to re-solve hundreds of times per episode, exact-gradient enough to steer on. The twin is trustworthy only where the drone has sensed — so the right panel paints in along the path and fades elsewhere.

§ 06 · cite

The paper.

Full method, proofs, ablations, and 17-PDE benchmark are described in Sulc, A. FastLSQ: Solving PDEs in One Shot via Fourier Features with Exact Analytical Derivatives, arXiv:2602.10541.

Code, examples, and the inverse-problem demos are available in the reference implementation on GitHub.

bibtex 
@misc{sulc2026solvingpdesshotfourier,
  title         = {FastLSQ: Solving PDEs in One Shot 
                   via Fourier Features with 
                   Exact Analytical Derivatives},
  author        = {Antonin Sulc},
  year          = {2026},
  eprint        = {2602.10541},
  archivePrefix = {arXiv},
  primaryClass  = {math.NA},
  url           = {https://github.com/sulcantonin/FastLSQ},
}