1 Introduction↩︎

Many physical systems exhibit memory effects or spatial nonlocality, where the state variables depend on the entire history of the system’s evolution and/or on long-range spatial interactions, including viscoelastic polymer relaxation [1], [2], anomalous diffusion in biological membranes [3]–[5], and heat conduction with thermal memory [6], [7]. These responses are modeled by integro-differential equations (IDEs), where a differential operator, as in standard ordinary or partial differential equations (ODEs/PDEs), is coupled with an integral operator. The key element in IDEs is a kernel that weights contributions from past states or spatially distant regions. Often, the kernel is unknown, and identifying it is crucial for solving IDEs efficiently [8], [9], incorporating it into reduced-order models [10], [11], and constructing physically interpretable surrogates. However, kernel identification from data remains a challenging problem.

The inverse problem of identifying memory and nonlocal kernels from observational data has a long history in applied mathematics. Classical approaches typically reformulate the problem as a Volterra integral equation for the unknown kernel and use problem-specific information such as boundary flux measurements [12] or integral overdetermination constraints [13]. The resulting equations are then solved using analytical techniques, including eigenfunction expansions [14], Laplace transforms [15], resolvent estimates [16], and fixed-point arguments [13], or formulated as regularized optimization problems solved using adjoint-based optimization, Tikhonov regularization [17], iterative reconstruction algorithms, deconvolution [18], and Prony-series fitting [19]. Although rigorous uniqueness and stability results have been established under suitable structural assumptions, many classical reconstruction methods rely on problem-specific analytical derivations, specialized observation requirements, or restrictive assumptions on the kernel to derive inversion formulas and guarantee well-posedness. Moreover, these reconstruction procedures are often tailored to particular classes of IDEs and observation settings, limiting their applicability across diverse physical systems. These limitations motivate flexible data-driven approaches capable of identifying unknown kernel structures directly from observations.

Recent advances in scientific machine learning (SciML) have made it possible to learn governing equations, constitutive relations, and unknown operators directly from observational data across a wide range of ODEs, PDEs, and dynamical systems [20]–[26]. The same data-driven approach has increasingly been applied to nonlocal systems to learn solution dynamics, kernel functions, or both. Existing work can be broadly grouped into neural-network-based methods, sparse and symbolic regression methods, and neural differential equations. Neural-network-based approaches, often within physics-informed learning frameworks, parameterize the unknown kernel as a neural network and identify its parameters by enforcing the governing equation during optimization. For example, [27] used an RBF-based physics-informed neural network to recover peridynamic kernels in a nonlocal wave equation, while [28] studied memory-kernel identification in viscoelasticity using a physics-informed neural architecture with a differentiable Prony state-space memory layer. These methods demonstrate that neural representations can reconstruct kernels in nonlocal and history-dependent systems, but the learned kernels are represented implicitly in network parameters or state-space parameterizations rather than recovered as closed-form expressions. To obtain more interpretable representations, sparse and symbolic regression have also been developed. Building on sparse identification of nonlinear dynamics (SINDy), [29] recovered kernels in delay IDEs and renewal equations, and [30] estimated nonlocal interaction kernels using sparse identification with partial inversion. Symbolic-regression approaches, such as the method of [31], combine preprocessing or transform-based steps with symbolic search to identify candidate kernel forms. These approaches can produce compact, interpretable expressions when the target kernel is well represented by the chosen candidate space, but their performance may depend strongly on preprocessing, prescribed libraries, or transform-based assumptions. A third related class of methods is neural ODEs and universal differential equations, in which unknown components of a dynamical system are replaced by trainable neural networks. Extending this idea to nonlocal dynamics, [32] introduced neural integro-differential equations (NIDE) and showed that incorporating integral terms can improve the modeling of nonlocal brain dynamics compared with standard Neural ODEs. Their formulation follows an optimize-then-discretize [21] strategy that requires deriving continuous adjoint equations specific to the IDE for gradient computation. A subsequent extension introduced Neural Integral Equations (NIEs) and Attentional Neural Integral Equations (ANIEs), in which integral operators are approximated via Monte Carlo sampling and attention mechanisms to improve scalability in higher-dimensional settings, with applications to Navier–Stokes equations and brain dynamics [33]. However, these methods learn the solution dynamics of the IDE, rather than recover the memory kernel itself as an explicit object of scientific interest. They also remain tied to case-specific formulations, requiring either IDE-specific adjoint derivations or iterative neural integral solvers.

In this work, we instead target the identification of unknown memory and nonlocal kernels in IDEs directly from observations. Our framework replaces the unknown kernel with a trainable neural representation embedded within a differentiable IDE solver, and optimizes its parameters by minimizing the mismatch between predicted and observed spatiotemporal data. In contrast to adjoint-based methods, we use a discretize-then-optimize strategy [34], [35] where the gradients are computed directly through the numerical solver using automatic differentiation, removing the need for case-specific adjoint derivations. The framework is therefore compatible with standard time integrators and spatial discretizations. Kernel identification from sparse and noisy observations is generally ill-posed [36], so physically meaningful constraints are required to restrict the solution space. In a broad class of memory-dominated systems, including the viscoelastic relaxation, anomalous diffusion, and thermal memory examples noted above, the kernel encodes a fading memory in which the influence of past states decays smoothly with elapsed time. We therefore focus on kernels that are positive, monotonically decreasing, and convex, consistent with standard exponential and power-law relaxation models [4], [19], [37]–[39]. To impose these properties, we consider two constrained Kolmogorov–Arnold Network (KAN) parameterizations. The proposed Monotone–Convex KAN (MC-KAN) is a Bernstein-polynomial-based architecture that enforces positivity, monotonic decrease, and convexity by construction through simple constraints on its polynomial coefficients. The second is a Chebyshev-based KAN (Cheb-KAN) [40], in which the same physical properties are encouraged through soft penalty terms during optimization. After training, the learned neural kernel is converted into an interpretable closed-form expression using symbolic regression [41]. We validate the proposed framework on three benchmark problems spanning temporal and spatial kernel discovery: a Volterra IDE, a one-dimensional viscoelastic wave partial integro-differential equation (PIDE), and a two-dimensional nonlocal reaction-diffusion equation. Across these benchmarks, both KAN parameterizations successfully identify the underlying kernels, while the hard-constrained formulation becomes advantageous for sparse, noisy, and multidimensional inverse problems.

The remainder of this paper is organized as follows. Section 2 presents the methodology, including the problem formulation, memory kernel properties, and the MC-KAN framework. Section 3 presents the numerical results and discussion. Section 4 concludes with key findings and future directions. 8 presents the theoretical guarantee showing that MC-KAN preserves the prescribed kernel constraints, 9 describes the coefficient reparameterization, and 10 provides additional details on the numerical experiments.

2 Methodology↩︎

In this section, we present the proposed framework for kernel discovery, illustrated in Fig. 1. We first formulate the kernel discovery problem and define the observational data. We then introduce the two KAN parameterizations considered in this work, namely the hard-constrained MC-KAN and the soft-constrained Cheb-KAN. Finally, we describe the differentiable forward solver, the training procedure, and the symbolic regression step used to obtain an interpretable closed-form representation of the learned kernel.

2.1 Problem Setup↩︎

Let \(\mathcal{D}=\{(t_k,\mathbf{x}_i,u(t_k,\mathbf{x}_i))\}\), \(k=1,\ldots,N_t\), \(i=1,\ldots,N_x\), denote a dataset of spatiotemporal observations of a scalar field \(u:[0,T]\times\Omega\rightarrow\mathbb{R}\), where \(\Omega\subset\mathbb{R}^n\) is the spatial domain; see the problem setup panel in Fig. 1. We assume that \(u\) satisfies an IDE of the form, \[\mathcal{T}[u](t,\mathbf{x}) = \mathcal{N}[u](t,\mathbf{x}) + \int_0^t \int_{\Omega} K(t-s,\mathbf{x}-\mathbf{y})\, \mathcal{F}[u](s,\mathbf{y}) \, d\mathbf{y}\, ds + S(t,\mathbf{x}), \label{eq:general-IDE-method}\tag{1}\] where \(\mathcal{T}[u](t,\mathbf{x})\) is the temporal evolution operator, \(\mathcal{N}[u](t,\mathbf{x})\) is the local operator, \(S(t,\mathbf{x})\) is an external source term, and \(K(t-s,\mathbf{x}-\mathbf{y})\) is a spatiotemporal convolution kernel that weights the contribution of the field state at past time \(s\) or spatial location \(\mathbf{y}\). The operator \(\mathcal{F}\) is a (possibly differential) operator acting on \(u\) and evaluated at \((s,\mathbf{y})\). Equation 1 is supplemented with boundary and initial conditions, \[\mathcal{B}[u](t,\mathbf{x}) = 0, \qquad \mathbf{x} \in \partial\Omega, \quad t \in [0,T], \label{eq:bc-method}\tag{2}\] \[u(0,\mathbf{x}) = u_0(\mathbf{x}), \qquad \partial_t u(0,\mathbf{x}) = v_0(\mathbf{x}), \qquad \mathbf{x} \in \Omega, \label{eq:ic-method}\tag{3}\] where \(\mathcal{B}\) is a boundary operator, \(u_0(\mathbf{x})\) and \(v_0(\mathbf{x})\) are the prescribed initial field values and initial rates of change, respectively. The second condition is imposed only when the temporal operator \(\mathcal{T}\) is second order in time. The functional form of the kernel \(K(t-s,\mathbf{x}-\mathbf{y})\) encodes how strongly the current state is influenced by past states and/or spatial locations, and is assumed to be unknown. We seek to identify \(K(t-s,\mathbf{x}-\mathbf{y})\) directly from the spatiotemporal dataset \(\mathcal{D}\).

