Physics-Informed Functional Tucker Method with RKHS Factors for Sparse Spatiotemporal Reconstruction

Introduction

Reconstructing spatiotemporal fields aims to recover full physical fields from a small set of fixed sensor locations, which is a fundamental problem in scientific research [10, 14, 25, 43, 45, 47, 51]. Such sparsity leaves most locations and time instants unconstrained, making the inverse problem severely underdetermined: fine-scale structures such as vortices, fronts, and submesoscale processes may remain unresolved, and a model may fit the observed sensors while producing nonphysical off-sensor fields or unreliable missing-frame predictions [44]. This motivates reconstruction methods that are accurate at sensor locations while improving physical consistency on the missing-frame off-sensor locations.

Data-driven reconstruction methods have made substantial progress in sparse physical-field recovery [22, 48, 49], including sensor-to-field decoders, neural implicit representations, low-rank functional tensors, and generative latent models [6, 17, 18, 19]. A common paradigm is to first reconstruct fields or latent states from sparse measurements and then use interpolation, recurrent models, or generative priors to infer the temporal evolution [6]. Functional Tucker models offer compact representations of continuous fields by combining mode-wise factor functions with a low-dimensional core tensor [24, 38]. However, these methods remain largely driven by sparse observation fitting or learned priors. Continuous representability alone does not ensure physical consistency where observations are absent.

For systems with an available equation constraint, physical laws provide an additional source of supervision where sparse data provide none [44]. The observation loss is evaluated only at sensor locations, whereas the desired reconstruction is defined over a continuous domain. Finite-difference equation residuals can be imposed at off-sensor collocation points, providing an additional inductive bias in regions where no measurements are available [29]. Thus, PDE residuals are not merely post-hoc physical-consistency metrics or generic smoothing penalties; they reduce the ambiguity of sparse reconstruction by constraining the decoded field in the observation null space. To use such regularization, the representation must be continuously evaluable at sensors, collocation points, and neighboring stencil locations, while remaining compact enough for sparse inverse recovery. This motivates a structured continuous representation that separates spatial coordinate dependence from low-dimensional time-varying coefficients.

Motivated by this view, we propose PhysFTM, which represents the continuous dynamics as

where $D_{\Phi}$ is a spatially continuous Functional Tucker decoder and $G(t)$ is a continuous trajectory in Tucker-core space. The decoder can be evaluated at sensor locations, dense grids, and collocation points, while temporal evolution is modeled through low-dimensional core coefficients. To make the spatial decoder reliable under sparse supervision, PhysFTM parameterizes its mode-wise factors in reproducing kernel Hilbert spaces (RKHSs) and regularizes their RKHS norms [1, 39]. This replaces unconstrained neural mode functions with explicitly regularized functional factors. Together with finite-difference PDE residuals at collocation points, the resulting RKHS-FTM decoder provides both continuous spatial evaluability and physics-guided off-sensor guidance.

Our contributions are summarized as follows. 1) We propose PhysFTM, a functional Tucker framework that combines RKHS-based spatial factors with a continuous trajectory in Tucker-core space for sparse spatiotemporal reconstruction, paving the way for inquiries of the field at arbitrary times and positions. 2) We introduce a hybrid solve-learn optimization scheme that solves conditionally linear Tucker cores with regularized conjugate gradient method while learning RKHS factors under observation, RKHS, and PDE losses. 3) We provide theoretical guarantees for reconstruction via the proposed PhysFTM framework, and evaluate PhysFTM on Allen–Cahn, Navier–Stokes, Heat equations, and a gridded Copernicus ocean reanalysis, covering 2D/3D reconstruction, showing improved accuracy, off-sensor reconstruction accuracy, physical consistency, and missing frame reconstruction.

Related Work

Sparse reconstruction of physical fields. Early deep-learning approaches often transformed sparse measurements into structured grid representations, e.g., through Voronoi-assisted representations, followed by convolutional decoders for full-field recovery [13]. Senseiver [30] encodes unordered sparse sensor sets with attention and decodes field values at query coordinates. SHRED and SINDy-SHRED [16, 42] reconstruct high-dimensional spatiotemporal states from sensor histories using recurrent latent representations or identified latent dynamics. STRIDE [37] combines temporal encoders with modulated implicit decoders for continuous field reconstruction. Additionally, generative approaches were explored for physical field recovery, including conditional neural-field latent diffusion for spatiotemporal turbulence [9] and diffusion in Functional Tucker space for full-field dynamics from irregular sparse observations [6]. Existing approaches often rely on neural decoders that remain weakly constrained in off-sensor regions under sparse observations.

Functional tensor and low-rank continuous representations. Low-rank tensor decompositions provide compact multilinear representations for high-dimensional data, with Tucker decomposition being one of the classical formulations based on mode factors and a core tensor [15, 21, 38]. To extend such representations beyond fixed grids, one may consider continuous or functional variants that replace discrete factor matrices with coordinate-dependent factor functions, enabling field reconstruction and evaluation at arbitrary locations; see, e.g., [12, 24]. Functional Tucker models, where observations are reconstructed by shared latent functions and time-dependent Tucker cores [6], were adapted as latent spaces for sparse physical dynamics. Related functional tensor approaches include [8, 33, 40, 41], which offer alternative representations and priors for continuous tensor modeling.

Physics-informed reconstruction and PDE supervision. Physics-informed learning incorporates governing equations into neural approximation by penalizing PDE residuals at collocation points. PINNs are the canonical example, using residual losses together with data, boundary, or initial-condition terms to solve forward and inverse PDE problems [26]; broader physics-informed machine learning has been reviewed as a way to combine data with scientific constraints [20]. Physics-informed neural operators extend this idea to operator learning by combining data losses with PDE residuals or physical consistency terms [23, 46]. Discrete and spatiotemporal architectures such as PhyCRNet impose residual constraints in convolutional-recurrent models for PDE dynamics [29]. In reconstruction settings, Hidden Fluid Mechanics uses Navier–Stokes constraints to infer velocity and pressure fields from indirect observations [27], and physics-informed diffusion models incorporate physical conditions to improve high-fidelity flow-field reconstruction [31]. These works show that PDE residuals can provide supervision beyond labeled samples and improve physical consistency.

Methodology of PhysFTM

Problem setup and PhysFTM decomposition

Sparse observations. For clarity, we describe the formulation for a single trajectory in two spatial dimensions and omit the trajectory index. Let $u(\mathbf{x}, t, c)$ denote the $c$-th physical channel at time $t$ and spatial coordinate $\mathbf{x} = (x, y)$. The data are represented by $\mathbf{Y} \in\mathbb{R}^{H \times W \times T \times C}$, where $H \times W$ is the spatial grid resolution, $T$ is the number of frames, and $C$ is the number of physical channels. Only a sparse subset of entries is available as observations. We denote the observation mask by $\mathbf{M} \in\{0,1\}^{H \times W \times T \times C}$ and the observed index set by

All trajectories share a common nominal time grid. Sensor locations are fixed over time within each trajectory, and other spatial locations are referred to as off-sensor locations. Frames with sensor observations are referred to as support frames, while frames without sensor observations are referred to as missing frames. We aim to reconstruct the full field $\hat{u}(\mathbf{x}, t, c)$ on both support and missing frames from these sparse observations.

Spatial–temporal factorization. PhysFTM combines a spatially continuous Functional Tucker decoder with a time-dependent Tucker core:

where $\mathcal{D}_{\Phi}$ is a spatial decoder defined by RKHS-parameterized Functional Tucker factors, and $G(t)$ is the time-dependent Tucker core. Sparse sensor observations in support frames are used to infer the corresponding cores $G(t_j)$. A continuous-time model then estimates the cores of missing frames from these inferred support-frame cores, and $\mathcal{D}_{\Phi}$ decodes them into reconstructed fields. The spatial decoder provides a continuous spatial representation and can be regularized by PDE residuals, while temporal modeling operates in the compact Tucker-core space.

Physics-informed RKHS Functional Tucker decoder

RKHS mode factors and Tucker decoding. For each active spatial mode $d \in\{x,y\}$, PhysFTM defines $R_d$ factor functions $\{\phi_{d,r}\}_{r=1}^{R_d}$. Each factor belongs to a reproducing kernel Hilbert space induced by a one-dimensional positive definite kernel $k_d$ and is represented by a finite kernel expansion

where $\{z_{d,m}\}_{m=1}^{M_d}$ are inducing points and $A_d=[a_{d,m,r}]\in\mathbb{R}^{M_d\times R_d}$ are learnable RKHS coefficients. By using Eq. (2), we provide the following spatial decoders via Tucker decomposition. For a scalar field at a support frame $t$, the spatial decoder is

where $G$ is a core tensor. For multichannel fields, the spatial decoder is modeled by

