arXiv CSopen access

Stochastic-Dimension Frozen Sampled Neural Network for High-Dimensional Gross-Pitaevskii Equations on Unbounded Domains

Unknown · 2026 · arxiv_cs

arXiv CS · Papers · License: Open Access · 2026

machine learning, deep learning, neural networks

STOCHASTIC-DIMENSION FROZEN SAMPLED NEURAL NETWORK FOR HIGH-DIMENSIONAL GROSS-PITAEVSKII EQUATIONS ON UNBOUNDED DOMAINS ∗

arXiv:2604.09361v1 [cs.LG] 10 Apr 2026

ZHANGYONG LIANG† Abstract. In this paper, we propose a stochastic-dimension frozen sampled neural network (SD-FSNN) for solving a class of high-dimensional Gross–Pitaevskii equations (GPEs) on unbounded domains. SD-FSNN is unbiased across all dimensions, and its computational cost is independent of the dimension, avoiding the exponential growth in computational and memory costs associated with Hermite-basis discretizations. Additionally, we randomly sample the hidden weights and biases of the neural network, significantly outperforming iterative, gradient-based optimization methods in terms of training time and accuracy. Furthermore, we employ a space-time separation strategy, using adaptive ordinary differential equation (ODE) solvers to update the evolution coefficients and incorporate temporal causality. To preserve the structure of the GPEs, we integrate a Gaussianweighted ansatz into the neural network to enforce exponential decay at infinity, embed a normalization projection layer for mass normalization, and add an energy conservation constraint to mitigate long-time numerical dissipation. Comparative experiments with existing methods demonstrate the superior performance of SD-FSNN across a range of spatial dimensions and interaction parameters. Compared to existing random-feature methods, SD-FSNN reduces the complexity from linear to dimension-independent. Additionally, SD-FSNN achieves better accuracy and faster training compared to general high-dimensional solvers, while focusing specifically on high-dimensional GPEs on unbounded domains. Key words. High-dimensional Gross–Pitaevskii equations, unbounded domains, gradient-based optimization, mass normalization, energy conservation. AMS subject classifications. 65N22, 65N55, 68T07

1. Introduction. The Gross–Pitaevskii equation (GPE), a nonlinear Schrödinger equation describing the macroscopic wave function of a Bose–Einstein condensate in the low-temperature regime T ≪ Tc [31], successfully models ultracold-atom experiments with high fidelity [1, 21] and serves as a benchmark for emerging quantum simulation algorithms [11, 20], motivating the development of efficient and reliable numerical methods for its solution. The GPE is commonly defined on the whole space Rd , with solutions decaying at infinity. Early numerical methods tackled this by truncating the infinite domain and using Fourier or sine spectral discretizations [2, 12], which require careful selection of domain size and boundary conditions to avoid spurious reflections. Pseudo-spectral methods have been developed that naturally account for the decay at infinity by employing Hermite and Laguerre functions as basis functions [51, 42], enabling spectral accuracy without artificial boundary conditions. A fourth-order time-splitting Laguerre-Hermite method was introduced to reduce a three-dimensional cylindrically symmetric GPE to an effective two-dimensional problem [4], which was further extended by generalized Laguerre–Fourier–Hermite expansions to accommodate rotating and multi-component condensates while ensuring spectral accuracy and mass conservation [3]. More recently, adaptive Hermite spectral methods have further enhanced computational efficiency for unbounded domains by adjusting the scaling factor and expansion order according to the solution’s evolution [16, 18, 47]. More recently, and expansion order according to the solution’s evolution, enabling accurate and efficient ∗ Submitted to the editors DATE. † National Center for Applied Mathematics, Tianjin University, Tianjin, 300072, China.

(Corresponding author: Zhangyong Liang, E-mail: [email protected]). 1

Z. Y. LIANG

computations in high-dimensional and non-homogeneous settings [16, 18]. The adaptive schemes incorporate scaling adjustments and p-adaptivity, optimizing the spectral resolution and ensuring high accuracy across varying decay rates and solution behaviors [47, 16]. Despite these advances, Hermite-basis discretizations face a fundamental scalability barrier in high spatial dimensions: the number of tensor-product basis functions required for spectral accuracy grows exponentially with the ambient dimension d, leading to prohibitive computational and memory costs [16, 18, 47]. This exponential scaling renders traditional spectral methods intractable for high-dimensional GPEs arising in contemporary experimental contexts, such as dipolar condensates [30] and quantum simulation platforms with tens of spatial modes [52]. To mitigate the curse of dimensionality, domain truncation together with absorbing boundary conditions (ABCs) has been extensively adopted [54, 55, 17, 40]. However, the induced truncation error remains highly sensitive to the computational domain and boundary parameters. In parallel, deep-learning-based solvers for highdimensional PDEs have developed along three representative directions, depending on how differential operators are handled. The first direction reformulates PDEs through stochastic or probabilistic representations, as in DeepBSDE [23, 24] and its extensions [6, 13], deep splitting [5], forward-backward stochastic neural networks (FBSNNs) [46], multilevel Picard schemes [7, 8, 9], tensor neural networks [53, 34], and separable PINNs [15]. The second direction reduces operator cost via randomized estimators, including SDGD [26], HTE [25], RS-PINNs via Stein’s identity [27, 50], and STDEs [48]. The third direction avoids backpropagation and uses random features or analytic fitting, represented by RFM [44, 45], RaNNs [29], and ELMs [28]. It also includes related backpropagation-free variants, such as CAN-PINNs [14], HFDPINNs [36], fPINNs [41], the Deep Galerkin method [49], Taylor-mode automatic differentiation [10], and zeroth-order optimization [35, 37]. Despite these advances, unresolved limitations persist, including memory bottlenecks, bias or variance amplification in randomized estimators, and slow or unstable long-time training. Moreover, discretization-dependent methods often exhibit degraded accuracy in high dimensions relative to automatic differentiation. Despite the considerable progress outlined above, solving high-dimensional GPEs on unbounded domains remains challenging. Existing randomized methods for gradient estimation may introduce bias through nonlinear loss functions or suffer from variance amplification that grows with the ambient dimension, potentially causing spurious dissipation or energy pumping in long-time simulations. Gradient-based optimization methods for neural network training often suffer from slow convergence, training instability, and sensitivity to hyperparameters in high-dimensional spaces, resulting in suboptimal accuracy for long-time GPE evolution. Enforcing exponential decay at infinity through effective numerical representations without artificial boundary conditions remains challenging for neural network ansätze, since standard architectures lack natural mechanisms for spatial localization. Furthermore, global-in-time training or insufficient causal regularization leads to error accumulation, poor longterm accuracy, and failure to capture temporal dynamics faithfully in time-dependent GPEs. Finally, high-dimensional integration for enforcing conservation of mass and energy has not been adequately addressed in existing deep learning frameworks, leading to numerical dissipation and loss of structural properties over time. Motivated by the aforementioned challenges and existing methods, we develop the Stochastic-Dimension Frozen Sampled Neural Network, termed SD-FSNN, for solving the high-dimensional GPE on unbounded domains. Our main contributions can be summarized as follows:

STOCHASTIC-DIMENSION FROZEN SAMPLED NEURAL NETWORK

• Scalable weight and dimension sampling. We propose the SD-FSNN method by randomly sampling network weights and dimensions to avoid the exponential computational and memory costs associated with Hermite-basis discretizations. SD-FSNN provides an unbiased zeroth-order estimator over all dimensions with computational and memory costs independent of the ambient dimension, yielding a computing- and memory-efficient solver. • Space-time separation with causal marching. We employ a space-time separation strategy, in which time-dependent output-layer parameters are computed via least squares, coupled with adaptive ordinary differential equation (ODE) solvers. This enforces temporal causality by construction, avoids recomputing high-dimensional spatial features during time stepping, and enables accurate long-time prediction and temporal extrapolation. • Gradient-free randomized representation. We employ adaptively sampled weights to obtain a randomized representation and avoid gradient-based optimization, thereby alleviating slow training and low-accuracy issues that often arise in general-purpose solvers for high-dimensional PDEs. • Structure-preserving constraints. We incorporate a Gaussian-weighted ansatz into the neural network to enforce exponential decay at infinity. Then, we embed a normalization projection layer to enforce mass normalization and add an energy-conservation penalty to reduce long-time numerical dissipation. • Validation across dimensions and interaction. We evaluate the effectiveness and robustness of SD-FSNN across a range of spatial dimensions and interaction parameters. Numerical results demonstrate that SD-FSNN is efficient and accurate for GPEs, and outperforms existing general highdimensional PDE solvers on whole-space GPE benchmarks. The remainder of this paper is organized as follows. Section 2 presents the mathematical formulation of the high-dimensional GPE on unbounded domains, including the conserved quantities and numerical challenges arising from spectral discretization. Section 3 describes the SD-FSNN methodology in detail, including the neural network ansatz, the stochastic dimension sampling strategy, the space-time separation architecture, and the physical constraint enforcement mechanisms. Section 4 presents numerical experiments on benchmark problems to validate the efficiency and accuracy of SD-FSNN. Section 5 provides a detailed analysis of the advantages and limitations of our approach, including comparisons with existing methods. Section 6 summarizes the findings and discusses promising directions for future research. 2. Preliminaries. 2.1. Problem formulation. The dimensionless Gross-Pitaevskii equation (GPE) in d-dimensions [21, 32, 43] is given by (2.1)  1  i ∂t ψ(x, t) = − ∇2 ψ(x, t) + Vd (x)ψ(x, t) + βd |ψ(x, t)|2 ψ(x, t), x ∈ Rd , t > 0, 2   ψ(x, 0) = ψ0 (x), x ∈ Rd , where d

Vd (x) =

1X 2 2 γ x , 2 j=1 j j

βd = β

n Y j=d+1

γj . 2π

Z. Y. LIANG

with γj > 0 (j = 1, . . . , d) being constants. n is the original spatial dimension, and β is the interaction coefficient in the original n-dimensional equation. There are two extreme regimes of β: (1) β = o(1), (2.1) describes a weakly interacting condensation; (2) β ≫ 1, it corresponds to a strongly interacting condensation or to the semiclassical regime. The normalization condition for (2.1) is Z Z (2.2) N (ψ) = ∥ψ(·, t)∥2 = |ψ(x, t)|2 dx ≡ |ψ0 (x)|2 dx = 1, Rd

The GPE (2.1) conserves the Hamiltonian energy, defined as Z βd 1 |∇ψ(x, t)|2 + Vd (x)|ψ(x, t)|2 + |ψ(x, t)|4 dx, (2.3) E(ψ) = 2 Rd 2 During time evolution, the energy remains constant: E(ψ(·, t)) ≡ E(ψ0 ),

(2.4)

∀ t ≥ 0.

This conservation law follows from the fact that the GPE can be derived as the Hamiltonian flow associated with the energy functional (2.3) under the constraint (2.2). The energy consists of three terms: the kinetic energy ( 21 |∇ψ|2 ), the potential energy (Vd |ψ|2 ), and the interaction energy ( β2d |ψ|4 ). Since the above equation is posed on the whole line, it is natural to consider Hermite basis functions, which have been successfully applied to other equations [19, 22]. The numerical simulation is performed on a tensor product grid consisting of N + 1 points in each dimension. Consequently, the memory requirement for storing the solution is O(N d ). Regarding the computational complexity, the dominant cost per time step arises from the transform and its inverse (transferring data between physical space and phase space), as detailed in [4]. This process scales as O(dN d+1 ). 2.2. Unbiased randomized gradient estimation. Automatic differentiation (AD) [38, 39] evaluates derivatives of compositions of known primitives. Consider a network F = FL ◦ · · · ◦ F1 with hidden dimension h. Forward mode AD. Each primitive Fi is linearized as the Fréchet derivative ∂Fi , giving the Jacobian-vector-product (JVP): ∂F v = ∂F (x)(v) = [∂FL ◦ · · · ◦ ∂F1 ](x)(v). ∂x

(2.5)

Setting v to a standard basis of Rd computes one Jacobian column, requiring d JVPs total. Each call costs O(max(d, h)) memory and O dh + (L − 1)h2 operations for an MLP. Backward mode AD. Each primitive Fi is linearized as the adjoint ∂ ⊤ Fi , giving the vector-Jacobian-product (VJP): (2.6)

v⊤

∂F = ∂ ⊤ F (x)(v⊤ ) = [∂ ⊤ F1 (x) ◦ · · · ◦ ∂ ⊤ FL (yL−1 )](v⊤ ), ∂x

which requires a prior forward pass for the activation trace {yi }. Backward mode doubles sequential computation, stores the full trace (O(d + (L − 1)h) memory), but for scalar loss (d′ = 1) needs only one VJP, making it more efficient than forward mode in standard training.

