From Monte Carlo to neural networks approximations of boundary value problems

arXiv:2209.01432 · math.PR, cs.AI, cs.LG, cs.NA, math.AP, math.NA · Submitted 2026-08-17 · Read on arXiv

Listen

Radio episode about this paper

Transcript

Introduction to the show: ident: AI Radio. Generated commentary on the latest Artificial Intelligence papers.

Tom: Next we'll be talking about the paper "From Monte Carlo to neural networks approximations of boundary value problems".

Jane: The paper was written by Lucian Beznea, Iulian Cimpean, Oana Lupascu-Stamate, Ionel Popescu and Arghir Zarnescu from University of Bucharest and Simion Stoilow Institute of Mathematics of the Romanian Academy and Gheorghe Mihoc – Caius Iacob Institute of Mathematical Statistics and Applied Mathematics of the Romanian Academy and Basque Center for Applied Mathematics and IKERBASQUE.

Tom: Stay tuned as we take you through the paper and discuss its implications.

Title and Authors: Tom: Welcome back, everyone! Today we’re cracking open a fresh one from the arXiv — it’s called “From Monte Carlo to neural networks approximations of boundary value problems.” Jane, I’ve got to say, that title alone tells you we’re in for a ride.

Jane: Oh, absolutely, Tom. And the author list is a who’s who of Romanian and Spanish applied math — Lucian Beznea, Iulian Cîmpean, Oana Lupascu-Stamate, Ionel Popescu, and Arghir Zarnescu. These folks are coming at this from probability, analysis, and numerical methods all at once.

Tom: So what’s the big idea here? Because “boundary value problems” sounds like something from a math textbook, but I’m guessing this is way more than that.

Jane: Right. So imagine you have a region — like a room, or a metal plate — and you want to know how heat spreads through it, or how a wave bounces around inside it. That’s a boundary value problem. You know what’s happening at the edges, and you want to figure out what’s happening inside.

Tom: And the classic way to solve that is with something like finite elements, where you chop the room into tiny pieces and compute on each one. But that blows up when the room is in, say, a hundred dimensions.

Jane: Exactly. And that’s where this paper comes in. They’re saying, instead of chopping up space, let’s use randomness. Think of it like a drunkard’s walk — you start at a point, take random steps, and see where you end up. Do that millions of times, and the average tells you the solution.

Tom: So it’s a Monte Carlo method — throwing dice, basically.

Jane: Yes, but with a twist. They’re using something called the Walk-on-Spheres algorithm, which is a clever way to speed up those random walks. Instead of taking tiny steps, you jump in big circles — spheres — and that gets you to the boundary much faster.

Tom: And then they take that whole thing and turn it into a neural network. So you don’t just get a number at one point — you get a function that works everywhere.

Jane: That’s the real magic. The Monte Carlo part gives you a fast, high-dimensional solver. The neural network part gives you a global approximation. And they prove that this whole pipeline doesn’t suffer from the curse of dimensionality.

Tom: The curse of dimensionality — that’s when the amount of work grows exponentially with the number of dimensions. And they’re saying, nope, we’re going to keep it polynomial.

Jane: Polynomial in the dimension and in the inverse of the error. That’s a big deal. It means you can actually solve these problems in high dimensions, like d equals a hundred or more, without your computer melting.

Tom: And that’s not just a toy — that’s finance, that’s physics, that’s materials science. So this paper is basically giving us a new tool for the toolbox.

Jane: And the authors are careful to make it constructive. They’re not just saying “a neural network exists.” They’re showing you how to build it, step by step, with explicit sizes and error bounds.

Tom: So we’ve got the title, we’ve got the team, and we’ve got the big idea — randomness plus neural nets to beat the high-dimensional curse. Next up, we’re going to dig into the actual method and why it’s so clever.

Summary of the Paper: Tom: Alright, so we’ve set the stage. Now let’s get into the meat of “From Monte Carlo to neural networks approximations of boundary value problems.” Jane, what’s the core problem they’re solving?