Note that factors $\phi_{x,p}(x)$ and $\phi_{y,q}(y)$ are shared between channels, and the channel index is stored in the core tensor. The decoder can be evaluated at observed sensors, dense grid points, and off-grid PDE collocation points.

RKHS and physics-informed loss functions and alternating-minimization strategy. PhysFTM exploits the block structure of Eq. (1). In general, we propose an alternating-minimization strategy to optimize the factors and core tensor. Specifically, with fixed RKHS factors, support-time Tucker cores can be inferred by solving a least-squares problem. For RKHS factors, since the factors are into 143 the PDE residuals, they are learned by gradient descent under observation, RKHS, and PDE losses.

For a mini-batch and fixed factors, the support-time core sequence is obtained by

Here $T_s$ denotes the number of support frames, $\Phi(\mathbf{x})=\bigotimes_{d=1}^{d_s}\phi_d(x_d)$ is the tensor-product spatial feature, and $M_t$ is the sparse observation mask at time $t$. The subproblem is solved by conjugate gradients using implicit matrix-vector products. Note that the matrix $\Phi(\mathbf{x})$ is never explicitly computed.

After solving the cores, the RKHS basis coefficients are updated by AdamW by minimizing

with the loss terms defined below. The observation loss is

This is the RMSE over observed scalar entries of the channelwise normalized fields. For fixed factors, the core update uses the squared-error subproblem in Eq. (5). Note that in practice, for a mini-batch of $N_b$ trajectories, let $\Omega_b$ denote the observed index set of trajectory $b$, we evaluate the loss on each $\Omega_b$ separately and compute the weighted average by Eq. (5).

The RKHS regularization term controls the complexity of the spatial factor functions:

where $K_d$ is the kernel Gram matrix over inducing points of mode $d \in\{x,y\}$. Given a mini-batch of trajectories and spatial collocation points $Q$, the decoder evaluates $\hat{u}(x,t,c)$ for each $x \in Q$. For a dataset-specific finite-difference residual operator $\mathcal{R}_h$, the PDE loss is

where $\mathcal{T}'$ contains the timestamps of support frames with valid finite-difference stencils, and $\mathcal{T}' = \mathcal{T}'/\mathcal{T}$. The derivative stencils and valid spatial and temporal points depend on the dataset; see Appendix § 8.

The time-kernel term regularizes the inferred support-frame core sequence:

where $G_b^{\flat} \in\mathbb{R}^{T_s \times D_G}$ is the flattened core sequence of trajectory $b$, $T_s$ is the number of support frames, $D_G = C R_x R_y$ is the flattened-core dimension, and $K_t$ is an RBF kernel matrix over the timestamps of support frames. For three-dimensional fields, the RKHS regularization also includes the $z$-mode factors, and $D_G = C R_x R_y R_z$.

RKHS-Tucker reconstruction analysis. We provide the reconstruction error analysis of the RKHS-Tucker formulation in Appendix § 9.

Note that the first term in Eq. (11) is the truncation error induced by the finite Tucker rank, and the second term comes from the statistical error induced by sparse sampling.

Continuous-time core trajectory and missing frames reconstruction

After learning the spatial RKHS-FTM decoder, PhysFTM freezes the RKHS factor functions and represents temporal variation only through the time-dependent Tucker core. The Tucker core provides a compact coefficient representation for modeling temporal dynamics under the spatial decoder. We derive the temporal model from support-core estimates inferred with the RKHS decoder learned during the physics-constrained spatial stage.

Support-core inference with a mean-core reference. For a new trajectory, the learned RKHS factors are frozen and reused as a deterministic spatial decoder. The remaining task is to infer support-frame cores from sparse sensor measurements. Under severe sparsity, this inverse problem is highly underdetermined: the sensing operator observes only a small subset of entries, leaving a large null space in core space. If the core is solved directly with a ridge penalty, the solution can be biased toward a minimum-norm core; in practice, especially with finite-iteration CG and extremely sparse masks, this often leads to near-zero cores and low-energy reconstructions.

Therefore, we introduce a residual correction:

on the data-driven mean-core reference $G_t^{\mathrm{ref}}$ computed from the training set,

Instead of regularizing the absolute core, we penalize only the correction. The core-correction subproblem is defined on the complete native time grid:

Here $T$ is the number of native timestamps. The observation mask is zero at missing frames, while the last summand couples corrections at adjacent native timestamps. The inferred support-time cores are

This residual formulation keeps the solution close to a plausible training-set core trajectory when observations are insufficient, while still allowing sparse measurements to correct the trajectory through $\Delta G$.

Temporal core modeling. Let $\{(t_j, \bar{G}(t_j))\}_{j=1}^{T_s}$ denote the inferred support-frame cores and their timestamps. A continuous-time model $H_\psi$ estimates the cores of missing frames:

where $\mathcal{T}_{\mathrm{miss}}$ contains the timestamps of missing frames within $[t_1, t_{T_s}]$.

We implement $H_\psi$ as a learned residual correction to cubic Hermite interpolation in core space. Let $B(\tau)$ denote the Hermite interpolant through the support-frame cores. Between two adjacent support-frame timestamps $t_a$ and $t_b$, the estimated core is

where $\Delta_\psi(\tau)$ is generated by a lightweight ODE-RNN in flattened core space. The gate is zero at both endpoints, preserving the inferred support-frame cores.

Missing-frame decoding. The reconstructed field at spatial location $x$ and missing-frame timestamp $\tau$ is obtained through the learned spatial decoder:

Experiments

Experimental setup

Datasets. We use Allen–Cahn (AC2D), a scalar reaction–diffusion benchmark with 200 trajectories of 201 frames on a $256 \times256$ grid, split into 160/20/20 train/validation/test trajectories; Navier–Stokes (NS2D), an incompressible flow benchmark with $(u, v, \omega)$ channels and 800 trajectories of 201 frames on a $128 \times128$ grid, split into 640/80/80; Heat3D, a scalar three-dimensional heat-diffusion benchmark with 200 trajectories of 50 frames on a $64^3$ grid, split into 160/20/20 train/validation/test trajectories; and Ocean3D, a three-dimensional ocean-velocity reanalysis dataset derived from the Copernicus Marine GLORYS12V1 product [11], with $(u, v)$ channels and 100 sequences of 90 frames on a $32 \times64 \times64$ grid, split chronologically into 80/10/10 train/validation/test sequences. Dataset details are in Appendix § 10.

Observation and missing-frame protocols. During training, 15% of spatial locations are sampled to simulate sparse grid-based sensors. The sensor mask is fixed over time within each trajectory and shared by physical channels, while different trajectories use different masks. We refer to the sensor positions of each trajectory as sensor locations and all remaining spatial locations as off-sensor locations. We report spatial sensor ratios $s \in\{0.01, 0.03, 0.15\}$. We use two missing-frame protocols: random-drop-half independently retains internal frames as support frames with probability 0.5 and always retains both endpoints, whereas temporal-interval-missing withholds sensor observations from a contiguous internal block of 20% of the frames. Frame masks are trajectory specific; full protocol and implementation details are in Appendix § 10.§ 2.

Compared Methods. We compare PhysFTM with SDIFT [6], LRTFR [24], SIREN-MLP [32], PINNs [26], and RecFNO [50]. Their conditioning, training, and test-time protocols are summarized in Appendix § 5.

Metrics. We report root mean squared error (RMSE) at sensor locations in support frames (S-sup), off-sensor locations in support frames (Off-sup), sensor locations in missing frames (S-miss), and off-sensor locations in missing frames (Off-miss). The four disjoint subsets and aggregation rules are defined in Appendix § 10.§ 1. All metrics use inverse-normalized physical variables.

Reconstruction results

Overall reconstruction under missing-frame protocols. We first compare reconstruction at $s = 0.03$ across the four datasets and both missing-frame protocols. Table 1 reports metrics separately for support frames and missing frames, each evaluated at sensor locations and off-sensor locations.

Table 1: Reconstruction RMSE on four datasets at spatial sensor ratio $s = 0.03$. The two column groups report the four disjoint evaluation subsets, defined by sensor/off-sensor locations and support/missing frames, under random-drop-half and temporal-interval-missing. Lower is better. For Heat3D and Ocean3D, the PhysFTM interval-missing results use the random-drop-half-trained temporal checkpoint. (Bold: best, Underline: second best)

Method	Random-drop-half				Temporal-interval-missing
	S-sup	Off-sup	S-miss	Off-miss	S-sup	Off-sup	S-miss	Off-miss