STOCHASTIC-DIMENSION FROZEN SAMPLED NEURAL NETWORK

θ1

θ2

F1 x

∂F1

θ3 F3

F2 F2

θ1

y = F (x)

∂F1

∂F2

θ3 F3

F2 F2

y = F (x)

∂F3

∂F2

∂ ⊤ F1

θ2

∂F3

∂F ∂x v

∂ ⊤ F2 ∂ ⊤ F1

v⊤ ∂F ∂x

(a) Forward mode AD

v2⊤

∂ ⊤ F3 ∂ ⊤ F2

v1⊤

∂ ⊤ F3

v⊤

(b) Backward mode AD

Fig. 2.1: The computation graph of (a) forward mode AD and (b) backward mode AD of a function F (·) with 3 primitives Fi each parameterized by θi .

Repeating backward mode AD. With each repeated application of backward mode AD, the new evaluation trace will include the cotangents from the previous application of backward AD, so the length of sequential computation doubles. Furthermore, the size of the cotangent also increases by a factor of d. Therefore applying backward mode AD has additional memory cost of O(d + (L − 1)h) and additional computation cost of O 2dh + 2(L − 1)h2 , which is clear from Fig. 2.1. In general, with k repeated k−1 applications of backward mode AD will incur O 2 (d + (L − 1)h) memory cost and O 2k (dh + (L − 1)h2 ) computation cost. And O dk−1 calls are needed to evaluate the entire derivative tensor. So both memory and compute scale exponentially in derivative order k Repeating forward mode AD. Consider uθ : Rd → R. The input tangent dimension is d on the first application of forward mode AD, but on the second application, it will become d × d since we are now computing the forward mode AD for ∇uθ : Rd → Rd . So the size of the input tangent with k repeated application is O dk , so it grows exponentially. This is also inefficient. Building on this, we introduce the following unbiased randomized gradient estimation methods to avoid stacked back-propagation: • Stochastic dimension gradient descent (SDGD). SDGD [26] computes only a minibatch of derivatives in each iteration by replacing a differential operator D with a randomly sampled subset of additive terms, D :=

ND X j=1

Dj ≈

ND X Dj := D˜J , |J| j∈J

where Dj = ∂ 2 /∂x2j is the Laplacian and |J| is the batch size. By treating nonsampled dimensions as constants, D˜J reduces complexity from O(2k−1 (d + (L − 1)h)) to O(|J|2k (h + (L − 1)h2 )), though the exponential scaling with respect to the derivative order k persists. • Randomized smoothing PINNs (RS-PINNs). RS-PINNs [27] bypass automatic differentiation through Monte Carlo estimation via Stein’s identity,

Z. Y. LIANG

parameterizing the solution with a Gaussian-smoothed neural network: B

u(x) = Eδ∼N (0,σ2 I) [f (x + δ)] ≈

1 X f (x + δb ) := ũ, B b=1

where f is the base network, δb are i.i.d. Gaussian noise samples, and B is the batch size. Higher-order derivatives are similarly expressed as expectations, achieving O(BL) per-iteration complexity independent of d and k. However, nonlinearities in the PDE residual introduce bias; bias-correction techniques can recover unbiasedness but may increase variance in high dimensions. • Hutchinson trace estimation (HTE). HTE [25] approximates the Laplacian via random vector–vector products: B h i 1 X ⊤ 2 e vb ∇ u(x) vb := ∆, ∆u(x) = tr(∇2 u(x)) = Ev∼p v ⊤ ∇2 u(x) v ≈ B b=1

iid

where v ∼ p is drawn from an isotropic distribution (e.g., vi ∼ N (0, 1) or Rademacher ±1). The quadratic form v ⊤ ∇2 u(x)v is evaluated with two backward passes, achieving O(B·L) complexity, but does not naturally extend to higher-order derivatives. • Stochastic Taylor derivative estimator (STDE). STDE [48] leverages high-order Taylor-mode automatic differentiation on randomly constructed input k-jets: B 1 X Du(x) = EJk ∼p c · dk u(x; Jk )[k] ≈ cb · dk u(x; Jk,b )[k] := D̃, B b=1

where Jk = (x, v (1) , . . . , v (k) ) is a k-jet, dk u(x; Jk ) is the k-th order pushforward (Taylor-mode AD), [k] extracts the highest-order component, and c is a debiasing scalar. Sparse random jets (tangents supported on a small random subset of dimensions) reduce per-iteration complexity to O(Bk 2 L), making both memory and time independent of d and avoiding the exponential graph explosion O(2k−1 L) of repeated backward-mode AD. • Forward Laplacian. Forward Laplacian [33] removes redundancy in the AD pipeline for Laplacian calculation, achieving O(d) complexity. Since it is not randomized, the scaling deteriorates at high dimensions: when d > 1000, it becomes worse than even parallelized stacked backward-mode SDGD, whereas randomized estimators like STDE maintain O(|J|) complexity. 2.3. Randomized neural networks. Although unbiased randomized gradient estimation methods effectively address the computational challenges of automatic differentiation in high-order and high-dimensional settings, they still rely on gradientbased optimization. This requires thousands to tens of thousands of iterations, with training times typically ranging from minutes to hours. The highly non-convex loss landscape of PINNs often leads to slow convergence, training instability, and failure due to imbalanced loss terms or poor local minima. To eliminate these persistent challenges of gradient optimization, randomized neural networks (RaNNs) have been proposed as powerful gradient-free solvers for PDEs. Given a probability distribution ν on the space of square-integrable functions L2 (D; R) on the domain D ⊂ Rd , a general random feature model is defined as the

STOCHASTIC-DIMENSION FROZEN SAMPLED NEURAL NETWORK

weighted sum of several i.i.d. generated functions φ1 , . . . φN ∼ ν, that is UW : Rd → R : x 7→

(2.7)

N X

Wi φi (x),

i=1

where the weights W1 , . . . WN ∈ R can be chosen freely, and will be optimized such that UW ≈ u. Often, one chooses parametrized functions φi (·) := φ(·; θi ) as random features, such that only (finite-dimensional) parameters θ1 , . . . θN ∼ ν ∗ need to be drawn. RaNNs realize the random features as the activations of a single hidden-layer feedforward neural network with frozen parameters. In this setting, the definition of a RaNNs (2.7) can be specialized to A,B : Rd → R : x 7→ UW

(2.8)

N X

Wi σ(Ai · x + Bi ),

i=1

where we let σ : R → R is a nonlinear activation function, we let A1 , . . . AN ∈ Rd be random weight vectors and we let B1 , . . . BN ∈ R be random bias scalars; all of which are independent. The goal is to choose the output weights W = (W1 , . . . WN ) as a function of the random hidden weights A = (A1 , . . . AN ) and B = (B1 , . . . , BN ) such A,B that UW is a good approximation of u. We note that one can also add bias W0 to definitions (2.7) and (2.8). For ELMs, the weights of random single-hidden-layer forward neural networks (SLFNs) are determined analytically. Given access to a training set {(x1 , u(x1 )), . . . (xn , u(xn ))} the goal is to find the weight vector W ∈ RN such that ideally (2.9)

N X

Wi σ(Ai · xk + Bi ) = u(xk ),

∀k ∈ {1, . . . , n}.

i=1

If we define matrix H ∈ Rn×N with Hki = σ(Ai · xk + Bi ) and a vector T ∈ Rn with Tk = u(xk ) then the above n equations (2.9) can be written compactly as (2.10)

HW = T.

Rather than solving this with gradient-based learning algorithms, the ELMs method proposes to choose Ŵ = H † T , where H † is the Moore-Penrose generalized inverse of H. Hence, Ŵ is the smallest norm least-squares solution of the linear system (2.10). As such, the ELMs framework provides a simple and extremely fast alternative to training neural networks with gradient-based optimizers. However, both RaNNs and ELMs treat time as an additional spatial dimension, causing neural basis functions to span the entire spatiotemporal domain without accounting for temporal dependencies, which can lead to error accumulation and propagation failures in long-time simulations. 3. Stochastic-dimension frozen sampled neural network (SD-FSNN). In this section, we introduce a stochastic-dimension frozen sampled neural network (SD-FSNN) method to address the aforementioned challenges. 3.1. Neural network ansatz. We parameterize the approximation of a solution (2.1) with a neural network with one hidden layer, activation function σ = tanh, and M neurons so that (3.1)

û(x, t) = C(t)[Φ(x), 1]⊤ = c(t)σ(W x⊤ + b) + c0 (t),

Z. Y. LIANG

where c(t) ∈ C1×M and c0 (t) ∈ C are time-dependent parameters, W ∈ RM ×d and b ∈ RM ×1 are space-independent parameters, and C(t) := [c(t), c0 (t)] ∈ C1×(M +1) . The activation functions are stacked in Φ(x) = [ϕ1 (x), . . . , ϕM (x)], where ϕm (x) = σ(wm x⊤ + bm ). Note that our approach does not require the solution to be separable in space and time. We next discuss how to sample parameters W and b. We construct each weight and bias pair wk , bk by taking two points x(1) , x(2) ∈ Ω and construct the weight and bias as (3.2)

wk = s1

x(2) − x(1) , ∥x(2) − x(1) ∥2

bk = −⟨wk , x(1) ⟩ + s2 ,

where s1 , s2 ∈ R depend on the activation function. In the unsupervised setting, one can choose pairs of collocation points from a uniform distribution over all possible pairs of collocation points, which is the default setting in this paper, as we do not know the solution of the GPEs beforehand. In the supervised setting, collocation pairs (x(1) , x(2) ) are sampled with solution ∥f (x(2) ) − f (x(1) )∥/∥x(2) − x(1) ∥. This data-driven initialization, referred to as the SWIM (Sampling Weights by Interpolating Marginals) strategy, strategically places basis functions with steep gradients near regions with physical shocks or rapid variations. Neuron weights and biases are set so that the tanh output is −0.5 at x(1) and +0.5 at x(2) , ensuring centers of activations tanh lie inside the domain and are aligned with the direction x(1) → x(2) , unlike ELMs. The suitability of each of the proposed approaches depends on the ground truth solution’s gradient distribution. Data-driven sampling

Data-agnostic sampling

Fig. 3.1: Sampling in SD-FSNN: (Left): SWIM (data-driven, places bases with steep gradients near regions with shocks) vs. (Right): ELM (data-agnostic, no control over basis placement).

Using neural networks as an ansatz for the solution has proven a challenge in terms of training time and approximation accuracy. In this contribution, we discuss how sampling the hidden weights and biases of the ansatz network from data-agnostic and data-dependent probability distributions enables progress on both challenges. 3.2. Space-time separation. We now discuss the computation of time-dependent output layer parameters C(t). We insert the ansatz (3.1) into the GPEs (2.1), reformulating it as an ODE for C(t), preserving the inherent causal structure of timedependent PDEs, thereby enforcing temporal causality by design. We assemble Nc collocation points in X ∈ RNc ×d , sample weights and biases of M neurons, compute hidden layer output Φ(X), and obtain the complex-valued ODE + Ct (t) = R(X, C(t)) [Φ(X), 1Nc ]⊤ , where 1 (3.3) R(X, C(t)) = i C(t)∇2 [Φ(X), 1Nc ]⊤ − Vd (X)⊤ ⊙ C(t)[Φ(X), 1Nc ]⊤ 2 , −βd |C(t)[Φ(X), 1Nc ]⊤ |2 ⊙ C(t)[Φ(X), 1Nc ]⊤

STOCHASTIC-DIMENSION FROZEN SAMPLED NEURAL NETWORK

where 1Nc ∈ RNc ×1 is a column vector of ones, [Φ(X), 1Nc ]⊤ ∈ R(M +1)×Nc , C(t) ∈ C1×(M +1) , ⊙ denotes element-wise multiplication, and the pseudo-inverse is denoted by ·+ . Since standard ODE solvers are designed for real-valued systems, we decompose the complex-valued parameters as C(t) = Cre (t) + iCim (t), where Cre (t), Cim (t) ∈ R1×(M +1) , yielding the coupled real-valued ODE system, as follows  +  Re R(X, C(t)) [Φ(X), 1Nc ]⊤ d Cre (t) = (3.4) +  . dt Cim (t) Im R(X, C(t)) [Φ(X), 1N ]⊤ c