When \(\mathcal{D}\) is sparse or noisy, this identification problem is ill-posed in the sense of Hadamard [36], and prior information is required to obtain a stable, physically meaningful solution. We therefore restrict attention to a broad class of dissipative memory processes whose kernels are positive, monotonically decreasing, and convex along each lag dimension, rather than jointly convex over the full lag domain. Writing \(\boldsymbol{\tau} = (\tau_t,\boldsymbol{\tau}_X)\), with \(\tau_t=t-s, \boldsymbol{\tau}_X=|\mathbf{x}-\mathbf{y}|\) (componentwise absolute value), and \(\tau_j\) for its \(j\)-th component, we require, \[K(\boldsymbol{\tau}) > 0, \qquad \frac{\partial K}{\partial \tau_j} \le 0, \qquad \frac{\partial^2 K}{\partial \tau_j^2} \ge 0, \qquad j = 1,\ldots,n+1 , \label{eq:kernel-properties}\tag{4}\] on the lag domain, taken with \(\tau_t^{\min} > 0\) so that the weakly singular power-law and Mittag–Leffler kernels remain permissible while the derivatives in Eq. 4 stay well defined.

2.2 Neural Kernel Parameterization↩︎

To identify a kernel satisfying Eq. 4 , we replace the unknown kernel with a neural network \(K_{\boldsymbol{\theta}}\) (see neural kernel in Fig. 1), where \(\boldsymbol{\theta}\) denotes the trainable parameters, \[K_{\boldsymbol{\theta}} : \mathbb{R}_{\ge 0}\times\mathbb{R}^n \;\to\; \mathbb{R}, \qquad \boldsymbol{\tau} \mapsto K_{\boldsymbol{\theta}}(\boldsymbol{\tau}), \label{eq:kernel-nn}\tag{5}\] that maps a space–time lag \(\boldsymbol{\tau}\) to a scalar kernel value. In this work, we consider two KAN-based parameterizations. The first is a Chebyshev-based KAN (Cheb-KAN) [40], [42], in which positivity, monotonic decrease, and convexity are encouraged through soft penalty terms. The second is the proposed Monotone–Convex Kolmogorov–Arnold Network (MC-KAN), in which the same physical properties are enforced by construction. Our choice of KANs is motivated by their decomposition into learnable univariate edge functions [43]–[46], which makes shape constraints easier to impose. Standard multilayer perceptrons (MLPs) [47], on the other hand, often require separate specialized architectures such as input-convex [48] or monotone neural networks [49]. Recent work has shown that shape constraints can be imposed in KANs through conditions on spline coefficients, including input-convex KANs for polyconvex hyperelastic constitutive modeling [50] and convex approximation and optimal transport [51]. Building on this idea, MC-KAN adapts constrained KANs to the memory-kernel setting by using a Bernstein basis [52] whose coefficient constraints guarantee positivity, monotone decrease, and convexity along each lag dimension.

2.2.1 Soft-constrained Chebyshev KAN↩︎

For a Cheb-KAN architecture \([d,h_1,\ldots,h_{L-1},h_L]\), where \(d=n+1\) is the number of kernel inputs (one temporal lag and \(n\) spatial lags) and \(h_l\) is the number of nodes in layer \(l\), with \(h_0=d\) and \(h_L=1\), the neural kernel is constructed as a composition of univariate transformations [42], [53], \[K_{\text{Cheb-KAN}}(\boldsymbol{\tau}) = \left( \boldsymbol{\Phi}_L \circ \tanh \circ \boldsymbol{\Phi}_{L-1} \circ \cdots \circ \tanh \circ \boldsymbol{\Phi}_1 \right)\!\bigl(\boldsymbol{\chi}(\boldsymbol{\tau})\bigr), \label{eq:cheb-kernel}\tag{6}\] where \(\boldsymbol{\chi}(\boldsymbol{\tau})\) denotes the normalized lag input and \(\tanh\) acts componentwise. Since Chebyshev polynomials are defined on \([-1,1]\), each lag component \(\tau_i\) is linearly normalized as \[\chi_i(\boldsymbol{\tau}) = 2\,\frac{\tau_i-\tau_i^{\min}}{\tau_i^{\max}-\tau_i^{\min}} - 1, \qquad i=1,\ldots,d, \label{eq:cheb-normalization}\tag{7}\] so that \(\chi_i\in[-1,1]\). To ensure numerical stability during training and to respect the domain of the Chebyshev polynomials, a \(\tanh\) activation is applied after each hidden layer, constraining all intermediate values to \([-1,1]\). Each layer map \(\boldsymbol{\Phi}_l:\mathbb{R}^{h_{l-1}}\rightarrow\mathbb{R}^{h_l}\) is composed of learnable univariate edge functions. For \(\mathbf{z}\in\mathbb{R}^{h_{l-1}}\), \[\left(\boldsymbol{\Phi}_l(\mathbf{z})\right)_p = \sum_{q=1}^{h_{l-1}} \phi^{(l)}_{p,q}(z_q), \qquad p=1,\ldots,h_l, \label{eq:cheb-layer}\tag{8}\] and each edge function is parameterized in a Chebyshev polynomial basis [40], \[\phi^{(l)}_{p,q}(z) = \sum_{k=0}^{M} a^{(l)}_{p,q,k}\,T_k(z), \qquad z\in[-1,1], \label{eq:cheb-edge}\tag{9}\] where \(z\) denotes the scalar input to the edge, \(M\) is the maximum polynomial degree, \(T_k\) is the Chebyshev polynomial of the first kind of degree \(k\), and \(a^{(l)}_{p,q,k}\in\mathbb{R}\) are trainable coefficients, so that \(\boldsymbol{\theta}_{\text{Cheb-KAN}}=\{a^{(l)}_{p,q,k}\}\). The Chebyshev polynomials are defined recursively as [54], [55], \[T_0(z)=1, \qquad T_1(z)=z, \qquad T_{k+1}(z)=2z\,T_k(z)-T_{k-1}(z), \qquad k\ge1. \label{eq:cheb-recursion}\tag{10}\]

2.2.2 Hard-constrained Monotone–Convex KAN↩︎

For the MC-KAN architecture \([d,h_1,\ldots,h_{L-1},h_{L}]\), where \(d=n+1\) is the number of kernel inputs (one temporal lag and \(n\) spatial lags) and \(h_l\) is the number of nodes in layer \(l\), the neural kernel is written as [53], \[K_{\text{MC-KAN}}(\boldsymbol{\tau}) = \left( \boldsymbol{\Phi}_L \circ \boldsymbol{\Phi}_{L-1} \circ \cdots \circ \boldsymbol{\Phi}_1 \right)\!\bigl(\boldsymbol{\zeta}(\boldsymbol{\tau})\bigr), \label{eq:mc-kernel}\tag{11}\] where \(\boldsymbol{\zeta}(\boldsymbol{\tau})\) denotes the normalized lag input. Since Bernstein polynomials are defined on \([0,1]\), each lag component \(\tau_i\) is first normalized as, \[\zeta_i(\boldsymbol{\tau}) = \left( \frac{\tau_i-\tau_i^{\min}}{\tau_i^{\max}-\tau_i^{\min}} \right)^{\alpha_i}, \qquad i=1,\ldots,d. \label{eq:bernstein-normalization}\tag{12}\] Here, \(\alpha_i\in(0,1]\) is a scaling exponent that adjusts the resolution of the normalized coordinate along the \(i\)-th lag dimension. Depending on the application, \(\alpha_i\) may be specified a priori or treated as a trainable parameter learned jointly with the network weights. Linear normalization corresponds to \(\alpha_i=1\), while \(\alpha_i<1\) gives a power-law normalization that allocates more resolution near small lags. Restricting \(\alpha_i\le1\) makes each normalized coordinate \(\zeta_i\) monotone increasing and concave in \(\tau_i\), so the normalization preserves positivity, monotone decrease, and convexity under composition. When \(\alpha_i\) is learnable it is constrained to be \(\alpha_i\in(0,1)\). A formal argument is provided in 8.

