¹

²

1 Introduction↩︎

Originally introduced to handle nonsmooth solutions of elliptic equations, very weak formulations (VWF) trace back to the transposition method of Lions and Magenes [1] [2]. In the finite element method (FEM), for instance, they have also proven useful for the numerical treatment of problems with low regularity data. Beyond these well known applications they provide a natural alternative to avoid automatic differentiation in neural networks (NN) based numerical methods for partial differential equations. Moreover, they are suitable for the use of NN spaces with low regularity. Although these obvious particularities do not seem to have gone completely unnoticed, to the best of our knowledge, there are no studies exploring their potential applicability in detail. Several factors might explain this fact. Differentiation based methods are more natural and there is an obvious reluctance to employ very weak formulations in general contexts, particularly if smooth solutions are expected. Another disadvantage of these methods lies in the requirement of regular test functions -a computational challenge for some numerical methods such as the FEM - and the fact that they, in its simplest form, offer error estimates in \(L^2\) or even weaker norms. On the other hand, in recent years, deep learning has emerged as a powerful tool for approximating solutions to partial differential equations (PDEs) by minimizing loss functionals that enforce governing equations and boundary conditions [3]–[8], He20?. Although not always competitive with classical methods in low-dimensional regimes, neural networks exhibit significant promise in high-dimensional contexts and on complex geometries [9], [10]. Also, formulations based on first-order systems [4], [5], [11], [12] circumvent the need for second-order derivatives in the cost functional. However, these approaches impose still regularity requirements on the trial space, which precludes the use of discontinuous activation functions. Conversely, employing trial spaces with reduced regularity, such as those utilizing ReLU or step-function architectures, offers substantial computational advantages regarding energy consumption, processing speed, and memory efficiency. This is evidenced by the growing interest in quantized networks as a robust alternative for efficient computation. Ultimately, discontinuous and highly quantized neural representations are particularly well-suited for hardware-aware implementations, reducing computational complexity and aligning with the constraints of emerging low-precision computing paradigm [13], [14].

In this work, we propose a theoretical and computational framework grounded in a Neural Very Weak Formulation (NVWF), specifically tailored to this class of networks. Numerical results illustrate the potential of this approach; although not competitive on high precision architectures compared to more traditional NN methods, it enables a mathematically rigorous formulation for problems involving both smooth and singular data and solutions. More importantly, it demonstrates the viability of developing neural PDE solvers on binarized architectures.

This paper is organized as follows. In Section 2 we present the very weak residual least square formulations that are the basis of our algorithms. Section 3 introduce the implementation details and Section 4 is devoted to present numerical simulations.

2 VWF and NVWF↩︎

We begin this section by reviewing some elementary aspects of very weak formulations (VWF). Let \(\Omega \subset \mathbb{R}^n\) be a bounded domain with boundary \(\partial \Omega\). For simplicity, we consider the model problem: Find \(u\) such that

\[\label{eq:poisson} \left\lbrace \begin{align} -\Delta u&=f,\, \text{in }\, \Omega\\ u&=g,\,\text{on } \, \partial \Omega. \end{align} \right.\tag{1}\]

The regularity of solutions to 1 depends on both, the regularity of the data and the domain \(\Omega\) [2]. To simplify the exposition, we temporarily assume \(\Omega\) smooth (or a convex polygon in \(\mathbb{R}^2\)). Consequently, if \(f\in L^2(\Omega)\) and \(g\in H^{3/2}(\partial \Omega)\), then \(u\in H^2(\Omega)\) [2]. Multiplying 1 by a smooth test function \(w\in H^2(\Omega)\cap H_0^1(\Omega)\) vanishing on \(\partial \Omega\), and integrating by parts twice yields \[\label{eq:vweaknh} -\int_{\Omega}u\Delta w=-\int_{\Gamma}g\frac{\partial w}{\partial\eta}+\int_\Omega fw ,\tag{2}\] where \(\eta\) is the outward unit normal to \(\Gamma\). For homogeneous Dirichlet conditions (\(g=0\)), this reduces to \[\label{eq:vweak} -\int_{\Omega}u\Delta w=\int_\Omega fw .\tag{3}\] Notice that 3 is well defined for any \(u\in L^2(\Omega)\), motivating the VWF of 1 which, for homogeneous boundaries, reads:

(VWF) Find \(u\in L^2(\Omega)\) satisfying 3 for all \(w\in H^2(\Omega)\cap H_0^1(\Omega)\).

Well-posedness for this VWF follows from classical inf-sup arguments. Discrete approximations, like those given by the finite element method (FEM), similarly rely on discrete inf-sup conditions requiring tailor-made spaces that are often challenging to construct. The high regularity required for test functions further discourages a direct FEM approach.