The initial condition is computed via a least squares solution + (3.5) C(0) = ψ0 (X)⊤ [Φ(X), 1Nc ]⊤ , which is similarly decomposed into real and imaginary parts for initialization. We solve the coupled system (3.4) using standard ODE solvers with step-size control (e.g., DOPRI5 or RK45), treating the 2(M + 1)-dimensional real vector [Cre , Cim ]⊤ as the state variable, and reconstruct the complex solution via C(t) = Cre (t) + iCim (t). For the GPEs, we primarily require the Laplacian operator ∇2 . We describe how to compute the relevant spatial derivatives of the approximate solution ψ̂(x, t) = C(t)[Φ(x), 1]⊤ : • First-order spatial derivative of the approximate solution is computed as: Φ (x) ∇x ψ̂(x, t) = C(t) x 01×d (3.6) W ⊙ σ̃x (x) = C(t) ∈ C1×d , 01×d where ⊙ denotes the Hadamard product, 01×d is a zero row vector, and σ̃x (x) := [σ ′ (z)|z=W x⊤ +b , . . . , σ ′ (z)|z=W x⊤ +b ]⊤ ∈ RM ×d , with σ ′ (z) = 1 − tanh2 (z) being the first derivative of the tanh activation function. • Second-order spatial derivative of the approximate solution is computed as: Φxx (x) 2 ∇x ψ̂(x, t) = C(t) 01×d (3.7) W ⊙ W ⊙ σ̃xx (x) = C(t) ∈ C1×d , 01×d where σ̃xx (x) is defined analogously to σ̃x (x) with σ ′′ (z) = −2 tanh(z)(1 − tanh2 (z)) being the second derivative of the tanh activation function. • The Laplacian of the approximate solution, which is essential for the GPEs, is computed by summing across spatial dimensions: d X ∂ 2 ψ̂ Φxx (x) 2 = C(t) 1d ∇ ψ̂(x, t) = 01×d ∂x2j j=1 (3.8) W ⊙ W ⊙ σ̃xx (x) = C(t) 1d ∈ C, 01×d where 1d ∈ Rd×1 is a column vector of ones.

Z. Y. LIANG

Following the spatial derivative formulations above, we now describe how to compute the time derivatives of the approximate GPEs solution ψ̂(x, t). Since the GPEs (2.1) is first-order in time, we only require the first-order time derivative: • First-order time derivative of the approximate solution is computed as: ∂t ψ̂(x, t) = Ct (t)[Φ(x), 1]⊤ ∈ C,

(3.9)

1×(M +1) is obtained from the ODE system (3.4). where Ct (t) = dC(t) dt ∈ C For the coupled real-valued ODE formulation, the time derivative is decomposed

as: ∂t ψ̂(x, t) = ∂t ψ̂re (x, t) + i ∂t ψ̂im (x, t) (3.10) =

dCim (t) dCre (t) [Φ(x), 1]⊤ + i [Φ(x), 1]⊤ , dt dt

(t) re (t) where dCdt and dCim are given by the right-hand side of the coupled system (3.4). dt For multiple collocation points X ∈ RNc ×d , the time derivative is evaluated as:

(3.11)

∂t ψ̂(X, t) = Ct (t)[Φ(X), 1Nc ]⊤ ∈ C1×Nc .

Approximating the solution with the neural network ansatz (3.1) where C(t) ∈ C1×(M +1) , and substituting (3.11) and (3.8) into the GPEs, we get 1 Ct (t)[Φ(X), 1Nc ]⊤ = i C(t)∇2 [Φ(X), 1Nc ]⊤ − Vd (X)⊤ ⊙ C(t)[Φ(X), 1Nc ]⊤ 2 −βd |C(t)[Φ(X), 1Nc ]⊤ |2 ⊙ C(t)[Φ(X), 1Nc ]⊤ , 1 Ct (t) = i C(t)∇2 [Φ(X), 1Nc ]⊤ − Vd (X)⊤ ⊙ C(t)[Φ(X), 1Nc ]⊤ 2 + −βd |C(t)[Φ(X), 1Nc ]⊤ |2 ⊙ C(t)[Φ(X), 1Nc ]⊤ [Φ(X), 1Nc ]⊤ . Since the GPEs involves complex-valued solutions, we decompose the complex coefficients into real and imaginary parts as C(t) = Cre (t)+iCim (t), where Cre (t), Cim (t) ∈ R1×(M +1) . Let ψ̂ = C(t)[Φ(X), 1Nc ]⊤ = ψ̂re + iψ̂im . The GPE separates into coupled real-valued equations: 1 ∂t ψ̂re = − ∇2 ψ̂im + Vd ψ̂im + βd (|ψ̂|2 ψ̂im ), 2 1 2 ∂t ψ̂im = ∇ ψ̂re − Vd ψ̂re − βd (|ψ̂|2 ψ̂re ), 2 2 2 + ψ̂im . where |ψ̂|2 = ψ̂re This yields the coupled system of real-valued ODEs:

1 Cre,t (t)[Φ(X), 1Nc ]⊤ = − Cim (t)∇2 [Φ(X), 1Nc ]⊤ + Vd (X)⊤ ⊙ Cim (t)[Φ(X), 1Nc ]⊤ 2 + βd |C(t)[Φ(X), 1Nc ]⊤ |2 ⊙ Cim (t)[Φ(X), 1Nc ]⊤ , 1 Cim,t (t)[Φ(X), 1Nc ]⊤ = Cre (t)∇2 [Φ(X), 1Nc ]⊤ − Vd (X)⊤ ⊙ Cre (t)[Φ(X), 1Nc ]⊤ 2 − βd |C(t)[Φ(X), 1Nc ]⊤ |2 ⊙ Cre (t)[Φ(X), 1Nc ]⊤ ,

STOCHASTIC-DIMENSION FROZEN SAMPLED NEURAL NETWORK

which can be reformulated as 1 Cre,t (t) = − Cim (t)∇2 [Φ(X), 1Nc ]⊤ + Vd (X)⊤ ⊙ Cim (t)[Φ(X), 1Nc ]⊤ 2 + +βd |C(t)[Φ(X), 1Nc ]⊤ |2 ⊙ Cim (t)[Φ(X), 1Nc ]⊤ [Φ(X), 1Nc ]⊤ , 1 Cre (t)∇2 [Φ(X), 1Nc ]⊤ − Vd (X)⊤ ⊙ Cre (t)[Φ(X), 1Nc ]⊤ Cim,t (t) = 2 + −βd |C(t)[Φ(X), 1Nc ]⊤ |2 ⊙ Cre (t)[Φ(X), 1Nc ]⊤ [Φ(X), 1Nc ]⊤ . Note that the cubic nonlinearity |ψ|2 ψ is transferred to the right-hand side of the ODE system. The initial condition is computed via +

⊤ +

Cre (0) = Re(ψ0 (X))⊤ [Φ(X), 1Nc ]⊤ Cim (0) = Im(ψ0 (X))

⊤

[Φ(X), 1Nc ]

Remark 3.1. Unlike wave-type equations that require second-order time derivatives, the GPE is first-order in time due to the presence of the imaginary unit i in front of ∂t ψ. This structure is preserved in our ODE reformulation (3.3), where Ct (t) is directly determined by the spatial operators without requiring Ctt (t). 3.3. SVD decomposition layer. As the last step in the SD-FSNN framework, we add a linear layer to reduce the stiffness of the associated ODE (3.3) and the size of the ODE system. To achieve this, we propose orthogonalizing the basis functions using an SVD layer. Let the augmented feature matrix be Ψ(X) = [Φ(X), 1Nc ]⊤ ∈ R(M +1)×Nc . We compute its truncated singular value decomposition to obtain matrices Vr , Σr , and Ur with rank r ≤ M + 1 such that Vr Σr Ur⊤ ≈ Ψ(X). We then define the reduced orthogonal basis Ψr (X) := Vr⊤ Ψ(X) ∈ Rr×Nc and use it instead of Ψ(X), with the output coefficients adapted to C(t) ∈ C1×r . This ensures that the rows of Ψr (X) are orthogonal functions evaluated on the data X, and the corresponding pseudo-inverse operation possesses a bounded condition number. 3.4. Sampled structure-preserving weights. Now, we introduce structurepreserving features into the sampling weights, including exponential decay at infinity, mass normalization, and energy conservation. These three structure-preserving properties are essential for the physical space, with a focus on the sampled basis functions after space-time separation. For exponential decay at infinity. To enforce the required whole-space boundary condition |ψ(x, t)| → 0 as ∥x∥ → ∞ without artificial truncation of Ω = Rd , we introduce a decaying envelope into the frozen spatial representation. Let ρ : Rd → (0, ∞) be a smooth function with exponential or Gaussian decay, such as ρ(x) = exp(−κg(x)), where g(x) → ∞ as ∥x∥ → ∞. The wave function is factored as ψ(x, t) = ρ(x)u(x, t). We apply the proposed space-time separation in Subsection 3.2 to the de-enveloped field u, with û(x, t) = C(t)[Φ(x), 1]⊤ , and reconstruct the solution as ψ̂(x, t) = ρ(x)û(x, t) = ρ(x)C(t)[Φ(x), 1]⊤ . Since tanh features are bounded, the decay of ψ̂ at infinity is inherited from ρ, ensuring the decay property by construction, rather than through boundary conditions. Substituting ψ = ρu into (2.1) yields an equivalent evolution for u as follows h i (3.12) ut = i 12 ∇2 u + (∇ρ/ρ) · ∇u + 21 (∇2 ρ/ρ)u − Vd u − βd ρ2 |u|2 u ,

Z. Y. LIANG

and collocation at X with (3.6)–(3.8) leads to a space-time separated ODE that retains the algebraic form Ct (t) = Rρ (X, C(t)) [Φ(X), 1Nc ]⊤

(3.13) where aj (X) = by

∂

xj ρ

(X) and b(X) =

∇2 ρ ρ

(X). The operator Rρ (X, C) is given

d h X Rρ (X, C) = i 21 C ∇2 [Φ(X), 1Nc ]⊤ + aj (X)⊤ ⊙ C ∂xj [Φ(X), 1Nc ]⊤

(3.14)

j=1

+ 12 b(X)⊤ ⊙

C[Φ(X), 1Nc ]⊤ − Vd (X)⊤ ⊙ C[Φ(X), 1Nc ]⊤ i − βd (ρ(X)2 )⊤ ⊙ |C[Φ(X), 1Nc ]⊤ |2 ⊙ C[Φ(X), 1Nc ]⊤ ,

Φxj (X)⊤ . The initial coefficients are obtained from the 01×Nc de-enveloped initial data as follows with ∂xj [Φ(X), 1Nc ]⊤ =

(3.15)

C(0) = ρ(X)−1 ⊙ ψ0 (X)

⊤

[Φ(X), 1Nc ]⊤

which avoids forming a pseudo-inverse of features multiplied by a vanishing envelope and thus improves numerical conditioning on unbounded domains. Finally, to adapt weight sampling to the whole-space setting, collocation points (and hence collocation pairs used in (3.2)) are drawn from a probability measure induced by the envelope, e.g., with density proportional to ρ(x)q (typically q = 2), concentrating frozen activations where the condensate has non-negligible mass while retaining support over Rd . For R mass normalization. To enforce the physically required mass normalization Rd |ψ(x, t)|2 dx = 1 for all t ≥ 0 while fully preserving the space-time separated structure described in Subsection 3.2, we incorporate a normalization mapping inspired by constrained neural network architectures directly into the coefficient-level representation. The standard ansatz ψ̂(x, t) = C(t)[Φ(x), 1]⊤ expresses the solution as a time-dependent linear combination of the frozen spatial basis Ψ(x) = [Φ(x), 1]⊤ ∈ R(M +1)×1 . The continuous L2 -norm (mass) is then Z (3.16) ∥ψ̂(·, t)∥2L2 = |C(t)Ψ(x)|2 dx. Rd

This integral is approximated discretely using the same collocation points X = Nc c {xn }n=1 employed in the ODE reformulation (3.3), with quadrature weights {wn }N n=1 chosen consistently with the sampling distribution, yielding (3.17)

M (t) =

Nc X

wn |C(t)Ψ(xn )|2 = C(t) G C(t)∗ ,

n=1

where the Gram matrix G ∈ C(M +1)×(M +1) is time-independent and precomputable, (3.18)

Gkj =

Nc X n=1

wn Ψk (xn ) Ψj (xn ).

STOCHASTIC-DIMENSION FROZEN SAMPLED NEURAL NETWORK

To achieve exact discrete mass conservation M (t) = 1 by construction, we modify the ansatz to (3.19)

C(t)Ψ(x) , ψ̂(x, t) = p C(t)GC(t)∗

where C(t) ∈ C1×(M +1) now denotes unnormalized coefficients. This ensures (3.20)

M (t) =

Nc X

C(t)Ψ(xn ) wn p C(t)GC(t)∗ n=1