Let \(\mathbf{z}^{(l)}\in\mathbb{R}^{h_l}\) denote the output of layer \(l\), with \(\mathbf{z}^{(0)}=\boldsymbol{\zeta}(\boldsymbol{\tau})\in[0,1]^d\) and \(h_0=d\), \(h_L=1\). Each layer map \(\boldsymbol{\Phi}_l:[0,1]^{h_{l-1}}\rightarrow(0,1)^{h_l}\) is composed of learnable univariate edge functions, aggregated as an average [53], \[z_p^{(l)} = \left(\boldsymbol{\Phi}_l\!\left(\mathbf{z}^{(l-1)}\right)\right)_p = \frac{1}{h_{l-1}} \sum_{q=1}^{h_{l-1}} \phi_{p,q}^{(l)}\!\left(z_q^{(l-1)}\right), \qquad p=1,\ldots,h_l, \label{eq:mc-kan-layer}\tag{13}\] for \(l=1,\ldots,L\). Here, \(q\) indexes the input node from layer \(l-1\), and \(\phi_{p,q}^{(l)}:[0,1]\rightarrow(0,1)\) is a learnable univariate function on the edge from node \(q\) in layer \(l-1\) to node \(p\) in layer \(l\). The coefficient parameterization introduced below constrains each edge function to map into \((0,1)\). Since each layer computes the average of its incoming edge-function values, every layer output, including the final kernel value, remains in \((0,1)\). The factor \(1/h_{l-1}\) is therefore a deliberate normalization for bound control, in contrast to the unnormalized sum used in standard KANs. This boundedness is an architectural consequence rather than a physical requirement. If kernels with arbitrary amplitude are required, the final network output may be multiplied by a positive learnable scaling parameter without affecting the positivity, monotonicity, or convexity guarantees.

The key design choice in MC-KAN is to parameterize each edge function as a Bernstein polynomial of degree \(N\), \[\phi_{p,q}^{(l)}(z) = \sum_{k=0}^{N} c_k^{(l,p,q)}\,B_{k,N}(z), \qquad z\in[0,1], \label{eq:edge-bernstein}\tag{14}\] where \(z=z_q^{(l-1)}\) is the scalar input to the edge, and the Bernstein basis functions are \[B_{k,N}(z) = \binom{N}{k} z^{k}(1-z)^{N-k}, \qquad k=0,\ldots,N, \quad z\in[0,1]. \label{eq:bernstein-basis}\tag{15}\] Here, \(\binom{N}{k}=N!/[k!(N-k)!]\) is the binomial coefficient. Because the Bernstein basis satisfies \(B_{k,N}(z)\ge0\) for all \(z\in[0,1]\) [52], shape constraints can be imposed through conditions on the Bernstein coefficients. For a Bernstein polynomial, \[\phi(z) = \sum_{k=0}^{N} c_k\,B_{k,N}(z),\] the first and second derivatives are, \[\frac{d\phi}{dz} = N \sum_{k=0}^{N-1} (c_{k+1}-c_k)\,B_{k,N-1}(z), \label{eq:derivative}\tag{16}\] \[\frac{d^2\phi}{dz^2} = N(N-1) \sum_{k=0}^{N-2} (c_{k+2}-2c_{k+1}+c_k)\,B_{k,N-2}(z). \label{eq:second-derivative}\tag{17}\]

Since \(B_{k,N-1}(z)\ge0\) and \(B_{k,N-2}(z)\ge0\) on \([0,1]\), sufficient conditions for monotonicity and convexity are obtained by constraining the first and second finite differences of the coefficients. Thus, a Bernstein edge function satisfies the desired boundedness, positivity, monotone decrease, and convexity properties if its coefficients satisfy the conditions listed in Table 1. The coefficient conditions in Table 1 are imposed using a differentiable constrained reparameterization of the Bernstein coefficients. In the first layer, the coefficients are ordered as \(c_0\ge c_1\ge\cdots\ge c_N\) to obtain decreasing convex edge functions. In subsequent layers, the coefficient sequence is reversed, giving \(c_0\le c_1\le\cdots\le c_N\) and hence increasing convex edge functions. In all layers, the reparameterization keeps the coefficients in \((0,1)\), so each edge function maps \([0,1]\) into \((0,1)\). This construction allows the prescribed kernel constraints to propagate through the full network depth. The formal proof and coefficient reparameterization are provided in 8 and 9.

2.3 Training Procedures↩︎

Given a neural kernel \(K_{\boldsymbol{\theta}}\), we embed it in a differentiable numerical solver for Eq. 1 ; see the forward solver in Fig. 1. The integral term is discretized using a quadrature or convolution scheme appropriate to the governing equation. Depending on the structure of the kernel and the problem geometry, the resulting discrete integral operator is evaluated either directly or using Fast Fourier Transforms (FFTs). The temporal evolution operator is then advanced with a suitable time-integration scheme. The specific spatial discretizations, quadrature rules, convolution evaluations, and time integrators are chosen according to the structure of each benchmark problem and reported where necessary. During the forward solve, the neural kernel \(K_{\boldsymbol{\theta}}\) is evaluated at the required lag points to assemble the discrete integral operator and propagate the solution in time. Let \(\hat{u}(t_k,\mathbf{x}_i;K_{\boldsymbol{\theta}})\) denote the numerical solution obtained with the current neural kernel, evaluated at the observation points in \(\mathcal{D}\). The data loss is defined as, \[\mathcal{L}_{\mathrm{data}}(\boldsymbol{\theta}) = \frac{1}{|\mathcal{D}|} \sum_{(t_k,\mathbf{x}_i)\in\mathcal{D}} \left( \hat{u}(t_k,\mathbf{x}_i;K_{\boldsymbol{\theta}}) - u(t_k,\mathbf{x}_i) \right)^2 . \label{eq:loss-method}\tag{18}\]

For MC-KAN, the kernel constraints are enforced by construction through the Bernstein coefficient parameterization, so the training objective is simply \(\mathcal{L}_{\mathrm{data}}\). For the Cheb-KAN baseline, the same data loss is combined with soft constraint penalties evaluated at a finite set of sampled lag points. Let \(\{\boldsymbol{\tau}_j\}_{j=1}^{N_g}\) denote points in the lag domain used to evaluate constraint violations. Depending on the experiment, these points are selected either from a uniform grid or sampled randomly from the domain. We define, \[\mathcal{L}_{\mathrm{pos}} = \frac{1}{N_g} \sum_{j=1}^{N_g} \left[ \max \bigl( 0, -K_{\boldsymbol{\theta}}(\boldsymbol{\tau}_j) \bigr) \right]^2, \label{eq:cheb-pos-loss}\tag{19}\] \[\mathcal{L}_{\mathrm{mono}} = \frac{1}{N_g d} \sum_{j=1}^{N_g} \sum_{i=1}^{d} \left[ \max \left( 0, \frac{\partial K_{\boldsymbol{\theta}}}{\partial \tau_i} (\boldsymbol{\tau}_j) \right) \right]^2, \label{eq:cheb-mono-loss}\tag{20}\] \[\mathcal{L}_{\mathrm{conv}} = \frac{1}{N_g d} \sum_{j=1}^{N_g} \sum_{i=1}^{d} \left[ \max \left( 0, - \frac{\partial^2 K_{\boldsymbol{\theta}}}{\partial \tau_i^2} (\boldsymbol{\tau}_j) \right) \right]^2. \label{eq:cheb-conv-loss}\tag{21}\] Here, \(\mathcal{L}_{\mathrm{pos}}\) penalizes negative kernel values, \(\mathcal{L}_{\mathrm{mono}}\) penalizes positive first derivatives, and \(\mathcal{L}_{\mathrm{conv}}\) penalizes negative second derivatives along each lag dimension. The full Cheb-KAN training loss is then, \[\mathcal{L}_{\mathrm{Cheb-KAN}} = \mathcal{L}_{\mathrm{data}} + \lambda \left( w_{\mathrm{pos}} \mathcal{L}_{\mathrm{pos}} + w_{\mathrm{mono}} \mathcal{L}_{\mathrm{mono}} + w_{\mathrm{conv}} \mathcal{L}_{\mathrm{conv}} \right), \label{eq:cheb-total-loss}\tag{22}\] where \(\lambda\) controls the overall strength of the constraint regularization, and \(w_{\mathrm{pos}}\), \(w_{\mathrm{mono}}\), and \(w_{\mathrm{conv}}\) determine the relative weights of the positivity, monotonicity, and convexity penalties, respectively. In all experiments, these three weights are set equal. The value of \(\lambda\) is selected separately for each benchmark by a short calibration procedure: it is chosen so that, at the start of training, the total constraint penalty is of the same order of magnitude as the data loss. During training, the constraint loss is monitored to verify that it remains active and does not become negligible relative to the data loss. Gradients of the loss with respect to the trainable parameters \(\boldsymbol{\theta}\) are computed by automatic differentiation [56], [57] through the discrete solver. Thus, the method follows a discretize-then-optimize strategy. Unless otherwise specified, we train the parameters using Adam [58] followed by L-BFGS with a strong Wolfe line search [59], [60].