Here, we employ neural network (NN) trial spaces. Because the natural solution space for the VWF is \(L^2(\Omega)\), NN spaces with low regularity suffice. For the test spaces, we use scaled translations of a smooth, compactly supported generating function \(\varphi\) (in our case we take \(\varphi\in C^\infty_0(\mathbb{R}^n)\)). Such translation invariant spaces, under appropriate assumptions on \(\varphi\) are sufficiently rich [15]. In particular, they can approximate functions defined in \(\mathbb{R}^n\) up to a prescribed small saturation error, even for a function \(\varphi\) with a moderately large support. This crucial fact, advantageously attenuates numerical errors during the stochastic evaluation of oscillatory integrals arising from integration by parts. On bounded domains, these properties require careful handling. For smooth \(\Omega\), Sobolev extension theorems [16] allow functions on \(\Omega\) to be treated as restrictions of globally smooth functions. Building on this, we propose two distinct NN ansatzes, (A1) and (A2), using test functions not strictly supported in \(\Omega\). Both make use of a broadly supported generating test function within a residual least-squares formulation. Specifically, we define \(\varphi\) as a fixed scaling of the standard mollifier \[\label{eq:rho} \rho(x) := \left\lbrace \begin{align} e^{-\frac{1}{1 - \|x\|^2}} \quad & \text{if } \|x\| \leq 1,\\ 0 \quad & \text{if } \|x\| > 1.\\ \end{align} \right.\tag{4}\] For a fixed parameter \(h > 0\) (chosen to be relatively large), we set \(\varphi (x)=\rho_h(x) := \frac{1}{h^n}\rho\left(\frac{x}{h}\right)\). The support of \(\varphi\) is thus the closed ball \(\bar{B}_h(0)\). Other generating test functions enjoying similar approximation properties perform comparably.

Hereafter, \(v_\theta\) denotes a generic element of an NN space. Let \(D_{\partial \Omega}\) and \(D_{\partial \Omega}^S\) represent the distance and signed distance functions to \(\partial \Omega\), respectively, with \(D_{\partial \Omega}^S>0\) inside \(\Omega\). We introduce surrogates \(d_{\partial \Omega}\) and \(d_{\partial \Omega}^S\), which are equivalent to their exact counterparts up to multiplicative constants. These surrogates can be chosen to be as smooth as \(\Omega\) up to (or across) the boundary. For smooth \(\Omega\), we may assume \(\frac{\partial d_{\partial \Omega}}{\partial \eta}(x)=1\) on \(\partial\Omega\), and although \(d_{\partial \Omega}:\bar\Omega\to \mathbb{R}\) cannot be smooth if \(\Omega\) is non-smooth [17], certain Lipschitz domains (e.g., piecewise smooth) admit a surrogate satisfying \(\frac{\partial d_{\partial \Omega}}{\partial \eta}(x)=1\) almost everywhere on \(\partial\Omega\). We exploit this fact in ansatz (A2) below.

For simplicity, we restrict our attention to the homogeneous case \(g=0\) in our numerical experiments, the general case can be handled similarly. Here, Hardy’s inequality [18] guarantees that \(\frac{u}{d_{\partial \Omega}}\in L^2(\Omega)\), even for solutions with low regularity (e.g., for \(u\in H_0^1(\Omega)\) only). In both (A1) and (A2), we seek approximate solutions to 3 of the form \(u_\theta=d_0 v_\theta\), where \(d_0 \in \{d_{\partial \Omega}, d_{\partial \Omega}^S\}\). We see that shifting the unknown from \(u_\theta\) to \(v_\theta\) requires NN architectures with \(L^2(\Omega)\) approximation properties (a well-established feature for standard activation functions [19]). Remark [rem:density] confirms this also holds for the one-bit quantized networks used herein.

2.1 (A1): NVWF for smooth domains↩︎

Our first ansatz operates similarly to the fictitious domain method (FDM) in the sense that assumes the solution defined outside \(\Omega\) with a strong imposition of the Dirichlet boundary condition (FDM typically employs a weak imposition via a Lagrange multiplier [20] instead). Then, calling \(u_\theta=d_{\partial \Omega}^Sv_\theta\), we look for minimizers of \[\label{eq:fun1} \mathcal{L}(u_{\theta})=\int_{\Omega} \left[ \int_{B_h(x)} \tilde{f}(y) \cdot \rho_h(x-y) \, dy + \int_{B_h(x)} u_{\theta}(y) \cdot \Delta \rho_h(x-y) \, dy \right]^2 \, dx\tag{5}\] where \(\tilde{f}\) stands for an extension of \(f\), preserving the smoothness of \(f\) across \(\partial \Omega\) (note that \(B_h(x)\cap \Omega^c\) may be nonempty). We assume \(d^S_{\partial\Omega}\) is sufficiently smooth in the neighborhood of \(\Omega\) given by \(\tilde{\Omega}:=\{x\in \mathbb{R}^n: d^S_{\partial\Omega}(x)>-h \}\) (note that \(B_h(x)\subset \tilde{\Omega}\) for all \(x\in \Omega\)). Notice that equating the bracketed term in 5 to zero implies, assuming a smooth enough \(u_\theta\) and after integration by parts \[\int_{B_h(x)} \left(\tilde{f}(y) + \Delta u_{\theta}(y) \right)\cdot \rho_h(x-y) \, dy=0\] which for a fixed \(h\), and any \(x\in \Omega\), establishes the orthogonality of the residual w.r.t. the space generated by translations of \(\rho_h\).