for any nonzero C(t), while the continuous mass is preserved to quadrature accuracy owing to the concentration of weight-sampled collocation points and activation centers in regions of significant condensate density. The initial condition is enforced in the discrete sense by first solving the standard least-squares problem C(0) = ψ0 (X)⊤ Ψ(X)+T and then applying the normalization (3.19) at t = 0. To obtain the evolution equation for the unnormalized p coefficients, we substitute the normalized ansatz into the GPE. Let N (t) = 1/ C(t)GC(t)∗ and denote the unnormalized coefficients temporarily by C̃(t), so that ψ̂(x, t) = N (t) C̃(t) Ψ(x). Time differentiation yields (3.21)

∂t ψ̂(x, t) = Ṅ (t) C̃(t) Ψ(x) + N (t) C̃t (t) Ψ(x).

Because the normalization factor N (t) is a strictly real positive scalar, its time derivative Ṅ (t) is also purely real. When substituted into the GPE’s left-hand side (i∂t ψ̂), this term produces a purely imaginary contribution iṄ (t)C̃(t)Ψ(x) to the residual. In the physical context of the Schrödinger equation, purely imaginary components represent the non-unitary injection or dissipation of probability amplitude (mass), rather than a physical Hamiltonian phase rotation (which would arise strictly from real-valued energy terms). To prevent this non-physical mass fluctuation from continuously coupling into the projected ODE, we systematically decouple the normalization during the continuous time-step. Thus, we track the unnormalized coefficients C̃(t) using the standard space-time separated form (3.4) to produce candidate states C̃(t + ∆t), and subsequently enforce the invariant mass manifold via an explicit post-step orthogonal projection: (3.22)

C̃(t + ∆t) . C(t + ∆t) = q C̃(t + ∆t) G C̃(t + ∆t)∗

This projection enforces exact discrete mass conservation at each step, incurs negligible overhead of O(M 2 + M Nc ), preserves the causal ODE structure derived in Subsection 3.2, and remains fully compatible with the sampled frozen weights and the envelope-based whole-space treatment, as the Gram matrix G naturally accounts for the physical distribution of activation centers induced by the sampling strategy. For energy conservation. To enforce long-term conservation of the total energy Z βd 1 |∇ψ̂(x, t)|2 + Vd (x)|ψ̂(x, t)|2 + |ψ̂(x, t)|4 dx (3.23) E(t) = 2 Rd 2

Z. Y. LIANG

within the space-time separated framework of Subsection 3.2, we discretize E(t) using Nc c the same collocation points X = {xn }N n=1 and quadrature weights {wn }n=1 as in (3.3),

(3.24)

E(t) ≈

Nc X n=1

1 βd |∇x ψ̂(xn , t)|2 + Vd (xn )|ψ̂(xn , t)|2 + |ψ̂(xn , t)|4 , 2 2

W ⊙ σ̃x (xn ) follow directly 01×d from (3.6). With the ansatz ψ̂(x, t) = C(t)Ψ(x), the discrete energy becomes a scalar function of the coefficient vector alone, where ψ̂(xn , t) = C(t)Ψ(xn ) and ∇x ψ̂(xn , t) = C(t)

(3.25)

E(t) = Ekin (C(t)) + Epot (C(t)) + Eint (C(t)),

with explicit quadratic and quartic terms in C(t)C(t)∗ computable from pre-stored matrices of ∇2 Ψ(X), Ψ(X), and |Ψ(X)|2 ⊙ Ψ(X). The projected ODE system (3.3)–(3.4) is equivalent to a finite-dimensional Hamiltonian system (3.26)

δH dC(t) = , dt δC(t)∗

where the discrete Hamiltonian H(C, C ∗ ) = E(t) coincides with (3.24) up to quadrature error. The real-imaginary split system (3.4) admits the canonical Hamiltonian structure d Cre 0 I Cre (3.27) = J∇z H(z), J = , z= . −I 0 Cim dt Cim To preserve H(z(t)) up to bounded oscillations O(∆tp ) over arbitrarily long times, we employ implicit symplectic integrators or Hamiltonian operator splitting methods, since the nonlinear interaction term renders the Hamiltonian non-separable. Adaptive step-size control is retained via embedded error estimators. From the perspective of geometric numerical integration, applying the explicit mass-normalization projection (3.22) after a symplectic integration step inevitably breaks the strict mathematical symplectic structure. However, we adopt this projected symplectic integrator approach because explicitly confining the system to the physically correct mass manifold effectively suppresses long-term phase drifts and probability amplitude dissipation. In practice, combined with the discrete mass-normalization projection (3.22) applied after each accepted step, the resulting projected integrator demonstrates excellent long-term pseudo-symplectic stability, yielding (3.28)

|H(zn ) − H(z0 )| ≤ εmach + O(∆tp ),

∀n.

with observed relative energy errors typically bounded below 10−8 over t ∈ [0, 104 ] in benchmark simulations, while preserving the causal ODE structure, frozen sampled basis, and envelope-augmented whole-space representation derived in Subsection 3.2. 3.5. Stochastic dimension for sampling weights. We introduce a stochasticdimension frozen sampled neural network (SD-FSNN) in the sampled-weight formulation of the d-dimensional GPE. Recall the space–time separated ansatz (3.29)

ψ̂(x, t) = C(t)[Φ(x), 1]⊤ ,

STOCHASTIC-DIMENSION FROZEN SAMPLED NEURAL NETWORK

where the hidden-layer basis functions Φ(x) are fixed after sampling, and only the output coefficients C(t) evolve in time. After substituting this ansatz into (2.1), the time evolution is governed by the coefficient ODE (3.3), whose right-hand side repeatedly evaluates spatial operators on Φ at collocation points. In high dimensions, the dominant cost is typically the Laplacian operator term. The Laplacian operator admits a dimension-wise additive decomposition (3.30)

∇2 ψ̂(x, t) =

d X ∂ 2 [Φ(x), 1]⊤ (x, t) = C(t) , 2 ∂xj ∂x2j j=1 j=1

d X ∂ 2 ψ̂

so evaluating ∇2 ψ̂ amounts to aggregating contributions from all d spatial dimensions. For neural basis functions of the form ϕm (x) = σ(wm x⊤ + bm ), each second partial derivative satisfies (3.31)

∂ 2 ϕm (x) 2 = wm,j σ ′′ (wm x⊤ + bm ), ∂x2j

which highlights that the j-th dimension contribution is determined by the j-th weight component wm,j . Stochastic-dimension Laplacian estimation. To amortize the dimensionwise sum in (3.30), we sample a minibatch of dimensions. Let J ⊂ {1, . . . , d} be a random index set with |J| = m ≪ d, sampled uniformly (with or without replacement). We define the stochastic-dimension Laplacian estimator (3.32)

X ∂ 2 ψ̂ f2 J ψ̂(x, t) := d ∇ (x, t), |J| ∂x2j j∈J

where non-sampled dimensions are treated as constants during evaluation. By linearity and symmetry of uniform sampling, h i f2 J ψ̂(x, t) = ∇2 ψ̂(x, t). (3.33) EJ ∇ More generally, under non-uniform sampling with replacement where each draw selects j with probability pj > 0, a Horvitz–Thompson form remains unbiased: |J|

(3.34)

X 1 ∂ 2 ψ̂ f2 J ψ̂(x, t) := 1 (x, t), ∇ |J| pjℓ ∂x2jℓ

jℓ ∼ p(·).

ℓ=1

Since each sampled dimension j ∈ J only activates the corresponding weight components {wm,j }M m=1 , (3.32) effectively induces a stochastic sparsification of the hiddenlayer weights during residual evaluation, reducing the Laplacian evaluation cost from O(M d) to O(M |J|). SD-FSNN residual and induced stochastic ODE. Define the (complexvalued) GPE residual quantity at (x, t) as (3.35)

1 R(x, t) := i ∂t ψ̂(x, t) + ∇2 ψ̂(x, t) − Vd (x)ψ̂(x, t) − βd |ψ̂(x, t)|2 ψ̂(x, t), 2

f2 J : and its stochastic-dimension estimator by replacing ∇2 with ∇ (3.36)

f2 J ψ̂(x, t) − Vd (x)ψ̂(x, t) − βd |ψ̂(x, t)|2 ψ̂(x, t). e J (x, t) := i ∂t ψ̂(x, t) + 1 ∇ R 2

Z. Y. LIANG

By (3.33), this residual estimator is unbiased: h i e J (x, t) = R(x, t). (3.37) EJ R On a collocation set X ∈ RNc ×d , the deterministic coefficient ODE (3.3) can be f2 J [Φ(X), 1N ]⊤ : randomized by substituting ∇2 [Φ(X), 1Nc ]⊤ with ∇ c eJ (X, C(t)) [Φ(X), 1N ]⊤ + , (3.38) Ct (t) = R c h1 f2 J [Φ(X), 1N ]⊤ − Vd (X)⊤ ⊙ C(t)[Φ(X), 1N ]⊤ eJ (X, C(t)) := i C(t) ∇ R c c 2 i − βd |C(t)[Φ(X), 1Nc ]⊤ |2 ⊙ C(t)[Φ(X), 1Nc ]⊤ . For each fixed (X, C(t)), the stochastic right-hand side is unbiased: h i eJ (X, C(t)) = R(X, C(t)). (3.39) EJ R In practice, the variance induced by stochastic dimension sampling can be controlled by increasing |J|, averaging multiple independent samples per RHS evaluation, or holding J fixed within each solver step/stage and resampling at prescribed times. SD-FSNN is unbiased. We show that SD-FSNN provides an unbiased estimator for sampling the spatial Laplacian weights under the space-time separation. As a consequence, the induced estimators of the GPE right-hand side and the projected ODE vector field are also unbiased, and we provide a complete proof below. Let ψ̂(x, t) = C(t)[Φ(x), 1]⊤ be the Frozen-PINN ansatz. For a fixed time t and fixed frozen hidden-layer parameters (W, b), define the exact Laplacian ∇2 ψ̂(x, t) =

(3.40)

d X ∂ 2 ψ̂ j=1

∂x2j

(x, t).

Let J ⊂ {1, . . . , d} be a random index set of size |J| = m sampled uniformly without replacement. Define the SD-FSNN Laplacian estimator f2 J ψ̂(x, t) := d ∇ m

(3.41)

X ∂ 2 ψ̂ j∈J

∂x2j

(x, t).

Using the collocation batch X ∈ RNc ×d , define the exact right-hand side in (3.3) as h i (3.42) R(X, C(t)) := i 12 ∇2 ψ̂(X, t) − Vd (X)⊤ ⊙ ψ̂(X, t) − βd |ψ̂(X, t)|2 ⊙ ψ̂(X, t) , and the ODE vector field (3.43)

F (C(t)) := R(X, C(t)) [Φ(X), 1Nc ]⊤

f2 J ψ̂, as Define the SD-FSNN versions by replacing ∇2 ψ̂ with ∇ (3.44) h i f2 J ψ̂(X, t) − Vd (X)⊤ ⊙ ψ̂(X, t) − βd |ψ̂(X, t)|2 ⊙ ψ̂(X, t) , eJ (X, C(t)) := i 1 ∇ R 2

eJ (X, C(t)) [Φ(X), 1N ]⊤ + . (3.45) FeJ (C(t)) := R c

STOCHASTIC-DIMENSION FROZEN SAMPLED NEURAL NETWORK

Theorem 3.2 (Unbiasedness of SD-FSNN). Fix (W, b), the collocation batch X, and time t (hence C(t)). Then the SD-FSNN Laplacian estimator is unbiased in the sense that h i f2 J ψ̂(x, t) = ∇2 ψ̂(x, t). (3.46) EJ ∇ Consequently, the SD-FSNN estimators of the GPE right-hand side and the ODE vector field are unbiased as well (3.47)