Once training is complete, the neural kernel \(K_{\boldsymbol{\theta}}\) provides a pointwise approximation of the underlying kernel over its domain. To extract a closed-form expression, we sample the learned kernel and apply symbolic regression via PySR [41] (see kernel recovery in Fig. 1), which searches over algebraic expressions to fit the sampled values. Alternative sparse-regression approaches over a predefined function library, such as SINDy [20], could also be applied at this stage. The recovered expression yields a physically interpretable kernel that can be compared directly with known analytical forms when available.

3 Results and Discussion↩︎

In this section, we evaluate the proposed framework on three kernel identification problems involving representative memory and nonlocal kernels commonly encountered in applications of IDEs. For each problem, the objective is to recover the unknown kernel from observations of the system response using only the governing equation and the available data. To compare hard and soft constraint enforcement, both MC-KAN and Cheb-KAN use the same number of hidden layers and approximately the same number of trainable parameters. All experiments are repeated over five independent random seeds to account for sensitivity to initialization. Unless otherwise stated, the reported results correspond to the mean performance across all runs, with standard deviations reported where relevant. To quantify predictive accuracy, we use the relative \(L_2\) error, \[\mathcal{E}(\nu) = \frac{\|\hat{\nu}-\nu\|_2}{\|\nu\|_2}, \label{eq:rel95error}\tag{23}\] where \(\nu\) denotes the true quantity and \(\hat{\nu}\) denotes its learned approximation. This metric is applied to both the recovered kernel \(K\) and the predicted solution \(u\) throughout all experiments. When the training observations are noisy, the ground truth for \(\mathcal{E}(u)\) is the reference solution before the synthetic noise is added, evaluated at the observation points in \(\mathcal{D}\).

3.1 Experiment I: 1D Volterra IDE↩︎

As a first experiment, we consider a scalar Volterra IDE describing the temporal evolution of a state variable \(u(t)\) in a system with memory. Such equations arise in models of viscoelasticity, population dynamics with hereditary effects, and anomalous relaxation processes [61], [62]. The governing equation is,

\[\frac{du(t)}{dt} + u(t) = \kappa \int_0^t K(t-s)\, u(s)\,ds, \qquad t\in(0,T], \label{eq:volterra}\tag{24}\] subject to the initial condition \(u(0)=1\), and \(\kappa>0\) controls the memory strength. Throughout this experiment, we set \(\kappa=1.5\) and \(T=5\,\mathrm{s}\). To generate reference data, we prescribe the memory kernel, \[K(t-s) = e^{-(t-s)}. \label{eq:exp95kernel}\tag{25}\]

The reference data are generated by solving Eq. 24 using a forward Euler time-discretization scheme together with a trapezoidal approximation of the memory integral. A sufficiently small time step, \(\Delta t=5\times10^{-3} \mathrm{s}\), is used for time integration, for which the numerical solution was observed to be stable and converged; see Fig. 9 in Appendix. The resulting solution is sampled at \(N_t=1001\) uniformly spaced time points.

In this experiment, we investigate the effect of measurement noise on kernel discovery. We take \(501\) uniformly spaced samples from the reference solution and add independent Gaussian noise, \[u^{\mathrm{obs}}(t_i) = u(t_i) + \sigma \,\varepsilon_i, \qquad \varepsilon_i \sim \mathcal{N}(0,1),\] where \(\sigma\) is the noise standard deviation. The resulting training dataset is \(\mathcal{D} = \bigl\{\,\bigl(t_i,\,u^{\mathrm{obs}}(t_i)\bigr)\,\bigr\}_{i=1}^{501}.\)

The unknown kernel is represented by either MC-KAN or Cheb-KAN. In Cheb-KAN, positivity, monotonic decrease, and convexity are encouraged through equally weighted soft penalty terms with \(\lambda=10^{-1}\). For both models, the inputs are linearly normalized. Following the original Cheb-KAN formulation [40], we use cubic Chebyshev edge functions (\(M=3\)). For MC-KAN, the Bernstein degree is fixed at \(N=6\), which was selected empirically to provide sufficient approximation accuracy while keeping the number of trainable parameters moderate. During training, the differentiable solver uses the same time step and numerical discretization as the reference solver, advancing on the full temporal grid (\(N_t = 1001\)). At each optimization step, we compute the solution on this grid and evaluate the loss only at the \(501\) observation times in \(\mathcal{D}\). For evaluation, the learned kernel is sampled on the same 501-point temporal grid used for the training observations, and the kernel reconstruction error is computed on this common grid. Complete numerical, architectural, and optimization settings are summarized in Table 7. For each method, we evaluate three architectures of increasing depth (shallow, medium, and deep), with widths chosen to match parameter counts across the two methods, so that the comparison is not tied to a single network choice.

Table 2 reports the solution and kernel reconstruction errors for all six network configurations at each noise level. Both MC-KAN and Cheb-KAN successfully identify the unknown kernel across all architectures and noise levels. Even at \(\sigma=0.15\), both methods recover the kernel with relative errors below \(13\%\). The solution reconstruction errors are nearly identical for the two methods, increasing from almost zero in the noise-free case to below \(2\%\) at the highest noise level. Among the six network configurations, however, the lowest kernel reconstruction error at every noise level is achieved by MC-KAN. For clarity, Fig. 2 compares the two methods using the common architecture depth that produces the lowest mean kernel reconstruction error among all six configurations at each noise level. In other words, the shallow, medium, or deep architecture associated with the lowest kernel error is identified from Table 2, and both MC-KAN and Cheb-KAN are then evaluated using that same depth. Relative to Cheb-KAN, MC-KAN reduces \(\mathcal{E}(K)\) by approximately \(54\%\) in the noise-free case (\(0.85\%\) versus \(1.83\%\)), \(32\%\) at \(\sigma=0.05\) (\(4.47\%\) versus \(6.62\%\)), \(32\%\) at \(\sigma=0.10\) (\(6.83\%\) versus \(10.05\%\)), and \(13\%\) at \(\sigma=0.15\) (\(8.78\%\) versus \(10.09\%\)). Figures 2 (b) and (c) show that both methods accurately reconstruct the solution at low and high noise levels.

5pt

Figure 3 presents the recovered kernels as the pointwise mean over the five random seeds, with the shaded band indicating \(\pm\) one standard deviation. The figure corresponds to the same architecture selection used in Fig. 2. At all noise levels, both methods closely follow the reference exponential kernel. The hard-constrained MC-KAN remains positive, monotonically decreasing, and convex by construction. In the pointwise kernel bands, Cheb-KAN shows a wider spread near the origin and through the middle of the domain, particularly at higher noise levels. We note that the seed-to-seed standard deviation of the error \(\mathcal{E}(K)\) in Table 2 is comparable for the two methods in this experiment, so the reduced variability of the hard-constrained parameterization is visible in the pointwise reconstructions rather than in the \(\mathcal{E}(K)\) error. The corresponding predicted solution trajectories and optimization loss histories are provided in Appendix (Figures 10–11).

To obtain an interpretable closed-form representation, we apply PySR to the mean recovered kernel of each configuration, restricting the operator set to primitives expected for relaxation kernels. The resulting expressions are annotated in Fig. 3. In all cases, both MC-KAN and Cheb-KAN recover the exponential functional form \(K(\tau)=e^{-r\tau}\), with the learned decay rate remaining close to the true value. At low noise, the recovered rate is approximately \(0.99\) for both methods, decreasing to roughly \(0.95\)–\(0.96\) at \(\sigma=0.15\). Thus, although MC-KAN provides more accurate numerical kernel reconstructions, both KAN parameterizations successfully identify the correct analytical constitutive law after symbolic regression.

3.2 Experiment II: 1D Viscoelastic Wave PIDE↩︎

