Introduction
Background and Motivation
High-dimensional partial differential equations (PDEs) arise naturally in many scientific and engineering problems, including uncertainty quantification, kinetic equations, stochastic dynamics, molecular simulation, and high-dimensional control [4 , 12 ] . In these settings, the solution often depends on a large number of spatial, parametric, or latent variables, making accurate numerical approximation substantially more challenging than in low-dimensional regimes [7 , 13 ] . Developing scalable solvers for such problems is therefore important not only for extending classical numerical methods to higher dimensions, but also for enabling reliable simulation and prediction in systems whose intrinsic state spaces are inherently high-dimensional. However, the computational cost of conventional discretization-based methods typically grows rapidly with dimension, making the efficient representation and solution of high-dimensional PDEs a persistent challenge in scientific computing [20 ] .
The above difficulty suggests that the scalability of high-dimensional PDE solvers depends critically on how the solution is represented. Spectral coordinates provide a natural description for smooth solutions and expose differential operators through explicit algebraic structure [3 , 36 ] . However, the resulting multidimensional coefficient field is itself a high-order tensor, whose dense representation still grows exponentially with dimension [2 ] . A large ambient tensor does not necessarily imply equally high intrinsic complexity, since spectral coefficients may exhibit strong cross-dimensional separability and admit a compact low-rank representation [1 ] . This observation motivates a representation-centric question: Can high-dimensional PDEs be solved directly in spectral space while replacing the exponentially large coefficient tensor with a structured low-rank representation?
Related Work
Physics-informed neural PDE solvers. Physics-informed neural PDE solvers have evolved well beyond the original PINN formulation. Strong-form methods such as PINN [32 ] and DGM [34 ] minimize pointwise PDE residuals, while DRM [8 ] , VPINN [18 ] , and h p hp h p [19 ] introduce energy- or weak-form objectives; cPINN [15 ] and XPINNs [14 ] further exploit domain decomposition. Recent work has also improved neural representations through adaptive activations, trainable Gaussian features such as PIG [17 ] , and residual-adaptive architectures such as PirateNets [37 ] . Nevertheless, optimization remains challenging because differential operators can induce ill-conditioned objectives and imbalance among PDE, initial, and boundary losses [38 ] , motivating gradient-based balancing [23 ] , NTK-guided weighting [39 ] , and self-adaptive weighting [27 ] . Despite these advances in formulation, architecture, and optimization, most methods still rely on generic coordinate-to-solution neural representations whose internal structure is learned implicitly, leaving the explicit and dimension-scalable organization of high-dimensional PDE solutions largely unresolved.
Spectral methods and neural spectral solvers. Spectral methods approximate PDE solutions using global orthogonal bases, with Fourier bases naturally suited to periodic problems and Jacobi families, including Legendre and Chebyshev polynomials, providing flexible representations for non-periodic domains [3 , 33 , 36 ] . For sufficiently smooth solutions, they offer rapid convergence and structured algebraic representations of differential operators in coefficient space [9 ] . The spectral bias of neural networks, which favors low-frequency components over high-frequency and multiscale structures [31 , 35 ] , has further motivated frequency-aware neural PDE solvers. Representative developments include FNN [28 ] , FNO [24 ] , and PINO [25 ] , followed by spectral operator frameworks such as OPNO [26 ] and SCLON [6 ] . More recently, SINN formulate PDE solving directly in coefficient space using algebraic spectral differentiation [43 ] , with subsequent extensions incorporating coefficient scaling and basis-dependent embeddings for high-dimensional problems [42 ] . Despite these advances, the multidimensional spectral coefficient field itself is generally not endowed with explicit cross-dimensional low-rank structure, leaving its representation complexity in high dimensions only partially addressed.
Tensor decompositions and separable representations. Low-rank tensor methods, including canonical polyadic (CP), Tucker, and tensor-train decompositions, reduce high-dimensional complexity by exploiting separability across coordinates [10 , 11 , 22 , 29 , 30 ] . This principle has also been adopted in neural PDE solvers. SPINN [5 ] uses coordinate-wise subnetworks, while TNN [41 ] represents multivariate solutions through tensor products of one-dimensional neural factors. More closely related, LoRFS [16 ] employs low-rank separable Fourier expansions to enable closed-form, sampling-free evaluation of physics-informed objectives and gradients. However, a compact representation alone does not guarantee efficient PDE computation, as differential operators, nonlinear interactions, and limited representation capacity introduce additional challenges in preserving low-rank structure and solution accuracy. This motivates a tighter integration of tensor decomposition with spectral discretization, where the algebraic structure of spectral operators enables factorwise computation, while low-rank representations control the dimensional complexity of the coefficient space.
Main Contributions
Motivated by the above difficulties, we propose TensorSpectra , whose framework could see Figure 1 , a low-rank spectral framework for scalable high-dimensional PDE solving. Our central perspective is that the main obstacle in high-dimensional spectral computation is not necessarily the spectral representation itself, but the exponentially large dense coefficient tensor induced by conventional tensor-product discretization. TensorSpectra therefore represents the multidimensional spectral coefficient field directly in a structured low-rank form and performs PDE computation within this compressed coefficient space, rather than treating spectral coefficients as independent quantities or reconstructing the full tensor. This formulation brings together the high-accuracy and algebraic operator structure of spectral methods with the dimensional scalability of low-rank tensor representations, while retaining a unified framework for representation, operator evaluation, and progressive refinement. The main contributions are summarized as follows:
Figure 1. Tensorspectra framework.
Structure-preserving spectral operator calculus. We apply differential operators directly to one-dimensional spectral factors, avoiding dense tensor reconstruction and repeated spatial automatic differentiation. The same structured formulation extends to nonlinear PDEs through de-aliased spectral interactions.
Spectral–rank convergence analysis. We establish an error framework separating spectral truncation, low-rank approximation, and optimization errors, together with stability and convergence results that characterize the roles of spectral resolution, tensor rank, and progressive correction.
Organization
The remainder of this paper is organized as follows: Section 2 introduces the PDE setting and spectral notation. Section 3 presents Tensorspectra, and Section 4 analyzes its approximation and residual bounds. Section 5 and 6 report numerical results and ablation studies respectively. Section 7 closes the paper. Besides, the appendix provides technical details and proofs.
Preliminaries
In this section, we first introduce the general PDE formulation considered in this work and then establish the spectral representation and notation used throughout the paper, providing the necessary foundation for the Tensorspectra framework developed in the next section.
Problem Setup
We consider a general class of high-dimensional partial differential equations defined on a spatial domain Ω ⊂ R d \Omega\subset\mathbb{R}^d Ω ⊂ R d [ 0 , T ] [0,T] [ 0 , T ] u ( x , t ) u(x,t) u ( x , t )
N [ u ] ( x , t ) = f ( x , t ) , ( x , t ) ∈ Ω × [ 0 , T ] , u ( x , t ) = g ( x , t ) , ( x , t ) ∈ ∂ Ω × [ 0 , T ] , u ( x , 0 ) = h ( x ) , x ∈ Ω , (1) \begin{aligned}
\mathcal{N}[u](x,t) &= f(x,t), \quad(x,t) \in\Omega\times[0,T],\\
u(x,t) &= g(x,t), \quad(x,t) \in\partial\Omega\times[0,T],\\
u(x,0) &= h(x), \quad x \in\Omega, \tag*{(1)}
\end{aligned} N [ u ] ( x , t ) u ( x , t ) u ( x , 0 ) = f ( x , t ) , ( x , t ) ∈ Ω × [ 0 , T ] , = g ( x , t ) , ( x , t ) ∈ ∂ Ω × [ 0 , T ] , = h ( x ) , x ∈ Ω , (1) where N \mathcal{N} N f f f g g g h h h
For stationary problems, the temporal dependence and initial condition are omitted, and the governing equation reduces to N [ u ] ( x ) = f ( x ) \mathcal{N}[u](x) = f(x) N [ u ] ( x ) = f ( x ) u u u x x x d d d
To obtain a representation that is compatible with high-dimensional computation, we express the solution in a tensor-product spectral basis rather than parameterizing it directly in physical space. This transforms the original PDE into an equivalent problem for the corresponding multidimensional spectral coefficients, where differentiation and other PDE operators can be treated through their algebraic action in coefficient space. We introduce this spectral representation and the associated notation in the next subsection.
Spectral Representation
For each coordinate direction j = 1 , … , d j=1,\ldots,d j = 1 , … , d { ϕ k j ( j ) } \{\phi_{k_j}^{(j)}\} { ϕ k j ( j ) } x j x_j x j Φ k ( x ) = ∏ j = 1 d ϕ k j ( j ) ( x j ) \Phi_{\mathbf{k}}(x)=\prod_{j=1}^{d}\phi_{k_j}^{(j)}(x_j) Φ k ( x ) = ∏ j = 1 d ϕ k j ( j ) ( x j ) k = ( k 1 , … , k d ) \mathbf{k}=(k_1,\ldots,k_d) k = ( k 1 , … , k d )
u ( x , t ) ≈ ∑ k ∈ K u ^ k ( t ) Φ k ( x ) = ∑ k ∈ K u ^ k ( t ) ∏ j = 1 d ϕ k j ( j ) ( x j ) , (2) u(x,t)\approx\sum_{\mathbf{k}\in\mathcal{K}}\widehat{u}_{\mathbf{k}}(t)\Phi_{\mathbf{k}}(x)=\sum_{\mathbf{k}\in\mathcal{K}}\widehat{u}_{\mathbf{k}}(t)\prod_{j=1}^{d}\phi_{k_j}^{(j)}(x_j), \tag*{(2)} u ( x , t ) ≈ k ∈ K ∑ u k ( t ) Φ k ( x ) = k ∈ K ∑ u k ( t ) j = 1 ∏ d ϕ k j ( j ) ( x j ) , (2) where K \mathcal{K} K u ^ k ( t ) \widehat{u}_{\mathbf{k}}(t) u k ( t ) k \mathbf{k} k
Collecting the coefficients associated with all retained multidimensional modes yields the spectral coefficient tensor
U ^ = [ u ^ k 1 , … , k d ] ∈ C N 1 × ⋯ × N d , (3) \widehat{\mathbf{U}}=\left[\widehat{u}_{k_1,\ldots,k_d}\right]\in\mathbb{C}^{N_1\times\cdots\times N_d}, \tag*{(3)} U = [ u k 1 , … , k d ] ∈ C N 1 × ⋯ × N d , (3) where N j N_j N j j = 1 , … , d j=1,\ldots,d j = 1 , … , d j j j N N N U ^ ∈ C N × ⋯ × N \widehat{\mathbf{U}}\in\mathbb{C}^{N\times\cdots\times N} U ∈ C N × ⋯ × N U ^ \widehat{\mathbf{U}} U
Despite the approximation power of spectral methods, the multidimensional coefficient tensor becomes prohibitively large in high dimensions. If N N N U ^ \widehat{\mathbf{U}} U N d N^d N d d d d U ^ \widehat{\mathbf{U}} U
Methodology
In this section, we present TensorSpectra, a framework for solving high-dimensional PDEs by optimizing a structured spectral coefficient field. We introduce its CP representation with trainable temporal factors, factorwise spectral operator calculus, operator-aware residual objective, and additive low-rank refinement. Together, these components enable solution representation and residual evaluation through one-dimensional factors.
Low-Rank Spectral Tensor Representation
For time-dependent problems, we represent the spatial coefficient tensor for Equation (1) as
U ^ ( t ) = ∑ α = 1 r λ α T α ( t ) a α ( 1 ) ∘ ⋯ ∘ a α ( d ) , u ^ k ( t ) = ∑ α = 1 r λ α T α ( t ) ∏ j = 1 d a α , k j ( j ) . (4) \widehat{\mathbf{U}}(t)=\sum_{\alpha=1}^{r}\lambda_{\alpha}T_{\alpha}(t)\mathbf{a}_{\alpha}^{(1)}\circ\cdots\circ\mathbf{a}_{\alpha}^{(d)},\qquad\widehat{u}_{\mathbf{k}}(t)=\sum_{\alpha=1}^{r}\lambda_{\alpha}T_{\alpha}(t)\prod_{j=1}^{d}a_{\alpha,k_j}^{(j)}. \tag*{(4)} U ( t ) = α = 1 ∑ r λ α T α ( t ) a α ( 1 ) ∘ ⋯ ∘ a α ( d ) , u k ( t ) = α = 1 ∑ r λ α T α ( t ) j = 1 ∏ d a α , k j ( j ) . (4) Here, ∘ \circ ∘ r r r λ α \lambda_{\alpha} λ α a α ( j ) ∈ C N j \mathbf{a}_{\alpha}^{(j)}\in\mathbb{C}^{N_j} a α ( j ) ∈ C N j T α ( t ) T_{\alpha}(t) T α ( t ) T α = 1 T_{\alpha}=1 T α = 1 (4) , without assembling the complete tensor.
We expand each temporal factor in a finite basis,
T α ( t ) = ∑ q = 0 p t c q , α ψ q ( t ) . (5) T_{\alpha}(t) = \sum_{q=0}^{p_t} c_{q,\alpha}\psi_q(t). \tag*{(5)} T α ( t ) = q = 0 ∑ p t c q , α ψ q ( t ) . (5) Low-degree polynomials in normalized time, ψ q ( t ) = ( t / T ) q \psi_q(t) = (t/T)^q ψ q ( t ) = ( t / T ) q c α = ( c 0 , α , … , c p t , α ) T \mathbf{c}_{\alpha} = (c_{0,\alpha}, \ldots, c_{p_t,\alpha})^{\mathsf{T}} c α = ( c 0 , α , … , c p t , α ) T
u ^ s t = ∑ α = 1 r λ α a α ( 1 ) ∘ ⋯ ∘ a α ( d ) ∘ c α . (6) \hat{\mathbf{u}}_{\mathrm{st}} = \sum_{\alpha=1}^{r} \lambda_{\alpha}\mathbf{a}_{\alpha}^{(1)} \circ\cdots\circ\mathbf{a}_{\alpha}^{(d)} \circ\mathbf{c}_{\alpha}. \tag*{(6)} u ^ st = α = 1 ∑ r λ α a α ( 1 ) ∘ ⋯ ∘ a α ( d ) ∘ c α . (6) Equation (4) and Equation (6) describe the same separated approximation that contracting the final coefficient axis with the temporal basis recovers the spatial tensor at any time. The first form supports residual evaluation at sampled times, while the second permits projection of the residual onto a temporal spectral basis. Both forms use r ( 1 + ∑ j N j + p t + 1 ) r(1+\sum_j N_j+p_t+1) r ( 1 + ∑ j N j + p t + 1 ) § 8 .
Structure-Preserving Spectral Calculus
Let D j \mathbf{D}_j D j x j x_j x j
∂ x j u ^ ( t ) = ∑ α = 1 r λ α T α ( t ) a α ( 1 ) ∘ ⋯ ∘ ( D j a α ( j ) ) ∘ ⋯ ∘ a α ( d ) . (7) \partial_{x_j}\hat{\mathbf{u}}(t) = \sum_{\alpha=1}^{r} \lambda_{\alpha}T_{\alpha}(t)\mathbf{a}_{\alpha}^{(1)} \circ\cdots\circ(\mathbf{D}_j\mathbf{a}_{\alpha}^{(j)}) \circ\cdots\circ\mathbf{a}_{\alpha}^{(d)}. \tag*{(7)} ∂ x j u ^ ( t ) = α = 1 ∑ r λ α T α ( t ) a α ( 1 ) ∘ ⋯ ∘ ( D j a α ( j ) ) ∘ ⋯ ∘ a α ( d ) . (7) Higher derivatives apply the corresponding higher-order differentiation map to one factor, while mixed derivatives transform each participating factor. Thus, spatial and temporal derivatives can be combined by acting independently on their respective axes. Each separable operator term preserves the rank bound; sums of terms can increase it. For Fourier bases, D j = diag ( i k j ) \mathbf{D}_j = \operatorname{diag}(i k_j) D j = diag ( i k j ) Δ u ^ k = − ∥ k ∥ 2 2 u ^ k \widehat{\Delta u}_k = -\lVert k\rVert_2^2\hat{u}_k Δ u k = − ∥ k ∥ 2 2 u ^ k
Time differentiation acts on the temporal factors according to
∂ t s u ^ ( t ) = ∑ α = 1 r λ α T α ( s ) ( t ) a α ( 1 ) ∘ ⋯ ∘ a α ( d ) , s = 1 , 2. (8) \partial_t^s\hat{\mathbf{u}}(t) = \sum_{\alpha=1}^{r} \lambda_{\alpha}T_{\alpha}^{(s)}(t)\mathbf{a}_{\alpha}^{(1)} \circ\cdots\circ\mathbf{a}_{\alpha}^{(d)}, \qquad s=1,2. \tag*{(8)} ∂ t s u ^ ( t ) = α = 1 ∑ r λ α T α ( s ) ( t ) a α ( 1 ) ∘ ⋯ ∘ a α ( d ) , s = 1 , 2. (8) The derivatives are obtained analytically from Equation (5) , including the physical-time scaling. In the space–time coefficient form, the equivalent operation replaces c α \mathbf{c}_{\alpha} c α D t s c α \mathbf{D}_t^s\mathbf{c}_{\alpha} D t s c α D t \mathbf{D}_t D t [ 0 , T ] [0,T] [ 0 , T ]
Polynomial nonlinearities are evaluated by expanding interactions between CP components and projecting their one-dimensional factor products. Fourier products use padding and truncation; Jacobi products use over-integration. Space–time products also include the temporal factors. A product of two rank-one components remains separable, allowing each coordinate product to be evaluated and projected independently before recombining the factors. The enlarged Fourier bandwidth or polynomial quadrature order prevents unresolved product modes from aliasing into the retained coefficients. These operations avoid dense tensor assembly, although nonlinear interactions can increase the rank. De-aliasing and interaction counts are detailed in Appendix § 9 .
PDE-Aware Spectral Objective
Let R k ( t ) = N [ u ^ ] − f k ( t ) \mathcal{R}_k(t) = \mathcal{N}[\hat{u}] - f_k(t) R k ( t ) = N [ u ^ ] − f k ( t )
To compensate for operator-induced differences in residual scale, we apply a positive weight w k w_k w k L \mathcal{L} L w k = ρ ( σ ( L , k ) ) w_k = \rho(\sigma(\mathcal{L},k)) w k = ρ ( σ ( L , k )) ρ \rho ρ S S S
L P D E p r e c = 1 ∣ S ∣ ∑ ( k , t ) ∈ S ∣ w k R k ( t ) ∣ 2 . (9) \mathcal{L}_{\mathrm{PDE}}^{\mathrm{prec}}=\frac{1}{\lvert S\rvert}\sum_{(k,t)\in S}\left\lvert w_kR_k(t)\right\rvert^2. \tag*{(9)} L PDE prec = ∣ S ∣ 1 ( k , t ) ∈ S ∑ ∣ w k R k ( t ) ∣ 2 . (9) When time is projected onto a spectral basis, training instead uses the space–time residual coefficients in the objective
L P D E p r e c = 1 ∣ K ∣ ∣ Q ∣ ∑ k ∈ K ∑ q ∈ Q ∣ w k , q R k , q ∣ 2 . (10) \mathcal{L}_{\mathrm{PDE}}^{\mathrm{prec}}=\frac{1}{\lvert K\rvert\lvert Q\rvert}\sum_{k\in K}\sum_{q\in Q}\left\lvert w_{k,q}R_{k,q}\right\rvert^2. \tag*{(10)} L PDE prec = ∣ K ∣ ∣ Q ∣ 1 k ∈ K ∑ q ∈ Q ∑ ∣ w k , q R k , q ∣ 2 . (10) where Q Q Q R k , q R_{k,q} R k , q
For a PDE of temporal order s t s_t s t h ^ k ( ℓ ) \widehat{h}^{(\ell)}_k h k ( ℓ )
L I C = ∑ ℓ = 0 s t − 1 γ ℓ ∣ K 0 ∣ ∑ k ∈ K 0 ∣ ∂ t ℓ u ^ k ( 0 ) − h ^ k ( ℓ ) ∣ 2 , L = L P D E p r e c + λ I C L I C + L c o n . (11) \mathcal{L}_{\mathrm{IC}}=\sum_{\ell=0}^{s_t-1}\frac{\gamma_\ell}{\lvert K_0\rvert}\sum_{k\in K_0}\left\lvert\partial_t^\ell\widehat{u}_k(0)-\widehat{h}^{(\ell)}_k\right\rvert^2,\qquad\mathcal{L}=\mathcal{L}_{\mathrm{PDE}}^{\mathrm{prec}}+\lambda_{\mathrm{IC}}\mathcal{L}_{\mathrm{IC}}+\mathcal{L}_{\mathrm{con}}. \tag*{(11)} L IC = ℓ = 0 ∑ s t − 1 ∣ K 0 ∣ γ ℓ k ∈ K 0 ∑ ∂ t ℓ u k ( 0 ) − h k ( ℓ ) 2 , L = L PDE prec + λ IC L IC + L con . (11) Here, γ ℓ > 0 \gamma_\ell>0 γ ℓ > 0 L c o n \mathcal{L}_{\mathrm{con}} L con u = h K + t v u=h_K+tv u = h K + t v h K h_K h K v v v L I C \mathcal{L}_{\mathrm{IC}} L IC
Low-Rank Spectral Residual Correction
A fixed small rank model may reach an accuracy plateau because its available separated components cannot resolve the remaining error. TensorSpectra introduces additional capacity in stages. After training the base model, it freezes the learned parameters and adds a new low-rank component of the form
Δ U ^ ( m ) ( t ) = ∑ α = 1 r m δ λ α ( m ) δ T α ( m ) ( t ) δ a α ( m , 1 ) ∘ ⋯ ∘ δ a α ( m , d ) , U ^ [ m ] ( t ) = U ^ [ m − 1 ] ( t ) + Δ U ^ ( m ) ( t ) . (12) \Delta\widehat{U}^{(m)}(t)=\sum_{\alpha=1}^{r_m}\delta\lambda^{(m)}_\alpha\delta T^{(m)}_\alpha(t)\,\delta a^{(m,1)}_\alpha\circ\cdots\circ\delta a^{(m,d)}_\alpha,\qquad\widehat{U}^{[m]}(t)=\widehat{U}^{[m-1]}(t)+\Delta\widehat{U}^{(m)}(t). \tag*{(12)} Δ U ( m ) ( t ) = α = 1 ∑ r m δ λ α ( m ) δ T α ( m ) ( t ) δ a α ( m , 1 ) ∘ ⋯ ∘ δ a α ( m , d ) , U [ m ] ( t ) = U [ m − 1 ] ( t ) + Δ U ( m ) ( t ) . (12) After m m m r + ∑ ℓ = 1 m r ℓ r+\sum_{\ell=1}^{m}r_\ell r + ∑ ℓ = 1 m r ℓ U ^ [ 0 ] \widehat{U}^{[0]} U [ 0 ] θ m \theta_m θ m
θ m ∗ = arg min θ m L [ U ^ [ m − 1 ] + Δ U ^ ( m ) ( θ m ) ] . (13) \theta_m^*=\underset{\theta_m}{\operatorname{arg\,min}}\ \mathcal{L}\left[\widehat{U}^{[m-1]}+\Delta\widehat{U}^{(m)}(\theta_m)\right]. \tag*{(13)} θ m ∗ = θ m arg min L [ U [ m − 1 ] + Δ U ( m ) ( θ m ) ] . (13) The same full objective evaluates the cumulative solution, including any soft initial or boundary constraints. For linear PDEs, training seeks a correction whose operator image cancels the remaining residual; for nonlinear PDEs, the full nonlinear residual is recomputed without linearization. Hard-constrained models use increments satisfying the corresponding homogeneous constraints, such as t δ v t\,\delta v t δ v
Theoretical Analysis
In this section, we analyze the coefficient-space method described in § 3 , which combines a CP spectral representation, factorwise linear operations, weighted spectral residual training, and additive corrections whose earlier stages are frozen. Let H = L 2 ( Ω , μ ) H=L^2(\Omega,\mu) H = L 2 ( Ω , μ ) μ \mu μ V N = ⨂ j = 1 d V N j ( j ) V_N=\bigotimes_{j=1}^d V_{N_j}^{(j)} V N = ⨂ j = 1 d V N j ( j ) P N P_N P N H H H u N = P N u u_N=P_Nu u N = P N u u ^ N \widehat{u}_N u N
E s p e c ( N ) = ∥ u − u N ∥ H , E r a n k ( N , r ) = inf rank C P ( V ~ ) ≤ r ∥ u ^ N − V ~ ∥ F . (14) E_{\mathrm{spec}}(N)=\lVert u-u_N\rVert_H,\qquad E_{\mathrm{rank}}(N,r)=\inf_{\operatorname{rank}_{\mathrm{CP}}(\widetilde{V})\le r}\lVert\widehat{u}_N-\widetilde{V}\rVert_F. \tag*{(14)} E spec ( N ) = ∥ u − u N ∥ H , E rank ( N , r ) = rank CP ( V ) ≤ r inf ∥ u N − V ∥ F . (14) Here ⊗ \otimes ⊗ ∘ \circ ∘ § 3 . We use an infimum because the rank-r r r r r r u ^ θ \widehat{u}_\theta u θ V ^ η \widehat{V}_\eta V η r r r ∥ u ^ N − V ^ η ∥ F ≤ E r a n k ( N , r ) + η \lVert\widehat{u}_N-\widehat{V}_\eta\rVert_F\le E_{\mathrm{rank}}(N,r)+\eta ∥ u N − V η ∥ F ≤ E rank ( N , r ) + η E o p t ( θ ; V ^ η ) = ∥ V ^ η − u ^ θ ∥ F E_{\mathrm{opt}}(\theta;\widehat{V}_\eta)=\lVert\widehat{V}_\eta-\widehat{u}_\theta\rVert_F E opt ( θ ; V η ) = ∥ V η − u θ ∥ F
Representation and Computational Cost
We first quantify storage and factorwise operation costs for a prescribed CP rank.
Proposition 4.1 (Storage and factorwise linear work). A dense coefficient tensor for V N V_N V N ∏ j N j \prod_j N_j ∏ j N j r r r r + r ∑ j N j r+r\sum_j N_j r + r ∑ j N j N j = N N_j=N N j = N N d N^d N d r + d N r r+dNr r + d N r d d d O ( r ∑ j C s p ( N j ) ) O\left(r\sum_j C_{\mathrm{sp}}(N_j)\right) O ( r ∑ j C sp ( N j ) ) C s p ( N j ) C_{\mathrm{sp}}(N_j) C sp ( N j ) O ( d N r ) O(dNr) O ( d N r ) O ( d N 2 r ) O(dN^2r) O ( d N 2 r ) L L L L r Lr L r
See the proof in Appendix Appendix § 12 . Differentiation preserves each CP term, whereas nonlinear products create pairwise interactions. The next result counts those interactions after factorwise de-aliasing (Appendix § 9 ).
Proposition 4.2 (Polynomial nonlinearities). For a quadratic nonlinearity, the symmetric expansion of a rank-r r r P 2 = ( r + 1 2 ) P_2=\binom{r+1}{2} P 2 = ( 2 r + 1 ) N = 2 K + 1 N=2K+1 N = 2 K + 1 M > 3 K M>3K M > 3 K O ( d ( r + P 2 ) M log M + d P 2 N ) O\left(d(r+P_2)M\log M+dP_2N\right) O ( d ( r + P 2 ) M log M + d P 2 N ) p p p Q Q Q O ( d ( r + P 2 ) Q ( p + 1 ) ) O\left(d(r+P_2)Q(p+1)\right) O ( d ( r + P 2 ) Q ( p + 1 ) ) O ( d P 2 N ) O(dP_2N) O ( d P 2 N ) O ( d P 2 ( p + 1 ) ) O(dP_2(p+1)) O ( d P 2 ( p + 1 )) P 3 = ( r + 2 3 ) P_3=\binom{r+2}{3} P 3 = ( 3 r + 2 ) O ( d N R 2 ) O(dNR^2) O ( d N R 2 ) R R R
See the proof in Appendix Appendix § 12 . We next analyze spectral truncation and the rank needed to approximate the retained coefficient tensor.
Spectral and CP Approximation
The projection estimate isolates the error caused by retaining finitely many modes. Its proof uses Fourier coefficient tails and one-dimensional Jacobi projection estimates respectively.
Lemma 4.3 (Spectral projection). For the Fourier box K K = { k ∈ Z d : ∣ k j ∣ ≤ K } \mathcal{K}_K=\{k\in\mathbb{Z}^d: |k_j|\le K\} K K = { k ∈ Z d : ∣ k j ∣ ≤ K } u ∈ H s ( T d ) u\in H^s(\mathbb{T}^d) u ∈ H s ( T d ) s > m ≥ 0 s>m\ge0 s > m ≥ 0
∥ u − P K u ∥ H m ≤ ( 1 + K 2 ) ( m − s ) / 2 ∥ u ∥ H s . (15) \lVert u-P_Ku\rVert_{H^m}\le(1+K^2)^{(m-s)/2}\lVert u\rVert_{H^s}. \tag*{(15)} ∥ u − P K u ∥ H m ≤ ( 1 + K 2 ) ( m − s ) /2 ∥ u ∥ H s . (15) If ∣ u ^ k ∣ ≤ B e − a ∣ k ∣ 1 |\widehat{u}_k|\le Be^{-a|k|_1} ∣ u k ∣ ≤ B e − a ∣ k ∣ 1 ∥ u − P K u ∥ L 2 ≤ C ( d , a ) B e − a K \lVert u-P_Ku\rVert_{L^2}\le C(d,a)Be^{-aK} ∥ u − P K u ∥ L 2 ≤ C ( d , a ) B e − a K w w w ∥ v ∥ w 2 = ∫ Ω ∣ v ∣ 2 w d x \lVert v\rVert_w^2=\int_\Omega|v|^2w\,dx ∥ v ∥ w 2 = ∫ Ω ∣ v ∣ 2 w d x ∥ v ∥ j , s 2 = ∫ Ω ∣ v ∣ 2 ( 1 − x j 2 ) s w d x \lVert v\rVert_{j,s}^2=\int_\Omega|v|^2(1-x_j^2)^s w\,dx ∥ v ∥ j , s 2 = ∫ Ω ∣ v ∣ 2 ( 1 − x j 2 ) s w d x p j p_j p j j j j s j ≥ 1 s_j\ge1 s j ≥ 1 p j ≥ s j − 1 p_j\ge s_j-1 p j ≥ s j − 1
∥ u − P p u ∥ w ≤ ∑ j = 1 d C α j , β j , s j ( p j + 1 ) − s j ∥ ∂ j s j u ∥ j , s j . (16) \lVert u-P_pu\rVert_w\le\sum_{j=1}^{d}C_{\alpha_j,\beta_j,s_j}(p_j+1)^{-s_j}\lVert\partial_j^{s_j}u\rVert_{j,s_j}. \tag*{(16)} ∥ u − P p u ∥ w ≤ j = 1 ∑ d C α j , β j , s j ( p j + 1 ) − s j ∥ ∂ j s j u ∥ j , s j . (16) The factor ( 1 − x j 2 ) s j (1-x_j^2)^{s_j} ( 1 − x j 2 ) s j j j j ( α j , β j ) (\alpha_j,\beta_j) ( α j , β j ) ( α j + s j , β j + s j ) (\alpha_j+s_j,\beta_j+s_j) ( α j + s j , β j + s j )
See the Fourier- and Jacobi-type proofs in Appendix Appendix § 12 . The following rank property identifies what increasing the rank guarantees without assuming a decay rate.
Lemma 4.4 (Rank error). E r a n k ( N , r + 1 ) ≤ E r a n k ( N , r ) E_{\mathrm{rank}}(N,r+1)\le E_{\mathrm{rank}}(N,r) E rank ( N , r + 1 ) ≤ E rank ( N , r ) r ≥ rank C P ( U ^ N ) r\ge\operatorname{rank}_{\mathrm{CP}}(\widehat{U}_N) r ≥ rank CP ( U N ) r r r
See the proof in Appendix Appendix § 12 . Combining the two approximation terms with the gap left by training gives the overall solution error. The proof applies the triangle inequality in physical space and then in coefficient space, using the synthesis isometry.
Theorem 4.5 (Spectral–rank–training decomposition). Under the definitions above, for every η > 0 \eta>0 η > 0 r r r u ^ θ \widehat{u}_\theta u θ u θ u_\theta u θ
∥ u − u θ ∥ H ≤ E s p e c ( N ) + E r a n k ( N , r ) + η + E o p t ( θ ; V ^ N ) . (17) \lVert u-u_\theta\rVert_H\le E_{\mathrm{spec}}(N)+E_{\mathrm{rank}}(N,r)+\eta+E_{\mathrm{opt}}(\theta;\widehat{V}_N). \tag*{(17)} ∥ u − u θ ∥ H ≤ E spec ( N ) + E rank ( N , r ) + η + E opt ( θ ; V N ) . (17) Thus increasing N N N r r r
See the proof in Appendix Appendix § 12 .
Residual Stability and Spectral Weighting
The preceding decomposition measures a training gap in solution coefficients. To relate the PDE residual used in § 3 to solution error, we next apply a stability inequality to the solution error and distinguish a complete residual from sampled training coefficients.
Theorem 4.6 (Full residual and sampled-loss distinction). Let A : X 0 → Y A : X_0 \to Y A : X 0 → Y ∥ v ∥ X ≤ C s t a b ∥ A v ∥ Y \lVert v\rVert_X \le C_{\mathrm{stab}}\lVert Av\rVert_Y ∥ v ∥ X ≤ C stab ∥ A v ∥ Y v ∈ X 0 v \in X_0 v ∈ X 0 A u = f Au = f A u = f u θ u_\theta u θ
∥ u − u θ ∥ X ≤ C s t a b ∥ A u θ − f ∥ Y . (18) \lVert u-u_\theta\rVert_X \le C_{\mathrm{stab}}\lVert Au_\theta-f\rVert_Y . \tag*{(18)} ∥ u − u θ ∥ X ≤ C stab ∥ A u θ − f ∥ Y . (18) For a stationary residual, write R = A u θ − f R = Au_\theta-f R = A u θ − f R k R_k R k 0 < w min ≤ w k ≤ w max 0 < w_{\min} \le w_k \le w_{\max} 0 < w m i n ≤ w k ≤ w m a x K \mathcal{K} K n n n § 3 , L P D E p r e c , f u l l = n − 1 ∑ k ∈ K ∣ w k R k ∣ 2 \mathcal{L}^{\mathrm{prec,full}}_{\mathrm{PDE}} = n^{-1}\sum_{k\in\mathcal{K}} \lvert w_kR_k\rvert^2 L PDE prec , full = n − 1 ∑ k ∈ K ∣ w k R k ∣ 2
∥ u − u θ ∥ X ≤ C s t a b ( w min − 1 n L P D E p r e c , f u l l + ∥ R K c ∥ Y ) , (19) \lVert u-u_\theta\rVert_X \le C_{\mathrm{stab}}\left(w_{\min}^{-1}\sqrt{n\mathcal{L}^{\mathrm{prec,full}}_{\mathrm{PDE}}}+\lVert R_{\mathcal{K}^c}\rVert_Y\right), \tag*{(19)} ∥ u − u θ ∥ X ≤ C stab ( w m i n − 1 n L PDE prec , full + ∥ R K c ∥ Y ) , (19) provided the coefficient split is orthogonal in Y Y Y S ⊊ K \mathcal{S}\subsetneq\mathcal{K} S ⊊ K L S = ∣ S ∣ − 1 ∑ k ∈ S ∣ w k R k ∣ 2 \mathcal{L}_{\mathcal{S}}=\lvert\mathcal{S}\rvert^{-1}\sum_{k\in\mathcal{S}}\lvert w_kR_k\rvert^2 L S = ∣ S ∣ − 1 ∑ k ∈ S ∣ w k R k ∣ 2 R S R_{\mathcal{S}} R S R S c R_{\mathcal{S}^c} R S c
See the proof in Appendix E. When the method measures only a projected residual, the relevant target is the Galerkin solution. A discrete stability inequality and the triangle inequality then separate the projected residual from the consistency error.
Proposition 4.7 (A projected-residual alternative). Suppose A : dom ( A ) → H A : \operatorname{dom}(A) \to H A : dom ( A ) → H V N ⊂ dom ( A ) V_N \subset\operatorname{dom}(A) V N ⊂ dom ( A ) f ∈ H f \in H f ∈ H A N = P N A ∣ V N A_N=P_NA|_{V_N} A N = P N A ∣ V N ∥ z ∥ H ≤ β N − 1 ∥ A N z ∥ H \lVert z\rVert_H \le\beta_N^{-1}\lVert A_Nz\rVert_H ∥ z ∥ H ≤ β N − 1 ∥ A N z ∥ H z ∈ V N z \in V_N z ∈ V N u N G ∈ V N u_N^G \in V_N u N G ∈ V N A N u N G = P N f A_Nu_N^G=P_Nf A N u N G = P N f v ∈ V N v \in V_N v ∈ V N
∥ u − v ∥ H ≤ ∥ u − u N G ∥ H + β N − 1 ∥ P N ( A v − f ) ∥ H . (20) \lVert u-v\rVert_H \le\lVert u-u_N^G\rVert_H+\beta_N^{-1}\lVert P_N(Av-f)\rVert_H . \tag*{(20)} ∥ u − v ∥ H ≤ ∥ u − u N G ∥ H + β N − 1 ∥ P N ( A v − f ) ∥ H . (20) Thus a complete weighted Galerkin residual controls error relative to the discrete solution, with a separate consistency error and a possibly N , d N,d N , d β N \beta_N β N
See the proof in Appendix E. For periodic Poisson, Fourier diagonalization makes the stability factor and the effect of spectral weights explicit. This example also shows when the weighted residual equals the solution error exactly.
Corollary 4.8 (Zero-mean periodic Poisson). Use normalized Haar measure on T d \mathbb{T}^d T d u u u u θ u_\theta u θ R = Δ u θ − f R=\Delta u_\theta-f R = Δ u θ − f Δ u = f \Delta u=f Δ u = f R k = R ^ k R_k=\widehat{R}_k R k = R k e = u θ − u e=u_\theta-u e = u θ − u k ≠ 0 k\ne0 k = 0
e ^ k = − R ^ k ∣ k ∣ 2 , ∥ e ∥ L 2 2 = ∑ k ≠ 0 ∣ R k ∣ 2 ∣ k ∣ 4 , ∥ ∇ e ∥ L 2 2 = ∑ k ≠ 0 ∣ R k ∣ 2 ∣ k ∣ 2 . (21) \widehat{e}_k=-\frac{\widehat{R}_k}{\lvert k\rvert^2},\qquad\lVert e\rVert_{L^2}^2=\sum_{k\ne0}\frac{\lvert R_k\rvert^2}{\lvert k\rvert^4},\qquad\lVert\nabla e\rVert_{L^2}^2=\sum_{k\ne0}\frac{\lvert R_k\rvert^2}{\lvert k\rvert^2}. \tag*{(21)} e k = − ∣ k ∣ 2 R k , ∥ e ∥ L 2 2 = k = 0 ∑ ∣ k ∣ 4 ∣ R k ∣ 2 , ∥ ∇ e ∥ L 2 2 = k = 0 ∑ ∣ k ∣ 2 ∣ R k ∣ 2 . (21) Consequently, ideal weights w k = ∣ k ∣ − 2 w_k=\lvert k\rvert^{-2} w k = ∣ k ∣ − 2 L 2 L^2 L 2
See the proof in Appendix E. The same residual principle gives bounds for heat, wave, and selected nonlinear equations once initial mismatches or nonlinear stability conditions are included; see Proposition 2, E.3 and E.4 in Appendix E.
Frozen Low-Rank Correction
We now analyze the frozen corrections introduced in § 3 . The first result tracks their total representational capacity by concatenating CP terms. Let U ^ [ 0 ] \widehat{\mathcal{U}}^{[0]} U [ 0 ] r 0 r_0 r 0 U ^ [ m ] = U ^ [ m − 1 ] + Δ U ^ ( m ) \widehat{\mathcal{U}}^{[m]}=\widehat{\mathcal{U}}^{[m-1]}+\Delta\widehat{\mathcal{U}}^{(m)} U [ m ] = U [ m − 1 ] + Δ U ( m ) rank C P ( Δ U ^ ( m ) ) ≤ r m \operatorname{rank}_{\mathrm{CP}}(\Delta\widehat{\mathcal{U}}^{(m)})\le r_m rank CP ( Δ U ( m ) ) ≤ r m
Lemma 4.9 (Rank accumulation). rank C P ( U ^ [ M ] ) ≤ r 0 + ∑ m = 1 M r m \operatorname{rank}_{\mathrm{CP}}(\widehat{\mathcal{U}}^{[M]})\le r_0+\sum_{m=1}^{M}r_m rank CP ( U [ M ] ) ≤ r 0 + ∑ m = 1 M r m
See the proof in Appendix E. The next theorem states the extra approximation and optimization conditions that make a correction recurrence possible; its proof iterates a one-stage inequality.
Theorem 4.10 (Conditional correction recurrence). Fix a target coefficient tensor U ^ N ∗ \widehat{\mathcal{U}}_N^* U N ∗ E ^ [ m ] = U ^ N ∗ − U ^ [ m ] \widehat{\mathcal{E}}^{[m]}=\widehat{\mathcal{U}}_N^*-\widehat{\mathcal{U}}^{[m]} E [ m ] = U N ∗ − U [ m ] δ m = inf rank C P ( V ^ ) ≤ r m ∥ E ^ [ m − 1 ] − V ^ ∥ F \delta_m=\inf_{\operatorname{rank}_{\mathrm{CP}}(\widehat{\mathcal{V}})\le r_m}\lVert\widehat{\mathcal{E}}^{[m-1]}-\widehat{\mathcal{V}}\rVert_F δ m = inf rank CP ( V ) ≤ r m ∥ E [ m − 1 ] − V ∥ F m m m ∥ E ^ [ m ] ∥ F ≤ δ m + ϵ m \lVert\widehat{\mathcal{E}}^{[m]}\rVert_F\le\delta_m+\epsilon_m ∥ E [ m ] ∥ F ≤ δ m + ϵ m
then ∥ E ^ [ m ] ∥ F ≤ ∥ E ^ [ m − 1 ] ∥ F + ε m \|\widehat{E}^{[m]}\|_{F} \le\|\widehat{E}^{[m-1]}\|_{F} + \varepsilon_m ∥ E [ m ] ∥ F ≤ ∥ E [ m − 1 ] ∥ F + ε m δ m ≤ q m ∥ E ^ [ m − 1 ] ∥ F \delta_m \le q_m\|\widehat{E}^{[m-1]}\|_{F} δ m ≤ q m ∥ E [ m − 1 ] ∥ F 0 ≤ q m < 1 0 \le q_m < 1 0 ≤ q m < 1
∥ E ^ [ M ] ∥ F ≤ ( ∏ j = 1 M q j ) ∥ E ^ [ 0 ] ∥ F + ∑ m = 1 M ε m ∏ j = m + 1 M q j . (22) \|\widehat{E}^{[M]}\|_{F} \le\left(\prod_{j=1}^{M}q_j\right)\|\widehat{E}^{[0]}\|_{F}+\sum_{m=1}^{M}\varepsilon_m\prod_{j=m+1}^{M}q_j. \tag*{(22)} ∥ E [ M ] ∥ F ≤ ( j = 1 ∏ M q j ) ∥ E [ 0 ] ∥ F + m = 1 ∑ M ε m j = m + 1 ∏ M q j . (22) For q m ≤ q < 1 q_m \le q < 1 q m ≤ q < 1 ε m ≤ ε \varepsilon_m \le\varepsilon ε m ≤ ε ∥ E ^ [ M ] ∥ F ≤ q M ∥ E ^ [ 0 ] ∥ F + ε ( 1 − q M ) / ( 1 − q ) \|\widehat{E}^{[M]}\|_{F} \le q^M\|\widehat{E}^{[0]}\|_{F}+\varepsilon(1-q^M)/(1-q) ∥ E [ M ] ∥ F ≤ q M ∥ E [ 0 ] ∥ F + ε ( 1 − q M ) / ( 1 − q ) Proposition E.5 in Appendix § 12 states what a decrease in the training objective alone guarantees.
See the proof in Appendix § 12 . Finally, we combine the correction recurrence with the spectral projection error. The same solution can also be assessed through full-residual stability, which supplies an independent a posteriori certificate.
Theorem 4.11 (Combined approximation and certification). Take U ^ N ∗ = U ^ N \widehat{U}^{*}_{N}=\widehat{U}_{N} U N ∗ = U N Theorem 4.10 , assume its stage hypotheses, and set u T S [ M ] = S N U ^ [ M ] u_{\mathrm{TS}}^{[M]}=S_N\widehat{U}^{[M]} u TS [ M ] = S N U [ M ] (107) . Then
∥ u − u T S [ M ] ∥ H ≤ E s p e c ( N ) + ( ∏ j = 1 M q j ) ∥ U ^ N − U ^ [ 0 ] ∥ F + ∑ m = 1 M ε m ∏ j = m + 1 M q j . (23) \|u-u_{\mathrm{TS}}^{[M]}\|_{H}\le E_{\mathrm{spec}}(N)+\left(\prod_{j=1}^{M}q_j\right)\|\widehat{U}_{N}-\widehat{U}^{[0]}\|_{F}+\sum_{m=1}^{M}\varepsilon_m\prod_{j=m+1}^{M}q_j. \tag*{(23)} ∥ u − u TS [ M ] ∥ H ≤ E spec ( N ) + ( j = 1 ∏ M q j ) ∥ U N − U [ 0 ] ∥ F + m = 1 ∑ M ε m j = m + 1 ∏ M q j . (23) If U ^ [ 0 ] = U ^ θ \widehat{U}^{[0]}=\widehat{U}_{\theta} U [ 0 ] = U θ E r a n k ( N , r 0 ) + η + E o p t ( θ ; V η ) E_{\mathrm{rank}}(N,r_0)+\eta+E_{\mathrm{opt}}(\theta;\mathcal{V}_{\eta}) E rank ( N , r 0 ) + η + E opt ( θ ; V η ) Theorem 4.5 . Separately, whenever the full-residual assumptions of Theorem 4.6 hold, its a posteriori bound applies to the same u T S [ M ] u_{\mathrm{TS}}^{[M]} u TS [ M ]
See the proof in Appendix § 12 . The Fourier and Jacobi estimates use classical spectral results [33 , 36 ] , and the CP properties follow standard tensor analysis [22 ] . The additional stability and de-aliasing results in Appendix § 12 specify the assumptions used for those cases.
Numerical Experiments
In this section, we systematically evaluate TensorSpectra on a range of linear and nonlinear PDE benchmarks, with particular emphasis on solution accuracy and scalability to high dimensions.
Benchmarks, baselines and implementation details. We evaluate TensorSpectra on five representative PDE benchmarks comprising the Poisson, heat, wave, steady-state Allen–Cahn, and nonlinear reaction equations, covering high-dimensional linear and middle-dimensional nonlinear problems. We compare against five baselines spanning different methodological paradigms, including PINN [32 ] , separable and tensorized methods (SPINN [5 ] and TNN [40 ] ), the architecture-enhanced PirateNet [37 ] , and the spectral-domain solver SINN [42 ] . Experiments with d < 50 d < 50 d < 50 d ≥ 50 d \ge50 d ≥ 50 [21 ] is adopted as the default optimizer. Benchmark definitions, baseline descriptions, and detailed experimental configurations are provided in § 5.1 and § 5.2 , Appendix § 11 , and Appendix § 14 respectively.
Evaluation metrics. We use the relative L 2 L^2 L 2
Rel. L 2 = ∥ u p r e d − u t r u e ∥ L 2 ( Ω ) ∥ u t r u e ∥ L 2 ( Ω ) , Rel. L 2 ( T ) = ∥ u p r e d ( ⋅ , T ) − u t r u e ( ⋅ , T ) ∥ L 2 ( Ω ) ∥ u t r u e ( ⋅ , T ) ∥ L 2 ( Ω ) . (24) \operatorname{Rel.}L^2=\frac{\|u_{\mathrm{pred}}-u_{\mathrm{true}}\|_{L^2(\Omega)}}{\|u_{\mathrm{true}}\|_{L^2(\Omega)}},\qquad\operatorname{Rel.}L^2(T)=\frac{\|u_{\mathrm{pred}}(\cdot,T)-u_{\mathrm{true}}(\cdot,T)\|_{L^2(\Omega)}}{\|u_{\mathrm{true}}(\cdot,T)\|_{L^2(\Omega)}}. \tag*{(24)} Rel. L 2 = ∥ u true ∥ L 2 ( Ω ) ∥ u pred − u true ∥ L 2 ( Ω ) , Rel. L 2 ( T ) = ∥ u true ( ⋅ , T ) ∥ L 2 ( Ω ) ∥ u pred ( ⋅ , T ) − u true ( ⋅ , T ) ∥ L 2 ( Ω ) . (24) For all time-dependent problems, we report the relative L 2 L^2 L 2 T T T
Main results. Table 1 shows that TensorSpectra consistently achieves the lowest relative L 2 L^2 L 2 d = 100 d=100 d = 100 9.06 × 10 − 8 9.06\times10^{-8} 9.06 × 1 0 − 8 3.82 × 10 − 3 3.82\times10^{-3} 3.82 × 1 0 − 3 6.31 × 10 − 3 6.31\times10^{-3} 6.31 × 1 0 − 3 Figure 2 shows the training loss curves on nonlinear problems. Furthermore, Appendix § 13 reports the computational cost per training update on the 30D wave benchmark, showing that TensorSpectra achieves high accuracy with substantially lower GFLOPs than most baselines.
Table 1. Comparison of relative L 2 L^2 L 2 ± \pm ± Bold values indicate the best performance.
Figure 2. Optimization trajectories on nonlinear PDE benchmarks. Evaluation loss versus optimizer updates for (a) 5D Allen–Cahn and (b) 10D Allen–Cahn. TensorSpectra consistently converges to lower evaluation loss than the competing baselines.
Interpretability of the learned CP representation. Beyond prediction accuracy, the learned CP representation offers a direct view of how TensorSpectra organizes information across spectral coordinates and rank components. Figure 3 visualizes the learned factors along representative coordinate dimensions. The real-valued factors are predominantly concentrated in low-frequency modes, while the imaginary factors exhibit highly consistent patterns across both ranks and dimensions. These factor-wise visualizations reveal the internal spectral organization of the CP approximation without requiring reconstruction of the full high-dimensional coefficient tensor. Figure 4 further summarizes the relative contribution of individual CP components through the normalized amplitudes of the stored real and imaginary parts. Most of the representation mass is concentrated in only a few dominant components, whereas the remaining components contribute marginally. This concentration provides empirical evidence that TensorSpectra recovers an effectively low-rank structure in the spectral coefficient space.
Figure 3. Magnitudes of the learned real and imaginary CP spectral factors over signed Fourier modes for representative dimensions j = 1 j = 1 j = 1 5 5 5 10 10 10
Figure 4. Relative importance of the learned CP components, measured by the normalized component amplitudes λ r / ∑ s λ s \lambda_r / \sum_s \lambda_s λ r / ∑ s λ s
High-Dimensional Linear PDE Benchmarks
We consider the Poisson, heat and wave equations as three canonical linear PDE benchmarks. Together, they cover elliptic, parabolic and hyperbolic PDEs respectively, and therefore provide a representative testbed for evaluating the accuracy and dimensional scalability of TensorSpectra across fundamentally different classes of differential operators and solution behaviors.
Poisson equation. We first consider the d d d
Δ u ( x ) = f ( x ) , x ∈ [ 0 , 2 π ] d . (25) \Delta u(\mathbf{x}) = f(\mathbf{x}), \quad\mathbf{x} \in[0,2\pi]^d. \tag*{(25)} Δ u ( x ) = f ( x ) , x ∈ [ 0 , 2 π ] d . (25) The forcing term is constructed from
f p ( x ) = 1 π 2 + 1 130 ∣ D ∣ ∑ i ∈ D cosh ( 2 ( x i − π ) ) , (26) f_p(\mathbf{x}) = \frac{1}{\pi^2} + \frac{1}{130\lvert D\rvert}\sum_{i\in D}\cosh(2(x_i-\pi)), \tag*{(26)} f p ( x ) = π 2 1 + 130 ∣ D ∣ 1 i ∈ D ∑ cosh ( 2 ( x i − π )) , (26) where D ⊆ { 1 , … , d } D \subseteq\{1,\ldots,d\} D ⊆ { 1 , … , d } f p f_p f p
f ( x ) = ∑ k ∈ K Re ( f ^ k e i k ⋅ x ) , u ( x ) = ∑ k ∈ K Re ( u ^ k e i k ⋅ x ) , u ^ k = − f ^ k ∥ k ∥ 2 2 . (27) f(\mathbf{x}) = \sum_{\mathbf{k}\in\mathcal{K}} \operatorname{Re}\left(\widehat{f}_{\mathbf{k}}e^{\mathrm{i}\mathbf{k}\cdot\mathbf{x}}\right), \quad u(\mathbf{x}) = \sum_{\mathbf{k}\in\mathcal{K}} \operatorname{Re}\left(\widehat{u}_{\mathbf{k}}e^{\mathrm{i}\mathbf{k}\cdot\mathbf{x}}\right), \quad\widehat{u}_{\mathbf{k}} = -\frac{\widehat{f}_{\mathbf{k}}}{\lVert\mathbf{k}\rVert_2^2}. \tag*{(27)} f ( x ) = k ∈ K ∑ Re ( f k e i k ⋅ x ) , u ( x ) = k ∈ K ∑ Re ( u k e i k ⋅ x ) , u k = − ∥ k ∥ 2 2 f k . (27) We evaluate all methods for d ∈ { 2 , 5 , 10 , 30 , 50 , 100 } d \in\{2,5,10,30,50,100\} d ∈ { 2 , 5 , 10 , 30 , 50 , 100 }
Heat equation. We next consider the d d d
∂ u ∂ t ( x , t ) = Δ u ( x , t ) , ( x , t ) ∈ [ 0 , 2 π ] d × [ 0 , 0.1 ] . (28) \frac{\partial u}{\partial t}(\mathbf{x},t) = \Delta u(\mathbf{x},t), \quad(\mathbf{x},t) \in[0,2\pi]^d \times[0,0.1]. \tag*{(28)} ∂ t ∂ u ( x , t ) = Δ u ( x , t ) , ( x , t ) ∈ [ 0 , 2 π ] d × [ 0 , 0.1 ] . (28) The initial condition and reference solution are constructed from the solution of the Poisson problem introduced above,
u ( x , 0 ) = ∑ k ∈ K Re ( u ^ k e i k ⋅ x ) , u ( x , t ) = ∑ k ∈ K Re ( u ^ k e − ∥ k ∥ 2 2 t e i k ⋅ x ) , (29) u(\mathbf{x},0) = \sum_{\mathbf{k}\in\mathcal{K}}\operatorname{Re}\left(\widehat{u}_{\mathbf{k}}e^{\mathrm{i}\mathbf{k}\cdot\mathbf{x}}\right), \quad u(\mathbf{x},t) = \sum_{\mathbf{k}\in\mathcal{K}}\operatorname{Re}\left(\widehat{u}_{\mathbf{k}}e^{-\lVert\mathbf{k}\rVert_2^2t}e^{\mathrm{i}\mathbf{k}\cdot\mathbf{x}}\right), \tag*{(29)} u ( x , 0 ) = k ∈ K ∑ Re ( u k e i k ⋅ x ) , u ( x , t ) = k ∈ K ∑ Re ( u k e − ∥ k ∥ 2 2 t e i k ⋅ x ) , (29) where u ^ k \widehat{u}_{\mathbf{k}} u k d ∈ { 2 , 5 , 10 , 30 , 50 , 100 } d \in\{2,5,10,30,50,100\} d ∈ { 2 , 5 , 10 , 30 , 50 , 100 }
Wave equation. Finally, we consider a d d d
∂ 2 u ∂ t 2 ( x , t ) = Δ u ( x , t ) , ( x , t ) ∈ [ 0 , 2 π ] d × [ 0 , 0.1 ] . (30) \frac{\partial^2u}{\partial t^2}(\mathbf{x},t) = \Delta u(\mathbf{x},t), \quad(\mathbf{x},t) \in[0,2\pi]^d \times[0,0.1]. \tag*{(30)} ∂ t 2 ∂ 2 u ( x , t ) = Δ u ( x , t ) , ( x , t ) ∈ [ 0 , 2 π ] d × [ 0 , 0.1 ] . (30) We impose periodic boundary conditions and use the solution of the Poisson benchmark as the initial displacement, while setting the initial velocity to zero,
u ( x , 0 ) = ∑ k ∈ K Re ( u ^ k e i k ⋅ x ) , ∂ u ∂ t ( x , 0 ) = 0. (31) u(\mathbf{x},0) = \sum_{\mathbf{k}\in\mathcal{K}}\operatorname{Re}\left(\widehat{u}_{\mathbf{k}}e^{\mathrm{i}\mathbf{k}\cdot\mathbf{x}}\right), \quad\frac{\partial u}{\partial t}(\mathbf{x},0) = 0. \tag*{(31)} u ( x , 0 ) = k ∈ K ∑ Re ( u k e i k ⋅ x ) , ∂ t ∂ u ( x , 0 ) = 0. (31) Since the Fourier modes diagonalize the Laplace operator, each spatial mode evolves independently according to a harmonic oscillator. The corresponding reference solution is therefore
u ( x , t ) = ∑ k ∈ K Re ( u ^ k cos ( ∥ k ∥ 2 t ) e i k ⋅ x ) . (32) u(\mathbf{x},t) = \sum_{\mathbf{k}\in\mathcal{K}}\operatorname{Re}\left(\widehat{u}_{\mathbf{k}}\cos(\lVert\mathbf{k}\rVert_2t)e^{\mathrm{i}\mathbf{k}\cdot\mathbf{x}}\right). \tag*{(32)} u ( x , t ) = k ∈ K ∑ Re ( u k cos (∥ k ∥ 2 t ) e i k ⋅ x ) . (32) As in the Poisson and heat benchmarks, we consider d ∈ { 2 , 5 , 10 , 30 , 50 , 100 } d \in\{2,5,10,30,50,100\} d ∈ { 2 , 5 , 10 , 30 , 50 , 100 }
Middle-Dimensional Nonlinear PDE Benchmarks
We further evaluate TensorSpectra on two middle-dimensional nonlinear PDE benchmarks, including a steady-state Allen–Cahn equation and a nonlinear reaction equation. These problems introduce nonlinear interactions that destroy the simple mode-wise separability available in linear PDEs, providing a more challenging setting for assessing the robustness of the proposed low-rank spectral formulation. To further demonstrate the spectral extensibility of TensorSpectra beyond Fourier representations, we employ Legendre polynomial basis mentioned in Appendix § 8 .§ 2 for all nonlinear benchmarks in this section.
Steady-state Allen–Cahn equation. We consider the d d d Ω = [ − 1 , 1 ] d \Omega= [-1,1]^d Ω = [ − 1 , 1 ] d
− ε 2 Δ u ( x ) + u 3 ( x ) − u ( x ) = f ( x ) , x ∈ Ω , (33) -\varepsilon^2 \Delta u(\boldsymbol{x}) + u^3(\boldsymbol{x}) - u(\boldsymbol{x}) = f(\boldsymbol{x}), \quad\boldsymbol{x} \in\Omega, \tag*{(33)} − ε 2 Δ u ( x ) + u 3 ( x ) − u ( x ) = f ( x ) , x ∈ Ω , (33) subject to the homogeneous Dirichlet boundary condition
u ( x ) = 0 , x ∈ ∂ Ω . (34) u(\boldsymbol{x}) = 0, \quad\boldsymbol{x} \in\partial\Omega. \tag*{(34)} u ( x ) = 0 , x ∈ ∂ Ω. (34) We adopt the manufactured solution
u ( x ) = ∏ j = 1 d ( 1 − x j 2 ) , (35) u(\boldsymbol{x}) = \prod_{j=1}^{d}(1-x_j^2), \tag*{(35)} u ( x ) = j = 1 ∏ d ( 1 − x j 2 ) , (35) from which the forcing term is obtained analytically as
f ( x ) = 2 ε 2 ∑ j = 1 d ∏ k ≠ j ( 1 − x k 2 ) + ∏ j = 1 d ( 1 − x j 2 ) 3 − ∏ j = 1 d ( 1 − x j 2 ) . (36) f(\boldsymbol{x}) = 2\varepsilon^2 \sum_{j=1}^{d} \prod_{k\ne j}(1-x_k^2) + \prod_{j=1}^{d}(1-x_j^2)^3 - \prod_{j=1}^{d}(1-x_j^2). \tag*{(36)} f ( x ) = 2 ε 2 j = 1 ∑ d k = j ∏ ( 1 − x k 2 ) + j = 1 ∏ d ( 1 − x j 2 ) 3 − j = 1 ∏ d ( 1 − x j 2 ) . (36) We set ε = 0.5 \varepsilon= 0.5 ε = 0.5 ( 1 − x j 2 ) q j ( x j ) (1-x_j^2)q_j(x_j) ( 1 − x j 2 ) q j ( x j )
Nonlinear reaction equation. We further consider the logistic reaction equation
∂ u ∂ t ( x , t ) − ρ u ( x , t ) ( 1 − u ( x , t ) ) = 0 , ( x , t ) ∈ [ 0 , 2 π ] d × [ 0 , 1 ] , (37) \frac{\partial u}{\partial t}(\boldsymbol{x},t) - \rho u(\boldsymbol{x},t)(1-u(\boldsymbol{x},t)) = 0,\quad(\boldsymbol{x},t) \in[0,2\pi]^d \times[0,1], \tag*{(37)} ∂ t ∂ u ( x , t ) − ρ u ( x , t ) ( 1 − u ( x , t )) = 0 , ( x , t ) ∈ [ 0 , 2 π ] d × [ 0 , 1 ] , (37) with the Gaussian initial condition
u ( x , 0 ) = h ( x ) = exp ( − ∥ x − π 1 ∥ 2 2 2 ( π / 4 ) 2 ) . (38) u(\boldsymbol{x},0) = h(\boldsymbol{x}) = \exp\left(-\frac{\lVert\boldsymbol{x}-\pi\mathbf{1}\rVert_2^2}{2(\pi/4)^2}\right). \tag*{(38)} u ( x , 0 ) = h ( x ) = exp ( − 2 ( π /4 ) 2 ∥ x − π 1 ∥ 2 2 ) . (38) The initial profile is separable across spatial dimensions,
h ( x ) = ∏ j = 1 d exp ( − ( x j − π ) 2 2 ( π / 4 ) 2 ) , (39) h(\boldsymbol{x}) = \prod_{j=1}^{d}\exp\left(-\frac{(x_j-\pi)^2}{2(\pi/4)^2}\right), \tag*{(39)} h ( x ) = j = 1 ∏ d exp ( − 2 ( π /4 ) 2 ( x j − π ) 2 ) , (39) and the corresponding exact solution is
u ( x , t ) = h ( x ) e t h ( x ) e t + 1 − h ( x ) . (40) u(\boldsymbol{x},t) = \frac{h(\boldsymbol{x})e^t}{h(\boldsymbol{x})e^t+1-h(\boldsymbol{x})}. \tag*{(40)} u ( x , t ) = h ( x ) e t + 1 − h ( x ) h ( x ) e t . (40) To evaluate TensorSpectra beyond Fourier representations, both the spatial coordinates and time are expanded in shifted orthonormal Legendre bases, with time treated as an additional tensor dimension. The Gaussian initial condition is deterministically projected onto the Legendre basis and imposed through the hard construction u ( x , t ) = h K ( x ) + t v ( x , t ) u(\boldsymbol{x},t)=h_K(\boldsymbol{x})+tv(\boldsymbol{x},t) u ( x , t ) = h K ( x ) + t v ( x , t )
Ablation Study and Hyperparameter Analysis
In this section, we conduct ablation studies and hyperparameter analyses to examine the contribution of the key components in TensorSpectra and to assess the sensitivity of the method to important design parameters. We progressively evaluate the effects of residual correction and spectral weighting, and further study how the tensor rank and temporal polynomial degree influence solution accuracy.
Ablation Study
TensorSpectra component ablation. We progressively add the proposed components to the rank 8 backbone. For a fair comparison, the backbone is trained for 50k epochs, while the correction-based variants use 25k epochs for the base stage followed by 25k epochs for the correction stage, resulting in the same total training budget. All results are averaged over three random seeds 100, 200 and 300. Figure 5 shows that the rank 8 backbone achieves a relative L 2 L^2 L 2 3.96 × 10 − 3 3.96 \times10^{-3} 3.96 × 1 0 − 3 3.68 × 10 − 3 3.68 \times10^{-3} 3.68 × 1 0 − 3 w w w 2.68 × 10 − 3 2.68 \times10^{-3} 2.68 × 1 0 − 3
Figure 5. Ablation study of TensorSpectra on the contribution of residual correction and spectral weighting. Starting from a rank 8 backbone, we progressively add a rank 4 correction stage and the spectral weighting w w w L 2 L^2 L 2
Compatibility with LoRFS. Table 2 compares TensorSpectra with LoRFS on the high-dimensional Reaction benchmark under the same problem setting described in § 5.2 . To demonstrate the flexibility of our spectral framework, we evaluate TensorSpectra using Legendre, Chebyshev, and Jacobi polynomial bases. All three variants outperform the Fourier-based LoRFS in both rMAE and rRMSE, with the Legendre basis achieving the best results of 3.62 × 10 − 3 3.62 \times10^{-3} 3.62 × 1 0 − 3 2.36 × 10 − 3 2.36 \times10^{-3} 2.36 × 1 0 − 3 91.15 % 91.15\% 91.15% 71.22 % 71.22\% 71.22%
Table 2. Comparison with LoRFS on the high-dimensional Reaction benchmark (d = 9 d = 9 d = 9 Bold values indicate the best performance and green cells report the percentage improvement over the underlined LoRFS value for each metric, using the best TensorSpectra basis. All the experiments are repeated over five random seeds and the mean ± standard deviation is reported.
Hyperparameter Analysis
Sensitivity to temporal polynomial degree. We evaluate the temporal polynomial degree under a fixed TensorSpectra configuration with a rank 8 backbone and a rank 4 correction stage. For all settings, the backbone and correction stages are trained for 50k epochs each, resulting in the same total budget of 100k epochs. Figure 6 (left) shows that the performance is not monotonic with respect to the temporal degree. Degrees 4 and 8 achieve comparable errors of approximately 5 × 10 − 4 5 \times10^{-4} 5 × 1 0 − 4 10 − 3 10^{-3} 1 0 − 3
Figure 6. Relative L 2 L^2 L 2 ± \pm ±
Sensitivity to tensor rank. We investigate the effect of the backbone and correction ranks while fixing the temporal polynomial degree to 8. For all configurations, the backbone and correction stages are trained for 50k epochs each, resulting in the same total training budget of 100k epochs. Figure 6 (right) shows that increasing the configuration from 4 + 2 4 + 2 4 + 2 8 + 4 8 + 4 8 + 4 L 2 L^2 L 2 8.84 × 10 − 4 8.84 \times10^{-4} 8.84 × 1 0 − 4 5.84 × 10 − 4 5.84 \times10^{-4} 5.84 × 1 0 − 4 3.54 × 10 − 4 3.54 \times10^{-4} 3.54 × 1 0 − 4 4.45 × 10 − 4 4.45 \times10^{-4} 4.45 × 1 0 − 4
Discussions and Conclusions
Our results suggest that high-dimensional spectral PDE solving is fundamentally constrained by the representation complexity of the multidimensional coefficient field. By exploiting low-rank structure directly in spectral space, TensorSpectra retains the accuracy and operator structure of spectral methods while substantially reducing this complexity. More broadly, these results indicate that explicitly encoding solution structure can provide an effective alternative to relying entirely on generic neural representations to discover such structure through optimization.
The efficiency of TensorSpectra depends on spectral compressibility. Strongly non-separable dynamics, singular solutions, and nonlinear interactions can increase the required rank. Future work will therefore focus on adaptive rank control, more flexible tensor formats, and extensions to more challenging nonlinear and inverse problems. These directions may further clarify when low-rank spectral structure provides a practical route to scalable high-dimensional scientific computing.
Data Availability
The datasets generated during and/or analysed during the current study are available from the corresponding author on request.
Conflict of interest
The authors declared that they have no conflict of interest.
Acknowledgments
Zhaoxiang Li is supported by the National Natural Science Foundation of China (12671434, 12571391, 12271366).
Spectral Basis Definitions
In this section, we provide the detailed definitions and basic properties of the spectral bases used in TensorSpectra. Specifically, we summarize Fourier basis for periodic problems and Jacobi polynomial bases for non-periodic settings, together with the normalization, orthogonality, and coefficient conventions adopted throughout the paper. We also clarify the relationships among these polynomial bases, with Legendre and Chebyshev polynomials viewed as important special cases within the broader Jacobi family. For time-dependent problems, the spatial factors are independent of time and each CP component is multiplied by the temporal factor T α ( t ) T_{\alpha}(t) T α ( t ) (5) . Stationary problems are recovered by setting T α ≡ 1 T_{\alpha} \equiv1 T α ≡ 1
Fourier Basis
For periodic problems, TensorSpectra adopts the complex Fourier basis. On the canonical interval [ 0 , 2 π ] [0,2\pi] [ 0 , 2 π ]
ϕ k ( x ) = e i k x , k ∈ Z , (41) \phi_k(x)=e^{ikx},\qquad k\in\mathbb{Z}, \tag*{(41)} ϕ k ( x ) = e ik x , k ∈ Z , (41) with i = − 1 i=\sqrt{-1} i = − 1 j j j α \alpha α
a α ( j ) = [ a α , k j ( j ) ] k j ∈ K j , (42) \mathbf{a}_{\alpha}^{(j)}=\left[a_{\alpha,k_j}^{(j)}\right]_{k_j\in\mathcal{K}_j}, \tag*{(42)} a α ( j ) = [ a α , k j ( j ) ] k j ∈ K j , (42) which defines the corresponding one-dimensional Fourier expansion
ψ α ( j ) ( x j ) = ∑ k j ∈ K j a α , k j ( j ) e i k j x j . (43) \psi_{\alpha}^{(j)}(x_j)=\sum_{k_j\in\mathcal{K}_j}a_{\alpha,k_j}^{(j)}e^{ik_jx_j}. \tag*{(43)} ψ α ( j ) ( x j ) = k j ∈ K j ∑ a α , k j ( j ) e i k j x j . (43) Combining the one-dimensional factors across dimensions gives the separable representation
u ( x , t ) = ∑ α = 1 r λ α T α ( t ) ∏ j = 1 d ψ α ( j ) ( x j ) , (44) u(\mathbf{x},t)=\sum_{\alpha=1}^{r}\lambda_{\alpha}T_{\alpha}(t)\prod_{j=1}^{d}\psi_{\alpha}^{(j)}(x_j), \tag*{(44)} u ( x , t ) = α = 1 ∑ r λ α T α ( t ) j = 1 ∏ d ψ α ( j ) ( x j ) , (44) which is equivalent to the CP representation of the multidimensional Fourier coefficient tensor. In particular, the coefficient associated with the multidimensional mode k = ( k 1 , … , k d ) \mathbf{k}=(k_1,\ldots,k_d) k = ( k 1 , … , k d )
u ^ k ( t ) = ∑ α = 1 r λ α T α ( t ) ∏ j = 1 d a α , k j ( j ) . (45) \widehat{u}_{\mathbf{k}}(t)=\sum_{\alpha=1}^{r}\lambda_{\alpha}T_{\alpha}(t)\prod_{j=1}^{d}a_{\alpha,k_j}^{(j)}. \tag*{(45)} u k ( t ) = α = 1 ∑ r λ α T α ( t ) j = 1 ∏ d a α , k j ( j ) . (45) Jacobi Basis
For non-periodic problems, TensorSpectra employs Jacobi polynomials as a general polynomial spectral basis on the canonical interval [ − 1 , 1 ] [-1,1] [ − 1 , 1 ] n n n J n ( α , β ) ( x ) J_n^{(\alpha,\beta)}(x) J n ( α , β ) ( x ) α , β > − 1 \alpha,\beta>-1 α , β > − 1
w ( α , β ) ( x ) = ( 1 − x ) α ( 1 + x ) β , x ∈ [ − 1 , 1 ] . (46) w^{(\alpha,\beta)}(x)=(1-x)^{\alpha}(1+x)^{\beta},\qquad x\in[-1,1]. \tag*{(46)} w ( α , β ) ( x ) = ( 1 − x ) α ( 1 + x ) β , x ∈ [ − 1 , 1 ] . (46) The Jacobi polynomials form an orthogonal basis with respect to this weighted inner product,
∫ − 1 1 J m ( α , β ) ( x ) J n ( α , β ) ( x ) w ( α , β ) ( x ) d x = h n ( α , β ) δ m n , (47) \int_{-1}^{1}J_m^{(\alpha,\beta)}(x)J_n^{(\alpha,\beta)}(x)w^{(\alpha,\beta)}(x)\,\mathrm{d}x=h_n^{(\alpha,\beta)}\delta_{mn}, \tag*{(47)} ∫ − 1 1 J m ( α , β ) ( x ) J n ( α , β ) ( x ) w ( α , β ) ( x ) d x = h n ( α , β ) δ mn , (47) where h n ( α , β ) h_n^{(\alpha,\beta)} h n ( α , β ) n n n δ m n \delta_{mn} δ mn
Within TensorSpectra, the j j j ℓ \ell ℓ
a ℓ ( j ) = [ a ℓ , k j ( j ) ] k j = 0 N j − 1 , (48) \mathbf{a}_{\ell}^{(j)}=\left[a_{\ell,k_j}^{(j)}\right]_{k_j=0}^{N_j-1}, \tag*{(48)} a ℓ ( j ) = [ a ℓ , k j ( j ) ] k j = 0 N j − 1 , (48) which defines the corresponding one-dimensional Jacobi expansion
ψ ℓ ( j ) ( x j ) = ∑ k j = 0 N j − 1 a ℓ , k j ( j ) J k j ( α , β ) ( x j ) . (49) \psi_{\ell}^{(j)}(x_j)=\sum_{k_j=0}^{N_j-1}a_{\ell,k_j}^{(j)}J_{k_j}^{(\alpha,\beta)}(x_j). \tag*{(49)} ψ ℓ ( j ) ( x j ) = k j = 0 ∑ N j − 1 a ℓ , k j ( j ) J k j ( α , β ) ( x j ) . (49) Combining these one-dimensional expansions across coordinates gives the multidimensional separable representation
u ( x , t ) = ∑ ℓ = 1 r λ ℓ T ℓ ( t ) ∏ j = 1 d ψ ℓ ( j ) ( x j ) , (50) u(\mathbf{x},t)=\sum_{\ell=1}^{r}\lambda_{\ell}T_{\ell}(t)\prod_{j=1}^{d}\psi_{\ell}^{(j)}(x_j), \tag*{(50)} u ( x , t ) = ℓ = 1 ∑ r λ ℓ T ℓ ( t ) j = 1 ∏ d ψ ℓ ( j ) ( x j ) , (50) and the corresponding Jacobi coefficient associated with the multidimensional index k = ( k 1 , … , k d ) \mathbf{k}=(k_1,\ldots,k_d) k = ( k 1 , … , k d )
u ^ k ( t ) = ∑ ℓ = 1 r λ ℓ T ℓ ( t ) ∏ j = 1 d a ℓ , k j ( j ) . (51) \widehat{u}_{\mathbf{k}}(t)=\sum_{\ell=1}^{r}\lambda_{\ell}T_{\ell}(t)\prod_{j=1}^{d}a_{\ell,k_j}^{(j)}. \tag*{(51)} u k ( t ) = ℓ = 1 ∑ r λ ℓ T ℓ ( t ) j = 1 ∏ d a ℓ , k j ( j ) . (51) Thus, the same CP coefficient structure used in the Fourier setting is retained, while only the underlying one-dimensional spectral basis is changed.
Legendre basis. The Legendre basis is obtained from the Jacobi family by setting α = β = 0 \alpha= \beta= 0 α = β = 0 n n n
L n ( x ) = J n ( 0 , 0 ) ( x ) , w ( x ) = 1. (52) L_n(x) = J_n^{(0,0)}(x), \qquad w(x) = 1. \tag*{(52)} L n ( x ) = J n ( 0 , 0 ) ( x ) , w ( x ) = 1. (52) The Legendre polynomials are orthogonal on [ − 1 , 1 ] [-1,1] [ − 1 , 1 ] L 2 L^2 L 2
∫ − 1 1 L m ( x ) L n ( x ) d x = 2 2 n + 1 δ m n . (53) \int_{-1}^{1} L_m(x)L_n(x)\,\mathrm{d}x = \frac{2}{2n+1}\delta_{mn}. \tag*{(53)} ∫ − 1 1 L m ( x ) L n ( x ) d x = 2 n + 1 2 δ mn . (53) Accordingly, when TensorSpectra uses a Legendre backend, the one-dimensional factor associated with the j j j
ψ ℓ ( j ) ( x j ) = ∑ k j = 0 N j − 1 a ℓ , k j ( j ) L k j ( x j ) . (54) \psi_\ell^{(j)}(x_j) = \sum_{k_j=0}^{N_j-1} a_{\ell,k_j}^{(j)}L_{k_j}(x_j). \tag*{(54)} ψ ℓ ( j ) ( x j ) = k j = 0 ∑ N j − 1 a ℓ , k j ( j ) L k j ( x j ) . (54) Thus, the CP structure remains unchanged, while the underlying one-dimensional spectral factors are interpreted as Legendre coefficients.
Chebyshev basis. The Chebyshev polynomials of the first kind provide another important special case of the Jacobi family. They are associated with α = β = − 1 / 2 \alpha= \beta= -1/2 α = β = − 1/2
T n ( x ) = cos ( n arccos x ) , x ∈ [ − 1 , 1 ] . (55) T_n(x) = \cos(n\arccos x), \qquad x \in[-1,1]. \tag*{(55)} T n ( x ) = cos ( n arccos x ) , x ∈ [ − 1 , 1 ] . (55) Under the standard Jacobi normalization, the relation between Chebyshev and Jacobi polynomials is
T n ( x ) = n ! π Γ ( n + 1 2 ) J n ( − 1 / 2 , − 1 / 2 ) ( x ) . (56) T_n(x) = \frac{n!\sqrt{\pi}}{\Gamma\left(n+\frac{1}{2}\right)}J_n^{(-1/2,-1/2)}(x). \tag*{(56)} T n ( x ) = Γ ( n + 2 1 ) n ! π J n ( − 1/2 , − 1/2 ) ( x ) . (56) The corresponding weight function is w ( x ) = ( 1 − x 2 ) − 1 / 2 w(x) = (1-x^2)^{-1/2} w ( x ) = ( 1 − x 2 ) − 1/2
∫ − 1 1 T m ( x ) T n ( x ) 1 − x 2 d x = { π , m = n = 0 , π 2 , m = n ≥ 1 , 0 , m ≠ n . (57) \int_{-1}^{1} \frac{T_m(x)T_n(x)}{\sqrt{1-x^2}}\,\mathrm{d}x =
\begin{cases}
\pi, & m=n=0, \\
\frac{\pi}{2}, & m=n\ge1, \\
0, & m\ne n.
\end{cases} \tag*{(57)} ∫ − 1 1 1 − x 2 T m ( x ) T n ( x ) d x = ⎩ ⎨ ⎧ π , 2 π , 0 , m = n = 0 , m = n ≥ 1 , m = n . (57) When Chebyshev polynomials are used in TensorSpectra, the corresponding one-dimensional factor is represented as
ψ ℓ ( j ) ( x j ) = ∑ k j = 0 N j − 1 a ℓ , k j ( j ) T k j ( x j ) . (58) \psi_\ell^{(j)}(x_j) = \sum_{k_j=0}^{N_j-1} a_{\ell,k_j}^{(j)}T_{k_j}(x_j). \tag*{(58)} ψ ℓ ( j ) ( x j ) = k j = 0 ∑ N j − 1 a ℓ , k j ( j ) T k j ( x j ) . (58) This basis is particularly convenient for representing non-periodic coordinates and, in our implementation, temporal dependence when Chebyshev-in-time discretization is employed.
De-Aliased Nonlinear Spectral Calculus
In this section, we provide the detailed de-aliasing procedures used by TensorSpectra to evaluate nonlinear spectral interactions while preserving the low-rank tensor structure, covering both the Fourier and Jacobi spectral backends.
Fourier Backend
For the Fourier backend, nonlinear interactions are evaluated using a standard padding-and-truncation procedure. Let P N → M P_{N\to M} P N → M N N N M > N M > N M > N T M → N T_{M\to N} T M → N
a ~ α ( j ) = P N → M a α ( j ) . (59) \widetilde{a}_{\alpha}^{(j)} = P_{N\to M}a_{\alpha}^{(j)}. \tag*{(59)} a α ( j ) = P N → M a α ( j ) . (59) The corresponding one-dimensional physical-space factor on the enlarged grid is then obtained through the inverse Fourier transform,
ψ ~ α ( j ) = F M − 1 a ~ α ( j ) . (60) \widetilde{\psi}_{\alpha}^{(j)} = F_M^{-1}\widetilde{a}_{\alpha}^{(j)}. \tag*{(60)} ψ α ( j ) = F M − 1 a α ( j ) . (60) Combining the transformed factors gives the padded physical-space approximation
u M ( x ) = ∑ α = 1 r λ α ∏ j = 1 d ψ ~ α ( j ) ( x j ) . (61) u_M(x) = \sum_{\alpha=1}^{r} \lambda_\alpha\prod_{j=1}^{d} \widetilde{\psi}_{\alpha}^{(j)}(x_j). \tag*{(61)} u M ( x ) = α = 1 ∑ r λ α j = 1 ∏ d ψ α ( j ) ( x j ) . (61) The nonlinear operator is evaluated on this enlarged physical grid and then transformed back to spectral space. The resulting nonlinear spectral term is truncated to the original retained bandwidth,
N ^ ( u ) = T M → N F M [ N ( F M − 1 P N → M u ^ ) ] . (62) \widehat{\mathcal{N}}(u) = \mathcal{T}_{M\to N}\mathcal{F}_M\left[\mathcal{N}\left(\mathcal{F}_M^{-1}\mathcal{P}_{N\to M}\widehat{u}\right)\right]. \tag*{(62)} N ( u ) = T M → N F M [ N ( F M − 1 P N → M u ) ] . (62) This enlarged-grid evaluation prevents unresolved nonlinear interactions from folding back into the retained modes. For quadratic nonlinearities, TensorSpectra adopts the standard 3 / 2 3/2 3/2
M ≥ 3 2 N . (63) M \ge\frac{3}{2}N. \tag*{(63)} M ≥ 2 3 N . (63) More generally, for a polynomial nonlinearity of degree s s s
M ≥ s + 1 2 N , (64) M \ge\frac{s+1}{2}N, \tag*{(64)} M ≥ 2 s + 1 N , (64) which reduces to the 3 / 2 3/2 3/2 s = 2 s=2 s = 2
Importantly, the padding and inverse transforms are applied to the one-dimensional spectral factors rather than to an explicitly assembled N d N^d N d
Jacobi Backend
For the Jacobi backend, nonlinear interactions are de-aliased through over-integration rather than Fourier zero-padding. Let J n ( α , β ) ( x ) J_n^{(\alpha,\beta)}(x) J n ( α , β ) ( x ) n n n [ − 1 , 1 ] [-1,1] [ − 1 , 1 ] w ( α , β ) ( x ) = ( 1 − x ) α ( 1 + x ) β w^{(\alpha,\beta)}(x) = (1-x)^\alpha(1+x)^\beta w ( α , β ) ( x ) = ( 1 − x ) α ( 1 + x ) β g ( x ) g(x) g ( x ) p p p
[ Π p ( α , β ) g ] n = 1 h n ( α , β ) ∫ − 1 1 g ( x ) J n ( α , β ) ( x ) w ( α , β ) ( x ) d x , 0 ≤ n ≤ p , (65) \left[\Pi_p^{(\alpha,\beta)}g\right]_n = \frac{1}{h_n^{(\alpha,\beta)}}\int_{-1}^{1} g(x)J_n^{(\alpha,\beta)}(x)w^{(\alpha,\beta)}(x)\,\mathrm{d}x,\qquad0\le n\le p, \tag*{(65)} [ Π p ( α , β ) g ] n = h n ( α , β ) 1 ∫ − 1 1 g ( x ) J n ( α , β ) ( x ) w ( α , β ) ( x ) d x , 0 ≤ n ≤ p , (65) where h n ( α , β ) h_n^{(\alpha,\beta)} h n ( α , β ) J n ( α , β ) J_n^{(\alpha,\beta)} J n ( α , β ) Q Q Q { x q } q = 1 Q \{x_q\}_{q=1}^{Q} { x q } q = 1 Q { ω q ( α , β ) } q = 1 Q \{\omega_q^{(\alpha,\beta)}\}_{q=1}^{Q} { ω q ( α , β ) } q = 1 Q
[ Π p , Q ( α , β ) g ] n = 1 h n ( α , β ) ∑ q = 1 Q ω q ( α , β ) g ( x q ) J n ( α , β ) ( x q ) . (66) \left[\Pi_{p,Q}^{(\alpha,\beta)}g\right]_n = \frac{1}{h_n^{(\alpha,\beta)}}\sum_{q=1}^{Q}\omega_q^{(\alpha,\beta)}g(x_q)J_n^{(\alpha,\beta)}(x_q). \tag*{(66)} [ Π p , Q ( α , β ) g ] n = h n ( α , β ) 1 q = 1 ∑ Q ω q ( α , β ) g ( x q ) J n ( α , β ) ( x q ) . (66) Thus, nonlinear terms are first evaluated at an over-integrated set of Gauss–Jacobi nodes and then projected back to the retained degree-p p p
The quadrature order Q Q Q Q Q Q 2 Q − 1 2Q-1 2 Q − 1 p p p s s s g ( x ) J n ( α , β ) ( x ) g(x)J_n^{(\alpha,\beta)}(x) g ( x ) J n ( α , β ) ( x ) ( s + 1 ) p (s+1)p ( s + 1 ) p
Q ≥ ⌈ ( s + 1 ) p + 1 2 ⌉ . (67) Q \ge\left\lceil\frac{(s+1)p+1}{2}\right\rceil. \tag*{(67)} Q ≥ ⌈ 2 ( s + 1 ) p + 1 ⌉ . (67) For quadratic nonlinearities, this reduces to
Q ≥ ⌈ 3 p + 1 2 ⌉ . (68) Q \ge\left\lceil\frac{3p+1}{2}\right\rceil. \tag*{(68)} Q ≥ ⌈ 2 3 p + 1 ⌉ . (68) Within TensorSpectra, the nonlinear products are formed from one-dimensional factors and projected independently along each coordinate. The resulting projected factors are then recombined through their outer products, so the nonlinear evaluation can be carried out without explicitly forming the full N d N^d N d Q Q Q
Low-Rank Nonlinear Interaction and Rank Control
Although the Fourier and Jacobi backends use different de-aliasing mechanisms, nonlinear interactions share the same low-rank structure at the level of CP factors. Consider two rank-one components of the physical-space representation,
u α ( x ) = λ α ∏ j = 1 d ψ α ( j ) ( x j ) , u β ( x ) = λ β ∏ j = 1 d ψ β ( j ) ( x j ) . (69) u_{\alpha}(x) = \lambda_{\alpha} \prod_{j=1}^{d} \psi_{\alpha}^{(j)}(x_j), \qquad u_{\beta}(x) = \lambda_{\beta} \prod_{j=1}^{d} \psi_{\beta}^{(j)}(x_j). \tag*{(69)} u α ( x ) = λ α j = 1 ∏ d ψ α ( j ) ( x j ) , u β ( x ) = λ β j = 1 ∏ d ψ β ( j ) ( x j ) . (69) Their product remains separable across spatial coordinates,
u α ( x ) u β ( x ) = λ α λ β ∏ j = 1 d [ ψ α ( j ) ( x j ) ψ β ( j ) ( x j ) ] . (70) u_{\alpha}(x)u_{\beta}(x) = \lambda_{\alpha}\lambda_{\beta} \prod_{j=1}^{d} \left[\psi_{\alpha}^{(j)}(x_j)\psi_{\beta}^{(j)}(x_j)\right]. \tag*{(70)} u α ( x ) u β ( x ) = λ α λ β j = 1 ∏ d [ ψ α ( j ) ( x j ) ψ β ( j ) ( x j ) ] . (70) Hence, nonlinear products can be evaluated through one-dimensional factor products without forming a dense multidimensional tensor. For the quadratic term u 2 u^2 u 2 r 2 r^2 r 2 u 2 u^2 u 2 ( α , β ) (\alpha,\beta) ( α , β ) α ≤ β \alpha\le\beta α ≤ β
r q u a d = r ( r + 1 ) 2 , (71) r_{\mathrm{quad}} = \frac{r(r+1)}{2}, \tag*{(71)} r quad = 2 r ( r + 1 ) , (71) rather than r 2 r^2 r 2
PDE-Specific Spectral Objectives
In this section, we provide the PDE-specific instantiations of the spectral training objective introduced in the main text. For each benchmark, we derive the corresponding spectral residual, identify the characteristic operator scale used to construct the weighting w k w_k w k
Poisson Equation
We first consider the periodic Poisson equation used in our linear benchmarks,
Δ u ( x ) = f ( x ) , x ∈ [ 0 , 2 π ] d . (72) \Delta u(x) = f(x), \qquad x \in[0,2\pi]^d. \tag*{(72)} Δ u ( x ) = f ( x ) , x ∈ [ 0 , 2 π ] d . (72) Under the Fourier representation, the Laplace operator acts diagonally on each spectral mode, with Δ u ^ k = − ∥ k ∥ 2 2 u ^ k \widehat{\Delta u}_k = -\lVert k\rVert_2^2 \widehat{u}_k Δ u k = − ∥ k ∥ 2 2 u k k k k
R k = − ∥ k ∥ 2 2 u ^ k − f ^ k . (73) R_k = -\lVert k\rVert_2^2 \widehat{u}_k - \widehat{f}_k. \tag*{(73)} R k = − ∥ k ∥ 2 2 u k − f k . (73) The characteristic spectral scale of the Poisson operator is thus determined by the Laplacian symbol,
σ Δ ( k ) = ∥ k ∥ 2 2 , w k = ρ ( ∥ k ∥ 2 2 ) , (74) \sigma_{\Delta}(k) = \lVert k\rVert_2^2, \qquad w_k = \rho(\lVert k\rVert_2^2), \tag*{(74)} σ Δ ( k ) = ∥ k ∥ 2 2 , w k = ρ (∥ k ∥ 2 2 ) , (74) where ρ ( ⋅ ) \rho(\cdot) ρ ( ⋅ )
Accordingly, TensorSpectra minimizes the operator-aware spectral objective
L P o i s s o n = 1 ∣ K ∗ ∣ ∑ k ∈ K ∗ w k 2 ∣ − ∥ k ∥ 2 2 u ^ k − f ^ k ∣ 2 , (75) \mathcal{L}_{\mathrm{Poisson}} = \frac{1}{|\mathcal{K}_*|} \sum_{k\in\mathcal{K}_*} w_k^2 \left|-\lVert k\rVert_2^2 \widehat{u}_k - \widehat{f}_k\right|^2, \tag*{(75)} L Poisson = ∣ K ∗ ∣ 1 k ∈ K ∗ ∑ w k 2 − ∥ k ∥ 2 2 u k − f k 2 , (75) where K ∗ \mathcal{K}_* K ∗
Heat Equation
We next consider the d d d
∂ u ∂ t ( x , t ) = Δ u ( x , t ) , ( x , t ) ∈ [ 0 , 2 π ] d × [ 0 , T ] . (76) \frac{\partial u}{\partial t}(x,t) = \Delta u(x,t), \qquad(x,t) \in[0,2\pi]^d \times[0,T]. \tag*{(76)} ∂ t ∂ u ( x , t ) = Δ u ( x , t ) , ( x , t ) ∈ [ 0 , 2 π ] d × [ 0 , T ] . (76) Under the Fourier representation, the Laplace operator satisfies Δ u ^ k = − ∥ k ∥ 2 2 u ^ k \widehat{\Delta u}_k=-\lVert k\rVert_2^2\widehat{u}_k Δ u k = − ∥ k ∥ 2 2 u k
∂ t u ^ k ( t ) + ∥ k ∥ 2 2 u ^ k ( t ) = 0 , (77) \partial_t\widehat{u}_k(t)+\lVert k\rVert_2^2\widehat{u}_k(t)=0, \tag*{(77)} ∂ t u k ( t ) + ∥ k ∥ 2 2 u k ( t ) = 0 , (77) and the corresponding spectral residual is
R k ( t ) = ∂ t u ^ k ( t ) + ∥ k ∥ 2 2 u ^ k ( t ) . (78) R_k(t)=\partial_t\widehat{u}_k(t)+\lVert k\rVert_2^2\widehat{u}_k(t). \tag*{(78)} R k ( t ) = ∂ t u k ( t ) + ∥ k ∥ 2 2 u k ( t ) . (78) Following the operator-aware construction in the main text, the spatial operator scale is determined by the Laplacian symbol,
σ H e a t ( k ) = ∥ k ∥ 2 2 , w k = ρ ( ∥ k ∥ 2 2 ) . (79) \sigma_{\mathrm{Heat}}(k)=\lVert k\rVert_2^2,\qquad w_k=\rho\left(\lVert k\rVert_2^2\right). \tag*{(79)} σ Heat ( k ) = ∥ k ∥ 2 2 , w k = ρ ( ∥ k ∥ 2 2 ) . (79) The resulting operator-aware PDE objective is therefore
L H e a t P D E = 1 ∣ S ∣ ∑ ( k , t ) ∈ S w k 2 ∣ ∂ t u ^ k ( t ) + ∥ k ∥ 2 2 u ^ k ( t ) ∣ 2 . (80) \mathcal{L}_{\mathrm{Heat}}^{\mathrm{PDE}}=\frac{1}{\lvert S\rvert}\sum_{(k,t)\in S}w_k^2\left\lvert\partial_t\widehat{u}_k(t)+\lVert k\rVert_2^2\widehat{u}_k(t)\right\rvert^2. \tag*{(80)} L Heat PDE = ∣ S ∣ 1 ( k , t ) ∈ S ∑ w k 2 ∂ t u k ( t ) + ∥ k ∥ 2 2 u k ( t ) 2 . (80) The initial condition is imposed in spectral space by matching the predicted coefficients at t = 0 t=0 t = 0 u ^ k 0 \widehat{u}_k^0 u k 0
L I C = 1 ∣ K ∣ ∑ k ∈ K ∣ u ^ k ( 0 ) − u ^ k 0 ∣ 2 , (81) \mathcal{L}_{\mathrm{IC}}=\frac{1}{\lvert\mathcal{K}\rvert}\sum_{k\in\mathcal{K}}\left\lvert\widehat{u}_k(0)-\widehat{u}_k^0\right\rvert^2, \tag*{(81)} L IC = ∣ K ∣ 1 k ∈ K ∑ u k ( 0 ) − u k 0 2 , (81) leading to the overall Heat-equation objective
L H e a t = L H e a t P D E + λ I C L I C . (82) \mathcal{L}_{\mathrm{Heat}}=\mathcal{L}_{\mathrm{Heat}}^{\mathrm{PDE}}+\lambda_{\mathrm{IC}}\mathcal{L}_{\mathrm{IC}}. \tag*{(82)} L Heat = L Heat PDE + λ IC L IC . (82) In the benchmark used in this work, the initial spectral coefficients are inherited from the Poisson solution introduced above, and each Fourier mode subsequently evolves independently under the Heat operator.
Wave Equation
We next consider the d d d
∂ 2 u ∂ t 2 ( x , t ) = Δ u ( x , t ) , ( x , t ) ∈ [ 0 , 2 π ] d × [ 0 , T ] . (83) \frac{\partial^2u}{\partial t^2}(x,t)=\Delta u(x,t),\qquad(x,t)\in[0,2\pi]^d\times[0,T]. \tag*{(83)} ∂ t 2 ∂ 2 u ( x , t ) = Δ u ( x , t ) , ( x , t ) ∈ [ 0 , 2 π ] d × [ 0 , T ] . (83) Under the Fourier representation, the Laplace operator acts diagonally on each spatial mode. Hence, every Fourier coefficient satisfies the decoupled second-order evolution equation
∂ t t u ^ k ( t ) + ∥ k ∥ 2 2 u ^ k ( t ) = 0 , (84) \partial_{tt}\widehat{u}_k(t)+\lVert k\rVert_2^2\widehat{u}_k(t)=0, \tag*{(84)} ∂ tt u k ( t ) + ∥ k ∥ 2 2 u k ( t ) = 0 , (84) and the corresponding spectral residual is
R k ( t ) = ∂ t t u ^ k ( t ) + ∥ k ∥ 2 2 u ^ k ( t ) . (85) R_k(t)=\partial_{tt}\widehat{u}_k(t)+\lVert k\rVert_2^2\widehat{u}_k(t). \tag*{(85)} R k ( t ) = ∂ tt u k ( t ) + ∥ k ∥ 2 2 u k ( t ) . (85) Following the operator-aware construction, the characteristic spatial scale is determined by the Fourier symbol of the Laplacian,
σ W a v e ( k ) = ∥ k ∥ 2 2 , w k = ρ ( ∥ k ∥ 2 2 ) . (86) \sigma_{\mathrm{Wave}}(k)=\lVert k\rVert_2^2,\qquad w_k=\rho\left(\lVert k\rVert_2^2\right). \tag*{(86)} σ Wave ( k ) = ∥ k ∥ 2 2 , w k = ρ ( ∥ k ∥ 2 2 ) . (86) The resulting operator-aware PDE objective is
L W a v e P D E = 1 ∣ S ∣ ∑ ( k , t ) ∈ S w k 2 ∣ ∂ t t u ^ k ( t ) + ∥ k ∥ 2 2 u ^ k ( t ) ∣ 2 . (87) \mathcal{L}_{\mathrm{Wave}}^{\mathrm{PDE}}=\frac{1}{\lvert S\rvert}\sum_{(k,t)\in S}w_k^2\left\lvert\partial_{tt}\widehat{u}_k(t)+\lVert k\rVert_2^2\widehat{u}_k(t)\right\rvert^2. \tag*{(87)} L Wave PDE = ∣ S ∣ 1 ( k , t ) ∈ S ∑ w k 2 ∂ tt u k ( t ) + ∥ k ∥ 2 2 u k ( t ) 2 . (87) Because the Wave equation is second order in time, both the initial displacement and initial velocity must be prescribed. In our benchmark, the initial displacement is inherited from the Poisson solution and the initial velocity is set to zero. Accordingly, we define
L I C = 1 ∣ K ∣ ∑ k ∈ K ∣ u ^ k ( 0 ) − u ^ k 0 ∣ 2 + 1 ∣ K ∣ ∑ k ∈ K ∣ ∂ t u ^ k ( 0 ) − v ^ k 0 ∣ 2 , (88) \mathcal{L}_{\mathrm{IC}}=\frac{1}{\lvert\mathcal{K}\rvert}\sum_{k\in\mathcal{K}}\left\lvert\widehat{u}_k(0)-\widehat{u}_k^0\right\rvert^2+\frac{1}{\lvert\mathcal{K}\rvert}\sum_{k\in\mathcal{K}}\left\lvert\partial_t\widehat{u}_k(0)-\widehat{v}_k^0\right\rvert^2, \tag*{(88)} L IC = ∣ K ∣ 1 k ∈ K ∑ u k ( 0 ) − u k 0 2 + ∣ K ∣ 1 k ∈ K ∑ ∂ t u k ( 0 ) − v k 0 2 , (88) where v ^ k 0 = 0 \widehat{v}_k^0=0 v k 0 = 0
L W a v e = L W a v e P D E + λ I C L I C . (89) \mathcal{L}_{\mathrm{Wave}}=\mathcal{L}_{\mathrm{Wave}}^{\mathrm{PDE}}+\lambda_{\mathrm{IC}}\mathcal{L}_{\mathrm{IC}}. \tag*{(89)} L Wave = L Wave PDE + λ IC L IC . (89) Steady-State Allen–Cahn Equation
We next consider the steady-state Allen–Cahn equation on the non-periodic domain [ − 1 , 1 ] d [-1,1]^d [ − 1 , 1 ] d
− ε 2 Δ u ( x ) + u 3 ( x ) − u ( x ) = f ( x ) , x ∈ [ − 1 , 1 ] d , (90) -\varepsilon^2 \Delta u(x) + u^3(x) - u(x) = f(x), \quad x \in[-1,1]^d, \tag*{(90)} − ε 2 Δ u ( x ) + u 3 ( x ) − u ( x ) = f ( x ) , x ∈ [ − 1 , 1 ] d , (90) subject to the homogeneous Dirichlet boundary condition u ( x ) = 0 u(x)=0 u ( x ) = 0 ∂ [ − 1 , 1 ] d \partial[-1,1]^d ∂ [ − 1 , 1 ] d
In coefficient space, the corresponding spectral residual is
R ^ A C = − ε 2 Δ u ^ + u 3 ^ − u ^ − F ^ , (91) \widehat{\mathcal{R}}_{\mathrm{AC}}=-\varepsilon^2\widehat{\Delta u}+\widehat{u^3}-\widehat{u}-\widehat{F}, \tag*{(91)} R AC = − ε 2 Δ u + u 3 − u − F , (91) where u 3 ^ \widehat{u^3} u 3 k \boldsymbol{k} k
R k A C = − ε 2 [ Δ u ^ ] k + [ u 3 ^ ] k − u ^ k − f ^ k . (92) \mathcal{R}^{\mathrm{AC}}_{\boldsymbol{k}}=-\varepsilon^2\left[\widehat{\Delta u}\right]_{\boldsymbol{k}}+\left[\widehat{u^3}\right]_{\boldsymbol{k}}-\widehat{u}_{\boldsymbol{k}}-\widehat{f}_{\boldsymbol{k}}. \tag*{(92)} R k AC = − ε 2 [ Δ u ] k + [ u 3 ] k − u k − f k . (92) Unlike the Fourier cases, the Legendre differentiation operator does not admit a diagonal mode-wise symbol. We therefore construct the operator scale directly from the one-dimensional coefficient-space differentiation matrix. Let D ( 2 ) \mathbf{D}^{(2)} D ( 2 ) e n \boldsymbol{e}_n e n n n n η n = ∥ D ( 2 ) e n ∥ 2 \eta_n=\lVert\mathbf{D}^{(2)}\boldsymbol{e}_n\rVert_2 η n = ∥ D ( 2 ) e n ∥ 2 k = ( k 1 , … , k d ) \boldsymbol{k}=(k_1,\ldots,k_d) k = ( k 1 , … , k d ) − ε 2 Δ − I -\varepsilon^2\Delta-I − ε 2 Δ − I
σ A C ( k ) = ( 1 + ε 4 ∑ j = 1 d η k j 2 ) 1 / 2 . (93) \sigma_{\mathrm{AC}}(\boldsymbol{k})=\left(1+\varepsilon^4\sum_{j=1}^{d}\eta_{k_j}^2\right)^{1/2}. \tag*{(93)} σ AC ( k ) = ( 1 + ε 4 j = 1 ∑ d η k j 2 ) 1/2 . (93) The corresponding operator-aware spectral weight is chosen as an inverse power of this scale,
w k = σ A C ( k ) − p = ( 1 + ε 4 ∑ j = 1 d η k j 2 ) − p / 2 , 0 < p ≤ 1. (94) w_{\boldsymbol{k}}=\sigma_{\mathrm{AC}}(\boldsymbol{k})^{-p}=\left(1+\varepsilon^4\sum_{j=1}^{d}\eta_{k_j}^2\right)^{-p/2},\quad0<p\leq1. \tag*{(94)} w k = σ AC ( k ) − p = ( 1 + ε 4 j = 1 ∑ d η k j 2 ) − p /2 , 0 < p ≤ 1. (94) The weight is constructed from the linear operator, while the full residual includes the cubic term. The resulting PDE-aware objective is
L A C = 1 ∣ K ∣ ∑ k ∈ K w k 2 ∣ R k A C ∣ 2 . (95) \mathcal{L}_{\mathrm{AC}}=\frac{1}{|\mathcal{K}|}\sum_{\boldsymbol{k}\in\mathcal{K}}w_{\boldsymbol{k}}^2\left|\mathcal{R}^{\mathrm{AC}}_{\boldsymbol{k}}\right|^2. \tag*{(95)} L AC = ∣ K ∣ 1 k ∈ K ∑ w k 2 R k AC 2 . (95) Because the homogeneous Dirichlet boundary condition is enforced directly through the factor-wise parameterization, no additional boundary penalty is required in this benchmark.
Nonlinear Reaction Equation
We finally consider the nonlinear logistic reaction equation
∂ u ∂ t ( x , t ) − ρ u ( x , t ) ( 1 − u ( x , t ) ) = 0 , ( x , t ) ∈ [ 0 , 2 π ] d × [ 0 , T ] , (96) \frac{\partial u}{\partial t}(x,t)-\rho u(x,t)(1-u(x,t))=0,\quad(x,t)\in[0,2\pi]^d\times[0,T], \tag*{(96)} ∂ t ∂ u ( x , t ) − ρ u ( x , t ) ( 1 − u ( x , t )) = 0 , ( x , t ) ∈ [ 0 , 2 π ] d × [ 0 , T ] , (96) with a prescribed initial condition u ( x , 0 ) = h ( x ) u(x,0)=h(x) u ( x , 0 ) = h ( x ) h K h_K h K u ( x , t ) = h K ( x ) + t v ( x , t ) u(x,t)=h_K(x)+tv(x,t) u ( x , t ) = h K ( x ) + t v ( x , t )
The hard initial-condition parameterization is
u ( x , t ) = h K ( x ) + t v ( x , t ) . (97) u(x,t)=h_K(x)+tv(x,t). \tag*{(97)} u ( x , t ) = h K ( x ) + t v ( x , t ) . (97) Expanding the nonlinear term gives ∂ t u − ρ u + ρ u 2 = 0 \partial_tu-\rho u+\rho u^2=0 ∂ t u − ρ u + ρ u 2 = 0
R ^ R e a c t = ∂ t u ^ − ρ u ^ + ρ u 2 ^ , (98) \widehat{\mathcal{R}}_{\mathrm{React}}=\partial_t\widehat{u}-\rho\widehat{u}+\rho\widehat{u^2}, \tag*{(98)} R React = ∂ t u − ρ u + ρ u 2 , (98) where u 2 ^ \widehat{u^2} u 2 k = ( k 1 , … , k d ) \boldsymbol{k}=(k_1,\ldots,k_d) k = ( k 1 , … , k d ) q q q
R k , q R e a c t = [ ∂ t u ^ ] k , q − ρ u ^ k , q + ρ [ u 2 ^ ] k , q . (99) \mathcal{R}^{\mathrm{React}}_{\boldsymbol{k},q}=[\partial_t\widehat{u}]_{\boldsymbol{k},q}-\rho\widehat{u}_{\boldsymbol{k},q}+\rho\left[\widehat{u^2}\right]_{\boldsymbol{k},q}. \tag*{(99)} R k , q React = [ ∂ t u ] k , q − ρ u k , q + ρ [ u 2 ] k , q . (99) To construct an explicit operator-aware scale, we consider the linear part of the reaction operator, ∂ t − ρ I \partial_t-\rho I ∂ t − ρ I D t \mathbf{D}_t D t e q \boldsymbol{e}_q e q q q q
η q = ∥ D t e q ∥ 2 . (100) \eta_q=\left\|\mathbf{D}_t\boldsymbol{e}_q\right\|_2. \tag*{(100)} η q = ∥ D t e q ∥ 2 . (100) For the shifted orthonormal Legendre basis on [ 0 , T ] [0,T] [ 0 , T ]
ψ q ( t ) ( t ) = 2 q + 1 T L q ( 2 t T − 1 ) , (101) \psi_q^{(t)}(t)=\sqrt{\frac{2q+1}{T}}L_q\left(\frac{2t}{T}-1\right), \tag*{(101)} ψ q ( t ) ( t ) = T 2 q + 1 L q ( T 2 t − 1 ) , (101) the Euclidean norm of D t e q \mathbf{D}_t\boldsymbol{e}_q D t e q L 2 ( 0 , T ) L^2(0,T) L 2 ( 0 , T ) ∫ − 1 1 [ L q ′ ( ξ ) ] 2 d ξ = q ( q + 1 ) \int_{-1}^{1}[L_q'(\xi)]^2\mathrm{d}\xi=q(q+1) ∫ − 1 1 [ L q ′ ( ξ ) ] 2 d ξ = q ( q + 1 )
η q 2 = ∥ ∂ t ψ q ( t ) ∥ L 2 ( 0 , T ) 2 = 2 q ( q + 1 ) ( 2 q + 1 ) T 2 . (102) \eta_q^2=\left\|\partial_t\psi_q^{(t)}\right\|_{L^2(0,T)}^2=\frac{2q(q+1)(2q+1)}{T^2}. \tag*{(102)} η q 2 = ∂ t ψ q ( t ) L 2 ( 0 , T ) 2 = T 2 2 q ( q + 1 ) ( 2 q + 1 ) . (102) Because differentiation lowers the Legendre polynomial degree, D t e q \mathbf{D}_t\boldsymbol{e}_q D t e q e q \boldsymbol{e}_q e q ∂ t − ρ I \partial_t-\rho I ∂ t − ρ I q q q
σ R e a c t ( q ) = ∥ ( D t − ρ I ) e q ∥ 2 = ( ρ 2 + η q 2 ) 1 / 2 . (103) \sigma_{\mathrm{React}}(q)=\left\|(\mathbf{D}_t-\rho\mathbf{I})\boldsymbol{e}_q\right\|_2=\left(\rho^2+\eta_q^2\right)^{1/2}. \tag*{(103)} σ React ( q ) = ∥ ( D t − ρ I ) e q ∥ 2 = ( ρ 2 + η q 2 ) 1/2 . (103) Substituting the explicit Legendre differentiation scale gives
σ R e a c t ( q ) = [ ρ 2 + 2 q ( q + 1 ) ( 2 q + 1 ) T 2 ] 1 / 2 . (104) \sigma_{\mathrm{React}}(q)=\left[\rho^2+\frac{2q(q+1)(2q+1)}{T^2}\right]^{1/2}. \tag*{(104)} σ React ( q ) = [ ρ 2 + T 2 2 q ( q + 1 ) ( 2 q + 1 ) ] 1/2 . (104) We normalize this scale by its lowest-order value ρ \rho ρ
w k , q = w q = ( σ R e a c t ( q ) ρ ) − p = [ 1 + 2 q ( q + 1 ) ( 2 q + 1 ) ρ 2 T 2 ] − p / 2 , 0 < p ≤ 1. (105) w_{\boldsymbol{k},q}=w_q=\left(\frac{\sigma_{\mathrm{React}}(q)}{\rho}\right)^{-p}=\left[1+\frac{2q(q+1)(2q+1)}{\rho^2T^2}\right]^{-p/2},\qquad0<p\leq1. \tag*{(105)} w k , q = w q = ( ρ σ React ( q ) ) − p = [ 1 + ρ 2 T 2 2 q ( q + 1 ) ( 2 q + 1 ) ] − p /2 , 0 < p ≤ 1. (105) Since the reaction equation contains no spatial differential operator, the preconditioning weight depends only on the temporal spectral index q q q ρ u 2 \rho u^2 ρ u 2
L R e a c t = 1 ∣ K ∣ ∣ Q ∣ ∑ k ∈ K ∑ q ∈ Q w q 2 ∣ R k , q R e a c t ∣ 2 , (106) \mathcal{L}_{\mathrm{React}}=\frac{1}{|\mathcal{K}||\mathcal{Q}|}\sum_{\boldsymbol{k}\in\mathcal{K}}\sum_{q\in\mathcal{Q}}w_q^2\left|\mathcal{R}^{\mathrm{React}}_{\boldsymbol{k},q}\right|^2, \tag*{(106)} L React = ∣ K ∣∣ Q ∣ 1 k ∈ K ∑ q ∈ Q ∑ w q 2 R k , q React 2 , (106) where K \mathcal{K} K Q \mathcal{Q} Q
Baseline Methods
In this section, we summarize the main ideas of the baseline methods used for comparison.
Physics-Informed Neural Network (PINN). PINN approximates the PDE solution with a neural network and trains it by minimizing the governing-equation residual together with boundary- and initial-condition errors. The required derivatives are computed through automatic differentiation at collocation points, providing a general physics-informed framework for both stationary and time-dependent problems.
Separable Physics-Informed Neural Network (SPINN). SPINN represents the solution as a low-rank sum of products of one-dimensional neural-network outputs, with a separate subnetwork for each coordinate. This separable architecture reduces network evaluations on tensor-product collocation grids and enables efficient forward-mode automatic differentiation, while retaining a physics-informed objective based on PDE residuals and auxiliary constraints.
Physics-Informed Residual Adaptive Network (PirateNet). PirateNet improves the trainability of deep PINNs through random Fourier embeddings, gated feature transformations, and adaptive residual connections. Its residual weights are initialized to zero so that the network starts from a shallow representation and learns to activate nonlinear blocks during training. The model retains the physics-informed objective while using this architecture and initialization to facilitate optimization of PDE residuals involving higher-order derivatives.
Spectral-Informed Neural Network (SINN). SINN learns spectral coefficients from frequency indices and, for time-dependent problems, time, rather than predicting solution values directly from spatial coordinates. It enforces the governing equation in the spectral domain, where spatial differentiation becomes algebraic operations on coefficients. For high-dimensional problems, we use the modified SINN variant with coefficient-decay scaling and learnable basis embeddings to incorporate prior spectral structure into the representation.
Tensor Neural Network (TNN). TNN approximates high-dimensional solutions using a low-rank tensor-product representation with normalized one-dimensional neural factors. This structure decomposes high-dimensional integrals into combinations of one-dimensional integrals, enabling efficient numerical integration. The PDE-solving framework combines this representation with variational or a posteriori error-estimator objectives, using the latter to control approximation error in the energy norm.
Low-Rank Fourier Sums (LoRFS). LoRFS represents the solution as a low-rank sum of products of finite one-dimensional Fourier expansions and directly optimizes their coefficients. It combines this representation with continuous integral-based physics objectives, exploiting Fourier orthogonality and separability to evaluate residual integrals deterministically without high-dimensional stochastic collocation. Its central idea is to pair an explicit low-rank spectral representation with efficient evaluation of the physics-informed loss.
Further Theoretical Results
In this section, we provide the complete proofs and technical details supporting the theoretical results in Section § 4 , including the complexity analysis, spectral and low-rank approximation bounds, residual stability estimates, staged correction analysis, and de-aliasing conditions.
Norms, Notation, and CP Rank
For Fourier results we equip T d = [ 0 , 2 π ] d \mathbb{T}^{d} = [0,2\pi]^{d} T d = [ 0 , 2 π ] d ( 2 π ) − d d x (2\pi)^{-d}\mathrm{d}x ( 2 π ) − d d x e i k ⋅ x e^{i k\cdot x} e ik ⋅ x L w 2 ( [ − 1 , 1 ] d ) L^{2}_{w}([-1,1]^{d}) L w 2 ([ − 1 , 1 ] d ) § 4 . Here K = ∏ j = 1 d K j \mathcal{K} = \prod_{j=1}^{d}\mathcal{K}_{j} K = ∏ j = 1 d K j ∣ K j ∣ = N j |\mathcal{K}_{j}| = N_{j} ∣ K j ∣ = N j S N S_{N} S N
∥ S N v ^ ∥ H 2 = ∑ k ∈ K ∣ v ^ k ∣ 2 = ∥ v ^ ∥ 2 2 . (107) \left\|S_{N}\widehat{v}\right\|_{H}^{2} = \sum_{k\in\mathcal{K}}\left|\widehat{v}_{k}\right|^{2} = \left\|\widehat{v}\right\|_{2}^{2}. \tag*{(107)} ∥ S N v ∥ H 2 = k ∈ K ∑ ∣ v k ∣ 2 = ∥ v ∥ 2 2 . (107) In another solution norm X X X ∥ S N v ^ ∥ X ≤ C N , X ∥ v ^ ∥ F \left\|S_{N}\widehat{v}\right\|_{X} \le C_{N,X}\left\|\widehat{v}\right\|_{F} ∥ S N v ∥ X ≤ C N , X ∥ v ∥ F C N , X C_{N,X} C N , X N N N d d d
Proposition E.1 (Why an infimum is necessary). For order-three tensors over R 2 \mathbb{R}^{2} R 2
Proof. Let e 1 , e 2 e_{1},e_{2} e 1 , e 2 v ε = e 1 + ε e 2 v_{\varepsilon}=e_{1}+\varepsilon e_{2} v ε = e 1 + ε e 2
T ε = ε − 1 ( v ε ∘ v ε ∘ v ε − e 1 ∘ e 1 ∘ e 1 ) . (108) T_{\varepsilon}=\varepsilon^{-1}\left(v_{\varepsilon}\circ v_{\varepsilon}\circ v_{\varepsilon}-e_{1}\circ e_{1}\circ e_{1}\right). \tag*{(108)} T ε = ε − 1 ( v ε ∘ v ε ∘ v ε − e 1 ∘ e 1 ∘ e 1 ) . (108) Each T ε T_{\varepsilon} T ε ε → 0 \varepsilon\to0 ε → 0
T = e 2 ∘ e 1 ∘ e 1 + e 1 ∘ e 2 ∘ e 1 + e 1 ∘ e 1 ∘ e 2 . (109) T=e_{2}\circ e_{1}\circ e_{1}+e_{1}\circ e_{2}\circ e_{1}+e_{1}\circ e_{1}\circ e_{2}. \tag*{(109)} T = e 2 ∘ e 1 ∘ e 1 + e 1 ∘ e 2 ∘ e 1 + e 1 ∘ e 1 ∘ e 2 . (109) Its two slices along the first mode are
A = ( 0 1 1 0 ) , B = ( 1 0 0 0 ) . (110) A=\begin{pmatrix}0 & 1\\1 & 0\end{pmatrix},\quad B=\begin{pmatrix}1 & 0\\0 & 0\end{pmatrix}. \tag*{(110)} A = ( 0 1 1 0 ) , B = ( 1 0 0 0 ) . (110) If T T T A A A A − 1 B A^{-1}B A − 1 B
A − 1 B = ( 0 0 1 0 ) (111) A^{-1}B=\begin{pmatrix}0 & 0\\1 & 0\end{pmatrix} \tag*{(111)} A − 1 B = ( 0 1 0 0 ) (111) is a nonzero nilpotent Jordan block. Hence T T T □ \square □
Complexity Proofs
Proof of Proposition 4.1. There is one dense coefficient for each multi-index in K = ∏ j = 1 d K j \mathcal{K} = \prod_{j=1}^{d} \mathcal{K}_{j} K = ∏ j = 1 d K j ∏ j N j \prod_{j} N_{j} ∏ j N j r r r r N j rN_{j} r N j j j j r + r ∑ j N j r + r\sum_{j}N_{j} r + r ∑ j N j r r r j j j r C s p ( N j ) rC_{\mathrm{sp}}(N_{j}) r C sp ( N j ) j j j i k j ik_{j} i k j C s p ( N ) = O ( N ) C_{\mathrm{sp}}(N) = O(N) C sp ( N ) = O ( N ) O ( N 2 ) O(N^{2}) O ( N 2 ) L L L L r Lr L r O ( L d r N ) O(LdrN) O ( L d r N )
Proof of Proposition 4.2. Writing
u = ∑ α = 1 r λ α ∏ j = 1 d ψ α ( j ) (112) u = \sum_{\alpha=1}^{r} \lambda_{\alpha}\prod_{j=1}^{d}\psi_{\alpha}^{(j)} \tag*{(112)} u = α = 1 ∑ r λ α j = 1 ∏ d ψ α ( j ) (112) gives
u 2 = ∑ 1 ≤ α ≤ β ≤ r ( 2 − δ α β ) λ α λ β ∏ j ( ψ α ( j ) ψ β ( j ) ) . (113) u^{2} = \sum_{1\leq\alpha\leq\beta\leq r}(2-\delta_{\alpha\beta})\lambda_{\alpha}\lambda_{\beta}\prod_{j}\left(\psi_{\alpha}^{(j)}\psi_{\beta}^{(j)}\right). \tag*{(113)} u 2 = 1 ≤ α ≤ β ≤ r ∑ ( 2 − δ α β ) λ α λ β j ∏ ( ψ α ( j ) ψ β ( j ) ) . (113) There are P 2 = r ( r + 1 ) / 2 P_{2}=r(r+1)/2 P 2 = r ( r + 1 ) /2 r r r O ( d r M log M ) O(drM\log M) O ( d r M log M ) O ( d P 2 M log M ) O(dP_{2}M\log M) O ( d P 2 M log M ) N N N O ( d P 2 N ) O(dP_{2}N) O ( d P 2 N ) P 2 P_{2} P 2 d P 2 N dP_{2}N d P 2 N Q × ( p + 1 ) Q\times(p+1) Q × ( p + 1 ) O ( Q ( p + 1 ) ) O(Q(p+1)) O ( Q ( p + 1 )) Q Q Q ( p + 1 ) × Q (p+1)\times Q ( p + 1 ) × Q O ( Q ( p + 1 ) ) O(Q(p+1)) O ( Q ( p + 1 )) ( r + 2 3 ) \binom{r+2}{3} ( 3 r + 2 ) R R R R 2 R^{2} R 2 d d d N N N O ( d N R 2 ) O(dNR^{2}) O ( d N R 2 )
Spectral and Low-Rank Approximation Proofs
Proof of the Fourier part of Lemma 4.3. Define
∥ u ∥ H 2 s 2 = ∑ k ∈ Z d ( 1 + ∣ k ∣ 2 2 ) s ∣ u ^ k ∣ 2 . (114) \lVert u\rVert_{H^{s}_{2}}^{2}=\sum_{k\in\mathbb{Z}^{d}}(1+\lvert k\rvert_{2}^{2})^{s}\lvert\widehat{u}_{k}\rvert^{2}. \tag*{(114)} ∥ u ∥ H 2 s 2 = k ∈ Z d ∑ ( 1 + ∣ k ∣ 2 2 ) s ∣ u k ∣ 2 . (114) For k ∉ K K k\notin\mathcal{K}_{K} k ∈ / K K ∣ k j ∣ ≥ K + 1 \lvert k_{j}\rvert\geq K+1 ∣ k j ∣ ≥ K + 1 1 + ∣ k ∣ 2 2 ≥ 1 + K 2 1+\lvert k\rvert_{2}^{2}\geq1+K^{2} 1 + ∣ k ∣ 2 2 ≥ 1 + K 2 m − s < 0 m-s<0 m − s < 0
∥ u − P K u ∥ H 2 m 2 = ∑ k ∉ K K ( 1 + ∣ k ∣ 2 2 ) m ∣ u ^ k ∣ 2 ≤ ( 1 + K 2 ) m − s ∑ k ∉ K K ( 1 + ∣ k ∣ 2 2 ) s ∣ u ^ k ∣ 2 ≤ ( 1 + K 2 ) m − s ∥ u ∥ H 2 s 2 . \begin{aligned}
\lVert u-P_{K}u\rVert_{H^{m}_{2}}^{2}=\sum_{k\notin\mathcal{K}_{K}}(1+\lvert k\rvert_{2}^{2})^{m}\lvert\widehat{u}_{k}\rvert^{2} \\
&\leq(1+K^{2})^{m-s}\sum_{k\notin\mathcal{K}_{K}}(1+\lvert k\rvert_{2}^{2})^{s}\lvert\widehat{u}_{k}\rvert^{2}\leq(1+K^{2})^{m-s}\lVert u\rVert_{H^{s}_{2}}^{2}.
\end{aligned} ∥ u − P K u ∥ H 2 m 2 = k ∈ / K K ∑ ( 1 + ∣ k ∣ 2 2 ) m ∣ u k ∣ 2 ≤ ( 1 + K 2 ) m − s k ∈ / K K ∑ ( 1 + ∣ k ∣ 2 2 ) s ∣ u k ∣ 2 ≤ ( 1 + K 2 ) m − s ∥ u ∥ H 2 s 2 . Taking square roots proves (15) . Under the analytic assumption, applying the union bound over the d d d ∣ k j ∣ > K \lvert k_{j}\rvert>K ∣ k j ∣ > K
∥ u − P K u ∥ L 2 2 ≤ B 2 ∑ k ∉ K K e − 2 a ∣ k ∣ 1 ≤ d B 2 ( 2 ∑ n = K + 1 ∞ e − 2 a n ) ( ∑ n ∈ Z e − 2 a ∣ n ∣ ) d − 1 . \begin{aligned}
\lVert u-P_{K}u\rVert_{L^{2}}^{2}\leq B^{2}\sum_{k\notin\mathcal{K}_{K}}e^{-2a\lvert k\rvert_{1}} \\
&\leq dB^{2}\left(2\sum_{n=K+1}^{\infty}e^{-2an}\right)\left(\sum_{n\in\mathbb{Z}}e^{-2a\lvert n\rvert}\right)^{d-1}.
\end{aligned} ∥ u − P K u ∥ L 2 2 ≤ B 2 k ∈ / K K ∑ e − 2 a ∣ k ∣ 1 ≤ d B 2 ( 2 n = K + 1 ∑ ∞ e − 2 an ) ( n ∈ Z ∑ e − 2 a ∣ n ∣ ) d − 1 . The geometric series yields the stated C ( d , a ) B e − a K C(d,a)Be^{-aK} C ( d , a ) B e − a K
Proof of the Jacobi part of Lemma 4.3. Let ϕ n ( α , β ) = J n ( α , β ) / h n ( α , β ) \phi_{n}^{(\alpha,\beta)}=J_{n}^{(\alpha,\beta)}/\sqrt{h_{n}^{(\alpha,\beta)}} ϕ n ( α , β ) = J n ( α , β ) / h n ( α , β ) ∥ v ∥ α , β 2 = ∫ − 1 1 ∣ v ( x ) ∣ 2 ( 1 − x ) α ( 1 + x ) β d x \lVert v\rVert_{\alpha,\beta}^{2}=\int_{-1}^{1}\lvert v(x)\rvert^{2}(1-x)^{\alpha}(1+x)^{\beta}\,dx ∥ v ∥ α , β 2 = ∫ − 1 1 ∣ v ( x ) ∣ 2 ( 1 − x ) α ( 1 + x ) β d x s ≥ 1 s\geq1 s ≥ 1
∥ ∂ x s ∑ n ≥ 0 c n ϕ n ( α , β ) ∥ α + s , β + s 2 = ∑ n ≥ s Λ n , s ( α , β ) ∣ c n ∣ 2 . (115) \left\lVert\partial_{x}^{s}\sum_{n\geq0}c_{n}\phi_{n}^{(\alpha,\beta)}\right\rVert_{\alpha+s,\beta+s}^{2}
=\sum_{n\geq s}\Lambda_{n,s}^{(\alpha,\beta)}\lvert c_{n}\rvert^{2}. \tag*{(115)} ∂ x s n ≥ 0 ∑ c n ϕ n ( α , β ) α + s , β + s 2 = n ≥ s ∑ Λ n , s ( α , β ) ∣ c n ∣ 2 . (115) Here
Λ n , s ( α , β ) = ∏ ℓ = 0 s − 1 ( n − ℓ ) ( n + α + β + 1 + ℓ ) . (116) \Lambda_{n,s}^{(\alpha,\beta)}=\prod_{\ell=0}^{s-1}(n-\ell)(n+\alpha+\beta+1+\ell). \tag*{(116)} Λ n , s ( α , β ) = ℓ = 0 ∏ s − 1 ( n − ℓ ) ( n + α + β + 1 + ℓ ) . (116) For α , β > − 1 \alpha,\beta>-1 α , β > − 1 n > p ≥ s − 1 n>p\ge s-1 n > p ≥ s − 1 Λ n , s ( α , β ) ≥ c α , β , s ( p + 1 ) 2 s \Lambda_{n,s}^{(\alpha,\beta)}\ge c_{\alpha,\beta,s}(p+1)^{2s} Λ n , s ( α , β ) ≥ c α , β , s ( p + 1 ) 2 s p p p
∥ v − P p v ∥ α , β 2 = ∑ n > p ∣ c n ∣ 2 ≤ c α , β , s − 1 ( p + 1 ) − 2 s ∥ ∂ x s v ∥ α + s , β + s 2 . \lVert v-P_pv\rVert_{\alpha,\beta}^{2}=\sum_{n>p}|c_n|^2
\le c_{\alpha,\beta,s}^{-1}(p+1)^{-2s}\lVert\partial_x^s v\rVert_{\alpha+s,\beta+s}^{2}. ∥ v − P p v ∥ α , β 2 = n > p ∑ ∣ c n ∣ 2 ≤ c α , β , s − 1 ( p + 1 ) − 2 s ∥ ∂ x s v ∥ α + s , β + s 2 . Tensor projections commute and each is an L w 2 L_w^2 L w 2
I − ∏ j = 1 d P p j ( j ) = ∑ j = 1 d ( ∏ ℓ < j P p ℓ ( ℓ ) ) ( I − P p j ( j ) ) . (117) I-\prod_{j=1}^{d}P_{p_j}^{(j)}
=\sum_{j=1}^{d}\left(\prod_{\ell<j}P_{p_\ell}^{(\ell)}\right)\left(I-P_{p_j}^{(j)}\right). \tag*{(117)} I − j = 1 ∏ d P p j ( j ) = j = 1 ∑ d ℓ < j ∏ P p ℓ ( ℓ ) ( I − P p j ( j ) ) . (117) Apply the triangle inequality, contraction in the other axes, and the one-dimensional estimate to each summand. Fubini gives Equation (16) . □
Proof of Lemma 4.4. The feasible tensors of rank at most r r r r + 1 r+1 r + 1 r r r U ^ N \widehat{U}_N U N
Proof of Theorem 4.5. Using u N = P N u u_N=P_Nu u N = P N u (107) , we obtain
∥ u − u θ ∥ H ≤ ∥ u − u N ∥ H + ∥ u N − S N u ^ θ ∥ H = E s p e c ( N ) + ∥ U ^ N − u ^ θ ∥ F ≤ E s p e c ( N ) + ∥ U ^ N − V ^ η ∥ F + ∥ V ^ η − u ^ θ ∥ F ≤ E s p e c ( N ) + E r a n k ( N , r ) + η + E o p t ( θ ; V ^ η ) . (118) \begin{aligned}
\lVert u-u_\theta\rVert_H
&\le\lVert u-u_N\rVert_H+\lVert u_N-S_N\widehat{u}_\theta\rVert_H\\
&=E_{\mathrm{spec}}(N)+\lVert\widehat{U}_N-\widehat{u}_\theta\rVert_F\\
&\le E_{\mathrm{spec}}(N)+\lVert\widehat{U}_N-\widehat{V}_\eta\rVert_F+\lVert\widehat{V}_\eta-\widehat{u}_\theta\rVert_F\\
&\le E_{\mathrm{spec}}(N)+E_{\mathrm{rank}}(N,r)+\eta+E_{\mathrm{opt}}(\theta;\widehat{V}_\eta). \tag*{(118)}
\end{aligned} ∥ u − u θ ∥ H ≤ ∥ u − u N ∥ H + ∥ u N − S N u θ ∥ H = E spec ( N ) + ∥ U N − u θ ∥ F ≤ E spec ( N ) + ∥ U N − V η ∥ F + ∥ V η − u θ ∥ F ≤ E spec ( N ) + E rank ( N , r ) + η + E opt ( θ ; V η ) . (118) The second inequality applies the triangle inequality in coefficient space. □
Residual, Weighting, and PDE Stability Proofs
Proof of Theorem 4.6. For e = u θ − u ∈ X 0 e=u_\theta-u\in X_0 e = u θ − u ∈ X 0 A u = f Au=f A u = f
A e = A u θ − f = R . (119) Ae=Au_\theta-f=R. \tag*{(119)} A e = A u θ − f = R . (119) The assumed stability inequality yields Equation (18) . For the retained coefficient vector,
∣ R k ∣ ≤ w min − 1 ∣ w k R k ∣ . (120) |R_k|\le w_{\min}^{-1}|w_kR_k|. \tag*{(120)} ∣ R k ∣ ≤ w m i n − 1 ∣ w k R k ∣. (120) Consequently,
∥ R K ∥ Y ≤ w min − 1 ( ∑ k ∈ K ∣ w k R k ∣ 2 ) 1 / 2 = w min − 1 n L P D E p r e c , f u l l . (121) \begin{aligned}
\lVert R_K\rVert_Y
&\le w_{\min}^{-1}\left(\sum_{k\in K}|w_kR_k|^2\right)^{1/2}\\
&=w_{\min}^{-1}\sqrt{n\mathcal{L}_{\mathrm{PDE}}^{\mathrm{prec,full}}}. \tag*{(121)}
\end{aligned} ∥ R K ∥ Y ≤ w m i n − 1 ( k ∈ K ∑ ∣ w k R k ∣ 2 ) 1/2 = w m i n − 1 n L PDE prec , full . (121) Orthogonality and the triangle inequality give
∥ R ∥ Y ≤ ∥ R K ∥ Y + ∥ R K c ∥ Y . (122) \lVert R\rVert_Y\le\lVert R_K\rVert_Y+\lVert R_{K^c}\rVert_Y. \tag*{(122)} ∥ R ∥ Y ≤ ∥ R K ∥ Y + ∥ R K c ∥ Y . (122) Substitution proves the loss bound. If training uses S ⊆ K S\subseteq K S ⊆ K R S R_S R S p k > 0 p_k>0 p k > 0
E p [ ∣ w k R k ∣ 2 p k ] = ∑ k ∈ K ∣ w k R k ∣ 2 . (123) \mathbb{E}_p\left[\frac{|w_kR_k|^2}{p_k}\right]=\sum_{k\in K}|w_kR_k|^2. \tag*{(123)} E p [ p k ∣ w k R k ∣ 2 ] = k ∈ K ∑ ∣ w k R k ∣ 2 . (123) An uncorrected expectation gives only ∑ k ∈ K p k ∣ w k R k ∣ 2 \sum_{k\in K}p_k|w_kR_k|^2 ∑ k ∈ K p k ∣ w k R k ∣ 2 p min > 0 p_{\min}>0 p m i n > 0
Proof of Proposition 4.7. Set z = v − u N G ∈ V N z=v-u_N^G\in V_N z = v − u N G ∈ V N A N u N G = P N f A_Nu_N^G=P_Nf A N u N G = P N f A N z = P N ( A v − f ) A_Nz=P_N(Av-f) A N z = P N ( A v − f )
∥ z ∥ H ≤ β N − 1 ∥ P N ( A v − f ) ∥ H . (124) \|z\|_H\le\beta_N^{-1}\|P_N(Av-f)\|_H. \tag*{(124)} ∥ z ∥ H ≤ β N − 1 ∥ P N ( A v − f ) ∥ H . (124) The triangle inequality
∥ u − v ∥ H ≤ ∥ u − u N G ∥ H + ∥ u N G − v ∥ H (125) \|u-v\|_H\le\|u-u_N^G\|_H+\|u_N^G-v\|_H \tag*{(125)} ∥ u − v ∥ H ≤ ∥ u − u N G ∥ H + ∥ u N G − v ∥ H (125) then proves Equation (20) . A full weighted average bounds the projected residual using the weight lower bound as in Theorem 4.6.
Proof of Corollary 4.8. Since Δ e = R \Delta e=R Δ e = R
− ∣ k ∣ 2 2 e ^ k = R k , k ≠ 0. (126) -|k|_2^2\widehat{e}_k=R_k,\qquad k\ne0. \tag*{(126)} − ∣ k ∣ 2 2 e k = R k , k = 0. (126) Parseval on normalized T d \mathbb{T}^d T d
∥ e ∥ L 2 2 = ∑ k ≠ 0 ∣ R k ∣ 2 ∣ k ∣ 2 4 , (127) \|e\|_{L^2}^2=\sum_{k\ne0}\frac{|R_k|^2}{|k|_2^4}, \tag*{(127)} ∥ e ∥ L 2 2 = k = 0 ∑ ∣ k ∣ 2 4 ∣ R k ∣ 2 , (127) ∥ ∇ e ∥ L 2 2 = ∑ k ≠ 0 ∣ k ∣ 2 2 ∣ e ^ k ∣ 2 = ∑ k ≠ 0 ∣ R k ∣ 2 ∣ k ∣ 2 2 . \|\nabla e\|_{L^2}^2=\sum_{k\ne0}|k|_2^2|\widehat{e}_k|^2=\sum_{k\ne0}\frac{|R_k|^2}{|k|_2^2}. ∥∇ e ∥ L 2 2 = k = 0 ∑ ∣ k ∣ 2 2 ∣ e k ∣ 2 = k = 0 ∑ ∣ k ∣ 2 2 ∣ R k ∣ 2 . The latter is the squared homogeneous H − 1 H^{-1} H − 1 R R R ∣ k ∣ 2 − 2 |k|_2^{-2} ∣ k ∣ 2 − 2 − e ^ k -\widehat{e}_k − e k ℓ 2 \ell^2 ℓ 2 L 2 L^2 L 2
min k ∈ K ∣ w k ∣ ∣ k ∣ 2 2 ∥ e ∥ L 2 ≤ ∥ W R ∥ 2 ≤ max k ∈ K ∣ w k ∣ ∣ k ∣ 2 2 ∥ e ∥ L 2 (128) \begin{aligned}
\min_{k\in\mathcal{K}}|w_k||k|_2^2\|e\|_{L^2}\le\|WR\|_2 \\
&\le\max_{k\in\mathcal{K}}|w_k||k|_2^2\|e\|_{L^2} \tag*{(128)}
\end{aligned} k ∈ K min ∣ w k ∣∣ k ∣ 2 2 ∥ e ∥ L 2 ≤ ∥ W R ∥ 2 ≤ k ∈ K max ∣ w k ∣∣ k ∣ 2 2 ∥ e ∥ L 2 (128) on a fixed finite retained set, provided all relevant coefficients are included.
Proposition E.2 (Heat equation with initial mismatch). Let u t − Δ u = 0 u_t-\Delta u=0 u t − Δ u = 0 v t − Δ v = R v_t-\Delta v=R v t − Δ v = R T d × ( 0 , T ) \mathbb{T}^d\times(0,T) T d × ( 0 , T ) R ∈ L 2 ( 0 , T ; L 2 ) R\in L^2(0,T;L^2) R ∈ L 2 ( 0 , T ; L 2 ) e 0 = v ( 0 ) − u ( 0 ) ∈ L 2 e_0=v(0)-u(0)\in L^2 e 0 = v ( 0 ) − u ( 0 ) ∈ L 2
sup 0 ≤ t ≤ T ∥ v ( t ) − u ( t ) ∥ L 2 ≤ ∥ e 0 ∥ L 2 + T ∥ R ∥ L 2 ( 0 , T ; L 2 ) . (129) \sup_{0\le t\le T}\|v(t)-u(t)\|_{L^2}\le\|e_0\|_{L^2}+\sqrt{T}\|R\|_{L^2(0,T;L^2)}. \tag*{(129)} 0 ≤ t ≤ T sup ∥ v ( t ) − u ( t ) ∥ L 2 ≤ ∥ e 0 ∥ L 2 + T ∥ R ∥ L 2 ( 0 , T ; L 2 ) . (129) Proof. The heat semigroup is an L 2 L^2 L 2 e − ∣ k ∣ 2 2 t e^{-|k|_2^2t} e − ∣ k ∣ 2 2 t e = v − u e=v-u e = v − u
e ( t ) = e t Δ e 0 + ∫ 0 t e ( t − s ) Δ R ( s ) d s . (130) e(t)=e^{t\Delta}e_0+\int_0^t e^{(t-s)\Delta}R(s)\,ds. \tag*{(130)} e ( t ) = e t Δ e 0 + ∫ 0 t e ( t − s ) Δ R ( s ) d s . (130) Taking norms, applying the contraction, and then Cauchy–Schwarz gives
∥ e ( t ) ∥ L 2 ≤ ∥ e 0 ∥ L 2 + ∫ 0 t ∥ R ( s ) ∥ L 2 d s ≤ ∥ e 0 ∥ L 2 + t ∥ R ∥ L 2 ( 0 , t ; L 2 ) . (131) \begin{aligned}
\|e(t)\|_{L^2}\le\|e_0\|_{L^2}+\int_0^t\|R(s)\|_{L^2}\,ds \\
&\le\|e_0\|_{L^2}+\sqrt{t}\|R\|_{L^2(0,t;L^2)}. \tag*{(131)}
\end{aligned} ∥ e ( t ) ∥ L 2 ≤ ∥ e 0 ∥ L 2 + ∫ 0 t ∥ R ( s ) ∥ L 2 d s ≤ ∥ e 0 ∥ L 2 + t ∥ R ∥ L 2 ( 0 , t ; L 2 ) . (131) The estimate follows upon taking the supremum.
Proposition E.3 (Wave equation with both initial mismatches). Let u t t − Δ u = 0 u_{tt}-\Delta u=0 u tt − Δ u = 0 v t t − Δ v = R v_{tt}-\Delta v=R v tt − Δ v = R e 0 = v ( 0 ) − u ( 0 ) e_0=v(0)-u(0) e 0 = v ( 0 ) − u ( 0 ) e 1 = v t ( 0 ) − u t ( 0 ) e_1=v_t(0)-u_t(0) e 1 = v t ( 0 ) − u t ( 0 ) R ∈ L 2 ( 0 , T ; L 2 ) R\in L^2(0,T;L^2) R ∈ L 2 ( 0 , T ; L 2 )
∥ v ( t ) − u ( t ) ∥ L 2 ≤ ∥ e 0 ∥ L 2 + t ∥ e 1 ∥ L 2 + t 3 / 2 3 ∥ R ∥ L 2 ( 0 , t ; L 2 ) . (132) \|v(t)-u(t)\|_{L^2}\le\|e_0\|_{L^2}+t\|e_1\|_{L^2}+\frac{t^{3/2}}{\sqrt{3}}\|R\|_{L^2(0,t;L^2)}. \tag*{(132)} ∥ v ( t ) − u ( t ) ∥ L 2 ≤ ∥ e 0 ∥ L 2 + t ∥ e 1 ∥ L 2 + 3 t 3/2 ∥ R ∥ L 2 ( 0 , t ; L 2 ) . (132) The energy also satisfies E ( t ) 1 / 2 ≤ E ( 0 ) 1 / 2 + 2 − 1 / 2 ∫ 0 t ∥ R ( s ) ∥ L 2 d s E(t)^{1/2}\le E(0)^{1/2}+2^{-1/2}\int_0^t\|R(s)\|_{L^2}\,ds E ( t ) 1/2 ≤ E ( 0 ) 1/2 + 2 − 1/2 ∫ 0 t ∥ R ( s ) ∥ L 2 d s E = 1 2 ( ∥ e t ∥ L 2 2 + ∥ ∇ e ∥ L 2 2 ) E=\frac{1}{2}(\|e_t\|_{L^2}^2+\|\nabla e\|_{L^2}^2) E = 2 1 ( ∥ e t ∥ L 2 2 + ∥∇ e ∥ L 2 2 )
Proof. For L = − Δ ≥ 0 L=-\Delta\ge0 L = − Δ ≥ 0
e ( t ) = cos ( t L ) e 0 + sin ( t L ) L e 1 + ∫ 0 t sin ( ( t − s ) L ) L R ( s ) d s , (133) e(t)=\cos(t\sqrt{L})e_0+\frac{\sin(t\sqrt{L})}{\sqrt{L}}e_1+\int_0^t\frac{\sin((t-s)\sqrt{L})}{\sqrt{L}}R(s)\,ds, \tag*{(133)} e ( t ) = cos ( t L ) e 0 + L sin ( t L ) e 1 + ∫ 0 t L sin (( t − s ) L ) R ( s ) d s , (133) where the quotient equals t t t ∣ cos ( t λ ) ∣ ≤ 1 |\cos(t\sqrt{\lambda})| \le1 ∣ cos ( t λ ) ∣ ≤ 1 ∣ sin ( t λ ) ∣ / λ ≤ t |\sin(t\sqrt{\lambda})|/\sqrt{\lambda} \le t ∣ sin ( t λ ) ∣/ λ ≤ t λ ≥ 0 \lambda\ge0 λ ≥ 0
∥ e ( t ) ∥ L 2 ≤ ∥ e 0 ∥ L 2 + t ∥ e 1 ∥ L 2 + ∫ 0 t ( t − s ) ∥ R ( s ) ∥ L 2 d s . (134) \|e(t)\|_{L^2} \le\|e_0\|_{L^2} + t\|e_1\|_{L^2} + \int_0^t (t-s)\|R(s)\|_{L^2}\,ds. \tag*{(134)} ∥ e ( t ) ∥ L 2 ≤ ∥ e 0 ∥ L 2 + t ∥ e 1 ∥ L 2 + ∫ 0 t ( t − s ) ∥ R ( s ) ∥ L 2 d s . (134) Cauchy–Schwarz and ∫ 0 t ( t − s ) 2 d s = t 3 / 3 \int_0^t (t-s)^2ds = t^3/3 ∫ 0 t ( t − s ) 2 d s = t 3 /3 (132) . For energy, test e t t − Δ e = R e_{tt}-\Delta e=R e tt − Δ e = R e t e_t e t
E ′ ( t ) = ⟨ R , e t ⟩ ≤ 2 E ( t ) ∥ R ( t ) ∥ L 2 . (135) E'(t) = \langle R,e_t\rangle\le\sqrt{2E(t)}\|R(t)\|_{L^2}. \tag*{(135)} E ′ ( t ) = ⟨ R , e t ⟩ ≤ 2 E ( t ) ∥ R ( t ) ∥ L 2 . (135) Where E ( t ) > 0 E(t)>0 E ( t ) > 0 2 E 2\sqrt{E} 2 E
Proposition E.4 (Conditional nonlinear stability). (i) For F ( z ) = ρ z ( 1 − z ) F(z)=\rho z(1-z) F ( z ) = ρ z ( 1 − z ) ρ > 0 \rho>0 ρ > 0 B B B [ 0 , T ] [0,T] [ 0 , T ] R = v t − F ( v ) ∈ L 1 ( 0 , T ; L 2 ) R=v_t-F(v)\in L^1(0,T;L^2) R = v t − F ( v ) ∈ L 1 ( 0 , T ; L 2 )
∥ v ( t ) − u ( t ) ∥ L 2 ≤ e ρ ( 1 + 2 B ) t ( ∥ v ( 0 ) − u ( 0 ) ∥ L 2 + ∫ 0 t ∥ R ( s ) ∥ L 2 d s ) . (136) \|v(t)-u(t)\|_{L^2} \le e^{\rho(1+2B)t}\left(\|v(0)-u(0)\|_{L^2}+\int_0^t\|R(s)\|_{L^2}\,ds\right). \tag*{(136)} ∥ v ( t ) − u ( t ) ∥ L 2 ≤ e ρ ( 1 + 2 B ) t ( ∥ v ( 0 ) − u ( 0 ) ∥ L 2 + ∫ 0 t ∥ R ( s ) ∥ L 2 d s ) . (136) (ii) On [ − 1 , 1 ] d [-1,1]^d [ − 1 , 1 ] d A ( z ) = − ε 2 Δ z + z 3 − z A(z)=-\varepsilon^2\Delta z+z^3-z A ( z ) = − ε 2 Δ z + z 3 − z λ 1 \lambda_1 λ 1 − Δ -\Delta − Δ ε 2 > λ 1 − 1 \varepsilon^2>\lambda_1^{-1} ε 2 > λ 1 − 1 R = A ( v ) − A ( u ) R=A(v)-A(u) R = A ( v ) − A ( u )
∥ v − u ∥ H 0 1 , ∇ ≤ ( ε 2 − λ 1 − 1 ) − 1 ∥ R ∥ H − 1 , (137) \|v-u\|_{H^1_0,\nabla}\le(\varepsilon^2-\lambda_1^{-1})^{-1}\|R\|_{H^{-1}}, \tag*{(137)} ∥ v − u ∥ H 0 1 , ∇ ≤ ( ε 2 − λ 1 − 1 ) − 1 ∥ R ∥ H − 1 , (137) where ∥ z ∥ H 0 1 , ∇ = ∥ ∇ z ∥ L 2 \|z\|_{H^1_0,\nabla}=\|\nabla z\|_{L^2} ∥ z ∥ H 0 1 , ∇ = ∥∇ z ∥ L 2
Proof. (i) On [ − B , B ] [-B,B] [ − B , B ] ∣ F ′ ( z ) ∣ = ρ ∣ 1 − 2 z ∣ ≤ ρ ( 1 + 2 B ) = L |F'(z)|=\rho|1-2z|\le\rho(1+2B)=L ∣ F ′ ( z ) ∣ = ρ ∣1 − 2 z ∣ ≤ ρ ( 1 + 2 B ) = L e = v − u e=v-u e = v − u e t = F ( v ) − F ( u ) + R e_t=F(v)-F(u)+R e t = F ( v ) − F ( u ) + R ∥ e ∥ L 2 \|e\|_{L^2} ∥ e ∥ L 2
D + ∥ e ∥ L 2 ≤ L ∥ e ∥ L 2 + ∥ R ∥ L 2 . (138) D^+\|e\|_{L^2}\le L\|e\|_{L^2}+\|R\|_{L^2}. \tag*{(138)} D + ∥ e ∥ L 2 ≤ L ∥ e ∥ L 2 + ∥ R ∥ L 2 . (138) Multiplying by e − L t e^{-Lt} e − L t e − L s ≤ 1 e^{-Ls}\le1 e − L s ≤ 1 (136) . (ii) Put e = v − u e=v-u e = v − u e e e
⟨ R , e ⟩ = ε 2 ∥ ∇ e ∥ L 2 2 + ∫ ( v 3 − u 3 ) ( v − u ) d x − ∥ e ∥ L 2 2 ≥ ( ε 2 − λ 1 − 1 ) ∥ ∇ e ∥ L 2 2 . \langle R,e\rangle=\varepsilon^2\|\nabla e\|_{L^2}^2+\int(v^3-u^3)(v-u)\,dx-\|e\|_{L^2}^2
\ge(\varepsilon^2-\lambda_1^{-1})\|\nabla e\|_{L^2}^2. ⟨ R , e ⟩ = ε 2 ∥∇ e ∥ L 2 2 + ∫ ( v 3 − u 3 ) ( v − u ) d x − ∥ e ∥ L 2 2 ≥ ( ε 2 − λ 1 − 1 ) ∥∇ e ∥ L 2 2 . The cubic integrand is nonnegative because
( v 3 − u 3 ) ( v − u ) = ( v − u ) 2 ( v 2 + u v + u 2 ) ≥ 0. (139) (v^3-u^3)(v-u)=(v-u)^2(v^2+uv+u^2)\ge0. \tag*{(139)} ( v 3 − u 3 ) ( v − u ) = ( v − u ) 2 ( v 2 + uv + u 2 ) ≥ 0. (139) Poincare’s inequality supplied the λ 1 − 1 \lambda_1^{-1} λ 1 − 1 ∥ R ∥ H − 1 ∥ ∇ e ∥ L 2 \|R\|_{H^{-1}}\|\nabla e\|_{L^2} ∥ R ∥ H − 1 ∥∇ e ∥ L 2 (137) .
Correction and Combined-Bound Proofs
Proof of Lemma 4.9. Concatenate the r 0 r_0 r 0 U ^ [ 0 ] \widehat{U}^{[0]} U [ 0 ] r m r_m r m Δ U ^ ( m ) \Delta\widehat{U}^{(m)} Δ U ( m ) r 0 + ∑ m r m r_0+\sum_m r_m r 0 + ∑ m r m
Proof of Theorem 4.10 . By the definition of the cumulative update,
E ^ [ m ] = U ^ N ∗ − ( U ^ [ m − 1 ] + Δ U ^ ( m ) ) = E ^ [ m − 1 ] − Δ U ^ ( m ) . (140) \begin{aligned}
\widehat{\mathcal{E}}^{[m]}=\widehat{U}_N^*-\left(\widehat{U}^{[m-1]}+\Delta\widehat{U}^{(m)}\right) \\
&=\widehat{\mathcal{E}}^{[m-1]}-\Delta\widehat{U}^{(m)}. \tag*{(140)}
\end{aligned} E [ m ] = U N ∗ − ( U [ m − 1 ] + Δ U ( m ) ) = E [ m − 1 ] − Δ U ( m ) . (140) An attained exact minimizer therefore has ∥ E ^ [ m ] ∥ F = δ m \|\widehat{\mathcal{E}}^{[m]}\|_F=\delta_m ∥ E [ m ] ∥ F = δ m δ m ≤ ∥ E ^ [ m − 1 ] ∥ F \delta_m\le\|\widehat{\mathcal{E}}^{[m-1]}\|_F δ m ≤ ∥ E [ m − 1 ] ∥ F
∥ E ^ [ m ] ∥ F ≤ δ m + ε m ≤ ∥ E ^ [ m − 1 ] ∥ F + ε m . (141) \|\widehat{\mathcal{E}}^{[m]}\|_F\le\delta_m+\varepsilon_m\le\|\widehat{\mathcal{E}}^{[m-1]}\|_F+\varepsilon_m. \tag*{(141)} ∥ E [ m ] ∥ F ≤ δ m + ε m ≤ ∥ E [ m − 1 ] ∥ F + ε m . (141) Under the extra compressibility inequality, a m : = ∥ E ^ [ m ] ∥ F a_m:=\|\widehat{\mathcal{E}}^{[m]}\|_F a m := ∥ E [ m ] ∥ F a m ≤ q m a m − 1 + ε m a_m\le q_ma_{m-1}+\varepsilon_m a m ≤ q m a m − 1 + ε m
a M ≤ ( ∏ j = 1 M q j ) a 0 + ∑ m = 1 M ε m ∏ j = m + 1 M q j . (142) a_M\le\left(\prod_{j=1}^{M}q_j\right)a_0+\sum_{m=1}^{M}\varepsilon_m\prod_{j=m+1}^{M}q_j. \tag*{(142)} a M ≤ ( j = 1 ∏ M q j ) a 0 + m = 1 ∑ M ε m j = m + 1 ∏ M q j . (142) The formula is true for M = 1 M=1 M = 1 M − 1 M-1 M − 1 q M q_M q M ε M \varepsilon_M ε M M M M q m ≤ q < 1 q_m\le q<1 q m ≤ q < 1 ε m ≤ ε \varepsilon_m\le\varepsilon ε m ≤ ε q q q
Proposition E.5 (Residual Reduction and Solution Error). Fix a complete coefficient residual objective
J ( v ) = ∥ W P N ( A v − f ) ∥ 2 , (143) J(v)=\lVert W P_N(\mathcal{A}v-f)\rVert_2, \tag*{(143)} J ( v ) = ∥ W P N ( A v − f ) ∥ 2 , (143) where W W W C m \mathcal{C}_m C m m m m c m ∈ C m c_m \in\mathcal{C}_m c m ∈ C m
J ( v m − 1 + c m ) ≤ inf c ∈ C m J ( v m − 1 + c ) + ζ m , (144) J(v_{m-1}+c_m)\le\inf_{c\in\mathcal{C}_m}J(v_{m-1}+c)+\zeta_m, \tag*{(144)} J ( v m − 1 + c m ) ≤ c ∈ C m inf J ( v m − 1 + c ) + ζ m , (144) then J ( v m ) ≤ J ( v m − 1 ) + ζ m J(v_m)\le J(v_{m-1})+\zeta_m J ( v m ) ≤ J ( v m − 1 ) + ζ m
a ∥ e ∥ H ≤ J ( u + e ) ≤ b ∥ e ∥ H ( a > 0 ) , (145) a\lVert e\rVert_H\le J(u+e)\le b\lVert e\rVert_H \qquad(a>0), \tag*{(145)} a ∥ e ∥ H ≤ J ( u + e ) ≤ b ∥ e ∥ H ( a > 0 ) , (145) then
∥ e m ∥ H ≤ b a ∥ e m − 1 ∥ H + ζ m a . (146) \lVert e_m\rVert_H\le\frac{b}{a}\lVert e_{m-1}\rVert_H+\frac{\zeta_m}{a}. \tag*{(146)} ∥ e m ∥ H ≤ a b ∥ e m − 1 ∥ H + a ζ m . (146) Proof. Zero is an admissible correction, so the infimum is at most J ( v m − 1 ) J(v_{m-1}) J ( v m − 1 )
a ∥ e m ∥ H ≤ J ( v m ) ≤ J ( v m − 1 ) + ζ m ≤ b ∥ e m − 1 ∥ H + ζ m . (147) \begin{aligned}
a\lVert e_m\rVert_H &\le J(v_m)\\
&\le J(v_{m-1})+\zeta_m\\
&\le b\lVert e_{m-1}\rVert_H+\zeta_m. \tag*{(147)}
\end{aligned} a ∥ e m ∥ H ≤ J ( v m ) ≤ J ( v m − 1 ) + ζ m ≤ b ∥ e m − 1 ∥ H + ζ m . (147) The factor b / a b/a b / a a = b a=b a = b w k = ∣ k ∣ − 2 w_k=\lvert k\rvert^{-2} w k = ∣ k ∣ − 2 a = b = 1 a=b=1 a = b = 1
Proof of Theorem 4.11. The triangle inequality and Equation (107) give
∥ u − u T S [ M ] ∥ H ≤ ∥ u − P N u ∥ H + ∥ U ^ N − U ^ [ M ] ∥ F = E s p e c ( N ) + ∥ E ^ [ M ] ∥ F . (148) \lVert u-u_{\mathrm{TS}}^{[M]}\rVert_H\le\lVert u-P_Nu\rVert_H+\lVert\widehat{U}_N-\widehat{U}^{[M]}\rVert_F
=E_{\mathrm{spec}}(N)+\lVert\widehat{E}^{[M]}\rVert_F. \tag*{(148)} ∥ u − u TS [ M ] ∥ H ≤ ∥ u − P N u ∥ H + ∥ U N − U [ M ] ∥ F = E spec ( N ) + ∥ E [ M ] ∥ F . (148) Insert Equation (22) with U ^ N ∗ = U ^ N \widehat{U}_N^*=\widehat{U}_N U N ∗ = U N (23) . The same two triangle inequalities as in Theorem 4.5 bound ∥ E ^ [ 0 ] ∥ F \lVert\widehat{E}^{[0]}\rVert_F ∥ E [ 0 ] ∥ F E r a n k ( N , r 0 ) + η + E o p t ( θ ; V ^ η ) E_{\mathrm{rank}}(N,r_0)+\eta+E_{\mathrm{opt}}(\theta;\widehat{\mathcal{V}}_\eta) E rank ( N , r 0 ) + η + E opt ( θ ; V η ) U ^ [ 0 ] = U ^ θ \widehat{U}^{[0]}=\widehat{U}_\theta U [ 0 ] = U θ
De-aliasing and Dimension Dependence
Proposition E.6 (Exact retained nonlinear coefficients). (i) Let every Fourier factor be supported in [ − K , K ] [-K,K] [ − K , K ] N = 2 K + 1 N=2K+1 N = 2 K + 1 s s s M M M [ − K , K ] [-K,K] [ − K , K ] M > ( s + 1 ) K M>(s+1)K M > ( s + 1 ) K M > 3 K M>3K M > 3 K M = ⌈ 3 N / 2 ⌉ M=\lceil3N/2\rceil M = ⌈ 3 N /2 ⌉ (ii) If each Jacobi factor has polynomial degree at most p p p s s s p p p Q Q Q
2 Q − 1 ≥ ( s + 1 ) p , Q ≥ ⌈ ( s + 1 ) p + 1 2 ⌉ . (149) 2Q-1\ge(s+1)p,\qquad Q\ge\left\lceil\frac{(s+1)p+1}{2}\right\rceil. \tag*{(149)} 2 Q − 1 ≥ ( s + 1 ) p , Q ≥ ⌈ 2 ( s + 1 ) p + 1 ⌉ . (149) For a cubic product the bound is Q ≥ 2 p + 1 Q\ge2p+1 Q ≥ 2 p + 1
Proof. (i) The exact s s s [ − s K , s K ] [-sK,sK] [ − sK , sK ] M M M n ∈ [ − K , K ] n\in[-K,K] n ∈ [ − K , K ] n + ℓ M n+\ell M n + ℓ M ℓ \ell ℓ M > ( s + 1 ) K M>(s+1)K M > ( s + 1 ) K
n + M ≥ − K + ( s + 1 ) K = s K , n+M\ge-K+(s+1)K=sK, n + M ≥ − K + ( s + 1 ) K = sK , n − M ≤ K − ( s + 1 ) K = − s K . (150) n-M\le K-(s+1)K=-sK. \tag*{(150)} n − M ≤ K − ( s + 1 ) K = − sK . (150) Every nonzero alias is outside the product support, so only the true coefficient at n n n s = 2 s=2 s = 2 N = 2 K + 1 N=2K+1 N = 2 K + 1 ⌈ 3 N / 2 ⌉ > 3 K \lceil3N/2\rceil>3K ⌈ 3 N /2 ⌉ > 3 K (ii) A degree-s s s p p p s p sp s p p p p ( s + 1 ) p (s+1)p ( s + 1 ) p 2 Q − 1 2Q-1 2 Q − 1 p p p p + 1 , … , 2 p p+1,\ldots,2p p + 1 , … , 2 p
Computational Cost
Figure 7 compares the computational cost per training update on the 30-dimensional wave equation. TensorSpectra requires only 0.529 GFLOPs for the base stage and 0.468 GFLOPs for the correction stage, which is substantially lower than PINN (47.34 GFLOPs), SPINN (33.16 GFLOPs), and PirateNet (242.22 GFLOPs), and is also slightly lower than SINN (0.716 GFLOPs). TNN has the lowest per-update cost at 0.048 GFLOPs, but its accuracy on this benchmark is notably lower than that of TensorSpectra. These results indicate that the low-rank spectral formulation provides a favorable balance between computational cost and solution accuracy in high-dimensional PDE solving.
Figure 7. Computational cost comparison on the 30-dimensional Wave benchmark. We report the GFLOPs required per training update for each method on a logarithmic scale. TensorSpectra is evaluated separately for the base and residual-correction stages. Lower values indicate lower computational cost.
Experimental Configurations
This section summarizes the main experimental configurations used throughout our numerical studies. Table 3 reports the spectral basis, temporal polynomial degree, tensor rank, training epochs, and learning rates adopted for TensorSpectra across different PDE benchmarks.
Table 3. Main hyperparameter settings of TensorSpectra for different PDE benchmarks.