Jane: So they’re looking at the Poisson equation — that’s the classic “what’s the steady state of heat or charge” equation. You’ve got a domain, a source term inside, and a boundary condition. And they want to approximate the solution uniformly — meaning everywhere in the domain, not just at a few points.

Tom: And the old-school way to do that is to simulate Brownian motion — that’s the random walk I mentioned — and stop it when it hits the boundary. But that’s slow, and the error estimates are usually point-dependent. You get a good answer at one point, but you can’t reuse the same samples for another point.

Jane: Right. And that’s the first big contribution here. They show that you can use the same Monte Carlo samples to approximate the solution everywhere at once, with a uniform error bound. That’s a game-changer for parallel computing.

Tom: How do they pull that off?

Jane: They use a modified Walk-on-Spheres algorithm. Instead of stopping at a random time — which depends on where you start — they stop after a fixed number of steps, M. That’s uniform across all starting points. And they use a modified distance function, called r-tilde, which is a bit smaller than the true distance to the boundary.

Tom: Why make it smaller? That sounds like you’re being conservative for no reason.

Jane: Because it keeps the walk inside the domain, and it makes the whole thing compatible with neural networks. You see, the true distance function is hard to approximate with a ReLU network. But r-tilde — which is just a clipped version of an approximation — is easy. And they prove that using r-tilde instead of the true distance doesn’t hurt the accuracy much.

Tom: So they’re trading a tiny bit of precision for a huge gain in tractability.

Jane: Exactly. And then they get a beautiful error bound. The error splits into three parts: one from stopping early, one from the Monte Carlo averaging, and one from the grid discretization they use to control the sup-norm. And they show that all three can be made small with high probability.

Tom: And the key number — the complexity — what does that look like?

Jane: For a convex domain, they need M on the order of d log(d) steps, and N — the number of Monte Carlo samples — on the order of d squared log cubed of d. That’s polynomial in the dimension. No exponential blow-up.

Tom: So for d equals a hundred, that’s totally feasible on a GPU.

Jane: Totally. And they even test it numerically — they solve a Poisson problem in a hundred dimensions on a domain that’s a hypercube with a hole in it. And the errors go down as you increase N, just like the theory says.

Tom: That’s the kind of paper that makes you want to go run the code yourself.

Jane: And the best part is, they give you all the constants. You don’t have to guess. You can plug in your domain, your data, your error tolerance, and you know exactly how many samples and how many steps you need.

Tom: So we’ve got the method and we’ve got the bounds. But what about the neural network part? That’s where things get really interesting.

Improvements Suggested by the Paper: Tom: So we’ve seen the Monte Carlo side. Now let’s talk about the neural network side of “From Monte Carlo to neural networks approximations of boundary value problems.” Jane, what’s the big improvement here?

Jane: The big improvement is that they make the neural network construction fully explicit. You don’t have to train anything. You take the Monte Carlo estimator, you replace the distance function, the source, and the boundary data with their ReLU network approximations, and boom — you get a neural network that approximates the solution.

Tom: So it’s like a Lego set. You snap together the pieces and you get a solver.

Jane: Exactly. And the size of that network — the number of parameters — grows polynomially in the dimension and in the inverse of the error. That’s the curse-of-dimensionality breaker.

Tom: And they’re not just hand-waving. They give you the exact size formula.

Jane: Right. For a general domain satisfying the uniform exterior ball condition, the size is on the order of d to the seventh, times gamma to the minus sixteen over alpha minus four, times a bunch of log factors. That looks scary, but the point is it’s polynomial.

Tom: And if the domain is what they call “delta-defective convex” — which is a fancy way of saying it’s almost convex — they get an even better bound. The size drops to d cubed over gamma squared, times log factors.

Jane: And that’s a huge improvement over previous work. In an earlier paper by Grohs and Herrmann, the size depended on the volume of the domain, which can blow up exponentially in high dimensions. Here, the size only depends on the diameter and the geometry of the boundary.