AC2D
SDIFT	0.35804	0.38080	0.38030	0.38151	0.35699	0.36199	0.39071	0.39259
LRTFR	0.31170	0.32220	0.31340	0.32038	0.31006	0.31676	0.32059	0.32781
PhysFTM	0.01735	0.08710	0.05035	0.09988	0.01633	0.09397	0.03953	0.07924
RecFNO	0.01848	0.45283	0.02178	0.45523	0.01811	0.45875	0.02116	0.43680
SIREN-MLP	0.34868	0.34932	0.35013	0.35076	0.34605	0.34664	0.35819	0.35898
PINNs	0.35166	0.35243	0.35268	0.35331	0.34874	0.34943	0.36508	0.36622
NS2D
SDIFT	0.52867	0.50970	0.51003	0.51027	0.53019	0.53052	0.50434	0.50275
LRTFR	0.13245	0.14720	0.13447	0.14576	0.13515	0.14682	0.19025	0.19661
PhysFTM	0.00669	0.06210	0.03240	0.06498	0.00435	0.06282	0.14949	0.15056
RecFNO	0.06810	0.09757	0.07820	0.10441	0.06906	0.09910	0.35803	0.35943
SIREN-MLP	0.24555	0.24526	0.24465	0.24441	0.24797	0.24777	0.21369	0.21354
PINNs	0.34254	0.34177	0.34101	0.34055	0.34722	0.34657	0.31892	0.31848
Heat3D
SDIFT	0.05811	0.05808	0.06303	0.06272	0.06284	0.06268	0.05016	0.04989
LRTFR	0.01212	0.01222	0.01309	0.01315	0.01210	0.01225	0.00577	0.00577
PhysFTM	0.00108	0.00125	0.00575	0.00574	0.00122	0.00137	0.00310	0.00305
RecFNO	0.28736	0.77847	0.26372	0.72034	0.28629	0.77526	0.22597	0.64187
SIREN-MLP	0.04080	0.04087	0.04211	0.04205	0.04342	0.04344	0.02749	0.02742
PINNs	0.04000	0.04011	0.04271	0.04264	0.04314	0.04319	0.03353	0.03334
Ocean3D
SDIFT	0.27863	0.28065	0.27585	0.27838	0.27494	0.27725	0.28712	0.28941
LRTFR	0.11038	0.11450	0.11767	0.12155	0.10506	0.10829	0.18062	0.18377
PhysFTM	0.00632	0.06574	0.07133	0.08858	0.00595	0.06382	0.18985	0.18565
RecFNO	0.73827	0.97074	0.77688	0.96302	0.71841	0.94379	0.99191	1.05383
SIREN-MLP	0.20158	0.20168	0.20024	0.20094	0.19773	0.19833	0.21199	0.21151
PINNs	0.30163	0.30100	0.30232	0.30194	0.29765	0.29739	0.31870	0.31725

PhysFTM obtains the lowest RMSE in most reported dataset–protocol–subset comparisons. The gains are especially clear at off-sensor locations in support frames on AC2D, NS2D, and Heat3D. These results support the effectiveness of the proposed representation for fixed-sensor reconstruction with missing frames across both 2D and 3D grids.

Complementary physical diagnostics. Table 2 reports NS2D missing-frame accuracy together with divergence, curl-derived vorticity, and kinetic-energy-spectrum errors. The diagnostics are computed on inverse-normalized full grids of missing frames; their formulas and cross-frame reducers are given in Appendix § 8.§ 1.

Table 2: NS2D missing-frame reconstruction and complementary physical diagnostics at $s = 0.03$. Divergence, vorticity, and spectrum errors are evaluated on full grids of missing frames; Off-miss is evaluated at off-sensor locations in missing frames. Lower is better.

Protocol	Method	Off-miss	Divergence	Vort. RMSE	Spectrum
	SDIFT	0.510274	0.352407	0.916343	0.877309
	LRTFR	0.145757	0.459044	0.237923	0.049954
Random	PhysFTM	0.064976	0.134170	0.109304	0.033599
	RecFNO	0.104413	0.459014	0.118608	0.067430
	SIREN-MLP	0.244405	1.887604	0.623477	0.408029
	PINNs	0.340554	0.120243	0.586433	0.398297
	SDIFT	0.502748	0.344546	0.931635	0.892291
	LRTFR	0.196613	0.616958	0.272195	0.055034
Interval	PhysFTM	0.150558	0.129873	0.150436	0.063646
	RecFNO	0.359434	0.342736	0.211528	0.055423
	SIREN-MLP	0.213537	1.863564	0.546894	0.368662
	PINNs	0.318479	0.120243	0.586433	0.349142

PhysFTM has the lowest Off-miss, vorticity, and spectrum errors under both protocols, while PINNs has the lowest divergence residual. The result indicates a favorable accuracy–consistency trade-off. Complete three-ratio diagnostics are reported in Appendix § 8.§ 1.

Model analysis

The effects of individual components and block-update strategies are evaluated through ablations. Table 3 compares RKHS regularization, PDE regularization, and core/factor block-update strategies with their corresponding effective objectives, without temporal missingness.

Data	Variant	PDE	Full	Sensor	Off-sensor
AC2D	PhysFTM	0.053250	0.088078	0.015185	0.089395
AC2D	w/o physics	0.064986	<u>0.091135</u>	<u>0.010852</u>	0.092520
AC2D	w/o RKHS reg.	<u>0.048491</u>	0.096622	0.027551	0.097985
AC2D	solve-only	1.155543	0.171035	0.008657	0.173653
AC2D	gradient-only	1.76 × 10⁻⁵	0.352111	0.351441	0.352132
NS2D	PhysFTM	<u>0.662575</u>	0.060270	<u>0.004189</u>	0.061196
NS2D	w/o physics	1.045838	<u>0.085478</u>	0.009674	<u>0.086781</u>
NS2D	w/o RKHS reg.	1.703260	0.113713	0.016850	0.115421
NS2D	solve-only	1.625978	0.105585	0.000653	0.107207
NS2D	gradient-only	0.450487	0.090639	0.073760	0.091112

The complete model gives the lowest Full and Off-sensor RMSE on both datasets. Although solve-only achieves the lowest Sensor RMSE, its substantially higher Off-sensor RMSE indicates that accurate fitting at sensor locations does not necessarily translate into accurate reconstruction elsewhere. The gradient-only variant obtains a low PDE residual on AC2D together with a high reconstruction RMSE. These comparisons assess the combined choices of block updates and their effective objectives.

Dataset	$s$	Learned	Linear	Cubic
AC2D	0.01	0.180029	0.189510	0.189556
AC2D	0.03	0.085552	0.090391	0.090388
NS2D	0.01	0.130469	0.137667	0.137453
NS2D	0.03	0.055546	0.064605	0.063790

Overall, PhysFTM offers the most favorable trade-off among all evaluated variants between sensor fitting, off-sensor reconstruction, and physical consistency. The block updates and active objective terms are specified in Appendix § 12. Additional spatial and collocation ablations are in Appendix § 7.

Reconstruction at non-grid-aligned times. We evaluate the learned core trajectory at 100 non-grid-aligned timestamps within the interval spanned by the timestamps of support frames. These evaluation timestamps do not coincide with the original frame timestamps and are therefore distinct from missing frames. For both AC2D and NS2D, the ground truth is recorded directly by rerunning the PDE solver at the perturbed physical timestamps rather than interpolating stored frames. The timestamp and split rules are given in Appendix § 10.§ 3.

The learned trajectory reduces reconstruction RMSE at non-grid-aligned times relative to both interpolation baselines in all four settings. This experiment supports continuous temporal reconstruction within the interval spanned by the timestamps of support frames on AC2D and NS2D.

Conclusion

We introduce PhysFTM, a functional Tucker method for reconstructing continuous spatiotemporal dynamics from sparse sensors and missing frames. By parameterizing Tucker factors by RKHS and enforcing governing equations via finite-difference residuals, the proposed framework is able to improve the physical consistency and reconstruction accuracy. Experiments on Allen–Cahn, Navier–Stokes, Heat equations, and a gridded Copernicus ocean reanalysis demonstrate accurate and physically consistent reconstruction across 2D and 3D grids.

The main limitations are as follows. First, the current framework considers sparse observations collected from fixed spatial sensors, with missing temporal frames. Extending it to more general observation settings, where measurements are available at arbitrary space-time locations, remains an interesting direction. Second, PhysFTM relies on known governing equations to regularize the reconstruction, which may limit its applicability when the underlying physical models are inaccurate or incomplete. Third, sparse observations generally do not uniquely determine the underlying physical field. Current PhysFTM returns a single regularized field without quantifying the associated reconstruction uncertainty. Therefore, we may extend the method to observations at arbitrary space-time coordinates and incorporate uncertainty quantification in the future to characterize reconstruction ambiguity and provide reliable estimates for downstream analysis.

Acknowledgement

This work was supported by the National Natural Science Foundation of China (Grant No. 42450181) and by a grant from the Research Grants Council of the Hong Kong Special Administrative Region, China (Grant No. AoE/P-601/23-N).