As a second experiment, we consider the free vibration of a viscoelastic rod, a setting representative of polymer rods and viscoelastic structural elements undergoing free oscillation [63]. We take a homogeneous rod of length \(L\), constant mass density \(\rho\), and uniform cross-section, fixed at both ends. Neglecting body forces, the balance of linear momentum along the rod axis gives, \[\frac{\partial}{\partial x}\sigma(x,t) = \rho \frac{\partial^2}{\partial t^2}u(x,t), \label{eq:motion}\tag{26}\] where \(u(x,t)\) is the displacement field and \(\sigma(x,t)\) the axial stress. Deformation is measured by the axial strain, \[\epsilon(x,t) = \frac{\partial}{\partial x}u(x,t), \label{eq:kinematic}\tag{27}\] and the stress is related to the strain history through the Boltzmann superposition principle, which encodes the fading memory of the viscoelastic material, \[\sigma(x,t) = \epsilon(x,0^+)\,K(t) + \int_0^t K(t-s)\, \frac{\partial}{\partial s} \epsilon(x,s)\,ds , \label{eq:boltzmann}\tag{28}\] where \(K\) is the relaxation kernel and the first term accounts for the instantaneous elastic response at \(t=0^+\). Substituting the constitutive law Eq. 28 and the kinematic relation Eq. 27 into the momentum balance Eq. 26 , with the term \(\epsilon(x,0^+)K(t)\) omitted because of the zero-displacement initial condition, yields the following PIDE for the displacement, \[\rho \frac{\partial^2 u}{\partial t^2}(x,t) = \int_0^t K(t-s)\, \frac{\partial^3 u}{\partial x^2 \partial s}(x,s)\,ds , \qquad (x,t)\in(0,L)\times(0,T], \label{eq:pide}\tag{29}\] in which the memory kernel \(K(t-s)\) couples the current acceleration to the entire strain-rate history. In this experiment, we set the rod length to \(L=1~\mathrm{m}\) and the final time to \(T=5~\mathrm{s}\). The rod is fixed at both ends, \[u(0,t)=0, \qquad u(L,t)=0, \qquad t>0, \label{eq:rod-bc}\tag{30}\] and starts from zero displacement with an initial velocity prescribed in the fundamental mode, \[u(x,0)=0, \qquad \frac{\partial u}{\partial t}(x,0) = \sin\!\left(\frac{\pi x}{L}\right). \label{eq:rod-ic}\tag{31}\]

To generate reference data, we prescribe the relaxation modulus using the Kohlrausch–Williams–Watts (KWW) stretched exponential kernel,

\[K(t-s) = \exp\!\left[ -\left( \frac{t-s}{\tau_0} \right)^{\beta} \right], \label{eq:kww}\tag{32}\] where \(\tau_0>0\) is the characteristic relaxation time and \(0<\beta\le1\) is the stretching exponent. In this experiment we set \(\tau_0=1.5\,\mathrm{s}\) and \(\beta=0.5\). The KWW kernel is widely used to model relaxation in glassy polymers, amorphous solids, and other disordered materials, where memory decay cannot be captured by a single exponential time scale [4], [64], [65]. The exponent \(\beta\) sets the degree of non-exponential decay: \(\beta=1\) recovers the classical Maxwell exponential kernel, while \(\beta<1\) corresponds to a broad distribution of relaxation times characteristic of heterogeneous materials. The choice \(\beta=0.5\) therefore represents a more challenging identification problem than the exponential kernel of Experiment I.

The reference data are generated by solving Eq. 29 using an implicit Newmark–\(\beta\) time-integration scheme together with central finite differences in space and a trapezoidal approximation of the memory integral. The average-acceleration Newmark–\(\beta\) scheme (\(\gamma=1/2\), \(\beta_N=1/4\)) is used for time integration. This scheme is classically unconditionally stable for linear elastodynamic systems, and the chosen spatial and temporal resolutions provide adequate accuracy for the present viscoelastic problem (see Fig. 12 in Appendix). The resulting displacement field is computed on a grid of \(1001\) time steps and \(200\) spatial points.

In this experiment, we aim to investigate the effect of data sparsity on kernel discovery. To do so, we construct training datasets by uniformly subsampling both the spatial and temporal dimensions. In all cases the spatial observations are restricted to \(50\) uniformly spaced locations (\(25\%\) of the spatial grid), while the number of temporal snapshots is varied according to \(N_s \in \{101,\,51,\,21,\,11\}\). The resulting dataset is \(\mathcal{D} = \left\{ \bigl( t_k,\, \mathbf{x}_i,\, u(t_k,\mathbf{x}_i) \bigr) \right\},\) where \(\{\mathbf{x}_i\}_{i=1}^{50}\) denotes the retained spatial locations and \(\{t_k\}_{k=1}^{N_s}\) denotes the retained temporal snapshots. The final training datasets contain \(10.1\%\), \(5.1\%\), \(2.1\%\), and \(1.1\%\) of the \(1001\) reference time steps, corresponding to roughly \(2.5\%\), \(1.3\%\), \(0.5\%\), and \(0.3\%\) of the full spatiotemporal grid. No measurement noise is added in this experiment to isolate the effect of data sparsity from observational uncertainty.

As in Experiment I, the unknown kernel is represented by either MC-KAN or Cheb-KAN. In Cheb-KAN, positivity, monotonic decrease, and convexity are encouraged through equally weighted soft penalty terms with \(\lambda=3\times10^{-1}\). The two models use depth-matched architectures with approximately the same number of trainable parameters, and all inputs are linearly normalized. During training, the differentiable solver uses the same spatial grid, time step, and numerical discretization as the reference solver. At each optimization step, the predicted displacement field is computed on the full spatiotemporal grid, and the loss is evaluated only at the observation points in \(\mathcal{D}\). Inside the solver, the relaxation kernel is normalized so that \(K(0^+)=1\), enforced by dividing the network output by its value at \(\tau_{\min}=10^{-3}\). Fixing the instantaneous modulus keeps the implicit term of the memory operator independent of the kernel, so the effective system matrix is assembled and factored once and reused throughout the forward solve. This normalization is also physically motivated: \(K(0^+)\) is the instantaneous glassy modulus, and since only the shape of the relaxation kernel is identifiable from displacement data, fixing it to unity removes the amplitude degeneracy without loss of generality. For evaluation, the learned kernel is sampled at 500 uniformly spaced temporal points, and the kernel reconstruction error is computed with respect to the reference KWW kernel on the same grid. Complete numerical, architectural, and optimization settings are summarized in Table 8.

5pt

5pt

Our results show that sparsity has little effect on kernel identification in this experiment. For each method the kernel error is nearly constant as the number of snapshots is reduced from \(101\) to \(11\) (Table 3), and the solution errors remain near \(0.4\%\) throughout. The recovered KWW parameters are likewise stable, staying within \(1\%\) of the true values at every sparsity level, as demonstrated in Table 4. The stretched-exponential kernel thus remains identifiable even when the temporal data barely resolve the dynamics. The corresponding sparse solution reconstruction errors and optimization loss histories are provided in Appendix (Figs. 13 and14).

Although both methods are insensitive to temporal sparsity, MC-KAN produces lower kernel reconstruction errors than Cheb-KAN at every sparsity level, approximately \(1.2\%\) compared with \(2.0\%\), corresponding to a reduction of roughly \(40\%\). In addition, the variability across random initializations is consistently smaller for MC-KAN. These trends are reflected in Fig. 4, shown as the pointwise mean over seeds with a \(\pm\) one standard deviation band. The MC-KAN reconstructions remain positive, monotone, and convex by construction, and they track the true kernel with a narrow band across all sparsity levels. The Cheb-KAN reconstructions also visually fit the true kernel well, but slight traces of oscillatory behavior are visible between \(\tau=0\) and \(\tau=1\).

Despite the differences in reconstruction accuracy, symbolic regression discovers the same constitutive law for both KAN parameterizations. As shown in Fig. 4, the recovered expression of the form \(0.443^{\,\tau^{\beta}}\) can be rewritten as \(e^{-(\tau/\tau_0)^{\beta}}\), and recalculating its parameters gives values of \(\tau_0\) and \(\beta\) reported in Table 4, in close agreement with the true \(\tau_0=1.5\) and \(\beta=0.5\). Across all sparsity levels, the recovered \(\tau_0\) and \(\beta\) deviate from their true values by less than \(1\%\) for both methods.

3.3 Experiment III: 2D Nonlocal Reaction-Diffusion Equation↩︎

As a third experiment, we consider a two-dimensional nonlocal reaction-diffusion equation. Unlike Experiments I and II, which focus on temporal memory kernels, this problem involves the discovery of a spatial interaction kernel governing nonlocal interactions between points in the domain. This model arises in phase separation, pattern formation, and collective dynamics, where the evolution at a given location depends on the surrounding state through a spatial convolution operator [66], [67]. We study an attractive, distance-decaying interaction [68], [69], for which the shape constraints are exact properties of the kernel. The experiment tests whether enforcing them by construction is advantageous when the true kernel belongs to that class.

The governing equation is,