Tom: So they’ve essentially removed the volume dependence. That’s the key improvement.

Jane: Yes. And they also handle much rougher data. The previous work assumed the source and boundary data were twice continuously differentiable. Here, they only need Hölder continuity — which is a much weaker condition. That covers a lot of real-world data, like piecewise smooth functions.

Tom: And they even show how to extend the boundary data into the domain in a way that’s compatible with neural networks. So you don’t need to know the data everywhere — just on the boundary.

Jane: Right. That’s the Corollary three point eight and three point one one part. They construct an extension using the nearest-point projection, and then they approximate that with a ReLU network. So the whole pipeline is end-to-end constructive.

Tom: And what about the practical side? Lu, you’re the AI researcher — what do you make of this?

Lu: I think the most exciting part is that this gives you a way to initialize a neural network with a provably good starting point. You don’t have to train from scratch. You can take this explicit construction, and then fine-tune it with gradient descent to get even better accuracy.

Tom: So it’s not just a theoretical existence result — it’s a practical starting point.

Lu: Exactly. And the fact that the network is random — because it depends on the Monte Carlo samples — means you get a distribution of networks, all of which are good with high probability. That’s a very robust setup.

Meng: As the engineer, I want to know: how does this actually run? Because those size formulas look big.

Jane: They’re big, but they’re polynomial. And the Monte Carlo part is embarrassingly parallel. You can run millions of random walks on a GPU simultaneously. The neural network part is just a feed-forward pass, which is also highly parallel.

Meng: So the bottleneck is memory, not compute?

Jane: Probably. But the authors show that the network width is bounded by something like two d plus the width of the distance approximation. So it’s not crazy wide.

Tom: So we’ve got the theory, we’ve got the construction, and we’ve got the practical implementation. What’s the takeaway for the world?

Conclusion: Tom: Alright, we’re wrapping up our look at “From Monte Carlo to neural networks approximations of boundary value problems.” Jane, give us the one-sentence summary.

Jane: It’s a paper that shows you can solve high-dimensional boundary value problems — like the Poisson equation — using a combination of Monte Carlo sampling and neural networks, with provable error bounds and polynomial complexity in the dimension.

Tom: And that’s not just a theoretical curiosity. That’s a practical tool for finance, physics, materials science — anywhere you need to solve a PDE in a high-dimensional space.

Jane: And the beauty is that it’s fully constructive. You don’t have to guess the architecture. You don’t have to train for days. You just plug in your domain, your data, and your error tolerance, and the paper tells you exactly how many samples, how many steps, and how big the network needs to be.

Lu: I’d add that this is a stepping stone. The same ideas — using a stochastic representation, accelerating it with a walk-type algorithm, and then converting it to a neural network — could apply to other equations. Parabolic equations, fractional Laplacians, even some nonlinear problems.

Meng: And from an engineering standpoint, the fact that it’s all parallelizable means it’s ready for modern hardware. You could implement this today on a GPU cluster.

Tom: So what’s the impact? Is this going to change how people solve PDEs?

Jane: I think it will change how people think about neural network solvers. Instead of treating them as black boxes that you train and hope for the best, this paper shows you can build them with guarantees. That’s a big deal for trust and reliability.

Tom: And it’s a bridge between two communities — the Monte Carlo people and the deep learning people. They’re speaking the same language now.

Jane: Exactly. And the authors even mention that the method is “grid-independent” — the grid is just a proof tool, not part of the algorithm. That’s elegant.

Tom: Alright, we’ve covered the title, the method, the improvements, and the implications. Time to say goodbye to this paper and get ready for the next one.

Jane: Thanks for joining us, everyone. We’ll see you on the next episode, where we’ll crack open another fresh paper from the arXiv.

Tom: Take care, and keep solving those boundary value problems!

Lucian Beznea, Iulian Cimpean, Oana Lupascu-Stamate, Ionel Popescu, Arghir Zarnescu