Method	Test input	Training / fitting	Supervision / equation use	Test optimization	Spatial scope	Time interface / selection
SDIFT	Sparse spatiotemporal observations	Cross-trajectory diffusion training	Dense training trajectories; no test missing-frame targets	Posterior sampling	2D/3D	Native-grid protocols; executed checkpoint
LRTFR	Sparse coordinate–value tuples	Per-trajectory representation fitting	Sensor measurements; no equation loss	Per-trajectory fit	2D/3D	Native-grid evaluation; fitted objective
SIREN-MLP	Coordinates; no test-mask conditioning	Train once on training split	Training fields; no equation loss	None	2D/3D	Coordinate output; executed checkpoint
PINNs	Coordinates; no test-mask conditioning	Train once on training split	Training fields and implemented equation residual	None	2D/3D	Coordinate output; executed checkpoint
RecFNO	Gridded sparse field and mask	Cross-trajectory operator training	Dense training fields; no equation loss	None	2D/3D	Prescribed-grid support/missing frames; executed checkpoint
PhysFTM	Sparse support-frame sensor values	Cross-trajectory decoder/head training	Observations, missing-frame targets, and equation loss where enabled	Support-core CG solve	2D/3D	Native-grid and AC/NS in-range non-grid-aligned reconstruction; spatial: training dense RMSE; temporal: validation missing-frame relative $\ell_2$

Comparison methods and spatial analysis

Compared methods. SDIFT [6] is a cross-trajectory diffusion model conditioned on sparse spatiotemporal observations. LRTFR [24] fits a functional tensor representation to the sensor coordinate–value measurements of each target trajectory. RecFNO [50] is trained across trajectories and maps gridded sparse-field/mask inputs to prescribed-grid outputs. SIREN-MLP [32] and PINNs [26] are coordinate-to-field controls trained once on the training split and evaluated directly at test coordinates without test-time adaptation; they share the same dense prediction across spatial ratios and temporal protocols, while only the evaluated subsets change. The earlier response used “SIREN-INR” for the same SIREN-MLP implementation.

PhysFTM trains the spatial stage for 500 epochs, where one epoch traverses the complete training loader, and selects the spatial checkpoint by the lowest dense-field training RMSE after inverse normalization, averaging evaluation mini-batch RMSEs. Its temporal head is trained for 20 epochs with AdamW and selects the checkpoint with the lowest validation relative $\ell_2$ error over missing frames, using decoded fields at all spatial locations after inverse normalization. RecFNO is evaluated on prescribed-grid support frames and missing frames.

The Panel A comparisons isolate basis parameterization and core optimization, whereas Panel B isolates collocation placement. These diagnostics are separate from the component/optimization ablation in Table 3, whose PDE column uses a disjoint held-out off-sensor set and internal centered-stencil times.

Finite-difference residuals and physical diagnostics

All derivatives are evaluated after inverse normalization. AC2D uses the Allen–Cahn residual

NS2D stores $(u,v,\omega)$ rather than pressure. Its held-out residual reuses the centered finite-difference vorticity-equation component of the training operator,

and excludes the first and last temporal points. Training collocation and held-out off-sensor evaluation points are disjoint. Heat3D uses $\mathcal{R}_{\mathrm{heat}} = D_t\hat{u} - 0.01(D_{xx}+D_{yy}+D_{zz})\hat{u}$ in the spatial stage with weight $10^{-4}$, while its temporal equation weight is zero. Ocean3D has no governing-equation loss in the reported configuration.

Method	AC2D			NS2D
	RMSE	SSIM	PDE	RMSE	SSIM	PDE
Linear interp.	0.1275	0.6379	0.5877	0.0897	0.9322	1.1605
Cubic interp.	0.1051	0.7632	0.7174	0.0653	0.9656	0.8339
SIREN-MLP	0.2907	0.1747	1.1703	0.2155	0.5836	2.1245
MLP-FTM	0.1049	0.7515	0.2282	0.0859	0.9116	4.0426
MLP-FTM+CG	0.1026	0.7719	0.1995	0.0832	0.9276	1.2277
FTM-RKHS	0.0925	0.7893	0.1532	0.0868	0.9468	0.7371
PhysFTM	0.0894	0.7938	0.1325	0.0612	0.9636	0.3500

Placement	AC2D		NS2D
	RMSE	PDE	RMSE	PDE
w/o equation	0.0925	0.1532	0.0868	0.7371
sensor locations	0.0898	0.1358	0.0622	0.3584
off-sensor locations	0.0884	0.1316	0.0609	0.3338
full domain	0.0897	0.1344	0.0618	0.3556

For NS training, the equation-based loss is the sum of three mean squared residuals,

multiplied by the overall PDE-loss weight. The reported NS equation-residual metric is the mean squared vorticity-equation residual, $\operatorname{mean}(R_{\mathrm{NS}}^2)$. The main prediction and ground-truth PDE exports use this same single-equation convention, while the unweighted sum of the three constraints is retained as a separate auxiliary export. The ground-truth export evaluates the numerical equation residual of the reference field.

Missing-frame physical diagnostics

The NS2D diagnostics are evaluated on inverse-normalized full grids of missing frames. Spatial derivatives use centered differences, and Fourier wavenumbers use the physical domain lengths. Divergence measures the predicted velocity field itself, $\lvert\nabla\cdot\hat{\mathbf{u}}\rvert$. Vorticity is derived independently from predicted and reference velocities and then compared:

Here the sums pool all test trajectories, their missing frames, and all full-grid points before taking the square root. For the spectrum diagnostic, the spatial means of $(u,v)$ are removed, a radial shell-sum kinetic-energy spectrum $E_{n,i}(k)$ is formed, and the frame error is

Frame errors are averaged within each trajectory, after which the evaluator computes the mean and standard deviation of the trajectory means. The tables report the mean component.

Diagnostic	Missing-frame/spatial reduction	Trajectory reduction
Divergence error	Global RMS over all missing frames and grid points	Included in the same global pool
Vorticity RMSE	Global RMSE over all missing frames and grid points	Included in the same global pool
Spectrum error	Per-frame relative $L_2$, then missing-frame mean	Mean/std of trajectory means

Protocol	$s$	Method	Off-miss	Divergence	Vort. RMSE	Spectrum
Random	0.01	SDIFT	0.509261	0.307569	0.876079	0.865337
		LRTFR	0.170679	0.475408	0.237857	0.048102
		PhysFTM	0.130659	0.274083	0.232358	0.104305
		RecFNO	0.413223	0.514580	0.396548	0.161395
		SIREN-MLP	0.244360	1.886754	0.623300	0.409266
		PINNs	0.341072	0.120243	0.586433	0.402015
	0.03	SDIFT	0.510274	0.352407	0.916343	0.877309
		LRTFR	0.145757	0.459044	0.237923	0.049954
		PhysFTM	0.064976	0.134170	0.109304	0.033599
		RecFNO	0.104413	0.459014	0.118608	0.067430
		SIREN-MLP	0.244405	1.887604	0.623477	0.408029
		PINNs	0.340554	0.120243	0.586433	0.398297
	0.15	SDIFT	0.510589	0.782450	1.288744	0.882962
		LRTFR	0.138756	0.467318	0.236704	0.050423
		PhysFTM	0.033357	0.029145	0.042455	0.029433
		RecFNO	0.040051	0.460688	0.729675	0.041657
		SIREN-MLP	0.244476	1.887873	0.623739	0.407336
		PINNs	0.340193	0.120243	0.586433	0.399324
Interval	0.01	SDIFT	0.492242	0.304216	0.881166	0.870480
		LRTFR	0.206649	0.889423	0.311479	0.069965
		PhysFTM	0.164821	0.270320	0.251761	0.128132
		RecFNO	0.478397	0.374108	0.347071	0.178211
		SIREN-MLP	0.206503	1.858348	0.528179	0.356109
		PINNs	0.306602	0.120243	0.586433	0.346441
	0.03	SDIFT	0.502748	0.344546	0.931635	0.892291
		LRTFR	0.196613	0.616958	0.272195	0.055034
		PhysFTM	0.150558	0.129873	0.150436	0.036346
		RecFNO	0.359434	0.342736	0.211528	0.055423
		SIREN-MLP	0.213537	1.863564	0.546894	0.368662
		PINNs	0.318479	0.120243	0.586433	0.349142
	0.15	SDIFT	0.485171	0.816494	1.383403	0.925256
		LRTFR	0.172046	0.550371	0.254336	0.061462
		PhysFTM	0.134567	0.026805	0.105232	0.010536
		RecFNO	0.318191	0.321772	0.172333	0.061395
		SIREN-MLP	0.203190	1.853913	0.520432	0.374622
		PINNs	0.300163	0.120243	0.586433	0.373544

Theoretical analysis on spatial reconstruction

In this section, we conduct an idealized spatial reconstruction error analysis for the RKHS-parameterized Functional Tucker Model (FTM-RKHS) under sparse observations. We aim to characterize how the modewise truncation level, kernel spectrum, and number of measurement samples jointly determine the reconstruction error. The i.i.d. spatial-sampling assumption is an analytical idealization and is not identical to the fixed-within-trajectory sensor protocol used in the experiments. Moreover, the spectral truncation level $r$ in this analysis is not identified with the implemented Tucker rank.