\[\frac{\partial u}{\partial t}(\mathbf{x},t) = D\nabla^2 u(\mathbf{x},t) + \alpha \int_{\Omega} K(x-\xi,y-\eta) \bigl( u(\xi,\eta,t)-u(x,y,t) \bigr) \,d\xi\,d\eta + u(\mathbf{x},t) - u^3(\mathbf{x},t), \label{eq:nonlocal-ac}\tag{33}\] for \((\mathbf{x},t)\in\Omega\times(0,T]\), where \(\Omega\subset\mathbb{R}^2\), \(\mathbf{x}=(x,y)\), \(u(\mathbf{x},t)\) is the state variable, \(D>0\) is the diffusion coefficient, \(\alpha>0\) controls the strength of the nonlocal interactions, and \(K(x-\xi,y-\eta)\) is the unknown spatial interaction kernel. Throughout this experiment, we set \(D=0.001\), \(\alpha=0.12\), \(L_x=L_y=10\), and the final time to \(T=6~\mathrm{s}\). Periodic boundary conditions are imposed in both spatial directions, and the initial condition consists of a smooth random perturbation of the homogeneous state (10.3.1). To generate reference data, we prescribe the interaction kernel, \[K(x-\xi,y-\eta) = \exp\!\Bigg( - \left| \frac{x-\xi}{\tau_x} \right|^{\beta_x} - \left| \frac{y-\eta}{\tau_y} \right|^{\beta_y} - \delta \left| \frac{x-\xi}{\tau_x} \right|^{\gamma_x} \left| \frac{y-\eta}{\tau_y} \right|^{\gamma_y} \Bigg), \label{eq:spatial95kernel}\tag{34}\] where \((x,y)\) denotes the target location and \((\xi,\eta)\) denotes the source location contributing to the nonlocal interaction. Since the problem is on a periodic domain, the kernel is represented numerically as a function of the periodic displacement (\(x-\xi,y-\eta\)) sampled on the periodic displacement grid. The resulting nonlocal convolution operator is evaluated as a circular convolution using the FFT. The parameters \(\tau_x\) and \(\tau_y\) determine the characteristic interaction lengths in the horizontal and vertical directions, respectively, while \(\beta_x\) and \(\beta_y\) control the decay rate of the kernel. The mixed term weighted by \(\delta\) introduces anisotropic coupling between the two spatial directions. Throughout this experiment, we use, \[\tau_x=3.0, \qquad \tau_y=1.0, \qquad \beta_x=0.4, \qquad \beta_y=0.8,\] \[\gamma_x=0.5, \qquad \gamma_y=0.5, \qquad \delta=0.5.\]

We selected these parameters deliberately to introduce anisotropy through differing interaction lengths (\(\tau_x \neq \tau_y\)) and decay rates (\(\beta_x \neq \beta_y\)). The resulting kernel is therefore substantially more challenging to identify than an isotropic one, providing a more demanding benchmark for multidimensional kernel recovery.

The reference data are generated by solving Eq. 33 on a periodic domain using a Fourier spectral discretization in space together with a semi-implicit time-stepping scheme (IMEX). The nonlocal interaction integral is evaluated as a circular convolution on the periodic domain using the FFT, reducing the computational complexity of the interaction term (see Fig. 15). The solution is computed on a \(256\times256\) spatial grid at \(121\) time points. In this experiment, we aim to mimic realistic measurements by introducing both sparsity and additive noise. The full-resolution solution is subsampled in space and time. Observations are restricted to a \(32\times32\) spatial grid and taken at the \(7\) time points, \(t =\{0,\,1,\,2,\,3,\,4,\,5,\,6\},\) giving the dataset \(\mathcal{D} = \left\{ \bigl( t_k,\, x_i,\, y_j,\, u(t_k,x_i,y_j) \bigr) \right\}, k=1,\ldots,7, i,j=1,\ldots,32 .\) Additive Gaussian noise is then applied, \[u^{\mathrm{obs}} = u + \sigma\varepsilon, \qquad \varepsilon \sim \mathcal{N}(0,1),\] where \(\sigma\in\{0.00,0.02,0.04,0.06,0.08,0.10,0.15\}\). As in the previous experiments, we consider two KAN parameterizations for the unknown kernel: the hard-constrained MC-KAN and the soft-constrained Cheb-KAN. In Cheb-KAN, positivity, monotonic decrease, and convexity are encouraged through equally weighted soft penalty terms with \(\lambda=10^{-2}\). The two models use depth-matched architectures with approximately the same number of trainable parameters. Here, both KANs take two inputs corresponding to the spatial displacements in the horizontal and vertical directions. Because the interaction kernel is even in each coordinate, the networks are parameterized in terms of the absolute displacements \(\tau_1=|x-\xi|\) and \(\tau_2=|y-\eta|\). For MC-KAN, the shape constraints are enforced independently along each lag direction. For both models, the inputs are linearly normalized. During training, the differentiable solver uses the same spatial grid, time step, and numerical discretization as the reference solver. At each optimization step, the full solution field is computed on the \(256\times256\) grid and evaluated against the observations in \(\mathcal{D}\). The loss is the MSE over the spatial observation points at each time snapshot, averaged over the snapshots. For evaluation, the learned spatial kernel is sampled on the same \(256\times256\) lag grid as the reference kernel, and the kernel reconstruction error is computed pointwise on this common grid. Complete numerical, architectural, and optimization settings are summarized in Table 9.

Table 5 reports the solution and kernel reconstruction errors across noise levels for the linear normalization. In the noise-free case, Cheb-KAN is marginally lower (\(2.64\%\) versus \(3.64\%\)). However, at every noisy level (\(\sigma \ge 0.02\)), MC-KAN achieves a lower kernel error \(\mathcal{E}(K)\) than Cheb-KAN, and the gap widens with noise, reaching nearly a factor of two at \(\sigma = 0.15\) (\(12.36\%\) versus \(21.95\%\)). The solution error \(\mathcal{E}(u)\) remains small for both methods and increases smoothly with noise. The two methods are comparable at low noise, with MC-KAN increasingly favored as noise grows. The solution error is an order of magnitude smaller than the kernel error, because the forward solution is less sensitive to small variations in the kernel than the inverse problem itself. The corresponding final-time solution reconstruction errors are shown in Appendix (Fig. 16).

Figure 5 shows the true kernel together with the pointwise absolute error of each method at \(\sigma=0.00\) and \(\sigma=0.08\). At zero noise both methods reconstruct the kernel accurately, with errors confined to a narrow region near the origin. Under noise the difference between the two becomes pronounced: the MC-KAN error remains concentrated near the axes, whereas the Cheb-KAN error spreads into broad off-axis regions where the soft penalty fails to enforce the correct decay. In both methods, however, and even at zero noise, the error is concentrated along the lines \(\tau_1=0\) and \(\tau_2=0\), where the kernel exhibits a sharp gradient near the origin. During training, the origin is excluded to represent general kernels that may exhibit singular behavior at \(\tau=0\).

To examine this behavior, in Fig. 6 we plot the recovered kernel along the two axis slices, \(K(\tau_1,\tau_2{\approx}0)\) and \(K(\tau_1{\approx}0,\tau_2)\). For consistency with the training setup, the one-dimensional slices are evaluated at \(\tau_2=10^{-3}\) and \(\tau_1=10^{-3}\), respectively. The true kernel rises steeply toward the origin, a consequence of the stretched exponents \(\beta_i<1\), and neither method fully resolves this sharp gradient with the linear input normalization. MC-KAN follows the true profile more closely than Cheb-KAN, but it too underestimates the peak near \(\tau_i\approx0\).

This difficulty motivates a change of input representation. The linear normalization maps the lag uniformly, allocating limited resolution to the near-origin region where the kernel varies most rapidly. We therefore replace it with a power-law normalization that concentrates resolution near the origin. Figure 7 (a) compares the kernel error across noise levels for the linear normalization, a learnable power-law exponent, and a fixed exponent. The learnable variant improves the reconstruction at low noise, where the steep gradient can be resolved, but its advantage diminishes as noise increases and the exponent becomes harder to identify reliably. We therefore first optimize the exponent in the noise-free setting, obtaining \(\alpha_1=\alpha_2=0.75\), and then fix this value for all subsequent noisy experiments. This fixed normalization retains most of the benefit of the learnable normalization while providing a single normalization that remains stable across all noise levels.

Finally, we increase the temporal sampling from \(7\) to \(25\) snapshots, holding the spatial sparsity fixed, to assess how much of the remaining error is data-limited. Table 6 reports the results for the fixed power-law MC-KAN against the Cheb-KAN baseline. The corresponding optimization loss histories are provided in Appendix (Fig. 17). The additional data reduce the kernel error substantially for both methods, and MC-KAN again achieves the lower error at every noise level. For example, at \(\sigma=0.15\) its kernel error falls to \(8.03\%\), roughly half that of Cheb-KAN (\(14.55\%\)) and well below the \(12.36\%\) obtained with seven snapshots. Despite the increased temporal resolution, a residual error remains near the sharp central peak (Fig. 8), likely because temporal information alone is insufficient to fully resolve this localized feature when the observations remain spatially sparse. To determine whether this residual error is intrinsic to the inverse problem or a consequence of the data sparsity, we repeat the experiment using dense \(256\times256\) spatial observations with \(25\) temporal snapshots. As shown in Fig. 7 (b)–(c), the reconstructed kernel now closely matches the true peak, meaning that the previous mismatch was primarily caused by the limited spatial sampling rather than by the representational capacity of the network or non-identifiability of the inverse problem. The dense observations provide sufficient information to resolve the steep spatial gradients near the origin.