University of Bucharest · Simion Stoilow Institute of Mathematics of the Romanian Academy · Gheorghe Mihoc – Caius Iacob Institute of Mathematical Statistics and Applied Mathematics of the Romanian Academy · Basque Center for Applied Mathematics · IKERBASQUE

math.PR, cs.AI, cs.LG, cs.NA, math.AP, math.NA

Submitted: 2026-08-17

Updated: 2026-08-18

License: http://arxiv.org/licenses/nonexclusive-distrib/1.0/

Importance score: 69/100

Key concepts

Boundary Value Problem
This involves finding a solution inside a defined region, such as a room or plate, where conditions are specified at its edges (the boundary). Examples include calculating how heat spreads through a material or how waves bounce inside it.
Monte Carlo Method
This technique uses randomness to approximate solutions. Instead of precise calculations, it involves taking many random steps, like a drunkard's walk, and averaging the results over millions of trials to estimate the solution in high dimensions.
Curse of Dimensionality
This is a problem where the computational work required to solve a problem grows exponentially as the number of dimensions increases. The paper aims to show that their method keeps this complexity polynomial, allowing solutions for very high dimensions without excessive computation.
Neural Network Approximation
The authors use neural networks to create a global approximation of the solution. They build this network by replacing parts of the Monte Carlo estimator with ReLU network approximations, resulting in a solver that works everywhere in the domain.

Terminology

Summary

Summary

This paper studies probabilistic and neural network approximations for solutions to the Poisson equation subject to Hölder data in general bounded domains of Rd. The authors aim at two fundamental goals. The first, and the most important, they show that the solution to Poisson equation can be numerically approximated in the sup-norm by Monte Carlo methods, and that this can be done highly efficiently if we use a modified version of the walk on spheres algorithm as an acceleration method. This provides estimates which are efficient with respect to the prescribed approximation error and with polynomial complexity in the dimension and the reciprocal of the error. A crucial feature is that the overall number of samples does not depend on the point at which the approximation is performed. As a second goal, they show that the obtained Monte Carlo solver renders in a constructive way ReLU deep neural network (DNN) solutions to Poisson problem, whose sizes depend at most polynomialy in the dimension d and in the desired error. In fact they show that the random DNN provides with high probability a small approximation error and low polynomial complexity in the dimension.

The paper considers the Poisson boundary value problem:


1/2 ∆u = -f in D

u∂D = g,

where D is a bounded domain from Rd, whilst f: D → R and g: D → R are given continuous functions. It is well known that there exists a unique solution u ∈ C(D) ∩ H1 loc(D) to this problem.

The authors introduce a generalized walk on spheres process defined as:


X0 x:= x ∈ D,

X n+1 x:= X n x + r̃(X n x)U n+1, n ≥ 0,

where x is the starting point in the domain D, the function r̃ denotes the replacement of the distance function r to the boundary ∂D, and Ui are drawn independent and identically on the unit sphere in Rd. Essentially they use a Lipschitz r̃ such that r̃ ≤ r on the whole domain and βr ≤ r̃ as long as r ≥ ε. They call such a candidate a (β, ε)-distance. In all cases they can take β ≥ 1/3.

With this process at hand, they introduce the Monte Carlo estimator u N M of the solution u to the Poisson problem by:


u N M(x):= 1/N Σ i=1 N [g(X M x,i) + 1/d Σ k=1 M r̃ 2(X k-1 x,i) f(X k-1 x,i + r̃(X k-1 x,i)Y i)], x ∈ D, M, N ≥ 1.

Here, the sequences U n,i n,i≥0 and Y i i≥0 are all independent, U n,i is drawn uniformly on the unit sphere, X n x,i is given by the WoS process with U n replaced by U n,i, whilst Y i is drawn on the unit ball in Rd from the distribution µ which has an explicit density proportional to y 2-d - 1, y < 1 if d ≥ 3, and which is in fact the (normalized) Green kernel of the Laplacian on the unit ball with pole at 0.