eJ (X, C(t))] = R(X, C(t)), E J [R

EJ [FeJ (C(t))] = F (C(t)).

Moreover, the same unbiasedness holds componentwise for the coupled real-valued system Cre (t), Cim (t) obtained by taking Re(·) and Im(·). 2

Proof. We first prove (3.46). For each j ∈ {1, . . . , d}, define aj (x, t) := ∂∂xψ̂2 (x, t). j Pd f2 J ψ̂(x, t) = d P Then ∇2 ψ̂(x, t) = j=1 aj (x, t) and ∇ a (x, t). Let δ := 1{j∈J} j j∈J j m be the inclusion indicator. Then     d d X X X (3.48) EJ  aj  = EJ  δj aj  = EJ [δj ] aj . j=1

j∈J

j=1

Since J is a uniform subset of size m, P(j ∈ J) = m/d, hence EJ [δj ] = m/d. Therefore,

(3.49)

d d h i X d Xm f2 J ψ̂(x, t) = d EJ [δj ] aj (x, t) = aj (x, t) EJ ∇ m j=1 m j=1 d

d X

aj (x, t) = ∇2 ψ̂(x, t),

j=1

which proves (3.46). eJ differs from R only by replacing ∇2 ψ̂ with its SD-FSNN estimator, and Next, R the mapping in (3.42) is affine in ∇2 ψ̂ (the other terms Vd ⊙ ψ̂ and |ψ̂|2 ⊙ ψ̂ do not involve dimension summation). Thus, applying linearity of expectation together with eJ (X, C(t))] = R(X, C(t)). (3.46) pointwise on X yields EJ [R + Finally, [Φ(X), 1Nc ]⊤ is deterministic given (W, b) and X, hence

(3.50)

eJ (X, C(t))] [Φ(X), 1N ]⊤ + EJ [FeJ (C(t))] = EJ [R c + = R(X, C(t)) [Φ(X), 1Nc ]⊤ = F (C(t)),

which proves (3.47). The real/imaginary coupled system follows immediately because Re(·) and Im(·) are linear maps, so unbiasedness is preserved componentwise. Remark 3.3 (What is unbiased estimation). The theorem establishes unbiasedness of the instantaneous spatial-operator estimate, the GPE right-hand side R, and the ODE vector field F (C(t)) at any fixed time t. It does not, in general, imply that the integrated ODE solution trajectory C(t) is an unbiased estimate of the trajectory obtained with the full Laplacian, due to nonlinearity of ODE flows.

Z. Y. LIANG

3.5.1. Well-posedness of the semi-discrete SD-FSNN system. To ensure that the semi-discrete approximation of the GPE precludes non-physical finite-time blow-up induced by the cubic nonlinearity |ψ|2 ψ, we establish the rigorous global well-posedness of the resulting ordinary differential equation (ODE) system. Let K = r denote the dimension of the reduced basis after the truncated SVD layer. For the collocation points X ∈ RNc ×d , let Ψr (X) ∈ RK×Nc be the orthogonalized feature matrix. For simplicity of notation, let Ψ = Ψr (X). The continuous-time semi-discrete ODE system for the coefficient vector C(t) ∈ C1×K can be abstracted as the initial value problem (IVP): (3.51)

dC(t) = F(C(t)) := L(C(t)) + N (C(t)), dt

C(0) = C0 ,

where the linear mapping L : C1×K → C1×K and the nonlinear mapping N : C1×K → C1×K are defined as: 1 2 (3.52) L(C) = i C∇ Ψ − CΨdiag(Vd (X)) Ψ+ , 2 N (C) = −iβd |CΨ|2 ⊙ (CΨ) Ψ+ , (3.53) with Ψ+ denoting the Moore-Penrose pseudo-inverse. Lemma 3.4 (Local Lipschitz Continuity). For any finite radius R > 0, the vector field F(C) is locally Lipschitz continuous on the open ball BR = {C ∈ C1×K : ∥C∥2 ≤ R}. Proof. Let ∥ · ∥2 denote the Euclidean norm for vectors and the induced spectral norm for matrices. Since the activation function σ = tanh and the Gaussian envelope ρ(x) possess uniformly bounded continuous derivatives up to the second order, the linear operator satisfies ∥L(C1 ) − L(C2 )∥2 ≤ ∥L∥2 ∥C1 − C2 ∥2 with a strictly finite constant ∥L∥2 . For the nonlinear operator, let Uk = Ck Ψ ∈ C1×Nc for k ∈ {1, 2}. For any C ∈ BR , the discrete wave function satisfies the norm bound ∥Uk ∥∞ ≤ ∥Uk ∥2 ≤ ∥Ck ∥2 ∥Ψ∥2 ≤ R∥Ψ∥2 . Utilizing the algebraic inequality |z1 |2 z1 − |z2 |2 z2 ≤ 3 max(|z1 |, |z2 |)2 |z1 − z2 | applied element-wise, we bound the difference: (3.54)

|U1 |2 ⊙ U1 − |U2 |2 ⊙ U2 2 ≤ 3 max(∥U1 ∥∞ , ∥U2 ∥∞ )2 ∥U1 − U2 ∥2 ≤ 3R2 ∥Ψ∥32 ∥C1 − C2 ∥2 ,

Consequently, substituting this into (3.53) yields: (3.55)

∥N (C1 ) − N (C2 )∥2 ≤ 3|βd |∥Ψ+ ∥2 ∥Ψ∥32 R2 ∥C1 − C2 ∥2 .

Thus, F (C) is locally Lipschitz continuous with the explicit constant LR = ∥L∥2 + 3|βd |∥Ψ+ ∥2 ∥Ψ∥32 R2 < ∞. By the Picard-Lindelöf theorem, (3.51) admits a unique local solution C(t) ∈ C 1 ([0, Tmax ); C1×K ). Theorem 3.5 (Global Stability of the Projected Numerical Scheme via A Priori Bounds). Let G = ΨW Ψ∗ ∈ CK×K be the mass Gram matrix with strictly positive quadrature weights W = diag(w1 , . . . , wNc ). Given a bounded initial condition C 0 , the numerical trajectory {C n }∞ n=0 generated by the semi-discrete SD-FSNN system, subject to the discrete mass normalization projection (3.22), is unconditionally bounded, strictly preventing numerical finite-time blow-up.

STOCHASTIC-DIMENSION FROZEN SAMPLED NEURAL NETWORK

Proof. To establish global stability for the numerical scheme, we must derive a uniform a priori bound for ∥C n ∥2 across all discrete time steps. The SVD truncation explicitly guarantees that the reduced basis matrix Ψ has full row rank. Because W is positive definite, the Hermitian Gram matrix G is strictly positive-definite, possessing a minimum eigenvalue λmin (G) > 0. The explicit mass normalization projection rigorously constrains the discrete-time trajectory to the invariant discrete mass manifold at the end of each integration step tn :

M n = C n G(C n )∗ ≤ M 0 ,

(3.56)

∀n ≥ 0.

Applying the Rayleigh quotient inequality to the Hermitian matrix G, we obtain:

λmin (G)∥C n ∥22 ≤ C n G(C n )∗ .

(3.57)

Combining this with (3.56) yields a strict, time-independent global a priori bound:

s (3.58)

∥C ∥2 ≤

M0 := Bmax < ∞, λmin (G)

∀n ≥ 0.

Because the projected numerical sequence {C n } is perpetually confined to the compact domain B Bmax , the trajectory cannot diverge to infinity (supn ∥C n ∥2 ≤ Bmax ). Consequently, the projected numerical scheme remains globally stable and uniformly bounded for arbitrarily long time integrations. Remark 3.6 (Stability of the Stochastic Estimator). The established well-posedness natively extends to the randomized SD-FSNN framework (subsection 3.5). When the f2 J Ψ, the exact Laplacian ∇2 Ψ is replaced by the unbiased stochastic estimator ∇ e randomized linear operator LJ remains almost surely bounded because the spatial derivatives σ ′′ are bounded and the subset J ⊂ {1, . . . , d} is finite. Crucially, the Gram matrix G and the invariant manifold bounding radius Bmax depend solely on the basis Ψ and weights W , which are entirely independent of the stochastic dimension sub-sampling. Therefore, the uniform a priori bound unconditionally holds path-wise for any sequence of subsets J, robustly preventing stochastic divergence.

3.6. Training procedures and implementation. Having established the network structure and the well-posedness of the semi-discrete system, we outline the initialization and step-by-step implementation of SD-FSNN. Because our method strictly circumvents gradient-based optimization, we do not require a traditional loss function or iterative back-propagation. The complete procedure is summarized in Algorithm 3.1.

Z. Y. LIANG

Algorithm 3.1 SD-FSNN Training Input: Test grid points Xtest , evaluation times Ttest , dimension batch size m ≪ d, Gaussian envelope ρ(x) Output: Predicted PDE solution ψ̂(Xtest , Ttest ) Parameters: Nc , M , ϵSVD 1: Sample Nc collocation points X ∈ RNc ×d according to the envelope density ρ(x)2 . 2: Initialize frozen parameters {W, b} via the SWIM strategy. 3: Compute augmented feature matrix: Ψ(X) = [Φ(X), 1Nc ]⊤ ∈ R(M +1)×Nc . 4: Compute truncated SVD: Vr Σr Ur⊤ ≈ Ψ(X) using threshold ϵSVD . 5: Obtain reduced orthogonal bases: Ψr (X) := Vr⊤ Ψ(X) ∈ Rr×Nc .

⊤ ρ(X)−1 ⊙ ψ0 (X) Ψr (X)+ . p 7: Apply initial mass normalization: C(0) ← C(0)/ C(0)GC(0)∗ . 8: for each ODE integrator step t → t + ∆t do 9: Uniformly sample a random dimension subset J ⊂ {1, . . . , d} with |J| = m. f2 J Ψr (X) using (3.32). 10: Compute the unbiased stochastic Laplacian estimator: ∇ 6: Initialize coefficients via least-squares: C(0) =

11: 12: 13:

e J (X, C(t)) ((3.38)). Evaluate the stochastic ODE right-hand side R Evolve unnormalized coefficients C̃(t + ∆t) using a suitable implicit symplectic integrator or splitting method ((3.4)). C̃(t+∆t) Apply mass-normalization projection: C(t + ∆t) = √ . ∗ C̃(t+∆t)GC̃(t+∆t)

14: end for 15: Return ψ̂(Xtest , Ttest ) = ρ(Xtest ) ⊙

C(Ttest )Ψr (Xtest ) .

4. Numerical experiments. In this section, we demonstrate the numerical performance of the proposed SD-FSNN method on high-dimensional GPEs over unbounded domains. We evaluate short-time dynamical behavior on benchmarks ranging from one to three spatial dimensions, assess the long-time accuracy with strict preservation of physical invariants, and investigate both in-distribution generalization and out-of-distribution transfer capability. All baseline experiments are conducted on an Nvidia A100-SXM4-80GB GPU. Code and data for the following experiments are available at https://github.com/liangzhangyong/SD-FSNN. As previously discussed, we consider the spatio-temporal domain Ω × [0, T ] ⊂ Rd × R, where d represents the spatial dimension. To evaluate the solution predicted by SD-FSNN, we use two error metrics: relative L1 error (rMAE) and relative Root Mean Square Error (rRMSE). We compute the test error on a uniform grid for all PDEs with 256 points in space and 100 points in time, unless otherwise specified. We use float64 numerical precision in all the experiments. For ground truth ψ and model prediction ψθ , these two metrics can be given by v v u PN u PN 2 u u |ψθ (xi , ti ) − ψ(xi , ti )| (ψθ (xi , ti ) − ψ(xi , ti )) rMAE:t i=1 PN , rRMSE:t i=1 PN . 2 i=1 |ψ(xi , ti )| i=1 (ψ(xi , ti )) where {xi }N i=1 are selected collocation points for evaluation. For each experiment, the mean and standard deviation of the rMAE and the rRMSE are computed with three seeds. 4.1. Short-time dynamical behavior.

STOCHASTIC-DIMENSION FROZEN SAMPLED NEURAL NETWORK

Example 4.1 (1D GPE). First, we consider the 1D GPE to validate the baseline accuracy of the SD-FSNN method. By naturally incorporating the Gaussian-weighted ansatz into the random feature network, SD-FSNN seamlessly enforces the exponential decay at infinity, thus avoiding the boundary reflection artifacts typical of truncated domains. The numerical results demonstrate that SD-FSNN successfully captures the short-time wave function evolution, achieving highly competitive rMAE and rRMSE compared to ground truth solutions obtained via refined Hermite spectral methods. Figure 4.1 presents the space-time evolution of the real part of the wave function Re(ψ) for the 1D GPE with β = 10. The left panel shows the reference solution obtained from a refined Hermite spectral method, the middle panel displays the SD-FSNN prediction, and the right panel depicts the pointwise absolute error. The SD-FSNN prediction is visually indistinguishable from the reference solution, accurately reproducing the spatial density profile concentrated near x = 0 and its temporal evolution over t ∈ [0, 1]. The absolute error remains at the level of 10−7 throughout the entire space-time domain, with slightly larger errors appearing near the center of the condensate at later times due to the accumulation of nonlinear interaction effects. These results confirm that SD-FSNN achieves high-fidelity approximation of the 1D GPE dynamics on the unbounded domain. We further investigate the accuracy of SD-FSNN across a wide range of interaction parameters β ∈ [10−2 , 102 ] and compare it with PINN, Causal PINN, and ELM. Figure 4.2 presents the relative L2 error as a function of β on a log-log scale. For moderate interaction strengths (β ≤ 10), PINN, Causal PINN, and ELM all exhibit relative errors in the range of 10−5 –10−4 , with ELM performing slightly better than the two PINN variants. As β increases beyond 10, the errors of all three baseline methods rise sharply, reaching 10−1 at β = 100, indicating a significant degradation in accuracy for strongly interacting condensates. In contrast, SD-FSNN consistently achieves relative L2 errors that are one to two orders of magnitude lower than all baselines across the entire range of β. Even at β = 100, SD-FSNN maintains an error below 10−5 , demonstrating robust convergence behavior regardless of the interaction strength. This superior performance can be attributed to the structure-preserving design of SD-FSNN: the Gaussian-weighted ansatz naturally captures the solution’s decay profile, while the mass normalization projection and energy conservation constraint stabilize the computation even in the strongly nonlinear regime. Example 4.2 (2D GPE). First, we extend our evaluation to the 2D GPE, incorporating an anisotropic harmonic potential to illustrate the method’s capability in handling multi-dimensional spatial dynamics. The 2D case introduces more complex interference patterns. By leveraging the space-time separation and randomized spatial basis, SD-FSNN efficiently approximates the 2D wave function. Unlike traditional grid-based methods, the dimension-independent computational structure effectively captures the anisotropic spread of the density profile across the unbounded domain, yielding stable and accurate short-time predictions. Figure 4.3 reports snapshots of the solution |ψ|2 on [−1, 1]2 at t ∈ {0, 0.3, 0.6, 1} for β = 10. Each row shows, from left to right, the reference solution, the SD-FSNN prediction, and the pointwise absolute error. The predicted fields are visually indistinguishable from the ground truth at all sampled times; the spatial error is negligible at t = 0 and remains on the order of 10−6 for t ≥ 0.3 without visible blow-up through t = 1. These panels corroborate that SD-FSNN resolves two-dimensional condensate dynamics with high fidelity over the short-time horizon and that the error growth in time stays mild, consistent with the structure-preserving time integration of the coef-

Z. Y. LIANG Reference

6 4

0.4

Prediction

6 0.6

0.2

0.4

0.2

1.2 1.0

0.8

-2

-0.2 -2

-4

-0.4 -4

-0.6

0.4

0.6

0.8

1.0

-6 0.0

0.0

0.2

0.4

0.6

0.8

1.0

1.4

0.0

0.2

1e-7

-6 0.0

Absolute Error

6 0.6

-6 0.0

0.6 0.4 0.2

0.2

0.4

0.6

0.8

1.0

0.0

(a) Re(ψ) Reference

0.6

0.4

Prediction

0.6

0.4

Absolute Error

1e-7 1.75

1.50

0.2

1.25

0.0

1.00

-2

-0.2 -2

0.75

-0.4

-4 -6 0.0

-0.6

0.2

0.4

0.6

0.8

1.0

-4 -6 0.0

-0.6

0.2

0.4

0.6

0.8

1.0

0.50

-4 -6 0.0

0.25

0.2

0.4

0.6

0.8

1.0

0.00

(b) Im(ψ) Reference

0.7

Prediction

0.7

Absolute Error

1e-7 1.75

0.6 4

0.5 2

1.25

0.4 0 0.3

1.00

0.2

0.50

-2 -4 -6 0.0

-2

0.1 -4

0.2

0.4

0.6

0.8

1.0

0.0 -6

0.0

1.50

0.75

-2

0.1 -4

0.2

0.4

0.6

0.8

1.0

0.0 -6

0.0

0.25

0.2

0.4

0.6

0.8

1.0

0.00

Fig. 4.1: SD-FSNN results for the 1D GPE with β = 10: reference solution (left), SD-FSNN prediction (middle), and absolute error (right) for (a) the real part, (b) the imaginary part, and (c) the norm part.

Relative

2 error

PINN

Causal PINN

ELM

SDRE

101

102

10 1 10 6 10 2

10 1

100

Fig. 4.2: Relative L2 error as a function of the interaction parameter β for PINN, Causal PINN, ELM, and SD-FSNN.

STOCHASTIC-DIMENSION FROZEN SAMPLED NEURAL NETWORK

ficient ODE.

t = 0.00

1.0

Ground Truth

0.560

1.0

0.5

0.480

0.0

0.400 0.320

-0.5 -1.0 -1.0 1.0

0.240 -0.5

0.0

0.5

1.0 0.480

t = 0.30

0.5

0.400

0.0

0.320

-0.5

t = 0.60

-1.0 -1.0 1.0

0.240

0.5

0.480

0.5

0.0

0.400

0.0

-1.0 -1.0 1.0

0.200

-1.0 -1.0

0.200 -0.5

0.0

0.5

1.0

0.240

-0.5

0.0

0.5

1.0 1.2e-06

0.5

9e-07

0.0

6e-07

-0.5

3e-07

4.5e-06

0.0

0.400

0.0

3e-06

-0.5

1.5e-06

-0.5

0.0

0.5

1.0

0.300

-0.5

0.400 0.300

0.320

-1.0 -1.0 1.0

3e-13

0.5

0.0 -0.5

0.400

6e-13

-0.5

0.500

1.0

0.480

9e-13

0.5

1.0

1.2e-12

0.600

0.0

0.5

Absolute Error

-1.0 -1.0 1.0

-0.5

0.0

-0.5

0.400 0.300

-0.5

0.0

1.0

0.240

0.5

0.500

0.5

0.320

-0.5

0.5

-1.0 -1.0 1.0

1.0

0.600

0.0

-0.5

t = 1.00

0.560

-1.0 -1.0 1.0

-0.5

Prediction

0.200 0.600

-1.0 -1.0 1.0

0.5

0.500

0.5

0.0

0.400

0.0

-0.5

0.0

0.5

1.0

0.300

-0.5 -1.0 -1.0

0.200 -0.5

0.0

0.5

1.0

-0.5

0.0

0.5

1.0 6e-06

-0.5

0.0

0.5

1.0 8e-06 6e-06 4e-06

-0.5 -1.0 -1.0

2e-06 -0.5

0.0

0.5

1.0

Fig. 4.3: SD-FSNN results for the 2D GPE with β = 10: at each time t ∈ {0, 0.3, 0.6, 1}, reference solution (left), SD-FSNN prediction (middle), and absolute error (right) for the solution |ψ|2 on [−1, 1]2 .

Example 4.3 (3D GPE). First, we assess the 3D GPE, a regime where traditional Hermite spectral methods face severe memory and computational bottlenecks due to the exponential scaling of basis functions. By employing the SD-FSNN framework, the stochastic dimension sampling mechanism bypasses this curse of dimensionality. The method accurately resolves the 3D condensate profile within the short-time regime, confirming its scalability and practical utility for physically realistic 3D whole-space settings. Figure 4.4 shows subsampled spatial points in [−4, 4]3 at t ∈ {0, 0.3, 0.6, 1} for β = 10. Each row displays the reference solution, the SD-FSNN prediction, and the pointwise absolute error. The predicted point clouds match the reference in both support and amplitude at every time slice, with a concentrated high-intensity core and decay toward the tails of the trap. The pointwise absolute error is largest in the central high-density region, where the field magnitude is greatest. At t = 0 the peak error is on the order of 10−8 , consistent with an accurate representation of the initial state on the subsampled set. For t ∈ {0.3, 0.6, 1} the peak error increases gradually to the order of 10−6 and 10−5 yet remains uniformly small through t = 1, without visible runaway

Z. Y. LIANG

growth. These results indicate that SD-FSNN captures three-dimensional short-time dynamics with high fidelity while avoiding the cost of full tensor-product spectral grids.

Fig. 4.4: SD-FSNN results for the 3D GPE with β = 10: at each time t ∈ {0, 0.3, 0.6, 1}, reference solution (left), SD-FSNN prediction (middle), and absolute error (right) for the solution |ψ|2 on subsampled points in [−4, 4]3 .

To further highlight the computational advantage of SD-FSNN, we compare its running time against the fourth-order time-splitting spectral (TSSP) method across dimensions d = 1, . . . , 8 for two time horizons T = 1 and T = 10. Figure 4.5 presents the elapsed time on a logarithmic scale. For TSSP, we test four temporal step sizes h ∈ {0.5, 0.25, 0.125, 0.0625}. As the dimension increases, the computational cost of TSSP grows exponentially due to the tensor-product grid structure, reaching 105 –106 s at d = 8 regardless of the step size. Even with the coarsest step size h = 0.5, TSSP requires over 104 s at d = 7. In contrast, SD-FSNN maintains a nearly flat cost curve across all dimensions, staying below 15 s for both T = 1 and T = 10. At d = 8 with T = 10, SD-FSNN is approximately four to five orders of magnitude faster than TSSP. This dramatic speedup stems from the dimension-independent computational structure of SD-FSNN: the randomized spatial basis avoids constructing high-dimensional tensor-product grids, and the space-time separation reduces time integration to a lowdimensional ODE system. These results confirm that SD-FSNN effectively overcomes

STOCHASTIC-DIMENSION FROZEN SAMPLED NEURAL NETWORK

the curse of dimensionality that fundamentally limits traditional spectral methods for GPEs on unbounded domains.

T=1 106

106

104 103 102 101 10

T=10, h=0.5 T=10, h=0.25 T=10, h=0.125 T=10, h=0.0625 SDRE

107

Elapsed Time (s)

105

Elapsed Time (s)

T=10

T=1, h=0.5 T=1, h=0.25 T=1, h=0.125 T=1, h=0.0625 SDRE

105 104 103 102 101

100 1

Dimension (d)

Fig. 4.5: Elapsed time comparison between TSSP (with various step sizes h) and SD-FSNN across spatial dimensions d = 1, . . . , 8 for T = 1 (left) and T = 10 (right). The y-axis is in logarithmic scale.

4.2. Long-time dynamical behavior. In this section, we investigate the performance of SD-FSNN over extended temporal domains. Standard global-in-time neural network solvers often fail to learn long-time dynamics due to error accumulation and a lack of causal structure. Our space-time separation approach, combined with symplectic ODE integration and mass normalization, natively addresses these issues, ensuring physical consistency and mitigating numerical dissipation over arbitrarily long periods. Figure 4.6 displays the reference solution, the SD-FSNN prediction, and the pointwise absolute error in the (t, x) plane for t ∈ [0, 100] and x ∈ [−6, 6]. The prediction matches the reference to visual accuracy across the full time window, preserving the localized oscillatory structure near the trap center. The error panel remains uniformly small at the level of about 10−7 without systematic growth toward t = 100, which indicates stable long-time integration without visible drift or blow-up in the spatial–temporal field. This qualitative behavior is consistent with projected symplectic evolution of the coefficients together with discrete mass normalization suppressing unbounded amplitude and phase errors over long horizons. Figure 4.7 summarizes the same experiment through the relative L2 error of the wave function predicted by SD-FSNN against the reference solution as a function of time T ∈ [0, 100]. After a brief initial transient in which the error rises from machineprecision levels to roughly 10−7 , the curve settles into a nearly flat regime. From T ≈ 5 onward the error grows only slowly and remains of order 10−6 at T = 100, without visible exponential blow-up. Taken together, Figures 4.6 and 4.7 support the claim that the dynamics predicted by SD-FSNN remain accurate and stable for long-time integration, in contrast to typical global-in-time neural solvers that exhibit rapid error growth.

Z. Y. LIANG Reference

0.7

Prediction

0.7

0.6 4

0.5 2

0.4 0 0.3

0.2

-2

0.1 -4

-6

100

1e-7 8 6 4

-2

-4 0

Absolute Error

100

Relative 2 error

Fig. 4.6: Long-time integration on t ∈ [0, 100] and x ∈ [−6, 6]. Reference solution (left), SD-FSNN prediction (middle), and pointwise absolute error (right).

10 6 10 7 10 8 10 9 10 10 0

50 Time T

100

Fig. 4.7: Relative L2 error of the solution predicted by SD-FSNN versus time T .

4.3. Combination coefficients across β. After space-time separation the wave function admits the discrete expansion (4.1)

ψ(x, t) =

n−1 X

ci (t) ϕi (x),

i=0 n−1 where {ϕi }n−1 i=0 are n fixed spatial features and {ci (t)}i=0 are the corresponding timedependent combination coefficients. For the 1D GPE we use n = 83 basis functions; the full set {ci (t)}82 i=0 is governed by the coupled nonlinear ODE system derived from the GPE residual. Figure 4.8 displays c3 (t), . . . , c8 (t) on t ∈ [0, 1] for β ∈ {1, 10, 30, 50, 70, 100} (a 2 × 3 panel). Modes i = 0, 1, 2 carry the dominant projection of the initial data and evolve smoothly, so indices 3–8 are selected to expose the mid-spectrum dynamics driven by nonlinear interactions. For β = 1 the coefficient curves are low-amplitude and slowly varying; the condensate behaves near the linear regime and the nonlinear term |ψ|2 ψ perturbs the basis weights only weakly. As β increases, the nonlinear self-focusing strengthens. The coefficients ci (t) transition from nearly monochromatic oscillations at β = 10 to more vigorous oscillations with larger peak amplitudes at β = 30 and β = 50, reflecting the growing contribution of higher spatial frequencies to the condensate profile. At β = 70

STOCHASTIC-DIMENSION FROZEN SAMPLED NEURAL NETWORK

and β = 100 the trajectories acquire irregular, multi-frequency structure: different ci (t) acquire distinct phase offsets, and some coefficients that are near-zero at small β become active, indicating that the nonlinear coupling coherently transfers energy from the leading modes into the higher-order spatial spectrum. These observations confirm that the coefficient ODE system captures the physically expected transition from linear harmonic evolution to strongly nonlinear condensate dynamics as β grows. β=1

0.0

0.4

1.0

0.2

0.5

0.0

0.2

0.4

0.6

0.8

1.0

-0.4

0.0

0.2

0.4

0.6

0.8

1.0

0.0

β = 70 2

1 0

c(t)

-1.0

β = 50

0.0 -0.5

-0.2

-1

-1 -2

β = 30

c(t)

0.2

c(t)

0.4

c(t)

β = 10

c3 c4 c5 c6 c7 c8

-2 0.0

0.2

0.4

0.6

0.8

1.0

0.0

0.2

0.4

0.6

0.8

1.0

3 2 1 0 -1 -2 -3

0.2

0.4

0.6

0.8

1.0

0.8

1.0

β = 100

0.0

0.2

0.4

0.6

Fig. 4.8: Time evolution of temporal basis coefficients under various β.

4.4. Ablation study. To ascertain the source performance gain of our method, we conduct a detailed ablation study on the core components of SD-FSNN: (1) the space-time separation, (2) the stochastic dimension subsampling, (3) the mass normalization projection, and (4) the energy-conserving symplectic integrator. Removing the space-time separation forces a standard global optimization that collapses during long-time evolution. Omitting the dimension subsampling drastically increases the computational cost per step in high dimensions. Without mass normalization and symplectic integration, the system exhibits severe numerical dissipation and amplitude drift over extended timeframes. These findings confirm that each structural prior is strictly necessary for the stability, accuracy, and efficiency of SD-FSNN. 5. Extended discussion of SD-FSNN. In this section, we provide an extended discussion on the theoretical and practical implications of the SD-FSNN framework. We detail its application to extremely high-dimensional nonlinear GPEs, provide comprehensive comparisons with existing general high-dimensional solvers, and address theoretical aspects such as the Kolmogorov n-width barrier and mesh dependency. 5.1. Extremely high-dimensional GPEs. In this section, we apply the Method of Manufactured Solutions (MMS) to construct exact solutions for high-dimensional nonlinear GPEs. This method introduces a source term to the original GPE, creating a modified equation whose exact solution is known by construction. Here, we aim to

Z. Y. LIANG

consider nonseparable and anisotropic solutions for nonlinear GPEs to form complicated and nontrivial high-dimensional GPEs. The resulting exact solution is of the following form  (5.1)

ψexact (x) = 

d Y γj

j=1

1/4 1



e− 2

2 j=1 γj xj

F(x),

where the function F(x) is given by " # d−1 X F(x) = 1 + ϵ ci sin (xi + cos(xi+1 ) + xi+1 cos(xi )) , i=1

with ϵ is a small perturbation term and ci ∼ N (0, 1) are random variables. In addition to the exact solution, the following GPEs defined within the unit ball Bd associated with zero boundary conditions on the unit sphere are under consideration, as follows (5.2)

1 2 ∇ ψ = Vd (x)ψ + βd |ψ|2 ψ + g(x), 2

x ∈ Bd .

where g(x) = 12 ∇2 ψexact − Vd (x)ψexact − βd |ψexact |2 ψexact . 5.2. Comparison with general high-dimensional solvers. We compare SDFSNN with five representative gradient estimation methods discussed in Section ??, highlighting the key distinctions. • SDGD [26] samples a minibatch of input dimensions but incurs O(2k−1 ) graph growth for derivative order k. SD-FSNN bypasses automatic differentiation via explicitly derived analytical derivatives of frozen spatial features. • RS-PINNs [27] introduce bias from nonlinearities in the PDE residual, and bias-correction increases variance in high dimensions. SD-FSNN evaluates exact second derivatives element-wise, yielding an exactly unbiased Laplacian estimate without smoothing-induced bias (Theorem ??). • HTE [25] requires two backward passes per iteration and does not generalize to higher-order derivatives. SD-FSNN constructs the Laplacian via direct linear and Hadamard products, eliminating graph-based autodiff overhead. • STDE [48] relies on algorithmic differentiation passes and iterative gradientbased optimization, which can suffer from slow convergence. SD-FSNN restricts the ansatz to a single hidden layer with explicitly known derivatives, reducing the problem to linear algebra operations and ODE integration. • Forward Laplacian [33] scales linearly with the ambient dimension d. SDFSNN exploits the frozen random feature structure to evaluate the Laplacian as a direct linear transformation ∆û = C(t) (W ⊙ W ) ⊙ σ ′′ (z) 1, enabling high-order ODE integration with per-sample cost controlled by a subsample size m independent of d. To provide a concrete quantitative comparison, we evaluate SD-FSNN against SDGD, STDE, HTE, and RS-PINN on the static high-dimensional GPE benchmark across dimensions d ∈ {10, 50, 100, 200, 500, 1000}. Figure 5.1 presents the running time and L2 error as functions of the spatial dimension d. In the left panel, the training times of all baseline methods (SDGD, HTE, RS-PINN, Score PINN, STDE) grow substantially with d, reaching 150–300 s

STOCHASTIC-DIMENSION FROZEN SAMPLED NEURAL NETWORK

at d = 10000. In contrast, SD-FSNN remains nearly flat across all dimensions, staying below 15 s even at d = 10000, confirming that its computational cost is essentially independent of the ambient dimension. In the right panel, the L2 errors of all baseline methods increase steadily with d, exceeding 40%–55% at d = 10000. SD-FSNN, however, maintains a consistently low error that grows only mildly, remaining below 5% across the entire range of dimensions tested. These results demonstrate that SD-FSNN achieves both dimension-independent computational cost and superior accuracy compared to existing general high-dimensional solvers. Running Time Comparison

250

L2 Error (%)

Time (s)

SDGD HTE RS_PINN Score_PINN STDE SDRE

200 150

100

Accuracy Comparison

SDGD HTE RS_PINN Score_PINN STDE SDRE

300

100

1000

10000

100

1000

10000

Fig. 5.1: Running time (left) and L2 error (right) comparison of SD-FSNN against baseline methods across spatial dimensions d ∈ {4, 10, 20, 50, 100, 1000, 10000}.

Several observations can be drawn from Table 5.1. First, SD-FSNN achieves the lowest L2 and L1 errors across all dimensions by a wide margin. For instance, at d = 200, SD-FSNN attains an L2 error of 8.80 × 10−5 , which is approximately 39× lower than SDGD/STDE (3.42 × 10−3 ), 63× lower than HTE (5.52 × 10−3 ), and over 400× lower than RS-PINN (3.71 × 10−2 ). Even at d = 1000, SD-FSNN maintains an L2 error of 3.45 × 10−4 , still more than one order of magnitude below the best competing method (SDGD/STDE at 7.43 × 10−3 ). Second, among the baseline methods, SDGD and STDE achieve nearly identical errors and consistently outperform HTE and RS-PINN. HTE ranks as the next most accurate, with L2 errors typically 1.5–2× larger than those of SDGD/STDE, while RS-PINN exhibits the largest errors, particularly in high dimensions where its L2 error can exceed that of SDGD by nearly an order of magnitude (e.g., 4.86 × 10−2 vs. 7.35 × 10−3 at d = 100). Third, regarding computational cost, SD-FSNN is by far the fastest method, requiring only 2.12 s at d = 10 and 12.18 s at d = 1000, which is orders of magnitude faster than all competing methods. By comparison, RS-PINN is the fastest among the baselines due to its avoidance of automatic differentiation, but still requires 30.80 s at d = 10 and 696.48 s at d = 1000. Among the AD-based methods, SDGD and STDE exhibit comparable training times that scale approximately linearly with dimension, while HTE is generally faster at moderate dimensions (d ≤ 100) but becomes comparable or slower at higher dimensions. Overall, SD-FSNN simultaneously achieves the highest accuracy and the lowest computational cost across all dimensions, outperforming existing gradient-based methods (SDGD, STDE, HTE) and randomized smoothing approaches (RS-PINN) by significant margins on the static GPE benchmark.

Z. Y. LIANG

Table 5.1: Comparison of methods on the static high-dimensional GPE (5.2). Dim

Method

L2 Error

L1 Error

−3

Time (s)

10 10 10 10 10

SDGD STDE HTE RS-PINN SD-FSNN

1000 1000 1000 1000 –

– – – – 128

4.42 × 10 4.42 × 10−3 6.52 × 10−3 1.20 × 10−2 2.41 × 10−6

3.14 × 10 3.14 × 10−3 4.54 × 10−3 8.16 × 10−3 1.58 × 10−6

199.53 198.11 617.58 30.80 2.12

50 50 50 50 50

SDGD STDE HTE RS-PINN SD-FSNN

1000 1000 1000 1000 –

– – – – 128

9.16 × 10−3 9.16 × 10−3 1.01 × 10−2 2.55 × 10−2 5.28 × 10−5

7.07 × 10−3 7.07 × 10−3 7.11 × 10−3 1.59 × 10−2 3.49 × 10−5

771.94 816.41 628.55 41.67 5.68

100 100 100 100 100

SDGD STDE HTE RS-PINN SD-FSNN

1000 1000 1000 1000 –

– – – – 196

7.35 × 10−3 7.35 × 10−3 1.03 × 10−2 4.86 × 10−2 8.05 × 10−5

5.31 × 10−3 5.31 × 10−3 8.25 × 10−3 4.09 × 10−2 6.98 × 10−5

836.11 848.18 647.88 48.26 6.21

200 200 200 200 200

SDGD STDE HTE RS-PINN SD-FSNN

1000 1000 1000 1000 –

– – – – 256

3.42 × 10−3 3.42 × 10−3 5.52 × 10−3 3.71 × 10−2 8.80 × 10−5

2.42 × 10−3 2.42 × 10−3 3.98 × 10−3 2.99 × 10−2 7.95 × 10−5

779.41 807.25 682.24 80.19 7.85

500 500 500 500 500

SDGD STDE HTE RS-PINN SD-FSNN

1000 1000 1000 1000 –

– – – – 324

3.45 × 10−3 3.45 × 10−3 8.25 × 10−3 2.73 × 10−2 1.05 × 10−4

2.65 × 10−3 2.65 × 10−3 7.48 × 10−3 2.20 × 10−2 9.98 × 10−5

1031.43 999.85 767.55 131.54 9.94

1000 1000 1000 1000 1000

SDGD STDE HTE RS-PINN SD-FSNN

1000 1000 1000 1000 –

– – – – 428

7.43 × 10−3 7.43 × 10−3 1.29 × 10−2 3.05 × 10−2 3.45 × 10−4

5.87 × 10−3 5.87 × 10−3 6.33 × 10−3 1.06 × 10−2 2.58 × 10−4

1573.35 2029.26 1638.89 696.48 12.18

5.3. Kolmogorov n-width barrier. The Kolmogorov n-width measures the optimal approximation error achievable by any linear subspace of a given dimension. Traditional spectral and reduced-order models often struggle to overcome the Kolmogorov n-width barrier for complex wave dynamics or highly oscillatory PDE solutions, as mapping the full spatio-temporal dynamics requires an exponentially growing number of parameters. Because SD-FSNN relies on a separation of variables ψ̂(x, t) = C(t)[Φ(x), I] with a frozen, randomized spatial basis Φ(x), its approximation capability is also inherently limited by the Kolmogorov n-width of the solution manifold. While the randomized feature networks effectively alleviate the curse of dimensionality for moderately complex wave functions, representing highly turbulent or strongly interacting condensates (βd ≫ 1) may require an exponentially large number

STOCHASTIC-DIMENSION FROZEN SAMPLED NEURAL NETWORK

of random basis functions M to effectively capture the complex solution manifold, representing a fundamental theoretical barrier compared to fully adaptive, iteratively trained deep neural networks. 6. Conclusions. We introduce SD-FSNN, an unbiased randomized estimation method for solving high-dimensional GPEs on unbounded domains. The computational cost of SD-FSNN is independent of dimension, avoiding the exponential growth associated with Hermite-basis discretizations. The framework enforces physical consistency through a Gaussian-weighted ansatz for asymptotic decay, a projection layer for mass conservation, and an energy constraint to curb numerical dissipation. SDFSNN with amortization achieves high accuracy while reducing computational and memory burdens compared to existing methods. Beyond the scope of this paper, SD-FSNN can be extended to multi-component BECs and other high-dimensional problems on unbounded domains. It establishes high-dimensional GPEs on unbounded domains as a novel benchmark for general unbiased high-dimensional solvers. The proposed space-time separation strategy improves performance for long-time dynamics. Potential limitations include the trade-off between Monte Carlo subset size, computational cost, and variance, and the diminishing representational capability of basis functions as the nonlinear interaction coefficient increases. Future work may explore variance reduction techniques and compensation strategies for large interaction coefficients. SD-FSNN represents a significant step toward practical neural PDE solvers on standard CPU architectures. Acknowledgments. We would like to thank the referees for their insightful comments, which greatly improved the paper. This work was carried out using computing resources at the High Performance Computing Platform of Wuhan University. REFERENCES [1] J. R. Anglin and W. Ketterle, Bose–einstein condensation of atomic gases, Nature, 416 (2002), pp. 211–218. [2] W. Bao, D. Jaksch, and P. A. Markowich, Numerical solution of the gross–pitaevskii equation for bose–einstein condensation, Journal of Computational Physics, 187 (2003), pp. 318–342. [3] W. Bao, H. Li, and J. Shen, A generalized-laguerre–fourier–hermite pseudospectral method for computing the dynamics of rotating bose–einstein condensates, SIAM Journal on Scientific Computing, 31 (2009), pp. 3685–3711. [4] W. Bao and J. Shen, A fourth-order time-splitting laguerre–hermite pseudospectral method for bose–einstein condensates, SIAM Journal on Scientific Computing, 26 (2005), pp. 2010– 2028. [5] C. Beck, S. Becker, P. Cheridito, A. Jentzen, and A. Neufeld, Deep splitting method for parabolic pdes, SIAM Journal on Scientific Computing, 43 (2021), pp. A3135–A3154. [6] C. Beck, W. E, and A. Jentzen, Machine learning approximation algorithms for highdimensional fully nonlinear partial differential equations and second-order backward stochastic differential equations, Journal of Nonlinear Science, 29 (2019), pp. 1563–1619. [7] C. Beck, L. Gonon, and A. Jentzen, Overcoming the curse of dimensionality in the numerical approximation of high-dimensional semilinear elliptic partial differential equations, Partial Differential Equations and Applications, 5 (2024), p. 31. [8] C. Beck, F. Hornung, M. Hutzenthaler, A. Jentzen, and T. Kruse, Overcoming the curse of dimensionality in the numerical approximation of allen–cahn partial differential equations via truncated full-history recursive multilevel picard approximations, Journal of Numerical Mathematics, 28 (2020), pp. 197–222. [9] S. Becker, R. Braunwarth, M. Hutzenthaler, A. Jentzen, and P. von Wurstemberger, Numerical simulations for full history recursive multilevel picard approximations for sys-

Z. Y. LIANG

tems of high-dimensional partial differential equations, arXiv preprint arXiv:2005.10206, (2020). [10] J. Bettencourt, M. J. Johnson, and D. Duvenaud, Taylor-mode automatic differentiation for higher-order derivatives in jax, in Program Transformations for ML Workshop at NeurIPS 2019, 2019. [11] I. Buluta and F. Nori, Quantum simulators, Science, 326 (2009), pp. 108–111. [12] M. M. Cerimele, M. L. Chiofalo, F. Pistella, S. Succi, and M. P. Tosi, Numerical solution of the gross-pitaevskii equation using an explicit finite-difference scheme: An application to trapped bose-einstein condensates, Physical Review E, 62 (2000), p. 1382. [13] Q. Chan-Wai-Nam, J. Mikael, and X. Warin, Machine learning for semi linear pdes, Journal of scientific computing, 79 (2019), pp. 1667–1712. [14] P.-H. Chiu, J. C. Wong, C. Ooi, M. H. Dao, and Y.-S. Ong, Can-pinn: A fast physicsinformed neural network based on coupled-automatic–numerical differentiation method, Computer Methods in Applied Mechanics and Engineering, 395 (2022), p. 114909. [15] J. Cho, S. Nam, H. Yang, S.-B. Yun, Y. Hong, and E. Park, Separable pinn: Mitigating the curse of dimensionality in physics-informed neural networks, arXiv preprint arXiv:2211.08761, (2022). [16] T. Chou, S. Shao, and M. Xia, Adaptive hermite spectral methods in unbounded domains, Applied Numerical Mathematics, 183 (2023), pp. 201–220. [17] T. Dohnal, Perfectly matched layers for coupled nonlinear schrödinger equations with mixed derivatives, Journal of Computational Physics, 228 (2009), pp. 8752–8765. [18] R. Doran and T. Billam, Numerical method for the projected gross–pitaevskii equation in an infinite rotating 2d bose gas, arXiv preprint arXiv:2002.09245, (2020). [19] D. Funaro, Polynomial approximation of differential equations, Springer, 1992. [20] I. M. Georgescu, S. Ashhab, and F. Nori, Quantum simulation, Reviews of Modern Physics, 86 (2014), pp. 153–185, https://doi.org/10.1103/RevModPhys.86.153. [21] M. Greiner, O. Mandel, T. Esslinger, T. W. Hänsch, and I. Bloch, Quantum phase transition from a superfluid to a mott insulator in a gas of ultracold atoms, nature, 415 (2002), pp. 39–44. [22] B.-y. Guo, J. Shen, and C.-l. Xu, Spectral and pseudospectral approximations using hermite functions: application to the dirac equation, Advances in Computational Mathematics, 19 (2003), pp. 35–55. [23] J. Han, A. Jentzen, and W. E, Solving high-dimensional partial differential equations using deep learning, Proceedings of the National Academy of Sciences, 115 (2018), pp. 8505–8510. [24] J. Han, A. Jentzen, et al., Deep learning-based numerical methods for high-dimensional parabolic partial differential equations and backward stochastic differential equations, Communications in mathematics and statistics, 5 (2017), pp. 349–380. [25] Z. Hu, Z. Shi, G. E. Karniadakis, and K. Kawaguchi, Hutchinson trace estimation for high-dimensional and high-order physics-informed neural networks, Computer Methods in Applied Mechanics and Engineering, 424 (2024), p. 116883. [26] Z. Hu, K. Shukla, G. E. Karniadakis, and K. Kawaguchi, Tackling the curse of dimensionality with physics-informed neural networks, Neural Networks, 176 (2024), p. 106369. [27] Z. Hu, Z. Yang, Y. Wang, G. E. Karniadakis, and K. Kawaguchi, Bias-variance tradeoff in physics-informed neural networks with randomized smoothing for high-dimensional pdes, SIAM Journal on Scientific Computing, 47 (2025), pp. C846–C872. [28] G.-B. Huang, H. Zhou, X. Ding, and R. Zhang, Extreme learning machine for regression and multiclass classification, IEEE Transactions on Systems, Man, and Cybernetics, Part B (Cybernetics), 42 (2012), pp. 513–529, https://doi.org/10.1109/TSMCB.2011.2168604. [29] G.-B. Huang, Q.-Y. Zhu, and C.-K. Siew, Extreme learning machine: Theory and applications, Neurocomputing, 70 (2006), pp. 489–501, https://doi.org/10.1016/j.neucom.2005. 12.126. [30] T. Lahaye, C. Menotti, L. Santos, M. Lewenstein, and T. Pfau, The physics of dipolar bosonic quantum gases, Reports on Progress in Physics, 72 (2009), p. 126401, https://doi. org/10.1088/0034-4885/72/12/126401. [31] L. D. Landau and E. M. Lifshitz, Quantum mechanics: non-relativistic theory, vol. 3, Elsevier, 2013. [32] L. D. Landau and E. M. Lifshitz, Quantum mechanics: non-relativistic theory, vol. 3, Elsevier, 2013. [33] R. Li, H. Ye, D. Jiang, X. Wen, C. Wang, Z. Li, X. Li, D. He, J. Chen, W. Ren, et al., A computational framework for neural network-based variational monte carlo with forward laplacian, Nature Machine Intelligence, 6 (2024), pp. 209–219. [34] Y. Liao, Z. Lin, J. Liu, Q. Sun, Y. Wang, T. Wu, and H. Xie, Solving schr\”{o} dinger

STOCHASTIC-DIMENSION FROZEN SAMPLED NEURAL NETWORK

equation using tensor neural network, arXiv preprint arXiv:2209.12572, (2022). [35] S. Liu, P.-Y. Chen, B. Kailkhura, G. Zhang, A. O. Hero III, and P. K. Varshney, A primer on zeroth-order optimization in signal processing and machine learning: Principals, recent advances, and applications, IEEE Signal Processing Magazine, 37 (2020), pp. 43–54. [36] C. Lv, L. Wang, and C. Xie, A hybrid physics-informed neural network for nonlinear partial differential equation, International Journal of Modern Physics C, 34 (2023), p. 2350082. [37] S. Malladi, T. Gao, E. Nichani, A. Damian, J. D. Lee, D. Chen, and S. Arora, Fine-tuning language models with just forward passes, Advances in Neural Information Processing Systems, 36 (2023), pp. 53038–53075. [38] P.-G. Martinsson and J. Tropp, Randomized numerical linear algebra: foundations & algorithms (2020), arXiv preprint arXiv:2002.01387, (2020). [39] R. Murray, J. Demmel, M. W. Mahoney, N. B. Erichson, M. Melnichenko, O. A. Malik, L. Grigori, P. Luszczek, M. Dereziński, M. E. Lopes, et al., Randomized numerical linear algebra: A perspective on the field with an eye to software, arXiv preprint arXiv:2302.11474, (2023). [40] A. Nissen and G. Kreiss, An optimized perfectly matched layer for the schrödinger equation, Communications in Computational Physics, 9 (2011), pp. 147–179. [41] G. Pang, L. Lu, and G. E. Karniadakis, fpinns: Fractional physics-informed neural networks, SIAM Journal on Scientific Computing, 41 (2019), pp. A2603–A2626. [42] A. S. Parkins and D. F. Walls, The physics of trapped dilute-gas bose–einstein condensates, Physics Reports, 303 (1998), pp. 1–80. [43] L. P. Pitaevskii, Vortex lines in an imperfect bose gas, Sov. Phys. JETP, 13 (1961), pp. 451– 454. [44] A. Rahimi and B. Recht, Random features for large-scale kernel machines, Advances in neural information processing systems, 20 (2007). [45] A. Rahimi and B. Recht, Uniform approximation of functions with random bases, in 2008 46th annual allerton conference on communication, control, and computing, IEEE, 2008, pp. 555–561. [46] M. Raissi, Forward–backward stochastic neural networks: deep learning of high-dimensional partial differential equations, in Peter Carr Gedenkschrift: Research Advances in Mathematical Finance, World Scientific, 2024, pp. 637–655. [47] S. Shao, Y. Wang, and J. Wu, An adaptive hermite spectral method for the boltzmann equation, arXiv preprint arXiv:2509.17981, (2025). [48] Z. Shi, Z. Hu, M. Lin, and K. Kawaguchi, Stochastic taylor derivative estimator: Efficient amortization for arbitrary differential operators, Advances in Neural Information Processing Systems, 37 (2024), pp. 122316–122353. [49] J. Sirignano and K. Spiliopoulos, Dgm: A deep learning algorithm for solving partial differential equations, Journal of computational physics, 375 (2018), pp. 1339–1364. [50] C. M. Stein, Estimation of the mean of a multivariate normal distribution, The annals of Statistics, (1981), pp. 1135–1151. [51] T. Tang, The hermite spectral method for gaussian-type functions, SIAM journal on scientific computing, 14 (1993), pp. 594–606. [52] J. Wang, S. Paesani, R. Santagati, et al., Experimental quantum hamiltonian learning, Nature Physics, 13 (2017), pp. 551–555, https://doi.org/10.1038/nphys4074. [53] Y. Wang, P. Jin, and H. Xie, Tensor neural network and its numerical integration, arXiv preprint arXiv:2207.02754, (2022). [54] Z. Xu and H. Han, Absorbing boundary conditions for nonlinear schrödinger equations, Physical Review E—Statistical, Nonlinear, and Soft Matter Physics, 74 (2006), p. 037704. [55] J. Zhang, Z. Xu, and X. Wu, Unified approach to split absorbing boundary conditions for nonlinear schrödinger equations, Physical Review E—Statistical, Nonlinear, and Soft Matter Physics, 78 (2008), p. 026709.

Record · ID 5968 · SHA-256 5e55c3c876870977
Conceptio Open Knowledge Archive — every document is proof-bundled with source, license, and retrieval metadata.