The next step is to apply symbolic regression to recover a closed-form expression for the kernel from the optimized network. Unlike Experiments I and II, where the true kernel is one-dimensional, here it is a two-dimensional anisotropic function with a coupling term between the two spatial lags. Directly recovering such a multivariate coupled form is beyond the reach of symbolic regression in this setting, as the search space grows rapidly with input dimension. We therefore apply PySR to one-dimensional slices of the learned kernel, obtained by fixing one lag coordinate, and report the recovered expressions for Figs. 5, 6, and 8 in Table ¿tbl:tab:p3-tab3?. The table lists the expressions returned directly by PySR and, where possible, their rewriting in the \(\exp(-c\,\tau^{p})\) form in the last column. The \(K(\tau_1,\tau_2\approx0)\) slice is recovered more accurately than \(K(\tau_1\approx 0,\tau_2)\) in all cases, probably because the fast-decaying \(K(\tau_1\approx 0,\tau_2)\) slice is near-zero over most of the lag range, concentrating the fitting information near the origin. In contrast, the slowly-decaying \(K(\tau_1,\tau_2\approx0)\) slice varies across the whole domain. At zero noise with \(7\) snapshots, both models recover a physically correct stretched-exponential kernel. At \(\sigma=0.08\), MC-KAN still recovers the \(K(\tau_1,\tau_2\approx0)\) slice in the correct exponential form, while Cheb-KAN returns a non-admissible power-law form for both slices despite comparable \(L_2\) error. Increasing the number of snapshots to \(25\) results in exponential recoveries in more cases and lower errors overall. Notably, the MC-KAN \(K(\tau_1,\tau_2\approx0)\) error is almost unchanged between \(\sigma=0\) and \(\sigma=0.08\), indicating that the hard-constrained recovery is robust to noise for the well-resolved slice.

The two slices capture only the separable part of the kernel. To recover the coupling term, we compute the residual between the full learned kernel and its separable slices, then apply symbolic regression to it. Because the recovered slices are exponential in form, we analyze the kernel in log space, where \(\log K\) decomposes into the two marginal contributions and a residual, \[\log K(\tau_1,\tau_2) = \log K(\tau_1,\tau_2\approx0) + \log K(\tau_1\approx0,\tau_2) + K_{\mathrm{res}}(\tau_1,\tau_2),\] with \(K_{\mathrm{res}} = -\delta\,C\). For our best case, MC-KAN with \(25\) snapshots at zero noise, PySR recovers \(\delta\,C_{\mathrm{rec}} = 0.156\,\tau_1^{0.481}\,\tau_2\), compared with the true coupling \(\delta\,C = 0.289\,\tau_1^{0.5}\,\tau_2^{0.5}\), with relative \(L_2\) error \(0.20\). The recovered term reproduces the product form and the \(\tau_1\) exponent (\(0.48\) versus \(0.5\)), but underestimates the coupling coefficient (\(0.156\) versus \(0.289\)) and recovers the \(\tau_2\) dependence as approximately linear rather than \(\tau_2^{0.5}\).

5pt

4 Conclusion↩︎

In this study, we developed a differentiable-solver-based framework for discovering memory and nonlocal kernels in integro-differential equations (IDEs) directly from spatiotemporal observations. Within this framework, we investigated two KAN parameterizations for the unknown kernel: a hard-constrained Bernstein-based Monotone–Convex KAN (MC-KAN), in which positivity, monotonic decrease, and convexity are enforced by construction through constraints on the Bernstein coefficients, and a soft-constrained Chebyshev-based KAN (Cheb-KAN), in which the same physical properties are encouraged through penalty terms in the loss function. We systematically evaluated both approaches across benchmark IDEs of increasing complexity, from a Volterra equation to a viscoelastic wave PIDE and a two-dimensional nonlocal reaction–diffusion equation, under sparse and noisy observation settings. Our experiments demonstrate that both KAN parameterizations provide an effective framework for kernel identification. When the kernel family is relatively simple, as in the one-dimensional Volterra equation and the viscoelastic wave PIDE, both approaches recover accurate kernel representations and symbolic closed-form expressions, with MC-KAN achieving the lowest kernel reconstruction error in most configurations of Experiment I and at every sparsity level of Experiment II. The advantage of hard constraint enforcement becomes more pronounced for the sparse, noisy, and multidimensional nonlocal reaction-diffusion problem. In this setting, the hard-constrained MC-KAN produces more accurate kernel reconstructions and suppresses the nonphysical oscillations that remain possible under soft penalty enforcement. For the most challenging configuration considered, consisting of \(7\) temporal snapshots on a \(32\times32\) spatial grid with \(\sigma=0.15\) observation noise, MC-KAN with linear input normalization reduced the kernel reconstruction error by roughly \(44\%\) relative to Cheb-KAN. These results indicate that the main benefit of hard architectural constraints is not necessarily improved performance for simpler, better-observed inverse problems, but a more robust and physically consistent kernel representation as the observations become increasingly sparse and noisy.

The proposed framework has several limitations. First, the current hard-constrained MC-KAN is designed for kernels that are positive, monotonically decreasing, and convex, and therefore is not directly applicable to oscillatory, sign-changing, or non-monotone kernels. Second, although symbolic regression recovered accurate expressions for the one-dimensional benchmarks, the two-dimensional example required a staged procedure in which the principal-axis contributions and the coupling residual were recovered separately. The resulting coupled expression reproduced the product structure only approximately and did not fully recover the coefficient and exponent of the true interaction term. Future work will therefore focus on extending the framework to broader classes of physical kernels, developing symbolic discovery methods capable of recovering fully coupled multidimensional kernel functions, and validating the proposed approach on experimental datasets.

5 Acknowledgements↩︎

S.A.F. acknowledges support from the U.S. National Science Foundation under the Collaborations in Artificial Intelligence and Geosciences (CAIG, Award No. 2530611). This research was also supported in part by funding from the Wilkes Center for Climate Science & Policy (Award No. 01-00055-6000-44773) at the University of Utah.

6 Conflict of Interest↩︎

The authors declare no conflict of interests.

7 Code Availability↩︎

The source code will be made publicly available upon publication of this work.

8 Theoretical Guarantee↩︎

We now show that the kernel constraints introduced in Section 2.2.2 are preserved throughout the full depth of MC-KAN. Since the architecture enforces monotonicity and convexity independently along each input dimension, all results are stated and proved axis-wise.

Proof.  Positivity follows immediately from \(\phi>0\). For monotonicity, let \(x\le y\). Since \(t\) is monotone increasing, \(t(x)\le t(y)\). Since \(\phi\) is monotone decreasing, this implies \(\phi(t(x))\ge \phi(t(y))\). Hence \(\phi\circ t\) is monotone decreasing. To prove convexity, let \(x,y\in[a,b]\) and \(\lambda\in[0,1]\). Since \(t\) is concave, \[t(\lambda x+(1-\lambda)y) \ge \lambda t(x)+(1-\lambda)t(y).\] Because \(\phi\) is monotone decreasing, \[\phi(t(\lambda x+(1-\lambda)y)) \le \phi\!\left(\lambda t(x)+(1-\lambda)t(y)\right).\] Using the convexity of \(\phi\), \[\phi\!\left(\lambda t(x)+(1-\lambda)t(y)\right) \le \lambda \phi(t(x))+(1-\lambda)\phi(t(y)).\] Therefore, \[\phi(t(\lambda x+(1-\lambda)y)) \le \lambda \phi(t(x))+(1-\lambda)\phi(t(y)),\] which proves that \(\phi\circ t\) is convex. \(\square\)