The first part of the main result is the following theorem:

Theorem (Part I). Fix a small ε0 > 0, β ∈ (0, 1], r̃ a (β, ε0)-distance, and consider uM and u N M given by (2.20) and (2.29). Also, assume that f and g are α-Hölder on D for some α ∈ (0, 1]. Then, for any compact subset F ⊂ D, for all N, M, K ≥ 1, γ > 0 and ε ∈ [0, ε0], then, there are explicit quantities A(F, M, K, d, ε) and B(M, K, d) given in terms of the boundary regularity, the parameters N, M, K, d, ε, the set F and the data f, g such that


E sup x∈F u(x) - u N M(x) ≤ A(F, M, K, d, ε) + B(M, K, d)/√N.

Moreover, for an arbitrary domain D, for any compact subset F ⊂ D, we have that


lim M→∞ lim N→∞ E sup x∈F u(x) - u N M(x) = 0.

In addition, if the domain D satisfies the uniform ball condition, we can take F = D in the estimates.

The full power of the result is a little more technical and states that we have a tail estimate in the form:


P(sup x∈F u(x) - u N M(x) ≥ γ) ≤ 2 exp(C1(M, K, d) - ((γ - A(F, M, K, d, ε))+) squared / C2(M, d)),

where


C1(M, K, d):= d (⌈M/α⌉ log(2 + r̃ 1) + log(K)),

C2(M, d):= g∞ + M diam(D) squared f∞ / d,

A(F, M, K, d, ε):= 2 (g α + (diam(D) squared f α + 2diam(D)f∞)/d) (diam(D)/K) α

+ d α/2 g α v α/2(F, ε) + f∞ v(F, ε) + (4g∞ + diam(D) squared f∞) e-β squared ε squared M/(4diam(D) 2).

The authors point out that the estimates are true for any arbitrary domain, both in expectation and also in the tail. However, in this very general case they do not get any quantitative estimates, only asymptotic convergence guarantees on compact subsets.

The paper also provides a DNN counterpart. The essential part of this construction is the formula for the Monte Carlo estimator. The key fact is that they choose r̃ to be already a ReLU DNN. This is the main reason for their modification of the walk on spheres algorithm with r̃ instead of the usual distance to the boundary. In addition, having some ReLU DNN approximations of the data f respectively g, they can use these building blocks together with some basic facts about the ReLU DNNs in conjunction with the tail estimate to get the following result.

Theorem (Part II). Under the same context as Theorem (Part I), assume that D satisfies the uniform exterior ball condition, and that we are given ReLU DNNs ϕf: D → R, ϕg: D → R, ϕr: D → R such that


f - ϕf∞ ≤ ϵf ≤ f∞, g - ϕg∞ ≤ ϵg, r - ϕr∞ ≤ ϵr.

If γ > 0 is chosen such that a.1) - a.3) from Theorem 3.10 hold, and 0 < η < 1, then we can construct a (random) ReLU DNN U(x) such that


P(sup x∈D u(x) - U(·, x) ≤ γ) ≥ 1 - η,

with


size(U(ω, ·)) = O(d 7 γ-16/α-4 log(1/γ) [d 4 γ-4/α log(1/γ) + log(1/η)] S),

where


S:= [max(d, W(ϕr), L(ϕr)) + size(ϕr) + size(ϕg) + size(ϕf)].

Furthermore, if D is defective convex then they get the significant improvement


size(U(ω, ·)) = O(d 3/γ squared log(d/γ) [d squared log(d/γ) + log(1/η)] S).

Here size denotes the number of non-zero parameters in the neural network, W(ϕ) and L(ϕ) represent the width respectively the length of the neural network ϕ. The implicit constants depend on g α, g∞, f∞, diam(D), adiam(D), δ, α, log(2 + ϕr 1).