2.2 (A2): NVWF for more general domains↩︎

If we assume the existence of a surrogate function \(d_{\partial \Omega}\) for which \(\frac{ \partial{d_{\partial \Omega}}}{\partial \eta}(x) =1\) a.e. along \(\partial \Omega\) then, calling \(u_\theta=d_{\partial \Omega}v_\theta\), we can look for minimizers of \[\label{eq:fun2} \mathcal{J}(u_{\theta})=\int_{\Omega} \Bigg[ \int_{B_h(x) \cap \Omega} f(y) \cdot \rho_h(x-y) \, dy + \int_{B_h(x) \cap \Omega} u_{\theta}(y) \cdot \Delta \rho_h(x-y) \, dy\tag{6}\] \[+ \int_{B_h(x) \cap \partial \Omega} \rho_h(x-y) \cdot v_{\theta}(y) \, dS \Bigg]^2 \, dx\] Replacing the expression \(u_\theta\) by \(d_{\partial \Omega}v_\theta\), we get, at least formally and after integration by parts twice \[\int_{B_h(x) \cap \Omega} u_{\theta}(y) \cdot \Delta \rho_h(x-y) \, dy=\int_{B_h(x) \cap \Omega} \Delta u_{\theta}(y) \cdot \rho_h(x-y) \, dy-\int_{B_h(x) \cap \partial \Omega} \rho_h(x-y) \cdot v_{\theta}(y) \, dS ,\] thanks to the fact that \(d_{\partial \Omega}\) vanishes identically and \(\frac{ \partial{d_{\partial \Omega}}}{\partial \eta}(x) =1\) a.e. \(x\) on \(\partial \Omega\), while for \(y\in \partial B_h(x)\), \(\rho_h(x-y)=0=\frac{\partial\rho_h}{\partial \eta}(x-y)\). In particular, equating to zero the bracketed term in 6 implies \[\int_{B_h(x)\cap \Omega} \left( f(y) + \Delta u_{\theta}(y) \right)\cdot \rho_h(x-y) \, dy=0,\] for fixed \(h\) and any \(x\in \Omega\), giving again an orthogonality of the residual w.r.t. translations of \(\rho_h\).

3 Implementation details↩︎

As detailed above, our experiments employ the test space generated by translations of the function \(\rho \in C_0^{\infty}(\Omega)\) defined in 4 . Observe that the Laplacian of \(\rho\) is explicitly given by \[\label{eq:lap95rho} \Delta \rho(x) = \left\lbrace \begin{align} \frac{ (8 - 2d)\|x\|^4 + 4(d - 1)\|x\|^2 - 2d}{(1 - \|x\|^2)^4}e^{-\frac{1}{1 - \|x\|^2}} \quad & \text{if } \|x\| \leq 1,\\ 0 \quad & \text{if } \|x\| > 1,\\ \end{align} \right.\tag{7}\] which scales according to the relation \(\Delta \rho_h(x) = \frac{1}{h^{d+2}}\Delta\rho\left(\frac{x}{h}\right)\). Given the oscillatory behavior of \(\Delta \rho\), the numerical integration of terms involving this function requires careful treatment.

3.1 Sampling strategy↩︎

The discretization of the functionals 5 and 6 relies on approximating the involved integrals through Monte Carlo methods. While the exterior integral is computed using uniform sampling, the interior integrals are resolved using Monte Carlo integration augmented by importance sampling. Empirical results suggest that a balanced sampling strategy, specifically centered on the integral involving the Laplacian, improves the stability of the method. More precisely, setting \(x = 0\) and \(h = 1\) and letting \(B = B_1(0)\) for brevity, we employ the following approximation:

\[\int_{B} u_{\theta}(y) \cdot \Delta \rho(-y) \, dy \approx \frac{ \int_{B} | \Delta \rho | }{4N_y} \Big( \sum_{y^+ \in Y^+} u_{\theta}(y^+) + u_{\theta}(-y^+) - \sum_{y^- \in Y^-} u_{\theta}(y^-) + u_{\theta}(-y^-) \Big),\] where \(Y^+ \subset \{ \Delta \rho > 0 \}\) and \(Y^- \subset \{ \Delta \rho < 0 \}\) are sets of points sampled according to the probability measures induced by the positive and negative parts of the Laplacian, \([\Delta \rho]^+\) and \([\Delta \rho]^-\), respectively. Here, \(N_y\) denotes the cardinality of each set, being \(N_y = \# Y^+ = \# Y^-\). By sampling an equal number of points in both regions and including their symmetric counterparts, we effectively emulate a centered difference operator, as illustrated in Figure 1.

3.2 Discretization of the loss functional \(\mathcal{L}\)↩︎

Following Section 2.1, we call \(\tilde{f}\) a smooth extension of the source term \(f\). The trial solution is then constructed as \(u_\theta = d_{\partial \Omega}^S v_\theta\), where \(v_\theta\) represents a neural network parameterized by \(\theta\), and \(d_{\partial \Omega}^S\) stands for a surrogate signed distance function. Consequently, we define:

\[\label{eq:discrete95func95L} L(u_\theta) := \frac{|\Omega|}{N_x} \sum_{x \in X} \Bigg[ \frac{C_f}{N_f} \sum_{y \in Y^f_x } f(y) + \frac{ C_{\Delta} }{4N_y} \Big( \sum_{y^+ \in Y^+_x} u_{\theta}(y^+) + u_{\theta}(-y^+) - \sum_{y^- \in Y^-_x} u_{\theta}(y^-) + u_{\theta}(-y^-) \Big) \Bigg]^2.\tag{8}\] As before, \(X\) represents a set of \(N_x\) points sampled uniformly over \(\Omega\), and \(Y^\Gamma\) denotes a set of \(N_\Gamma\) points sampled uniformly over the boundary \(\Gamma = \partial \Omega\). For each \(x \in X\), the set \(Y^f_x \subset B_h(x)\) is composed of \(N_f\) points sampled according to the probability measure induced by \(|\rho_h(x - \cdot)|\). Analogously, the sets \(Y^+_x, Y^-_x \subset B_h(x)\) are generated following the strategy outlined in Section 3.1 by sampling from the measure induced by \(|\Delta \rho_h(x - \cdot)|\), where \(N_y = \#Y^+_x = \#Y^-_x\). The constants \(C_f\) and \(C_\Delta\) are consistent with the definitions provided in 10 .

3.3 Discretization of the loss functional \(\mathcal{J}\)↩︎

The evaluation of the functional \(\mathcal{J}\) requires computing integrals over the intersections \(B_h(x) \cap \partial \Omega\) and \(B_h(x) \cap \Omega\). To simplify the numerical treatment of these terms, we perform the integration over \(\partial \Omega\) and \(B_h(x)\) instead. This can be achieved extending by zero the involved integrands. Specifically, let \(f_0: \mathbb{R}^n \to \mathbb{R}\) be the extension of \(f\) defined as \(f_0(x) = f(x)\) for \(x \in \Omega\) and \(f_0(x) = 0\) for \(x \in \Omega^c\). We represent the trial solution as \(u_\theta = d_0 v_\theta\), where \(d_0\) is a surrogate distance function to the boundary \(d_{\partial \Omega}\) extended by zero in \(\Omega^c\) and \(v_\theta\) is a neural network with parameters \(\theta\). Since \(\text{supp } \rho_h(x - \cdot) = B_h(x)\), the functional \(\mathcal{J}\) in 6 can be reformulated as:

\[\mathcal{J}(u_{\theta})=\int_{\Omega} \Bigg[ \int_{B_h(x) } f_0(y) \cdot \rho_h(x-y) \, dy + \int_{B_h(x) } u_{\theta}(y) \cdot \Delta \rho_h(x-y) \, dy + \int_{ \partial \Omega} \rho_h(x-y) \cdot v_{\theta}(y) \, dS \Bigg]^2 \, dx.\]

Building upon these observations, we define the discrete functional \(J\), which approximates \(\mathcal{J}\), as follows:

\[\label{eq:discrete95func95J} J(u_\theta) := \frac{|\Omega|}{N_x} \sum_{x \in X} \Bigg[ \frac{C_f}{N_f} \sum_{y \in Y^f_x } f_0(y) + \frac{ C_{\Delta} }{4N_y} \Big( \sum_{y^+ \in Y^+_x} u_{\theta}(y^+) + u_{\theta}(-y^+) - \sum_{y^- \in Y^-_x} u_{\theta}(y^-) + u_{\theta}(-y^-) \Big)\tag{9}\] \[+ \frac{|\Gamma|}{N_{\Gamma}} \sum_{y \in Y^\Gamma} \rho_h(x-y) \cdot v_\theta(y) \Bigg]^2.\] In this formulation, \(X\) represents a set of \(N_x\) points sampled uniformly over \(\Omega\), while \(Y^\Gamma\) denotes a set of \(N_\Gamma\) points sampled uniformly over the boundary \(\Gamma = \partial \Omega\). For every \(x \in X\), the set \(Y^f_x \subset B_h(x)\) contains \(N_f\) points sampled according to the probability measure induced by \(|\rho_h(x - \cdot)|\). Analogously, the sets \(Y^+_x, Y^-_x \subset B_h(x)\) are generated following the strategy described in Section 3.1 by sampling from the measure induced by \(|\Delta \rho_h(x - \cdot)|\), where \(N_y = \#Y^+_x = \#Y^-_x\). The constants \(C_f\) and \(C_\Delta\) are given by:

\[\label{eq:constants} C_f := \int_{B_h(0)} |\rho_h|, \quad C_\Delta := \int_{B_h(0)} |\Delta \rho_h|.\tag{10}\]

3.4 Description of the method↩︎

The proposed methodology focuses on approximating the minimizers of the functionals \(J\) and \(L\). Following the theoretical framework established in Section 2, we utilize the functional \(L\) specifically when the domain \(\Omega\) is characterized by a smooth boundary. This approach enables us to bypass the numerical evaluation of the boundary flux term inherent in \(J\). The implementation details for both variants of the proposed method are synthesized in Algorithms 3 and 2.

3.5 Architectures with step activation functions↩︎

We use the notation \(\mathcal{N}(L,P)\), to denote a generic neural network, with depth \(L\in \mathbb{N}\) (including the input and output layer) and width vector \(P = (p_1, \dots, p_{L-2})\). Each \(p_i\in \mathbb{N}\), for \(1\le i\le L-2\), represent the width of layer \(i\). We are interested in functions defined in \(\mathbb{R}^d\), returning scalar values and therefore only the width of hidden layers is declared. Furthermore, unless otherwise is stated, the simplified notation \(\mathcal{N}(L,P)\) assumes that only Heaviside activations \(\sigma_s: \mathbb{R}^{p_l} \to \mathbb{R}^{p_{l}}\) are applied. In this context, it is well known that even architectures with two hidden layers, i.e. \(\mathcal{N}(4, (m,m))\), can approximate functions in \(L^p({\mathbb{R}^d})\), for \(1\le p<\infty\) if we let \(m\to \infty\) (see, for instance, [19]). Regardless this theoretical result, for our implementation, we adopt the augmented Deep Heaviside Network (DHN) architecture proposed by [21]. This model extends the classical threshold network by integrating a small number of neurons with linear activations in each hidden layer alongside the standard Heaviside units. These linear components function as parallel information channels, allowing the network to preserve and propagate continuous input features throughout the deeper layers of the model.

The primary advantage of this hybrid design lies in its ability to overcome the representational bottlenecks typical of purely discrete activation functions. By combining the efficiency of Heaviside units with the expressivity of identity mappings, the architecture not only guarantees a richer feature propagation through deeper layers but also inherently mitigates the gradient flatlining problem during optimization, preserving active error signals across the continuous channels.

Formally, we denote the class of such networks as \(\mathcal{N}(L, P, s)\), where \(s\) is the number of augmented linear units per layer. The network is defined by the recursion \(x_{l} = \sigma_s(W_l x_{l-1} + b_l)\), where the hybrid operator \(\sigma_s: \mathbb{R}^{p_l + s} \to \mathbb{R}^{p_l + s}\) applies the Heaviside activation to the first \(p_l\) units and the identity mapping to the remaining \(s\) neurons.

Beyond its approximation properties, the localized discrete nature of the Heaviside units makes this architecture highly amenable to deployment on specialized digital accelerators (e.g., FPGAs). While the forward pass involves evaluating standard MAC operations for the linear bypass, the step activations drastically reduce the complexity of the nonlinear mapping overhead during the dense quadrature evaluations required by the very weak formulation.

The parameter optimization in our numerical experiments is performed via a surrogate gradient technique. Specifically, we employ the straight-through estimator (STE) with a sigmoid surrogate activation.

4 Computational experiments↩︎

We evaluate the proposed methodology across several test cases. Numerical results demonstrate the robustness of the scheme under non-ideal regimes, successfully recovering the expected profiles even when the exact solutions exhibit low regularity.

Optimization relies on the Adam algorithm across all computational configurations. To accelerate convergence, we employ a linear learning rate scheduler that decays monotonically to zero at the final epoch, with initial rates in the \([10^{-3}, 10^{-2}]\) range. We test our framework using distinct neural network trial spaces, denoted by \(u_*\), where \(* \in \{ \text{Tanh}, \text{ReLU}, \text{Step} \}\) specifies the activation function. The \(\text{Step}\) variant designates a multi-layer feedforward network parameterized by discontinuous Heaviside operators. As formalized in Section 3.5, these architectures incorporate a fixed proportion of linear units within each hidden layer. Furthermore, to explore low-precision limits, some experiments employ a quantized architecture, \(u_{\text{Qstep}}\), where hidden layer weights are constrained to 1-bit precision (\(\{-1, 1\}\)), while biases and output weights retain 32-bit floating-point precision. Network depth and width are specified for each respective test case. Finally, approximation accuracy is quantified via the error \(e_* := u_* - u_{\text{GT}}\) in the \(L^2\) norm, where \(u_{\text{GT}}\) denotes the exact analytical solution.

The numerical simulations were implemented in PyTorch, and the source code is available on GitHub ³. The current software implementation of the architectures described in Section 3.5 and their quantized variants is executed on conventional hardware (CPUs/GPUs). Consequently, this work does not exploit low-level hardware-specific parallelism, and the present experiments are not intended to quantify execution efficiency. Nevertheless, our empirical findings demonstrate that these Heaviside-based models, in both single-precision and 1-bit quantized formats, achieve competitive approximation accuracy while maintaining training stability under gradient-based optimization. These results establish a numerical proof of concept for the proposed framework, validating its mathematical viability prior to deployment on dedicated low-precision hardware.

Example 2 (Smooth Solution - Smooth Domain).

We consider the unit ball \(\Omega = B_1(0) = \{ x \in \mathbb{R}^2 : \|x\| < 1 \}\) and seek a function \(u: \Omega \to \mathbb{R}\) such that

\[\label{eq:suave2} \left\lbrace \begin{align} -\Delta u & = -9\|x\| \quad & in\Omega,\\ u & = 0 \quad & on\partial\Omega.\\ \end{align} \right.\qquad{(2)}\] The solution to ?? is \[u_{\text{GT}} = \|x\|^3 - 1.\]

To approximate the solution to this problem, we implement the algorithmic procedure outlined in Algorithm 2. The functional spaces for all baseline trial functions are parameterized via NN consisting of two hidden layers with an internal dimensionality of \(128\) neurons per layer. Specifically, \(u_{\text{Step}}\) incorporates a network topology \(\mathcal{N}(4, (128,128), 4)\) (see Section 3.5). To explore low-precision operational regimes, this example further incorporates a 1-bit weight quantization network \(u_{\text{QStep}}\) \(\mathcal{N}(4, (256, 256), 4)\) configuration.

Numerical results and the selected hyperparameter configurations, are summarized in Table 2 and visualized in Figure 5.

Table 2: Hyperparameters and resulting relative \(L^2\) errors for Algorithm [alg:dominio95ficticio] applied to Example [sec:suave2] (Smooth Solution - Smooth Domain).

Architecture	\(h\)	\(N_x\)	\(N_y\)	\(\|e_*\|_2 / \|u_{\text{GT}}\|_2\)
\(u_{\text{Tanh}}\)	0.2	20	100	0.0038
\(u_{\text{ReLU}}\)	0.2	20	100	0.0022
\(u_{\text{Step}}\)	0.2	100	20	0.013
\(u_{\text{QStep}}\)	0.2	40	100	0.043

Example 3 (Smooth Solution - Non-smooth Domain). We consider the following 2-dimensional problem. Let \(\Omega = [-1,1]^2,\) and \(k \in \mathbb{N}\). We seek \(u \colon \Omega \to \mathbb{R}\) such that \[\label{eq:example1} \left\lbrace \begin{align} -\Delta u & = \sin( \pi x_1) \left( \pi^2 (1-x_2^2) + 2 \right) \quad & in\Omega,\\ u & = 0 \quad & on\partial\Omega.\\ \end{align} \right.\qquad{(3)}\] The solution to ?? is \[u_{\text{GT}} = \sin( \pi x_1)(1-x^2_2).\]

Due to the non-smooth character of the boundary \(\partial \Omega\), this example is treated using the framework described in Algorithm 3. In alignment with our previous setup, the approximation spaces for all trial functions are parameterized via multilayer neural networks featuring two hidden layers of \(128\) neurons each. For \(u_{\text{Step}}\) integrates a neural network \(\mathcal{N}(4, (128,128), 16)\) (see Section 3.5). This example is also tested evaluating a quantized variant of the step-activation architecture, denoted as \(u_{\text{QStep}}\), using a network configuration \(\mathcal{N}(4, (2048,2048), 4)\). It is worth noting that this network topology is substantially larger than the one deployed in our previous experiments. This architectural expansion compensates the representation capacity degradation induced by 1-bit quantization when dealing with more complex solutions.

Table 3 and Figure 6 summarize the obtained numerical results and hyperparameter configurations.

Table 3: Hyperparameters and resulting relative errors for Algorithm [alg:flujo] applied to Example [sec:suave] (Smooth Solution - Non-smooth Domain).

Architecture	\(h\)	\(N_x\)	\(N_y\)	\(N_b\)	\(\|e_*\|_2 / \|u_{\text{GT}}\|_2\)
\(u_{\text{Tanh}}\)	0.2	30	100	5,000	0.0071
\(u_{\text{ReLU}}\)	0.2	30	100	5,000	0.0047
\(u_{\text{Step}}\)	0.2	30	100	5,000	0.041
\(u_{\text{QStep}}\)	0.25	150	100	8,000	0.067

Example 4 (High-dimensional case).

We consider nonlinear second-order Sine-Gordon equations characterized by anisotropic and inseparable solutions. To evaluate the performance of the proposed method in high-dimensional regimes, we adopt the following exact solution established in the literature [22]–[24]:

\[u_{\text{GT}}(x) = \left(1 - \|x\|^2\right) \left( \sum_{i=1}^{d-1} c_i \sin(x_i + \cos(x_{i+1}) + x_{i+1} \cos(x_i)) \right),\] where the coefficients \(c_i\) are sampled from a standard normal distribution \(\mathcal{N}(0, 1)\). This solution exhibits correlation between the \(i\)-th and \((i+1)\)-th coordinates, frequently referred to as two-body interaction.

Given this exact solution, we consider the Sine-Gordon equation subject to homogeneous Dirichlet boundary conditions:

\[\label{eq:hdim} \left\lbrace \begin{align} \Delta u(x) + \sin (u(x)) &= g(x), \quad & in\Omega,\\ u & = 0 \quad & on\partial\Omega,\\ \end{align} \right.\qquad{(4)}\] where \(\Omega = \{x \in \mathbb{R}^d : \|x\| < 1\}\) and the source term is prescribed as \(g(x) = \Delta u_{\text{GT}}(x) + \sin (u_{\text{GT}}(x))\).

To approximate the solution in this test case, we implement the framework described in Algorithm 2. Following the architectural choices of our previous experiment, all trial functions are parameterized using multilayer neural networks equipped with two hidden layers of \(128\) neurons each. In particular, the formulation for \(u_{\text{Step}}\) involves a network mapping characterized by the configuration \(\mathcal{N}(4, (128,128), 4)\), as detailed in Section 3.5. The underlying system relies on a collection of randomly initialized coefficients \(c_i\), which are explicitly given by: \[\begin{align} c = [&-0.2793, -1.5054, \phantom{-}0.6111, \phantom{-}0.7158, -0.9995, \phantom{-}1.2583, -1.1201, \phantom{-}1.1827, \\ &\phantom{-}0.1987, -1.1433, -1.3718, -0.3441, -1.0641, \phantom{-}0.7060, -2.0442, \phantom{-}0.1833, \\ &\phantom{-}0.0275, -0.1868, \phantom{-}1.0890] \end{align}\] The resulting numerical performance, alongside the comprehensive hyperparameter configurations, is documented in Table 4 and visualized in Figure 7.

Table 4: Hyperparameters and resulting relative errors for Algorithm [alg:dominio95ficticio] applied to Example [sec:hdim] (High-dimensional case).

Architecture	\(h\)	\(N_x\)	\(N_y\)	\(\|e_*\|_2 / \|u_{\text{GT}}\|_2\)
\(u_{\text{Tanh}}\)	0.15	400	100	0.011
\(u_{\text{ReLU}}\)	0.15	400	100	0.013
\(u_{\text{Step}}\)	0.15	400	100	0.020

Example 5 (Point source term).

To further evaluate the proposed method, we consider the Poisson equation with a Dirac delta distribution supported at a single point. Let \(\Omega = B(0,1) \subset \mathbb{R}^2\) denote the unit disk centered at the origin. We seek a solution \(u\) satisfying \[\begin{cases} -\Delta u = \delta & \text{in } \Omega, \\ u = 0 & \text{on } \partial\Omega, \end{cases}\] where \(\delta\) denotes the Dirac measure supported at the origin. The exact solution is given by the fundamental solution of the Laplacian in two dimensions, modified to satisfy the homogeneous Dirichlet boundary condition: \[u_{\text{GT}}(x) = -\frac{1}{2\pi} \log(\|x\|).\]

We approximate the solution using Algorithm 2. As before, the trial function approximation spaces employ neural networks with two \(128\)-neuron hidden layers. Specifically, \(u_{\text{Step}}\) uses the topology \(\mathcal{N}(6, (64, 64, 64, 64), 16)\) (Section 3.5). Table 5 and Figure 8 report the hyperparameters and numerical results.

In this example, the non-smooth architectures \(u_{\text{ReLU}}\) and \(u_{\text{Step}}\) outperform \(u_{\text{Tanh}}\). This performance gap is directly attributable to the singularity in the solution, a regime where non-smooth architectures inherently offer superior approximation capabilities.

Table 5: Hyperparameters and resulting relative errors for Algorithm [alg:dominio95ficticio] applied to Example [sec:dirac] (Point source term).

Architecture	\(h\)	\(N_x\)	\(N_y\)	\(\|e_*\|_2 / \|u_{\text{GT}}\|_2\)
\(u_{\text{Tanh}}\)	0.2	25	200	0.023
\(u_{\text{ReLU}}\)	0.2	25	200	0.013
\(u_{\text{Step}}\)	0.2	200	20	0.047

Example 6 (Segment source term).

We consider a Poisson equation featuring a singular source term supported on a one-dimensional segment. Let \(\xi\) denote the segment defined by \[\xi = \{(0, 0, z) \in \mathbb{R}^3 : -1 \le z \le 1\}.\] We prescribe the source term \(\gamma = \delta_\xi\). The corresponding free-space fundamental solution is expressed by the integral \[\begin{align} G(x,y,z) &= \int_{-1}^{1} \frac{1}{4\pi\sqrt{x^2+y^2+(z-t)^2}} \, dt \\ &= \frac{1}{4\pi} \log \left( \frac{\sqrt{x^2+y^2+(z-1)^2} + 1 - z}{\sqrt{x^2+y^2+(z+1)^2} - 1 - z} \right). \end{align}\]

For our numerical experiments, we define the domain \(\Omega\) as the solid ellipsoid bounded by the level set where the argument of the logarithm in \(G\) equals \(3\): \[\Omega = \left\{ (x, y, z) \in \mathbb{R}^3 : \frac{x^2}{3} + \frac{y^2}{3} + \frac{z^2}{4} \le 1 \right\}.\]

The boundary value problem seeks a function \(u\) such that \[\begin{cases} -\Delta u = \delta_\xi & \text{in } \Omega, \\ u = 0 & \text{on } \partial\Omega. \end{cases}\]

The exact singular solution to this problem is given by \[u_{\text{GT}}(x,y,z) = \frac{1}{4\pi} \log \left( \frac{\sqrt{x^2+y^2+(z-1)^2} + 1 - z}{3(\sqrt{x^2+y^2+(z+1)^2} - 1 - z)} \right).\]

In this example, we approximate the solution using Algorithm 2 and two-layer neural networks with \(128\) neurons each. Specifically, the architecture for \(u_{\text{Step}}\) is \(\mathcal{N}(6, (64, 64, 64, 64), 16)\) (Section 3.5). Table 6 and Figure 9 report the hyperparameters and numerical results.

As in the previous example, the non-smooth architectures \(u_{\text{ReLU}}\) and \(u_{\text{Step}}\) exhibit superior overall performance compared to \(u_{\text{Tanh}}\). Consistently, non-smooth architectures prove more suitable for approximating singular solutions.

Table 6: Hyperparameters and resulting relative errors for Algorithm [alg:dominio95ficticio] applied to Example [sec:segment] (Segment source term).

Architecture	\(h\)	\(N_x\)	\(N_y\)	\(\|e_*\|_2 / \|u_{\text{GT}}\|_2\)
\(u_{\text{Tanh}}\)	0.4	200	40	0.059
\(u_{\text{ReLU}}\)	0.4	200	40	0.047
\(u_{\text{Step}}\)	0.4	400	40	0.066

Example 7 (L-shaped domain).

We now consider a non-smooth domain featuring a re-entrant corner, a configuration where the solution is expected to exhibit a singular derivative. Let \(\Omega = B_1(0) \setminus [0,1]^2\). We seek a function \(u: \Omega \to \mathbb{R}\) such that \[\begin{cases} -\Delta u = 1 & \text{in } \Omega, \\ u = 0 & \text{on } \partial\Omega, \end{cases}\] In the absence of a closed-form expression for the exact solution, we evaluate the performance of our method against a reference solution, denoted by \(u_{\text{FEM}}\), obtained via the finite element method. The reference solution is computed using linear finite elements on a triangular mesh of approximately 20,000 nodes, with a mesh grading toward the origin to resolve the local singularity.

To accommodate the non-smooth boundary \(\partial\Omega\), we apply Algorithm 3. All trial functions are modeled using neural networks with two \(128\)-neuron hidden layers. Specifically, the architecture for \(u_{\text{Step}}\) is \(\mathcal{N}(4, (64, 64, 64, 64), 16)\) (Section 3.5). Table 7 and Figure 10 report the hyperparameters and numerical results.

Unlike the two preceding examples, the solution in this case exhibits a singularity only in its derivatives. Nonetheless, the smooth architecture \(u_{\text{Tanh}}\) underperforms relative to its non-smooth counterparts, \(u_{\text{ReLU}}\) and \(u_{\text{Step}}\).

Table 7: Hyperparameters and resulting relative errors for Algorithm [alg:flujo] applied to Example [sec:pacman] (L-shaped domain).

Architecture	\(h\)	\(N_x\)	\(N_y\)	\(N_b\)	\(\|e_*\|_2 / \|u_{\text{GT}}\|_2\)
\(u_{\text{Tanh}}\)	0.2	40	100	3,000	0.031
\(u_{\text{ReLU}}\)	0.2	40	100	3,000	0.017
\(u_{\text{Step}}\)	0.2	100	40	3,000	0.030

5 Conclusions↩︎

This work establishes a theoretical framework for employing low regularity neural networks for the numerical approximation of PDEs. Emphasis is placed on hardware-oriented architectures, specifically those utilizing Heaviside activations and quantized weights. While our current implementations do not directly exploit the computational efficiency of binarized networks, this study provides a proof of concept for future deployment on specialized hardware. Numerical experiments demonstrate the efficacy of the method, even for problems characterized by low regularity solutions. Notably, binarized architectures maintain robust performance when utilizing large support test functions and surrogate gradient methods during optimization. This indicates that neither the test space saturation error nor the intrinsic error of the surrogate gradient approximation poses a significant limitation. Nevertheless, empirical results reveal a sensitivity to the hyperparameters \(h\), \(N_x\), \(N_y\), and \(N_b\). Deriving their optimal scaling remain as an open problem.

Declaration of AI-Assisted Technologies in Manuscript Preparation↩︎

The authors acknowledge the use of AI tools to improve the English phrasing and grammar, assuming complete responsibility for all the scientific content.

References↩︎