Proof.  Positivity follows immediately from \(\phi>0\). On the interior of the domain, the chain rule gives \[(\phi\circ g)'=\phi'(g)g'.\] Since \(\phi'\ge0\) and \(g'\le0\), it follows that \[(\phi\circ g)'\le0.\] Similarly, \[(\phi\circ g)''=\phi''(g)(g')^2+\phi'(g)g''.\] Since \(\phi''\ge0\), \((g')^2\ge0\), \(\phi'\ge0\), and \(g''\ge0\), both terms are nonnegative. Therefore, \[(\phi\circ g)''\ge0.\] Thus \(\phi\circ g\) is positive, monotone decreasing, and convex on the interior. Since \(g\) and \(\phi\) are continuous on \([0,1]\), the composition is continuous on the closed domain, and the monotonicity and convexity properties extend to the boundary. \(\square\)

Layer 1 edge functions satisfy the decreasing coefficient conditions of Table 1, while all subsequent layers satisfy the corresponding increasing coefficient conditions. In addition, the reparameterization in  9 ensures that all Bernstein coefficients lie strictly in \((0,1)\). Therefore, every edge function maps \([0,1]\) into \((0,1)\), and all layer outputs, being averages of edge functions, also remain in \((0,1)\).

Proof.  The proof proceeds by induction on the layer index.

Base case. Let \[\mathbf{t}(\boldsymbol{\tau})=(t_1,\ldots,t_d)\] denote the normalized inputs. For each input dimension, \[t_i(\boldsymbol{\tau})=\left(\frac{\tau_i-\tau_i^{\min}}{\tau_i^{\max}-\tau_i^{\min}}\right)^{\alpha_i}, \qquad \alpha_i\in(0,1].\] The scalar map \(x\mapsto x^{\alpha_i}\) is monotone increasing and concave on \([0,1]\) for every \(\alpha_i\in(0,1]\). Therefore, each normalized input \(t_i\) is monotone increasing and concave with respect to \(\tau_i\) on \([\tau_i^{\min},\tau_i^{\max}]\). Each layer-1 edge function is positive, monotone decreasing, and convex as a function of its scalar argument by the coefficient conditions of Table 1. Let \[h_{p,q}^{(1)}=\phi_{p,q}^{(1)}(t_q).\] Fix an input dimension \(\tau_i\). If \(q\neq i\), then \(t_q\) does not depend on \(\tau_i\), so \(h_{p,q}^{(1)}\) is constant with respect to \(\tau_i\). Hence it is positive, monotone decreasing, and convex with respect to \(\tau_i\). If \(q=i\), then \(t_i\) is monotone increasing and concave with respect to \(\tau_i\), while \(\phi_{p,i}^{(1)}\) is positive, monotone decreasing, and convex. By Lemma 1, \[h_{p,i}^{(1)}=\phi_{p,i}^{(1)}(t_i)\] is positive, monotone decreasing, and convex with respect to \(\tau_i\) on the closed interval \([\tau_i^{\min},\tau_i^{\max}]\). Because layer outputs are averages of edge functions, \[z_p^{(1)}=\frac{1}{h_0}\sum_{q=1}^{h_0}h_{p,q}^{(1)}, \qquad h_0=d,\] the same properties hold for every layer-1 output. Moreover, since all Bernstein coefficients lie strictly in \((0,1)\), each edge function maps \([0,1]\) into \((0,1)\), and hence \[z_p^{(1)}\in(0,1).\] Inductive step. Assume that every output \(z_q^{(l)}\) is positive, monotone decreasing, and convex with respect to \(\tau_i\), and satisfies \[z_q^{(l)}\in(0,1).\] The output of layer \(l+1\) is \[z_p^{(l+1)}=\frac{1}{h_l}\sum_{q=1}^{h_l}\phi_{p,q}^{(l+1)}\!\left(z_q^{(l)}\right),\] where each \(\phi_{p,q}^{(l+1)}\) is positive, monotone increasing, and convex. By Lemma 2, every term \[\phi_{p,q}^{(l+1)}\!\left(z_q^{(l)}\right)\] is positive, monotone decreasing, and convex with respect to \(\tau_i\). Since sums and positive averages preserve positivity, monotonicity, and convexity, \(z_p^{(l+1)}\) inherits these properties. Furthermore, because \(z_q^{(l)}\in(0,1)\) and all Bernstein coefficients of \(\phi_{p,q}^{(l+1)}\) lie strictly in \((0,1)\), each edge output also lies in \((0,1)\). Therefore their average satisfies \[z_p^{(l+1)}\in(0,1).\] By induction, every layer output is positive, monotone decreasing, and convex with respect to \(\tau_i\), and remains in \((0,1)\). Since the kernel output is the final network output, \[K_{\boldsymbol{\theta}}=z^{(L)},\] the result follows. \(\square\)

9 Implementation↩︎

To enforce the coefficient conditions of Table 1, the Bernstein coefficients are generated from unconstrained trainable parameters through a differentiable reparameterization. For a Bernstein polynomial of degree \(n\), let \(\rho_k \in \mathbb{R}\), \(k=0,\ldots,n-2\), denote unconstrained parameters. We first define non-negative increments, \[s_k = \operatorname{ReLU}\bigl( \operatorname{softplus}(\rho_k)-\varepsilon \bigr),\] where \(\varepsilon>0\) is a small constant. For the first layer, a decreasing gap sequence is constructed as, \[g_{n-1} = -\operatorname{softplus}(\rho_{\mathrm{anchor}}),\] \[g_k = g_{k+1}-s_k, \qquad k=n-2,\ldots,0.\] The coefficient range is controlled through, \[c_{\mathrm{level}} = \operatorname{sigmoid}(\rho_{\mathrm{level}}),\] \[c_{\min} = c_{\mathrm{level}} \operatorname{sigmoid}(\rho_{\mathrm{floor}}),\] which satisfy \(0<c_{\min}<c_{\mathrm{level}}<1.\) The gap sequence is scaled according to, \[d = c_{\mathrm{level}}-c_{\min}, \qquad S = \sum_{k=0}^{n-1}|g_k|,\] \[\tilde{g}_k = g_k \frac{d}{S}.\] The Bernstein coefficients are then constructed recursively as, \[c_0 = c_{\mathrm{level}},\] \[c_{k+1} = c_k+\tilde{g}_k, \qquad k=0,\ldots,n-1.\] This construction yields, \[c_0 \ge c_1 \ge \cdots \ge c_n = c_{\min} > 0,\] with non-negative second differences. Consequently, the resulting Bernstein polynomial is positive, monotone decreasing, and convex.

For layers \(l\ge2\), monotone increasing edge functions are obtained by reversing the coefficient sequence, \[c_k^{\uparrow} = c_{n-k}, \qquad k=0,\ldots,n,\] which preserves convexity while reversing monotonicity. The learnable normalization exponent \(\alpha_i\) is parameterized through an unconstrained variable \(\tilde{\alpha}_i\) using the sigmoid transformation, \[\alpha_i=\sigma(\tilde{\alpha}_i)=\frac{1}{1+e^{-\tilde{\alpha}_i}},\] which guarantees that \(\alpha_i\in(0,1)\) throughout training. The network’s trainable parameters are, \[\boldsymbol{\theta} = \Bigl\{ \rho_k^{(l,p,q)}, \, \rho_{\mathrm{anchor}}^{(l,p,q)}, \, \rho_{\mathrm{level}}^{(l,p,q)}, \, \rho_{\mathrm{floor}}^{(l,p,q)}, \, \tilde{\alpha_i} \Bigr\},\] To improve numerical stability, Bernstein basis functions are evaluated in log-space, \[\log B_{k,n}(t) = \log\binom{n}{k} + k\log t + (n-k)\log(1-t),\] with inputs clamped to \(t \in [10^{-8},\,1-10^{-8}]\).

10 Numerical Solver Details↩︎

10.1 Experiment I: 1D Volterra IDE↩︎

10.2 Experiment II: 1D Viscoelastic Wave PIDE↩︎

10.3 Experiment III: 2D Nonlocal Reaction-Diffusion Equation↩︎

10.3.1 Initial Condition Construction↩︎

The initial condition is generated as a smooth random field in order to excite a broad range of spatial modes while avoiding grid-scale noise. Let \(\zeta(x,y)\) denote a realization of spatial white noise with independent standard normal samples. Its Fourier transform \(\widehat{\zeta}(\mathbf{k})\) is multiplied by a Gaussian low-pass filter, \[\widehat{u}_0(\mathbf{k}) = \widehat{\zeta}(\mathbf{k}) \exp\!\left( -\frac{\|\mathbf{k}\|^2}{2\sigma_s^2} \right),\] where \(\mathbf{k}=(k_x,k_y)\) is the wavevector and \(\sigma_s\) controls the smoothness of the resulting field. The initial condition is then obtained by inverse Fourier transform, \[u_0(x,y) = \mathcal{F}^{-1} \!\left[ \widehat{u}_0(\mathbf{k}) \right].\] Finally, the field is normalized to have prescribed standard deviation \(a\) and zero mean, \[u_0(x,y) \leftarrow a\, \frac{ u_0(x,y) }{ \operatorname{std}(u_0) },\] followed by \[u_0(x,y) \leftarrow u_0(x,y) - \frac{1}{|\Omega|} \int_{\Omega} u_0(x,y)\,dx\,dy.\] In this work, we use \[\sigma_s = 0.1, \qquad a = 0.5.\] The same initial condition is used for all simulations. This initialization serves two purposes. From a physical perspective, it represents a spatially heterogeneous state containing fluctuations across a range of length scales, which is characteristic of many phase separation and pattern-formation processes. From a numerical perspective, the broadband excitation ensures that the evolving field contains sufficient information about the underlying nonlocal interactions, producing gradients of adequate magnitude during backpropagation to enable reliable kernel identification.

References↩︎