The paper also discusses several key issues that are tackled throughout the paper:

  • Overcoming the curse of high dimensionality: A very important consequence of the results concerns the breaking of the curse of high dimensions in the sense that the size of the neural network approximating the solution u adds at most a (low degree) polynomial complexity to the overall complexity of the approximating networks for the distance function and the data. In terms of the dimension d, the main results state, in particular, that if the domain is sufficiently regular (e.g. convex) then the complexity of the Monte Carlo estimator of the exact solution scales at most like d cubed log 4(d), whilst the DNN estimator of the same solution scales at most like d 5 log 5(d)S.

  • Low dimensions improvements for general bounded domains with Dirichlet data: The approach also has interesting consequences in low dimensions. The key is that the samples for the Monte Carlo solver can be reused to simultaneously approximate the solution for all points in the domain, and furthermore the number of the steps required by the designed Walk-on-Spheres algorithm does not depend on the starting point; these two features make the proposed scheme highly parallelizable.

  • General bounded domains and Hölder continuous data: The curse of high dimensions can be overcome for a general class of domains, namely those that satisfy a uniform exterior ball condition. The results are even more general, covering the case of an arbitrary bounded domain in Rd, but then the estimates are given in terms of the behavior of the function vD in the proximity of the boundary of the domain. Concerning the regularity of the source and boundary data, the assumption is that they are merely Hölder continuous.

  • L∞(D) estimates: The errors are estimated in the uniform norm which gives much better results. The uniform norm of the error is small with large probability, whilst the approximation complexity depends on D merely through its (annular) diameter.

  • Walk-on-Spheres acceleration revisited: The walk on spheres (WoS) is modified in two respects. Firstly, the stopping rule for the walk on sphere is deterministic, namely, the walk-on-spheres chain is run for a given number of steps, uniformly for all trajectories and all points in the domain. Secondly, the walk-on-spheres scheme is performed with the maximal radius replaced by a more general radius, which is not necessarily maximal and is compatible with ReLU DNNs. Overall, a generalized walk-on-spheres algorithm is developed which is of self interest and which is much more compatible with parallel computing.

  • Universality with respect to given data: The estimator explicitly constructed essentially approximates the operator that maps the data (source and boundary) of the problem into the corresponding solution u. In particular, it means that the DNN solvers constructed consist of the composition of two separate neural networks: one which approximates the source and boundary data and one for the above-mentioned operator.

  • Explicit construction of the approximation: One key element of the approach is the explicit formulation of the approximation. This is reflected in the formula for the Monte Carlo estimator where all elements are fully determined. This structure can be exploited to initialize a DNN with significantly less complexity than the guaranteed Monte Carlo construction.

The paper also includes numerical results. The numerical tests are conducted for the domains Dc = [-1,1] d and Dac = [-1,1] d x ∈ Rd: x 1 ≤ 0.5. Test 1 numerically verifies the estimate for the probability that the WoS chain is still at distance at least ε from the boundary after M steps. Test 2 tests the approximation of the solution u by simulating its Monte Carlo estimator, confirming that the errors are decreasing to a small value as the number of WoS trajectories N increases. The numerical results are obtained for d = 10 and d = 100, and the tests turned out to be successful even on a worse domain geometry.

Improvements for AI systems

Based on the paper, here are specific improvements that can be made to AI systems, along with what the improved system can do:

Improvement: Implement the Monte Carlo estimator (equation 1.5) with the modified Walk-on-Spheres (WoS) algorithm as a standalone solver for Poisson equations with Dirichlet boundary conditions.

What the improved system can do:

  • Solve Poisson equations in dimensions d = 100+ without the curse of dimensionality

  • Provide explicit L∞ error bounds with high probability (Theorem 2.26)

  • Guarantee polynomial complexity in dimension and error (O(d3log4d) for convex domains)

  • Work with only Hölder continuous data (α ∈ (0,1]), not just C2 functions

  • Handle arbitrary bounded domains, not just convex ones

Sources

Related papers