Let $u : \mathcal{Z} \to\mathbb{R}$ denote the target field, where $\mathcal{Z}$ contains trajectory, time, spatial coordinates, and channel indices. The field is measured only on a sparse set $\Omega= \{z_i\}_{i=1}^n$, with $n = |\Omega|$, through measurements $y_i = u(z_i) + \epsilon_i$, where $\mathbb{E}[\epsilon_i \mid z_i] = 0$. The task is to recover the field $u$ via the measurements $y_i$.

For the sake of theoretical analysis, we consider the following assumptions. For high-probability bounds, we assume that the noise is conditionally sub-Gaussian, i.e., $\mathbb{E}[\exp(t\epsilon_i) \mid z_i] \le $\exp\left(\frac{\sigma^{2}t^{2}}{2}\right)$, $t \in\mathbb{R}$. We assume that the observation distribution echoes the target distribution $\rho$ over the full space $\mathcal{Z}$, i.e., the observations $z_i$ are sampled i.i.d. from $\rho$.

From one-dimensional RKHS to tensor-product RKHS. We first recall the one-dimensional RKHS construction for each mode, following the standard spectral characterization of RKHSs induced by compact positive kernel integral operators [1, 4, 35]. Let $(\mathcal{Z}_j,\mathcal{B}_j,\rho_j)$, $j=1,\ldots,d$, be the probability space associated with the $j$-th mode, where $\mathcal{Z}_j$ is the input domain, $\mathcal{B}_j$ is the corresponding sigma-algebra, and $\rho_j$ is the probability measure on $\mathcal{Z}_j$. Let $\mathcal{H}_j$ be an RKHS on $\mathcal{Z}_j$ with kernel $k_j$, the kernel integral operator is defined by $(T_jf)(z_j)=\int_{\mathcal{Z}_j} k_j(z_j,z'_j)f(z'_j)\,d\rho_j(z'_j)$. Assume that each $T_j:L^2(\rho_j)\to L^2(\rho_j)$ is compact, self-adjoint, and positive semidefinite, so that the Mercer-type spectral decomposition applies [7, 35]. Then, there exist eigenpairs $T_j\phi_{j,m}=\lambda_{j,m}\phi_{j,m}$, with $\lambda_{j,1}\ge\lambda_{j,2}\ge\cdots>0$, where $\{\phi_{j,m}\}_{m\ge1}$ is orthonormal in $L^2(\rho_j)$.

For $f_j=\sum_{m\ge1}a_{j,m}\phi_{j,m}$ and $g_j=\sum_{m\ge1}b_{j,m}\phi_{j,m}$, the one-dimensional $L^2$ and RKHS inner products are defined by

which induces the norms $\|f_j\|_{L^2(\rho_j)}^2=\sum_{m\ge1}a_{j,m}^2$ and $\|f_j\|_{\mathcal{H}_j}^2=\sum_{m\ge1}a_{j,m}^2/\lambda_{j,m}$. Equivalently, $\{\sqrt{\lambda_{j,m}}\phi_{j,m}\}_{m\ge1}$ forms an orthonormal system in $\mathcal{H}_j$.

We now extend the construction to the tensor setting using the standard tensor-product RKHS construction [1, 4]. Let $\mathcal{Z}=\mathcal{Z}_1\times\cdots\times\mathcal{Z}_d$, and assume that the full coordinate distribution is $\rho=\rho_1\otimes\cdots\otimes\rho_d$. For the sake of brevity, we assume that the observation distribution echoes the target distribution.

Define the tensor-product RKHS $\mathcal{H}=\bigotimes_{j=1}^{d}\mathcal{H}_j$ with tensor-product kernel $K(z,z')=\prod_{j=1}^{d}k_j(z_j,z'_j)$, which is positive definite by the closure properties of reproducing kernels under tensor products [1, 4]. For a multi-index $m=(m_1,\ldots,m_d)\in\mathbb{N}^d$, define

Because $\rho$ is a product measure, it holds that $\langle\Phi_m,\Phi_n\rangle_{L^2(\rho)}=\prod_{j=1}^{d}\delta_{m_jn_j}=\delta_{mn}$, i.e., $\{\Phi_m\}_{m\in\mathbb{N}^d}$ is orthonormal in $L^2(\rho)$.

The tensor-product integral operator is $T=\bigotimes_{j=1}^{d}T_j$, and the preceding definitions imply $T\Phi_m=\Lambda_m\Phi_m$. Thus, the tensor spectrum follows from the one-dimensional spectra and the product-measure structure.

For two functions $f=\sum_m a_m\Phi_m$ and $g=\sum_m b_m\Phi_m$, the tensor-product RKHS inner product is given by

which introduces the norm

Truncated RKHS-Tucker function class. For a fixed truncation level $r$, define the rectangular modewise truncation region

Its complement is $\mathbb{N}^d\setminus\mathcal{R}_r=\{m\in\mathbb{N}^d:\exists j,\ m_j>r\}$, which contains all index tuples where at least one mode enters the high-frequency region $\{r+1,r+2,\ldots\}$.

Using the tensor-product bases defined above, the truncated RKHS-Tucker function class is

Let $u_{r}^{\star}$ be the best approximation of $u$ within $\mathcal{F}_{r}$, namely

The empirical estimator is denoted by $\hat{u}_{r}$. The total reconstruction error satisfies

approximation error \hspace{4em} statistical error

Approximation error. We evaluate the approximation error by using a specific RKHS-Tucker function in $\mathcal{F}_{r}$. Assume that the target function satisfies the following tensor source condition, which is a standard smoothness assumption in spectral regularization and kernel learning theory [5, 7, 36]:

Using the spectral notation defined above, it holds that $g=\sum_{m}g_{m}\Phi_{m}$ and $u=\sum_{m}\Lambda_{m}^{v}g_{m}\Phi_{m}$. Let $P_{\mathcal{R}_{r}^{d}}$ denote the tensor-product projection onto the span of bases indexed by $\mathcal{R}_{r}$. By the RKHS norm, the truncation error can be given by

The exponent $(2v-1)$ appears because the RKHS norm divides the squared coefficient by the corresponding kernel eigenvalue.

Assume further that each mode has polynomial eigenvalue decay $\lambda_{j,m}\leq K_{0}m^{-\beta}$, a standard capacity condition for kernel integral operators [5, 36]. Let $\gamma=\beta(2v-1)$. If $\gamma>1$, then the one-dimensional spectral tail satisfies $\sum_{m>r}\lambda_{j,m}^{2v-1}\lesssim r^{1-\gamma}$, while the corresponding low-frequency sum is uniformly bounded in $r$. Therefore, the truncation error satisfies

or equivalently,

Here $C_{A}$ absorbs the source radius $C_{g}$, the spectral decay constant $K_{0}$, the tensor order $d$, and the bounded low-frequency spectral sums; in particular, $C_{A}$ is independent of $r$. Under the continuous embedding $\|h\|_{L^{2}(\rho)}\leq C_{\mathrm{emb}}\|h\|_{\mathcal{H}}$, the same rate gives

where $C_{A}$ additionally absorbs $C_{\mathrm{emb}}$.

Statistical error. For the statistical component, define the truncated tensor effective dimension as a finite-dimensional analogue of the effective dimension used in kernel ridge regression and regularized least squares [5, 36]:

Here $\lambda>0$ is the regularization parameter. It determines the learnability of each tensor spectral direction: directions with $\prod_{j=1}^{d}\lambda_{j,m_{j}}\gg\lambda$ are effectively retained, while directions with $\prod_{j=1}^{d}\lambda_{j,m_{j}}\ll\lambda$ are suppressed by regularization. Thus, $d_{\mathrm{eff},r}^{\mathrm{tensor}}(\lambda)$ measures the effective number of learnable tensor spectral directions within the truncated RKHS class. Equivalently, the same statistical term can be obtained through localized Rademacher complexity for squared-loss regression under boundedness and local curvature conditions [2, 3].

Under bounded-kernel, sub-Gaussian noise, and local curvature assumptions, standard kernel regression and localized Rademacher arguments yield that, with probability at least $1-\delta$ [3, 5],

Since $d_{\mathrm{eff},r}^{\mathrm{tensor}}(\lambda)\le r^d$, a conservative rank-based bound is

Balancing reconstruction rate. Combining the approximation and statistical bounds, we yield

Ignoring the lower-order concentration term, the optimal truncation level is obtained by balancing $C_A r^{-s}\simeq C_S r^{d/2}n^{-1/2}$, which gives

Substituting $r_\star$ into the total error, the right-hand side boils down to $Cn^{-\frac{s}{2s+d}}$.

Since $2s=\beta(2v-1)-1$, we obtain

This rate reflects two effects: the kernel spectrum controls how fast the rectangular truncation error decays, while the retained tensor-product spectral dimension controls the statistical term.

Sample complexity. Finally, the sample complexity can be inferred from the balanced rate. To ensure that $\|\hat{u}_r-u\|_{L^2(\rho)}\le\epsilon$ with probability at least $1-\delta$, it suffices to choose

Specialization to the Matérn-3/2 kernel. We now adopt the preceding bound to the Matérn-3/2 kernel. For a one-dimensional mode, the Matérn kernel with smoothness parameter $\nu_M=3/2$ has polynomial spectral decay $\lambda_{j,m}\lesssim m^{-\beta}$, with $\beta=2\nu_M+1=4$, under the standard one-dimensional Sobolev/Matérn spectral correspondence [28, 34]. The value of $\beta$ is determined by the kernel spectrum, whereas the source exponent $v$ describes the regularity of the target field relative to the kernel integral operator.

If the target field satisfies the second-order source condition $u=T^2g$ with $\|g\|_{L^2(\rho)}\le C_g$, then $v=2$, so $\gamma=\beta(2v-1)=12$ and $s=(\gamma-1)/2=11/2$. The approximation and total reconstruction bounds become

Balancing $r^{-11/2}\simeq r^{d/2}n^{-1/2}$ gives $r_\star\simeq n^{1/(11+d)}$. Substituting this choice yields

This rate shows that, under a Matérn-3/2 kernel and a second-order source condition, the approximation term decays rapidly with the truncation level, while the statistical term still depends on the effective tensor dimension induced by the retained modes.

Dataset Information

The experiments involve three synthetic benchmarks and a gridded ocean reanalysis. Table 9 summarizes the tensors from the datasets.

Dataset	Split	Stored tensor	Channels	Domain and time span
AC2D	160/20/20	$200 \times201 \times256 \times256$	scalar $u$	$[0,1)^2$, $t \in[0,2]$
NS2D	640/80/80	$800 \times201 \times128 \times128 \times3$	$(u,v,\omega)$	$[0,1)^2$, $t \in[0,20]$
Heat3D	160/20/20	$200 \times50 \times64 \times64 \times64 \times1$	temperature $u$	$[0,1]^3$, $t \in[0,4]$
Ocean3D	80/10/10	$100 \times90 \times32 \times64 \times64 \times2$	$(u,v)$	depth–latitude–longitude, 90 days

AC2D and NS2D. AC2D solves the periodic two-dimensional Allen–Cahn equation with $\lambda= 1$ and $\epsilon= 0.01$ and saves 201 frames at interval 0.01. NS2D solves the forced periodic vorticity equation with viscosity $10^{-4}$, saves 201 frames at interval 0.1, and stores $(u,v,\omega)$ after downsampling to $128^2$. The AC2D and NS2D training and held-out equation-residual evaluators use centered differences in both space and time and discard the first and last time points.

Heat3D. Heat3D contains 200 independent scalar-temperature trajectories. The trajectories are permuted with seed 123 and split 160/20/20. All trajectories share 50 physical save times in $[0,4]$; these times are slightly nonuniform because the explicit-Euler step and the save interval are rounded, while the logical model coordinate is a uniform 50-point grid in $[0,1]$. All $64^3$ voxels are valid, the spatial grid includes both endpoints, and the generator imposes homogeneous Dirichlet conditions on all six faces.

Ocean3D. Ocean3D is constructed from daily-mean horizontal currents $(u,v)$ in the CMEMS product GLOBAL_MULTIYEAR_PHY_001_030. The local source consists of annual NetCDF files, starting with uv_fulldepth_1993-01-01_to_1993-12-31.nc. Files are ordered chronologically, and each file is divided into nonoverlapping 90-day sequences beginning at its first record. Trailing records that do not form a complete sequence are discarded; sequences never cross file boundaries. We retain the first 100 sequences, comprising 9000 daily fields from 1993 to 2017 and split them sequentially into 80 training, 10 validation, and 10 test sequences. The first and last retained dates of these splits are 1993-01-01 to 2012-12-25, 2013-01-01 to 2015-06-29, and 2015-06-30 to 2017-12-26, respectively; these endpoint ranges do not imply uninterrupted daily coverage because 5 or 6 trailing days are omitted in each year. We directly retain the source-grid crop of $64 \times64$ horizontal points spanning latitudes $[-2.6666667461, 2.5833332539]$ and longitudes $[167.33333282471, 172.5833282471]$ degrees, without horizontal interpolation. We select zero-based depth indices 0 to 31; these 32 nonuniform levels range from 0.494024992 to 541.088928 m.Each sequence has shape $90 \times32 \times64 \times64 \times2$.

Evaluation subsets

For trajectory $b$, let $\mathcal{X}_b$ be its evaluation grid, $\mathcal{X}^{(b)}_{\mathrm{sen}} \subset\mathcal{X}_b$ its fixed set of sensor locations, and $\mathcal{T}^{(b,p)}_{\mathrm{sup}}$ and $\mathcal{T}^{(b,p)}_{\mathrm{miss}}$ the indices of support frames and missing frames under protocol $p$, respectively. The four disjoint subsets are

Support frames provide sparse measurements only at sensor locations, not complete spatial fields. Missing frames provide no sensor measurements to the model; their ground-truth fields, including values at sensor locations, are retained for evaluation. Let $C_D$ denote the number of channels. For each evaluation subset, we pool squared errors across trajectories, selected space–time points, and channels, divide by the total number of selected scalar values, and take the square root. An example is shown in Figure 3.

Missing-frame protocols

Let $T$ be the max timestamp of the trajectory. Let $M_{t}^{(b)}$ be the spatial mask fixed within the trajectory, and let $K_{t}^{(b,p)}\in\{0,1\}$ be a frame gate specific to the trajectory. The observation mask is $\widetilde{M}_{t}^{(b,p)}=K_{t}^{(b,p)}M_{t}^{(b)}$; when $K_{t}=1$, the frame is a support frame with sparse sensor measurements, and when $K_{t}=0$, it is a missing frame with all sensor measurements withheld across depths and channels. This mask affects the input of the model, not the availability of ground truth for evaluation.

Random-drop-half. Both endpoints are retained, while internal frames are sampled independently:

The implementation uses numpy.random.RandomState and split seed $10^{5}\cdot123+10g(s)+i_{\mathrm{split}}$, where $g(0.01)=1$, $g(0.03)=3$, $g(0.15)=15$, and $i_{\mathrm{split}}\in\{0,1,2\}$.

Temporal-interval-missing. A contiguous internal block is sampled independently for each trajectory:

and $K_{t}^{(b)}=0$ for $t\in[a_b,a_b+L_{\mathrm{miss}})$. Both endpoints remain support frames. An example is shown in Figure 4.

Reconstruction at non-grid-aligned times

This AC2D/NS2D ablation is separate from the native-grid missing-frame protocols. Its evaluation timestamps are non-grid-aligned relative to the original coarse time grid and lie within the interval spanned by the timestamps of support frames. They are not defined as missing frames from the original stored sequence. Starting from 201 coarse times $t_i^{(0)} = i\Delta t_c$, both endpoints are fixed and every internal time is perturbed by

AC2D uses $\Delta t_c = 0.01$ on $[0,2]$, and NS2D uses $\Delta t_c = 0.1$ on $[0,20]$. 100 evaluation indices are sampled without replacement from $\{1,\ldots,199\}$:

Sensor measurements are provided only at the timestamps indexed by $I_s$; the timestamps indexed by $I_e$ are used for reconstruction evaluation without providing their sensor measurements to the model. The endpoints therefore always belong to the support set. AC2D and NS2D use the same normalized offset vector and index-level partition between support frames and evaluation snapshots at non-grid-aligned times, shared across all test trajectories and sensors. Ground-truth states are saved directly by rerunning each PDE solver at the perturbed timestamps, not by interpolating uniformly saved snapshots. An example is shown in Figure 5.

Complete reconstruction results across sensor ratios

Tables 10–13 report the complete three-ratio four-subset comparison. Rankings, when discussed, use only available values within the same dataset, temporal protocol, spatial ratio, and subset.

Implementation Details

The spatial decoder uses a Matérn-3/2 kernel. Each active spatial mode has its own factor matrix; two-dimensional data use two modes, while Heat3D and Ocean3D use three. The number of inducing points along each mode equals that mode’s grid side length in the reported default configuration; the reported AC2D sensitivity baseline therefore uses 256 inducing points. Spatial factors are shared across channels, with the channel index stored in the Tucker core.

The spatial stage uses batch size 8 and is trained for 500 epochs, where one epoch traverses the complete training loader. For each mini-batch, the conditionally linear cores are updated by conjugate gradients and the RKHS factors are updated with AdamW. The learning rate is $2 \times10^{-4}$, the RKHS, core-ridge, core-TV, and temporal-GP weights are $10^{-4}$, $10^{-2}$, $10^{-2}$, and $10^{-1}$, respectively. Core updates use five CG iterations and tolerance $10^{-3}$.

Spatial checkpoints are selected by the lowest dense-field RMSE on the training trajectories after inverse normalization, averaging the RMSE values of evaluation mini-batches. Dense training fields are used for spatial checkpoint selection; spatial reconstruction gradients use only available sensor measurements, together with the stated regularizers and equation constraints.

The complete model samples equation-collocation centers from off-sensor locations on the spatial grid. For each trajectory, the set of sensor locations is the union of spatial locations with sensor measurements at any saved time and in any channel; the eligible collocation set is its spatial complement. AC2D and NS2D use two-dimensional grid centers, while Heat3D uses three-dimensional voxel centers. For Heat3D, centers whose finite-difference stencils would require moving the center inward from a boundary are excluded before sampling. The off-sensor restriction applies to collocation centers; neighboring stencil evaluations may lie at sensor locations. Ocean3D uses no governing-equation loss in the reported setting. AC2D and NS2D use centered spatial and temporal differences and exclude the temporal endpoints from residual training and evaluation.

Before each AdamW factor update, the gradient of each RKHS coefficient matrix is left-multiplied by the inverse Gram matrix of its spatial-mode kernel, $K_d^{-1}\nabla_{A_d}\mathcal{L}$. The spatial AdamW optimizer uses weight decay $10^{-2}$; the gradient-only core updater uses the same weight decay.

We compare block-update strategies with their corresponding effective objectives. The complete model alternates CG core solves with gradient-based RKHS-factor updates. The solve-only variant retains the same ridge-TV core subproblem, but updates each spatial factor by solving for target factor values with a value-ridge penalty and projecting them to RKHS coefficients with kernel ridge regression. These factor updates use observational least squares and kernel projection. The gradient-only variant updates both cores and factors with AdamW using the RMSE over available sensor measurements, RKHS penalty, equation-based loss, and time-kernel penalty. The RKHS penalty depends only on the factors, whereas the time-kernel penalty depends only on the cores. Core-ridge and temporal-TV penalties belong to the CG subproblem. The time-kernel term is defined in Eq. (10).

Fields are normalized with channel-wise train-set statistics. Validation and test splits reuse these statistics, and all reported RMSEs, residuals, and physical diagnostics are computed after inverse normalization. The Ocean3D RecFNO RMSE export originally aggregates both channels jointly; the values reported in Table 1 use the corresponding per-channel convention by dividing those four exported RMSE columns by $\sqrt{2}$.

At inference, the spatial decoder is frozen and the core-correction sequence is solved around the train-set mean-core reference on the complete native time grid with ridge $10^{-2}$, TV weight $10^{-2}$, 30 CG iterations, tolerance $10^{-4}$, and batch size 32. Missing frames have zero observation masks; after the joint solve, the cores at the timestamps of support frames are extracted as support anchors.

The residual ODE-RNN temporal head operates in Tucker-core space, uses hidden dimension 256, and is trained for 20 epochs with AdamW (batch size 4, learning rate $10^{-3}$, and weight decay $10^{-6}$, except that the learned residual-scale parameter has zero weight decay). Let $B_i(t)$ be the cubic Hermite baseline constructed from the inferred support cores, and let $C_\psi(t)$ denote the actual correction added to that baseline, including the interval gate and the learned bounded residual scale. At the timestamps of missing frames, we minimize the mean squared error between $C_\psi(t)$ and $G_i^{\mathrm{train}}(t)-B_i(t)$, where $G_i^{\mathrm{train}}$ is the frozen training-core sequence saved by the spatial stage. Equivalently, the objective is the missing-frame MSE between the predicted core $B_i(t)+C_\psi(t)$ and $G_i^{\mathrm{train}}(t)$. The core-loss weight is 1; decoded-field, spatial-gradient, and temporal-equation losses are disabled. We select the checkpoint with the lowest validation relative $\ell_2$ error over missing frames, computed from decoded fields after inverse normalization.

The MSE averages all scalar coefficients of the selected cores. Missing-frame validation uses all spatial locations at the selected times.

Protocol	$s$	Method	S-sup	Off-sup	S-miss	Off-miss
Random-drop-half	0.01	SDIFT	0.373381	0.381700	0.378977	0.381200
Random-drop-half	0.01	LRTFR	0.311380	0.329100	0.308521	0.328814
Random-drop-half	0.01	PhysFTM	0.009898	0.187900	0.042055	0.191673
Random-drop-half	0.01	RecFNO	<u>0.020782</u>	0.697782	0.021170	0.696267
Random-drop-half	0.01	SIREN-MLP	0.351270	0.350988	0.349308	0.349007
Random-drop-half	0.01	PINNs	0.354653	0.354871	0.350601	0.350840
Random-drop-half	0.03	SDIFT	0.358035	0.380800	0.380301	0.381514
Random-drop-half	0.03	LRTFR	0.311700	<u>0.322200</u>	0.313404	<u>0.320383</u>
Random-drop-half	0.03	PhysFTM	0.017354	0.087100	0.050349	0.099881
Random-drop-half	0.03	RecFNO	<u>0.018483</u>	0.452826	0.021782	0.455233
Random-drop-half	0.03	SIREN-MLP	0.348683	0.349317	0.350125	0.350761
Random-drop-half	0.03	PINNs	0.351663	0.352425	0.352679	0.353309
Random-drop-half	0.15	SDIFT	0.294393	0.374100	0.371273	0.371736
Random-drop-half	0.15	LRTFR	0.320074	0.317200	0.319075	0.320370
Random-drop-half	0.15	PhysFTM	<u>0.029960</u>	0.035900	<u>0.057787</u>	0.062454
Random-drop-half	0.15	RecFNO	0.011095	0.230709	0.011986	0.231022
Random-drop-half	0.15	SIREN-MLP	0.350102	<u>0.350055</u>	0.350040	0.349975
Random-drop-half	0.15	PINNs	0.353934	0.353904	0.351765	0.351766
Temporal-interval-missing	0.01	SDIFT	0.372720	0.379670	0.372029	0.373811
Temporal-interval-missing	0.01	LRTFR	0.309972	<u>0.330427</u>	0.306973	<u>0.326716</u>
Temporal-interval-missing	0.01	PhysFTM	0.006558	0.189538	0.032476	0.183397
Temporal-interval-missing	0.01	RecFNO	<u>0.020743</u>	0.701355	0.028592	0.681053
Temporal-interval-missing	0.01	SIREN-MLP	0.351543	0.351296	0.342732	0.342199
Temporal-interval-missing	0.01	PINNs	0.353767	0.354026	0.347975	0.348079
Temporal-interval-missing	0.03	SDIFT	0.356991	0.361986	0.390707	0.392586
Temporal-interval-missing	0.03	LRTFR	0.310055	<u>0.316757</u>	0.320590	<u>0.327806</u>
Temporal-interval-missing	0.03	PhysFTM	0.016328	0.093969	<u>0.039527</u>	0.079239
Temporal-interval-missing	0.03	RecFNO	<u>0.018107</u>	0.458752	0.021163	0.436798
Temporal-interval-missing	0.03	SIREN-MLP	0.346053	0.346642	0.358185	0.358978
Temporal-interval-missing	0.03	PINNs	0.348742	0.349428	0.365077	0.366217
Temporal-interval-missing	0.15	SDIFT	0.294740	0.299284	0.349127	0.350140
Temporal-interval-missing	0.15	LRTFR	0.320387	<u>0.321854</u>	0.305895	<u>0.307641</u>
Temporal-interval-missing	0.15	PhysFTM	<u>0.031414</u>	0.040888	<u>0.045196</u>	0.048007
Temporal-interval-missing	0.15	RecFNO	0.011196	0.234129	0.012961	0.218094
Temporal-interval-missing	0.15	SIREN-MLP	0.353258	0.353191	0.332458	0.332446
Temporal-interval-missing	0.15	PINNs	0.356292	0.356274	0.338695	0.338699

Protocol	$s$	Method	S-sup	Off-sup	S-miss	Off-miss
Random-drop-half	0.01	SDIFT	0.518934	0.512800	0.510791	0.509261
		LRTFR	0.123841	0.172000	0.126188	0.170679
		PhysFTM	0.007849	0.137700	0.030927	0.130659
		RecFNO	0.261520	0.418079	0.262461	0.413223
		SIREN-MLP	0.244818	0.245620	0.243905	0.244360
		PINNs	0.338968	0.341335	0.339094	0.341072
	0.03	SDIFT	0.528668	0.509700	0.510029	0.510274
		LRTFR	0.132450	0.147200	0.134470	0.145757
		PhysFTM	0.006687	0.062100	0.032396	0.064976
		RecFNO	<u>0.068098</u>	<u>0.097567</u>	<u>0.078196</u>	<u>0.104413</u>
		SIREN-MLP	0.245551	0.245255	0.244651	0.244405
		PINNs	0.342539	0.341771	0.341006	0.340554
	0.15	SDIFT	0.609226	0.510500	0.510615	0.510589
		LRTFR	0.135290	0.139400	0.136582	0.138756
		PhysFTM	0.005888	<u>0.020300</u>	0.032096	0.033357
		RecFNO	<u>0.007173</u>	0.007929	<u>0.039964</u>	<u>0.040051</u>
		SIREN-MLP	0.245407	0.245486	0.244408	0.244476
		PINNs	0.342011	0.342206	0.340038	0.340193
Temporal-interval-missing	0.01	SDIFT	0.522547	0.520767	0.494061	0.492242
		LRTFR	0.127082	0.173888	0.179695	0.206649
		PhysFTM	0.004291	0.141533	0.136380	0.164821
		RecFNO	0.267930	0.428534	0.405006	0.478397
		SIREN-MLP	0.249049	0.249814	0.206410	0.206503
		PINNs	0.346718	0.349269	0.306146	0.306602
	0.03	SDIFT	0.530194	0.530518	0.504341	0.502748
		LRTFR	0.135148	0.146821	0.190250	0.196613
		PhysFTM	0.004346	0.062824	0.149493	0.150558
		RecFNO	<u>0.069056</u>	<u>0.099100</u>	<u>0.358033</u>	<u>0.359434</u>
		SIREN-MLP	0.247968	0.247773	0.213687	0.213537
		PINNs	0.347219	0.346569	0.318917	0.318479
	0.15	SDIFT	0.610478	0.609820	0.485120	0.485171
		LRTFR	0.138750	0.140954	<u>0.170839</u>	<u>0.172046</u>
		PhysFTM	0.005690	<u>0.012067</u>	0.134833	0.134567
		RecFNO	<u>0.007474</u>	0.008253	0.318252	0.318191
		SIREN-MLP	0.250208	0.250287	0.203004	0.203190
		PINNs	0.350376	0.350661	0.300498	0.300163

Table 12: Complete four-subset reconstruction RMSE on Heat3D across spatial sensor ratios. Lower is better. (Bold: best, <u>Underline</u>: second best)

Protocol	$s$	Method	S-sup	Off-sup	S-miss	Off-miss
Random-drop-half	0.01	SDIFT	0.058457	0.058594	0.061825	0.062284
		LRTFR	<u>0.010659</u>	<u>0.010815</u>	<u>0.012014</u>	<u>0.012549</u>
		PhysFTM	0.001344	0.001907	0.006379	0.006824
		RecFNO	0.339625	0.888925	0.307609	0.813064
		SIREN-MLP	0.040199	0.040240	0.042215	0.042397
		PINNs	0.040309	0.040409	0.041993	0.042295
	0.03	SDIFT	0.058111	0.058081	0.063030	0.062716
		LRTFR	<u>0.012121</u>	<u>0.012216</u>	<u>0.013093</u>	<u>0.013146</u>
		PhysFTM	0.001083	0.001250	0.005745	0.005744
		RecFNO	0.287363	0.778472	0.263722	0.720340
		SIREN-MLP	0.040803	0.040865	0.042108	0.042053
		PINNs	0.040002	0.040113	0.042714	0.042640
	0.15	SDIFT	0.059938	0.059993	0.060820	0.060836
		LRTFR	<u>0.012518</u>	<u>0.012578</u>	<u>0.011607</u>	<u>0.011652</u>
		PhysFTM	0.001015	0.001040	0.005282	0.005295
		RecFNO	0.171750	0.322953	0.158704	0.258328
		SIREN-MLP	0.041383	0.041470	0.041042	0.041105
		PINNs	0.041778	0.041929	0.040826	0.040917
Temporal-interval-missing	0.01	SDIFT	0.059872	0.060202	0.060550	0.060789
		LRTFR	<u>0.011883</u>	<u>0.012311</u>	<u>0.012546</u>	<u>0.012876</u>
		PhysFTM	0.001293	0.002132	0.007193	0.007070
		RecFNO	0.320509	0.837414	0.321806	0.904206
		SIREN-MLP	0.040678	0.040735	0.038204	0.038412
		PINNs	0.040636	0.040735	0.043391	0.044007
	0.03	SDIFT	0.062840	0.062680	0.050158	0.049889
		LRTFR	<u>0.012100</u>	<u>0.012250</u>	<u>0.005766</u>	<u>0.005767</u>
		PhysFTM	0.001218	0.001372	0.003103	0.003046
		RecFNO	0.286290	0.775257	0.225969	0.641868
		SIREN-MLP	0.043423	0.043440	0.027492	0.027421
		PINNs	0.043136	0.043190	0.033532	0.033341
	0.15	SDIFT	0.063097	0.063136	0.047558	0.047574
		LRTFR	<u>0.011794</u>	<u>0.011881</u>	<u>0.009763</u>	<u>0.009595</u>
		PhysFTM	0.001105	0.001136	0.003756	0.003682
		RecFNO	0.175643	0.317359	0.122943	0.169169
		SIREN-MLP	0.043686	0.043784	0.025485	0.025478
		PINNs	0.044045	0.044195	0.027758	0.027711

Protocol	$s$	Method	S-sup	Off-sup	S-miss	Off-miss
Random-drop-half	0.01	SDIFT	0.281171	0.278458	0.280994	0.280682
Random-drop-half	0.01	LRTFR	0.113447	0.126022	0.122721	0.133850
Random-drop-half	0.01	PhysFTM	0.006599	0.115414	0.076277	0.125325
Random-drop-half	0.01	RecFNO	0.758963	1.005671	0.793739	1.002127
Random-drop-half	0.01	SIREN-MLP	0.200624	0.203220	0.196868	0.200081
Random-drop-half	0.01	PINNs	0.297127	0.301674	0.295973	0.301402
Random-drop-half	0.03	SDIFT	0.278627	0.280655	0.275848	0.278376
Random-drop-half	0.03	LRTFR	0.110377	0.114503	0.117666	0.121548
Random-drop-half	0.03	PhysFTM	0.006319	0.065738	0.071328	0.088581
Random-drop-half	0.03	RecFNO	0.738273	0.970741	0.776879	0.963019
Random-drop-half	0.03	SIREN-MLP	0.201582	0.201681	0.200241	0.200943
Random-drop-half	0.03	PINNs	0.301625	0.300997	0.302323	0.301937
Random-drop-half	0.15	SDIFT	0.275481	0.276918	0.281502	0.282907
Random-drop-half	0.15	LRTFR	0.109859	0.110456	0.127294	0.127726
Random-drop-half	0.15	PhysFTM	0.007103	0.017565	0.080431	0.080677
Random-drop-half	0.15	RecFNO	0.711010	0.854068	0.755037	0.861364
Random-drop-half	0.15	SIREN-MLP	0.199900	0.200691	0.201670	0.202316
Random-drop-half	0.15	PINNs	0.298829	0.300486	0.301210	0.302968
Temporal-interval-missing	0.01	SDIFT	0.283314	0.282061	0.272620	0.270666
Temporal-interval-missing	0.01	LRTFR	0.104156	0.119385	0.171227	0.182007
Temporal-interval-missing	0.01	PhysFTM	0.005952	0.115920	0.191076	0.194097
Temporal-interval-missing	0.01	RecFNO	0.763253	1.007166	0.944159	0.992394
Temporal-interval-missing	0.01	SIREN-MLP	0.201337	0.204394	0.182836	0.184601
Temporal-interval-missing	0.01	PINNs	0.299977	0.306116	0.282315	0.282459
Temporal-interval-missing	0.03	SDIFT	0.274940	0.277254	0.287116	0.289410
Temporal-interval-missing	0.03	LRTFR	0.105062	0.108289	0.180623	0.183770
Temporal-interval-missing	0.03	PhysFTM	0.005946	0.063821	0.189849	0.185648
Temporal-interval-missing	0.03	RecFNO	0.718410	0.943793	0.991909	1.053829
Temporal-interval-missing	0.03	SIREN-MLP	0.197733	0.198327	0.211994	0.211505
Temporal-interval-missing	0.03	PINNs	0.297650	0.297391	0.318698	0.317248
Temporal-interval-missing	0.15	SDIFT	0.277062	0.278427	0.284215	0.285884
Temporal-interval-missing	0.15	LRTFR	0.110311	0.110758	0.190074	0.190728
Temporal-interval-missing	0.15	PhysFTM	0.007105	0.017647	0.159473	0.158267
Temporal-interval-missing	0.15	RecFNO	0.712010	0.855126	0.930787	0.906187
Temporal-interval-missing	0.15	SIREN-MLP	0.202733	0.203547	0.189757	0.190244
Temporal-interval-missing	0.15	PINNs	0.300200	0.301877	0.299327	0.301165