Comparative Analysis of New Unbiased and Biased Monte Carlo Algorithms for the Fredholm Integral Equation of the Second Kind

1. Introduction

Integral equations [1,19] serve as one of the cornerstones of modern mathematical analysis and applied computation [9,10] . They represent one of the most powerful analytical tools in applied mathematics and theoretical physics [11,12]. They arise naturally in a wide range of disciplines [13,14,15] such as potential theory, elasticity, electromagnetism, geophysics, kinetic gas theory, quantum mechanics, and financial mathematics. In many of these applications, the Fredholm integral equations of the second kind serve as the fundamental mathematical structure describing steady-state or equilibrium processes. Their versatility allows them to model a wide variety of physical and engineering phenomena—ranging from potential theory and heat conduction to electromagnetic scattering, quantum mechanics, and stochastic processes. When formulated in multidimensional domains, these equations become indispensable tools for describing systems characterized by spatial or temporal interactions over continuous spaces. Consequently, the development of accurate and efficient numerical methods [20,21,22] for solving such equations remains a topic of enduring significance.

A particularly important class is represented by the Fredholm integral equations of the second kind, which can be expressed in operator form as

$$u\left(x\right)={\int }_{\Omega }k(x,x^{\prime })u\left(x^{\prime }\right)d{x}^{\prime }+f\left(x\right),$$

(1)

or equivalently

$$u=Ku+f,$$

(2)

where K denotes the integral operator defined by kernel $k(x,x^{\prime })$, both $k(x,x^{\prime })$ and $f\left(x\right)$ belong to $L^2(\Omega )$, and the solution $u\left(x\right)$ is sought in the same functional space. Here, $\Omega$ represents a bounded domain in ${\mathbb{R}}^n$.

In many practical problems, the evaluation of linear functionals of the solution,

$$J\left(u\right)={\int }_{\Omega }\phi \left(x\right)u\left(x\right)dx=(\phi ,u),$$

(3)

is of primary interest. The function $\phi \left(x\right)$ is given and belongs to $L^2(\Omega )$. For such functionals, stochastic formulations based on the Monte Carlo method [23,25] have proven to be highly advantageous, especially in high-dimensional domains where deterministic quadrature approaches become computationally prohibitive.

Analytical solutions to such equations exist only for limited cases, and even then, their computation becomes infeasible when the domain $\Omega \subset {\mathbb{R}}^n$ is multidimensional or when the kernel exhibits complex structure.

Classical deterministic numerical techniques—such as quadrature-based discretization, Galerkin methods, and collocation schemes—tend to suffer from rapid growth in computational cost with increasing dimension n. This phenomenon, known as the curse of dimensionality [4,26], severely limits their applicability to high-dimensional integral equations.

Stochastic methods [16,17,18] and particularly Monte Carlo (MC) and quasi-Monte Carlo (QMC) algorithms, overcome this limitation by providing a convergence rate that depends weakly—or not at all—on the number of dimensions. These methods allow for direct estimation of functionals[3] of the solution without requiring the full discretization of the integral equation. Moreover, stochastic formulations are naturally parallelizable and provide probabilistic error estimates, making them highly suitable for high-dimensional or complex-kernel problems.

Two principal stochastic paradigms are typically considered:

• Biased stochastic methods [24], based on probabilistic representations of truncated Liouville–Neumann series, which yield efficient approximations at the cost of a systematic (truncation) error.
• Unbiased stochastic methods [27], designed such that the mathematical expectation of the estimator equals the exact solution, thereby eliminating systematic bias but possibly increasing variance.

Recent works by Maire and Dimov et al. [6,7] have significantly advanced both approaches, leading to new unbiased stochastic algorithms that exhibit excellent accuracy and scalability in multidimensional settings. The present paper continues this line of research by developing and analysing refined unbiased algorithms and systematically comparing them with the most efficient biased stochastic methods available in the literature.

Mathematical Formulation of the Problem 

Consider the Fredholm integral equation of the second kind on a bounded domain $\Omega \subset {\mathbb{R}}^n$:

$$u\left(x\right)={\int }_{\Omega }k(x,x^{\prime })u\left(x^{\prime }\right)d{x}^{\prime }+f\left(x\right),$$

(4)

or equivalently, in operator notation,

$$u=Ku+f,$$

(5)

where K is the integral operator defined by the kernel $k(x,x^{\prime })$. We assume that $k(x,x^{\prime })\in L^2(\Omega \times \Omega )$ and $f\left(x\right)\in L^2(\Omega )$, which ensures that the solution $u\left(x\right)\in L^2(\Omega )$ exists and is unique when ${\parallel K\parallel }_{L^2}<1$.

In many applications, it is not necessary to compute $u\left(x\right)$ everywhere in $\Omega$; rather, one seeks to evaluate a particular linear functional of the solution:

$$J\left(u\right)={\int }_{\Omega }\phi \left(x\right)u\left(x\right)dx=(\phi ,u),$$

(6)

where $\phi \left(x\right)$ is a known weight or test function in $L^2(\Omega )$. This formulation enables the computation of quantities such as average values, fluxes, or solution evaluations at specific points (e.g., $J\left(u\right)=u\left(x_0\right)$ for $\phi =\delta (x-x_0)$).

Liouville–Neumann Series Representation. 

A natural approach to solving (4) involves successive approximations, leading to the Liouville–Neumann expansion:

$$u^{\left(i\right)}=\sum _{j=0}^i{K}^{\left(j\right)}f=f+Kf+\dots +K^{\left(i\right)}f,i=1,2,\dots ,$$

(7)

where $K^{\left(0\right)}$ denotes the identity operator and $K^{\left(j\right)}$ is the j-fold composition of K. Under the contraction condition ${\parallel K\parallel }_{L^2}<1$, the sequence $\left\{u^{\left(i\right)}\right\}$ converges to the true solution u in $L^2$, i.e.,

$$u^{\left(i\right)}\to u,\mathrm{as}i\to \infty .$$

Accordingly, the desired functional (6) can be expressed as

$$J\left(u\right)=\underset{i\to \infty }{lim}(\phi ,u^{\left(i\right)})=\underset{i\to \infty }{lim}\sum _{j=0}^i(\phi ,K^{\left(j\right)}f).$$

(8)

This representation serves as the analytical foundation for both deterministic and stochastic approximation schemes.

Integral Representation of Iterative Terms. 

Each term in (8) can be expressed as a multidimensional integral:

$$\begin{array}{cc}\hfill I^{\left(j\right)}& =(\phi ,K^{\left(j\right)}f)={\int }_{{\Omega }^{j+1}}\phi \left(t_0\right)k(t_0,t_1)k(t_1,t_2)\dots k(t_{j-1},t_j)f\left(t_j\right)d{t}_{0}d{t}_1\dots d{t}_j.\hfill \end{array}$$

(9)

Defining

$$F_j(t_0,\dots ,t_j)=\phi \left(t_0\right)k(t_0,t_1)k(t_1,t_2)\dots k(t_{j-1},t_j)f\left(t_j\right),$$

(10)

we have

$$I^{\left(j\right)}={\int }_G{F}_j\left(t\right)dt,G={\Omega }^{j+1}\subset {\mathbb{R}}^{n(j+1)}.$$

Hence, computing a truncated approximation

$$(\phi ,u^{\left(i\right)})=\sum _{j=0}^i{I}^{\left(j\right)}$$

(11)

reduces to evaluating a finite number of high-dimensional integrals. This observation directly motivates the application of stochastic integration methods.

Alternative Formulation via the Adjoint Problem.

When viewed in a Banach space framework, Eq. (4) admits a dual (adjoint) formulation:

$$v=K^{*}v+\phi ,$$

(12)

where $v,\phi \in X^{*}(\Omega )$ and $K^{*}:X^{*}\to X^{*}$ is the adjoint operator of K. In this case, the same functional can be expressed as

$$J\left(u\right)=(f,v)={\int }_{\Omega }f\left(x\right)v\left(x\right)dx,$$

providing an alternative computational pathway—often advantageous when solving the adjoint equation is simpler or more stable.

Error Structure and Balancing. 

Approximating (8) by truncating the Neumann series at a finite level i introduces two distinct error components:

• The systematic (truncation) error $R_{\mathrm{sys}}=|(\phi ,u)-(\phi ,u^{\left(i\right)})|$, and
• The stochastic (sampling) error $R_N$, resulting from the numerical approximation of integrals $I^{\left(j\right)}$.

Hence, the total error satisfies

$$|(\phi ,u)-(\phi ,\tilde{u})|\le R_{\mathrm{sys}}+R_N<\epsilon ,$$

(13)

where $\tilde{u}$ denotes the stochastic estimate of $u^{\left(i\right)}$. Properly balancing these two errors—by tuning both i and the number of samples N—is essential for minimizing computational effort while maintaining accuracy. Optimal computational efficiency is achieved by balancing the two error components through an appropriate choice of i and N.

This work aims to develop, analyse, and compare advanced stochastic algorithms for the efficient solution of Fredholm integral equations, focusing on high-dimensional and complex-kernel settings. To set the stage, we begin by formulating the mathematical problem and reviewing the stochastic interpretation that underpins both biased and unbiased approaches. The present study specially focuses on developing a refined unbiased stochastic approach for multidimensional Fredholm integral equations. We first revisit the classical biased and unbiased algorithm [6], and we propose a novel unbiased stochastic algorithm (NUSA) that generalizes the traditional Unbiased Stochastic Approach (USA). The new method employs modified probabilistic transition mechanisms, inspired by the dynamics of continuous-state Markov chains, ensuring convergence even in multidimensional case.

The primary contributions of this paper can be summarized as follows:

• We present a unified theoretical and algorithmic framework for solving multidimensional Fredholm integral equations of the second kind using stochastic techniques.
• A systematic comparison between several biased stochastic algorithms—including the Crude Monte Carlo method, Markov Chain Monte Carlo (MCM) approach, and quasi-Monte Carlo algorithms based on modified Sobol $\Lambda \Pi \tau$ sequences—is performed.
• We propose and analyse a novel unbiased stochastic algorithm (NUSA), designed to overcome the limitations of the traditional unbiased approach (USA) and to improve numerical stability in high-dimensional settings.
• Comprehensive numerical experiments are conducted for both one-dimensional and multidimensional cases, demonstrating the superior accuracy and robustness of the proposed NUSA method compared to state-of-the-art algorithms.

In the subsequent section, we begin by revisiting and detailing the class of biased stochastic algorithms that serve as benchmarks for evaluating the performance of the new unbiased approaches. The rest of this paper is organized as follows. Next Section 2 revisits the conventional biased and unbiased stochastic approach and outlines its algorithmic framework. Then it introduces the newly developed NUSA algorithm, discussing both its theoretical foundation and implementation details for basic and general kernel scenarios. Section 3 presents numerical experiments, comparing the proposed method with existing techniques for one- and multi-dimensional problems. Finally, Section ?? summarizes the main contributions and highlights potential directions for future research.

2. Methods

2.1. Biased Stochastic Algorithms for Solving Fredholm Integral Equations

The first part of our analysis focuses on biased stochastic algorithms, which form the baseline for the unbiased approaches developed later in this work. These algorithms approximate the functional $(\phi ,u)$ through a truncated Liouville–Neumann expansion and probabilistic representations of its terms. Although the truncation introduces a systematic bias, these methods are valuable benchmarks for the performance and variance characteristics of the novel unbiased techniques proposed in this paper (see Section 2.4).

Following the methodology described in Section 1, the truncated expansion

$$(\phi ,u^{\left(i\right)})=\sum _{j=0}^i(\phi ,K^{\left(j\right)}f)$$

(14)

is approximated stochastically. In the context of this work, we consider four representative biased algorithms:

• the Markov Chain Monte Carlo (MCM) method for integral equations;
• the Monte Carlo algorithm based on Modified Sobol $\Lambda \Pi \tau$ sequences (MCA-MSS-1);
• the symmetrised version MCA-MSS-2;
• and the Stratified Symmetrised Monte Carlo approach (MCA-MSS-2-S).

These algorithms differ mainly in how they generate and distribute random (or quasi-random) samples to approximate multidimensional integrals.

2.1.1. Crude Monte Carlo Method (CMCM) for Evaluating Integrals

The Crude Monte Carlo Method (CMCM) provides the simplest stochastic approach for estimating integrals of the form

$$I^{\left(j\right)}={\int }_G{F}_j\left(t\right)dt,t\in G\subset {\mathbb{R}}^{n(j+1)}.$$

(15)

Let ${\left\{x^{\left(l\right)}\right\}}_{l=1}^N$ be a sequence of points uniformly distributed over G with total volume $V\left(G\right)$. Then, by the probabilistic interpretation of the integral,

$$I^{\left(j\right)}\approx I_N^{\left(j\right)}={\displaystyle \frac{V\left(G\right)}{N}}\sum _{l=1}^N{F}_j\left(x^{\left(l\right)}\right),$$

(16)

where the random variable $F_j\left(x^{\left(l\right)}\right)$ has expectation $\mathbb{E}\left[F_j\right]=I^{\left(j\right)}/V\left(G\right)$. The corresponding mean-square error is

$${\epsilon }_N^{\left(j\right)}=|I^{\left(j\right)}-I_N^{\left(j\right)}|\approx \sqrt{{\displaystyle \frac{\mathrm{Var}\left(F_j\right)}{N}}},$$

which demonstrates the $O\left(N^{-1/2}\right)$ convergence rate characteristic of Monte Carlo methods. This approach, although basic, serves as the foundation for more advanced biased stochastic algorithms.

Algorithm 1: Crude Monte Carlo Method (CMCM) for Integral Evaluation

The CMCM algorithm establishes the conceptual foundation for biased stochastic techniques such as the Markov Chain Monte Carlo and quasi-Monte Carlo methods based on Sobol $\Lambda \Pi \tau$ sequences, which are discussed next.

2.1.2. Markov Chain Monte Carlo Method for Integral Equations

The biased Monte Carlo approach is founded on representing the functional $(\phi ,u)$ as a finite sum of linear functionals $(\phi ,K^{j}f)$ for $j=0,\dots ,i$. Each term is estimated using trajectories generated by a Markov chain

$$T_i:x_0\to x_1\to \dots \to x_i,$$

(17)

constructed with initial probability density $\pi \left(x\right)$ and transition probability $p(x,x^{\prime })$ satisfying

$${\int }_{\Omega }\pi \left(x\right)dx=1,{\int }_{\Omega }p(x,x^{\prime })d{x}^{\prime }=1,\forall x\in \Omega .$$

Typically, $p(x,x^{\prime })$ is chosen proportional to $\left|k\right(x,x^{\prime }\left)\right|$, that is $p(x,x^{\prime })=c\left|k(x,x^{\prime })\right|$, where

$$c={\left({\int }_{\Omega }\left|k(x,x^{\prime })\right|d{x}^{\prime }\right)}^{-1}.$$

The corresponding random variable

$${\theta }_i\left[\phi \right]={\displaystyle \frac{\phi \left(x_0\right)}{\pi \left(x_0\right)}}\sum _{j=0}^i{W}_{j}f\left(x_j\right),W_0=1,W_j=W_{j-1}{\displaystyle \frac{k(x_{j-1},x_j)}{p(x_{j-1},x_j)}},j\ge 1$$

provides a biased estimator for $(\phi ,u)$:

$$(\phi ,u^{\left(i\right)})\approx {\displaystyle \frac{1}{N}}\sum _{l=1}^N{\theta }_i^{\left(l\right)}\left[\phi \right].$$

Although the estimate converges to $(\phi ,u^{\left(i\right)})$ with $O\left(N^{-1/2}\right)$ variance, it remains biased due to the truncation of the Liouville–Neumann series.

Algorithm 2: Biased MCM Algorithm for Integral Equations

This MCM-based scheme is straightforward to implement and can reuse the same random trajectories for different test functions or right-hand sides, provided the kernel $k(x,x^{\prime })$ remains fixed.

2.1.3. Monte Carlo Algorithms Based on Modified Sobol $\Lambda \Pi \tau$ Sequences

To enhance uniformity and reduce variance, quasi-Monte Carlo (qMC) methods employ deterministic low-discrepancy point sets. The Sobol $\Lambda \Pi \tau$ sequences form a class of $(t,s)$-sequences in base 2 with excellent equidistribution properties. Two modified algorithms, denoted MCA-MSS-1 and MCA-MSS-2, introduce randomization techniques (“shaking”) to further improve convergence.

2.1.4. MCA-MSS-1 Algorithm

MCA-MSS-1 introduces a random perturbation around each Sobol point within a small spherical region of radius $\rho \ll 1$. For an integrand $F_j\left(t\right)$ defined over $U_d={[0,1]}^d$, the randomized points ${\xi }^{\left(l\right)}\left(\rho \right)=x^{\left(l\right)}+\rho {\omega }^{\left(l\right)}$, where ${\omega }^{\left(l\right)}$ is uniformly distributed on the unit sphere, are used to estimate

$$I^{\left(j\right)}={\int }_{U_d}{F}_j\left(t\right)dt.$$

It is proven that $\mathbb{E}\left[F_j\left({\xi }^{\left(l\right)}\left(\rho \right)\right)\right]=I^{\left(j\right)}$. This leads to a Monte Carlo estimator with reduced discrepancy and improved convergence rate for smooth integrands.

Algorithm 3: MCA-MSS-1 Algorithm

For integrands $F_j\in W^1(L;U_d)$, the error decays as $O\left(N^{-1/2-1/d}\right)$.

2.1.5. MCA-MSS-2 Algorithm

MCA-MSS-2 extends the previous method by employing symmetric point pairs with respect to subdomain centers. Each unit cube $U_d$ is subdivided into $m^d$ subdomains $K_l$, each containing one Sobol point $x^{\left(l\right)}$. A shaken point ${\xi }^{\left(l\right)}$ and its symmetric counterpart ${\xi }^{\left(l\right)\prime }$ relative to the center $s^{\left(l\right)}$ are generated, providing an enhanced estimate:

$$I\left(f\right)\approx {\displaystyle \frac{1}{2N}}\sum _{l=1}^N\left[f\left({\xi }^{\left(l\right)}\right)+f\left({\xi }^{\left(l\right)\prime }\right)\right].$$

Algorithm 4: MCA-MSS-2 Algorithm

For functions with continuous and bounded second derivatives ($F_j\in W^2(L;U_d)$), the MCA-MSS-2 method achieves an optimal convergence rate $O\left(N^{-1/2-2/d}\right)$.

2.2. Stratified Symmetrized Monte Carlo (MCA-MSS-2-S)

A simplified version, denoted MCA-MSS-2-S, combines the stratification and symmetry ideas without explicitly using the fixed radius $\rho$. Each random point ${\xi }^{\left(l\right)}$ is drawn uniformly inside a subdomain $K_l$, and its symmetric counterpart ${\xi }^{\left(l\right)\prime }$ is computed relative to the local center. Although the control parameter $\rho$ is no longer used, this algorithm reduces computational cost while maintaining near-optimal accuracy.

Algorithm 5: MCA-MSS-2-S Algorithm

This method is computationally cheaper than MCA-MSS-2 but lacks the radius-based fine-tuning mechanism. While MCA-MSS-2-S lacks a tunable shaking radius, it remains nearly optimal for functions in $W^2(L;U_d)$ and is computationally more efficient for large-scale multidimensional problems.

2.3. Error Balancing in Biased Stochastic Approaches

Each biased stochastic algorithm approximates the functional $(\phi ,u)$ with two distinct errors:

• A systematic error (truncation bias) $R_{\mathrm{sys}}$, due to terminating the Liouville–Neumann expansion at depth i;
• A stochastic error $R_N$, arising from the finite number of Monte Carlo or quasi-Monte Carlo samples.

The overall accuracy requirement

$$|(\phi ,u)-(\phi ,\tilde{u})|\le R_{\mathrm{sys}}+R_N<\epsilon$$

(18)

is met by balancing both components. Increasing i reduces $R_{\mathrm{sys}}$ exponentially under ${\parallel K\parallel }_{L^2}<1$, while increasing N decreases $R_N$ algebraically as $N^{-1/2}$ or faster for quasi-Monte Carlo schemes. This principle of error balancing guides the numerical experiments presented later in this paper.

The algorithms described in this section establish a well-defined framework for biased stochastic approximations of Fredholm integral equations. In the following subsection, we turn our attention to unbiased stochastic algorithms, which remove the systematic bias entirely while preserving the stochastic convergence properties.

2.4. Unbiased Stochastic Algorithms for Fredholm Integral Equations

The second stage of our study focuses on unbiased stochastic algorithms, whose key advantage lies in the complete elimination of the systematic truncation error inherent to biased Monte Carlo approaches. The main objective of these algorithms is to construct random variables whose expected values coincide exactly with the analytical solution of the Fredholm integral equation (4), without truncating the Liouville–Neumann series. This idea was first introduced in [6] and further extended here through a new unbiased stochastic algorithm (NUSA) that improves variance control and numerical stability, particularly for multidimensional problems.

The unbiased stochastic approach (USA) transforms the integral equation into a probabilistic representation involving random trajectories (Markov chains). In contrast to biased methods, where the number of iterations is fixed, here each trajectory can continue indefinitely but terminates probabilistically via an absorption event. The expected accumulated contribution along all possible trajectories equals the true solution $u\left(x\right)$.

2.5. Classical Unbiased Stochastic Approach (USA)

2.5.1. Mathematical Formulation

Consider the integral equation

$$u\left(x\right)={\int }_{\Omega }k(x,x^{\prime })u\left(x^{\prime }\right)d{x}^{\prime }+f\left(x\right),$$

(19)

where $k(x,x^{\prime })\in L^2(\Omega \times \Omega )$ and $f\left(x\right)\in L^2(\Omega )$. The goal is to compute $u\left(x\right)$ as the expectation of a random functional constructed through a Markov process with transition and absorption probabilities derived from $k(x,x^{\prime })$.

For kernels satisfying $0\le k(x,x^{\prime })\le 1$, the probability of continuation is given by $k(x,x^{\prime })$, while the probability of absorption is $1-k(x,x^{\prime })$. If the chain begins at $x_0=x$, it evolves as

$$T_i:x_0\to x_1\to \dots \to x_i,$$

with

$$\begin{array}{cc}\hfill P(x_{i+1}\in \Omega \mid x_i\ne \partial )& =k(x_i,U_i)P(U_i\in \Omega ),\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill P(x_{i+1}=\partial \mid x_i\ne \partial )& =1-k(x_i,U_i),\hfill \end{array}$$

(21)

$$\begin{array}{cc}\hfill P(x_{i+1}=\partial \mid x_i=\partial )& =1,\hfill \end{array}$$

(22)

where ∂ denotes an absorbing state and $f(\partial )=0$.

The stochastic representation of the solution is then

$$u\left(x\right)=\mathbb{E}\left(\sum _{i=0}^{\tau }f\left(x_i\right)|x_0=x\right),$$

(23)

where $\tau$ is the random stopping time corresponding to absorption. Equation (23) satisfies the original Fredholm equation in expectation:

$$u\left(x\right)={\int }_{\Omega }k(x,y)u\left(y\right)dy+f\left(x\right).$$

2.5.2. Classical USA Algorithm for Simple Kernels

When the kernel satisfies $0\le k(x,y)\le 1$ for all $x,y\in {[0,1]}^n$, the unbiased algorithm takes a particularly simple form. Here, the term $(1-k(x,y\left)\right)$ can be interpreted as the absorption probability, representing the likelihood that the random trajectory terminates at a given state.

Algorithm 6: Classical USA for $0\le k(x,y)\le 1$

Initialization: Provide the kernel $k(x,y)$, source term $f\left(x\right)$, and the total number of stochastic trajectories N. Execution: (a) Set score = 0, choose the initial point $x_0=x$, and initialize test = 1. (b) For $j=1,\dots ,N$: While test $\ne 0$: i. Update score ←score + $f\left(x_0\right)$. ii. Generate two random numbers: $U\sim \mathcal{U}{[0,1]}^n$, $V\sim \mathcal{U}[0,1]$. iii. If $k(x_0,U)<V$, set test = 0 (absorption occurs); otherwise, update $x_0=U$ and continue. Estimation: Compute the final approximation $$u\left(x\right)={\displaystyle \frac{\mathtt{score}}{N}}.$$ (24)

This procedure performs a sequence of random walks, accumulating the contributions of $f\left(x\right)$ along each path until absorption occurs. The algorithm yields an unbiased estimator for $u\left(x\right)$ since the expectation over all random trajectories reproduces the integral equation (1).

2.5.3. Algorithmic Implementation for $0\le k(x,y)\le 1$

This method provides an unbiased estimator of $u\left(x\right)$, since the expected value over infinitely many random walks equals the analytical solution. Its computational cost depends primarily on the expected trajectory length, which in turn is governed by the magnitude of $k(x,x^{\prime })$.

Algorithm 6: Classical USA for $0\le k(x,y)\le 1$

2.5.4. General Kernel Case

When the kernel can assume negative or large values, i.e. $k(x,y)\in L^2(\Omega \times \Omega )$ with $\left|k\right(x,y\left)\right|\ge 1$ possible, the continuation and absorption rules are modified to maintain unbiasedness. The product weight $\mathtt{prod}$ accumulates the sign and magnitude of the kernel along each path. This general algorithm, which can be regarded as a stochastic analogue of the Neumann series representation, is summarized below.

Algorithm 7: General USA for Arbitrary Kernels

Initialization: Specify $k(x,y)$, $f\left(x\right)$, and the number of trajectories N. Execution: (a) Set score = 0. (b) For $j=1,\dots ,N$: Initialize $x_0=x$, test = 1, and prod = 1. While test $\ne 0$: i. Update score ←score + $f\left(x_0\right)\times$prod. ii. Generate $U\sim \mathcal{U}{(0,1)}^n$. iii. If $\left|k\right(x_0,U\left)\right|>1$, set prod$\leftarrow \mathtt{prod}\times k(x_0,U)$. iv. Otherwise, draw $V\sim \mathcal{U}(0,1)$ and: – If $\left|k\right(x_0,U\left)\right|<V$, set test$=0$ (absorption); – Else, set prod $\leftarrow \mathtt{prod}\times \mathrm{sign}\left(k(x_0,U)\right)$ and $x_0\leftarrow U$. Estimation: The solution is computed as $$u\left(x\right)={\displaystyle \frac{\mathtt{score}}{N}}.$$ (25)

2.5.5. Algorithmic Implementation for Arbitrary Kernels

Although the USA guarantees unbiasedness, the resulting estimator may exhibit high variance, particularly for large ${\parallel K\parallel }_{L^2}$ or highly oscillatory kernels. This motivates the development of a refined method, the Novel Unbiased Stochastic Approach (NUSA), presented below.

Algorithm 7: General USA for Arbitrary Kernels

2.5.6. Novel Unbiased Stochastic Approach (NUSA)

The NUSA algorithm generalizes the USA by introducing an adaptive weighting and absorption control mechanism that mitigates variance growth and improves convergence in multidimensional contexts. Let $Y\sim \mathtt{Uniform}\left({[0,1]}^n\right)$ and suppose $0\le k(x,Y)\le 1$. The stochastic representation of the solution is

$$u\left(x\right)=\mathbb{E}\left\{k(x,Y)u\left(Y\right)+[1-k(x,Y\left)\right]\right\}+f\left(x\right),$$

(26)

which leads naturally to a random-walk process with probabilistic continuation and absorption events governed by $k(x,Y)$.

2.5.7. Motivation and Theoretical Formulation

Let $T_i:x_0\to x_1\to \dots \to x_i$ denote a random trajectory initiated at $x_0=x$. The probabilistic dynamics of this procedure can be outlined as:

$$P({\mathrm{x}}_{i+1}\in \Omega |{\mathrm{x}}_i\ne \partial )=k(\mathrm{x},U_i)P(U_i\in \Omega ),$$

$$P({\mathrm{x}}_{i+1}=\partial |{\mathrm{x}}_i\ne \partial )=1-k(\mathrm{x},U_i),$$

and

$$P({\mathrm{x}}_{i+1}=\partial |{\mathrm{x}}_i=\partial )=1,\mathrm{and}\mathrm{f}(\partial )=0.$$

From the above, $u\left(\mathrm{x}\right)$ can be expressed as:

$$u\left(\mathrm{x}\right)=E\left(\sum _{i=0}^{\infty }f\left({\mathrm{x}}_i\right)|{\mathrm{x}}_0=\mathrm{x}\right)=E\left(\sum _{i=0}^{\tau }f\left({\mathrm{x}}_i\right)|{\mathrm{x}}_0=\mathrm{x}\right),$$

with

$$\tau =\underset{i\ge 0}{inf}\left({\mathrm{x}}_i\right)=\partial .$$

Conceptually, this expression can be reformulated and extended as:

$$\begin{array}{ccc}{\displaystyle u\left(\mathrm{x}\right)=E\left(\sum _{i=0}^{\infty }f\left({\mathrm{x}}_{i+1}\right)|{\mathrm{x}}_0=\mathrm{x}\right)+f\left(\mathrm{x}\right)}& =& {\displaystyle E\left(E\left(\sum _{i=0}^{\infty }f\left({\mathrm{x}}_{i+1}\right)|{\mathrm{x}}_1\right)\right)+f\left(\mathrm{x}\right)}\hfill \\ & =& {\displaystyle {\int }_{\Omega }k(\mathrm{x},\mathrm{y})u\left(\mathrm{y}\right)d\mathrm{y}+f\left(\mathrm{x}\right).}\hfill \end{array}$$

(27)

The algorithm introduces an auxiliary control parameter $P_d\in (0,1)$, representing the probability of forced absorption. This artificial damping shortens excessively long trajectories and reduces estimator variance without introducing bias.

2.5.8. Algorithm for Basic Kernels ($0\le k(x,y)\le 1$)

The following algorithm operationalizes the stochastic formulation for kernels satisfying $0\le k(x,y)\le 1$.

Algorithm 8: NUSA for Basic Kernels

Initialization: Provide the kernel $k(x,y)$, source function $f\left(x\right)$, number of random trajectories N, and evaluation point $x_0$. Execution: (a) Set $S=0$. (b) For $k=1$ to N (Monte Carlo cycles): Initialize $\mathtt{weight}=1$, $x=x_0$, and define a fixed probability threshold $P_d=0.5$. Set $\mathtt{test}=0$. While $\mathtt{test}\ne 1$: i. Generate $U\sim \mathcal{U}(0,1)$. ii. If $U<P_d$, set $\mathtt{test}=1$ (absorption occurs). iii. Otherwise: A. Select V according to the probability density proportional to $\left|K\right(x,·\left)\right|$. B. If $K(x,V)<0$, update $\mathtt{weight}\leftarrow \mathtt{weight}\times (-1)/(1-P_d)$ and set $x=V$. End while. Update $\mathtt{weight}\leftarrow \mathtt{weight}\times f\left(x\right)/P_d$ and accumulate $S\leftarrow S+\mathtt{weight}$. Estimation: The unbiased estimator of the solution is $$u\left(x\right)={\displaystyle \frac{S}{N}}.$$ (28)

The algorithm simulates a collection of stochastic trajectories originating at $x_0$, each subjected to probabilistic absorption determined by $k(x,y)$. Unlike the traditional USA, NUSA introduces a normalization factor ${(1-P_d)}^{-1}$ and an adaptive weight mechanism that compensate for variable kernel amplitudes, thereby reducing estimator variance.

2.5.9. Algorithmic Implementation for Basic Kernels

This adaptive mechanism ensures the expected length of trajectories remains finite, thus improving numerical stability while preserving the unbiased nature of the estimator.

Algorithm 8: NUSA for Basic Kernels ($0\le k(x,y)\le 1$)

2.5.10. Algorithm for General Kernels

For general kernels that may assume negative values or magnitudes exceeding unity, NUSA employs an auxiliary sampling distribution $M(x,y)$ to normalize transition probabilities. The control parameter $P_d$ again governs the average trajectory length.

2.5.11. NUSA Algorithm for General Kernel Scenarios

When the kernel $k(x,y)$ may assume negative values or magnitudes exceeding unity, the above procedure is extended by redefining the sampling distribution and weight accumulation to preserve unbiasedness. The corresponding method is described below.

Algorithm 9: NUSA for General Kernels

Initialization: Specify the kernel $k(x,y)$, source term $f\left(x\right)$, number of trajectories N, and initial point $x_0$. Execution: (a) Initialize $S=0$. (b) For $k=1$ to N: Set $\mathtt{weight}=1$, $x=x_0$, and $P_d=0.5$. While $\mathtt{test}\ne 1$: i. Generate $U\sim \mathcal{U}(0,1)$. ii. If $U<P_d$, set $\mathtt{test}=1$ (absorption event). iii. Otherwise: A. Draw V using the probability distribution $M(x,V)$ associated with the kernel. B. Update $\mathtt{weight}\leftarrow \mathtt{weight}\times {\displaystyle \frac{K(x,V)}{M(x,V)(1-P_d)}}$ and set $x=V$. Estimation: After all trajectories are simulated, compute $$u\left(x\right)={\displaystyle \frac{S}{N}},$$ (29) where S is the accumulated weighted score.

This generalized NUSA framework accommodates a wide class of kernels, including those with sign variations and large magnitudes. The absorption probability $P_d$ controls the expected trajectory length: increasing $P_d$ leads to shorter walks but higher variance, while smaller values extend the paths at the expense of computation time. In practice, setting $P_d\approx 0.5$ achieves a favorable balance between computational efficiency and statistical accuracy.

2.5.12. Algorithmic Implementation for General Kernels

By decoupling the kernel sampling from the absorption probability, this formulation achieves enhanced stability and flexibility. Empirical studies (see Section 3) show that NUSA attains substantially lower relative errors compared to both USA and traditional biased methods for the same computational cost.

The USA and NUSA methods share the same theoretical foundation—representing the Fredholm equation as a stochastic expectation—but differ in their variance control strategies. While USA offers an exact theoretical framework, its variance can increase significantly for strongly coupled kernels. NUSA mitigates this issue by introducing adaptive absorption control and kernel-weight normalization, producing smoother convergence behavior, especially in higher dimensions.

Algorithm 9: NUSA for General Kernels

2.6. Extensions of the Unbiased Approach

The unbiased Monte Carlo algorithms described in the preceding sections (USA and NUSA) provide pointwise estimates of the solution $u\left({\mathbf{x}}_0\right)$ to the Fredholm integral equation of the second kind. However, in many scientific and engineering applications, interest extends beyond point evaluations to more global quantities of physical or statistical relevance. Two important generalizations of the unbiased approach address these broader objectives. The first aims at estimating weak or functional averages of the form ${\int }_{\Omega }\phi \left(\mathbf{x}\right)u\left(\mathbf{x}\right)d\mathbf{x}$, which frequently arise in sensitivity analysis, uncertainty quantification, and model reduction. The second extension reconstructs the full spatial distribution of the solution $u\left(\mathbf{x}\right)$ from data collected during the random walk simulations, enabling global field approximation across the entire computational domain. Both procedures retain the unbiased character of the original algorithms while exploiting the statistical information generated by multiple random trajectories. The following subsections detail the methodology and implementation of these two extensions.

The unbiased Monte Carlo framework can be extended beyond pointwise estimation of the solution to encompass the computation of linear functionals and global field approximations. Such extensions are of significant practical relevance in applications where the complete solution $u\left(\mathbf{x}\right)$ is not required explicitly, but rather its weak form or statistical averages involving a prescribed weighting function $\phi \left(\mathbf{x}\right)$. Typical examples include sensitivity analysis [5], uncertainty quantification, and mean-value estimation in systems with spatially localized observables.

2.6.1. Weak Approximation

In many cases, the quantity of interest is not $u\left({\mathbf{x}}_0\right)$ at a single point, but the functional

$$I_{\phi }={\int }_{\Omega }\phi \left(\mathbf{x}\right)u\left(\mathbf{x}\right)d\mathbf{x},$$

where $\phi \left(\mathbf{x}\right)$ is a known weight or test function. If $\phi$ is a probability density function on $\Omega$, then $I_{\phi }$ can be interpreted as an expected value

$$I_{\phi }=\mathbb{E}\left[u\left(Y\right)\right],Y\sim \phi \left(\mathbf{x}\right),$$

and thus estimated using the same unbiased randomization principle as in the USA or NUSA algorithm. In this case, the only modification to the algorithm is that the initial point ${\mathbf{x}}_0$ is sampled from the density $\phi$.

When direct sampling from $\phi$ is difficult, or when $\phi$ is not itself a valid density function, an alternative approach is used. The algorithm draws ${\mathbf{x}}_0$ uniformly from $\Omega$ and multiplies the resulting Monte Carlo score by the bounded weight $\phi \left({\mathbf{x}}_0\right)$. This modification retains the unbiasedness of the estimator while accounting for nonuniform weighting in the functional. Such an approach allows direct computation of weighted mean values without additional structural changes to the underlying random walk process.

Algorithm 10: Weak Approximation via the Unbiased Stochastic Algorithm

2.6.2. Global Approximation

Another natural extension of the unbiased approach is the recovery of the full solution field $u\left(\mathbf{x}\right)$ across the domain $\Omega$ based on a single ensemble of random trajectories. During each random walk, multiple points of the domain are visited. The corresponding accumulated scores contain valuable statistical information that can be reused to reconstruct an approximation of $u\left(\mathbf{x}\right)$ over $\Omega$. This concept parallels the idea proposed in [7] for solving large linear systems, where visited states of a discrete random walk were treated as new starting points. The current context, however, involves a continuous domain and thus requires additional statistical smoothing of the visited data.

The procedure can be summarized as follows. Each trajectory starts from a uniformly distributed random initial point ${\mathbf{x}}_0$. Because the process is spatially uniform, the collection of all visited points across all walks can be considered as a uniform sampling of $\Omega$. For each visited point ${\mathbf{x}}_i$, an associated score $\mathrm{score}\left({\mathbf{x}}_i\right)$ is obtained, whose expected value equals $u\left({\mathbf{x}}_i\right)$. Hence, the mean score over points within a small neighborhood (hypercube) of $\mathbf{x}$ provides an unbiased local estimate of $u\left(\mathbf{x}\right)$. To improve the estimate, one can assume local smoothness of u and apply a second-order Taylor expansion:

$$u\left({\mathbf{x}}_i\right)=u\left(\mathbf{x}\right)+\nabla u\left(\mathbf{x}\right)·{\mathbf{h}}_i+\frac{1}{2}{\mathbf{h}}_i^T\left[D^{2}u\left(\mathbf{x}\right)\right]{\mathbf{h}}_i+O\left(h^3\right),$$

where ${\mathbf{h}}_i={\mathbf{x}}_i-\mathbf{x}$ and $\left[D^{2}u\left(\mathbf{x}\right)\right]$ is the Hessian matrix. Averaging the scores over the $N_h$ visited points within the hypercube of size $2h$ centered at $\mathbf{x}$ yields

$${\displaystyle \frac{1}{N_h}}\sum _{i=1}^{N_h}u\left({\mathbf{x}}_i\right)=u\left(\mathbf{x}\right)+O\left(h^2\right)$$

since the first-order term vanishes in expectation. Thus, $u\left(\mathbf{x}\right)$ can be approximated by the mean of the recorded scores within a small local neighborhood.

Algorithm 11: Global Approximation via Unbiased Monte Carlo

Both weak and global extensions of the unbiased Monte Carlo framework rely on the same fundamental random walk process that underpins the USA and NUSA algorithms. By reinterpreting the information collected during stochastic trajectories, these approaches allow efficient estimation of aggregate or spatially distributed quantities without additional model evaluations. The weak approximation reformulates the pointwise estimator into a functional form, enabling direct computation of weighted averages or mean responses with respect to arbitrary functions $\phi \left(\mathbf{x}\right)$. The global approximation, in turn, leverages all visited states and their associated scores to reconstruct a continuous representation of the solution across the domain. Because both extensions reuse the same ensemble of random trajectories, they maintain the unbiasedness and scalability of the original algorithm while dramatically improving data efficiency. Together, these methods extend the applicability of the unbiased stochastic framework to a wider class of problems where full-field or expectation-based quantities are of primary interest.

In the next section, we validate these algorithms through extensive numerical experiments and provide a quantitative comparison between biased and unbiased stochastic methods in both one- and multi-dimensional settings.

3. Numerical Simulations and Results

This section reports numerical experiments that validate and compare the biased and unbiased stochastic algorithms presented earlier. We consider both one-dimensional and multidimensional test problems to examine the behaviour of the estimators in terms of accuracy, variance, and computational cost. All implementations were executed in double precision with independent random or quasi-random sequences.

3.1. Benchmark Problem Description

To evaluate the algorithms under controlled conditions, we adopt the model problem originally introduced in [6,8]:

$$u\left(x\right)={\int }_{\Omega }k(x,x^{\prime })u\left(x^{\prime }\right)d{x}^{\prime }+f\left(x\right),\Omega =[0,1],$$

(30)

where

$$\begin{array}{cc}\hfill k(x,x^{\prime })& =x^2{e}^{x^{\prime }(x-1)},\hfill \end{array}$$

(31)

$$\begin{array}{cc}\hfill f\left(x\right)& =x+(1-x)e^x,\hfill \end{array}$$

(32)

$$\begin{array}{cc}\hfill \phi \left(x\right)& =\delta (x-x_0).\hfill \end{array}$$

(33)

The exact solution is $u\left(x\right)=e^x$, which enables the computation of the absolute and relative errors. The experiments aim to estimate $u\left(x_0\right)$ for several positions $x_0\in [0,1]$.

3.2. One-Dimensional Numerical Experiments

In this subsection, we present the one-dimensional numerical results for the test problem defined in Eq. (30). The computations focus on evaluating $u\left(x_0\right)$ at $x_0=0.5$ (and, for the unbiased case, also at $x_0=0.99$) using the biased and unbiased stochastic algorithms described earlier. We systematically compare the performance of the Crude Monte Carlo (MCM), quasi-Monte Carlo (MCA–MSS–1, MCA–MSS–2, MCA–MSS–2–S), and unbiased (USA) methods. All computations were performed using double-precision arithmetic and independent random or quasi-random sequences.

3.2.1. Performance of Crude Monte Carlo Sampling

Table 1 shows the convergence of the Crude Monte Carlo (MC) and Sobol quasi-Monte Carlo (qMC) methods applied directly to the integral equation. Both methods exhibit the expected monotonic decrease of the relative error with increasing N and number of series terms i. At small sample sizes ($N=128$), both achieve about 6% relative error, which drops below 0.1% for $N>{10}^4$. The Sobol sequence yields slightly smaller errors and much lower variance than pure random sampling, confirming the benefits of low-discrepancy sequences. The runtime grows sublinearly with N, remaining below 0.5 s even for the largest test.

3.2.2. Monte Carlo Integration of Liouville–Neumann Series Terms

Table 2 compares the Crude and Sobol methods for computing a finite number of Liouville–Neumann integrals. The relative error decays geometrically with the number of series terms i, demonstrating the convergence of the Neumann expansion. For $i=4$ and $N=65,000$, both implementations achieve sub-0.1% relative error. The Sobol sequence slightly outperforms random sampling, while the Finite Number of Integrals (FNI) residuals $R_N$ confirm that the stochastic integration error becomes negligible at large N. Computation time remains below one second, illustrating the efficiency of the algorithm even at high precision.

3.2.3. Quasi–Monte Carlo Algorithms MCA–MSS–1, MCA–MSS–2, and MCA–MSS–2–S

Table 3 and Table 4 consolidates the results for all quasi–Monte Carlo variants. The convergence pattern mirrors that observed for the basic Monte Carlo methods but with significantly reduced variance. MCA–MSS–1 and MCA–MSS–2 deliver almost identical results, confirming that both random shaking and symmetry correction achieve the same theoretical order of convergence. The MCA–MSS–2–S variant maintains accuracy within 5% of the other two methods while dramatically reducing computation time (approximately one-tenth of the MSS–1 runtime at $i=4$). The small FNI residuals $R_N$ confirm the stability of all three approaches even for large N and deep series truncations.

The results in Table 3 and Table 4 demonstrate that

• For all three algorithms, the relative error decreases rapidly as the number of integrals i increases from 1 to 4. At $i=4$, all methods achieve relative errors on the order of ${10}^{-3}$ or better, confirming convergence of the Liouville–Neumann expansion.
• Smaller ${\rho }_{\mathrm{m}}$ values correspond to tighter random perturbations around Sobol points, leading to reduced variance but increased computational time. As ${\rho }_{\mathrm{m}}$ decreases from $2\times {10}^{-3}$ to $4\times {10}^{-6}$, the error decreases by nearly two orders of magnitude.
• Both algorithms MCA–MSS–1 vs. MCA–MSS–2 exhibit nearly identical accuracy for each i and N, but MCA–MSS–2 provides marginally smoother error decay and slightly better stability, particularly at larger N. The introduction of symmetry in MCA–MSS–2 helps cancel first-order integration bias.
• The stratified symmetrized variant achieves similar accuracy (within $5\%$) while drastically reducing runtime. For instance, at $i=4$ and $N=65000$, MCA–MSS–2–S computes results in 6 s compared to 55 s for MCA–MSS–1/2, making it substantially more efficient.

Overall, the quasi-Monte Carlo algorithms clearly outperform classical random Monte Carlo techniques in terms of both precision and efficiency. Among them, MCA–MSS–2 offers the best balance of accuracy and stability, while MCA–MSS–2–S provides a computationally lighter alternative with nearly equivalent results. These findings underscore the effectiveness of quasi-Monte Carlo stratification and symmetrization for high-precision stochastic evaluation of Fredholm integral equations.

The quasi–Monte Carlo (MCA–MSS) results in Table 3 and Table 4 demonstrate that deterministic low-discrepancy sampling, particularly when combined with stratification and symmetrisation, can dramatically reduce the stochastic error in biased estimators derived from the Liouville–Neumann expansion. However, despite their high precision, all such biased algorithms remain fundamentally limited by the truncation error inherent in the finite-series representation. Even at large N and deep iteration levels ($i\ge 4$), the residual bias does not vanish completely, especially for strongly coupled or oscillatory kernels.

To overcome this intrinsic limitation, we now turn to the class of unbiased stochastic algorithms—namely the USA and its refined version NUSA—where the expectation of the random estimator equals the exact solution of the Fredholm integral equation without series truncation. The following subsection presents their formulation and numerical performance, highlighting how these algorithms remove systematic bias while retaining the statistical efficiency achieved by the best quasi–Monte Carlo methods.

3.2.4. Unbiased USA and NUSA Algorithms

The results presented in Table 5 provide a detailed comparison between the classical unbiased USA algorithm and the proposed Novel Unbiased Stochastic Algorithm (NUSA) for estimating $u\left(x_0\right)$ at different spatial points and sampling levels N. Both algorithms demonstrate systematic error reduction as the number of trajectories increases, confirming the expected convergence of the unbiased stochastic estimators. However, the NUSA method consistently achieves smaller relative errors than USA across all N values and evaluation points, with an improvement factor typically between $1.3$ and $1.5$. This advantage becomes especially evident for small and moderate sample sizes, where NUSA’s adaptive variance-control mechanism suppresses estimator dispersion more effectively.

The convergence behaviour shows that the relative error for both algorithms decreases approximately as $\mathcal{O}\left(N^{-1/2}\right)$, characteristic of Monte Carlo estimators, while the computational time scales linearly with the number of simulated trajectories. Although NUSA requires slightly higher computational time—about $1.5$ times that of USA—its accuracy gains are substantial, especially at larger N, where it attains sub-${10}^{-4}$ precision. For example, at $x_0=0.5$ and $N=4,194,304$, the USA method yields a relative error of $3\times {10}^{-5}$, whereas NUSA reduces it to $2.2\times {10}^{-5}$, maintaining numerical stability throughout the simulation range. Similar improvements are observed at $x_0=0.99$, confirming that NUSA preserves its unbiasedness and accuracy even near boundary regions where variance tends to increase.

The results in Table 6 provide a clear comparison between biased and unbiased Monte Carlo algorithms for solving integral equations. The biased approaches—including the Crude and Sobol-sequence implementations of the finite-integral (FNI) method and the standard Markov Chain Monte Carlo (MCM)—exhibit consistent convergence as the number of trajectories N increases, but their accuracy plateaus once the truncation of the Liouville–Neumann series dominates the total error. In contrast, the unbiased algorithms (USA and NUSA) continue to improve steadily with increasing N, since they do not rely on a finite-series approximation.

Between the two unbiased schemes, NUSA consistently achieves smaller relative errors across all sampling levels, typically by a factor of $1.3$–$1.5$, while maintaining similar average numbers of moves before absorption. This improvement can be attributed to NUSA’s adaptive variance-control strategy and dynamic weighting, which reduce estimator dispersion without affecting the unbiased expectation. For example, at $N=65,536$, the USA method attains a relative error of $7.5\times {10}^{-4}$, whereas NUSA lowers it to $5.2\times {10}^{-4}$ with negligible additional cost. Moreover, NUSA exhibits smoother convergence and reduced fluctuation in relative errors, indicating enhanced numerical stability and robustness.

3.3. Systematic Error Analysis

To better understand the magnitude and behavior of deterministic bias in the iterative estimation of $u\left(x_0\right)$, we evaluated the systematic component of the total error, denoted by $R_{\mathrm{sys}}$. Table 7 reports the computed values of $R_{\mathrm{sys}}$ at two representative points, $x_0=0.1$ and $x_0=0.5$, for three iteration levels $i=1,2,3$. These results illustrate the typical order of magnitude of this component under the conditions used in our Monte Carlo simulations.

The results confirm that the systematic error decreases as the number of iterations i increases, indicating the expected convergence behavior of the iterative Monte Carlo scheme. At small iteration counts, however, the systematic component remains dominant in most cases, particularly for points closer to the center of the domain (e.g., $x_0=0.5$), where the kernel’s influence and accumulated truncation effects are more pronounced. This behavior arises because the balancing analysis between stochastic and systematic errors is typically derived under the assumption that the Crude Monte Carlo (CMCM) method is used, which tends to exhibit slower reduction of systematic bias.

When the results in Table 7 are compared with those in Table 6, a consistent trend emerges. The unbiased stochastic algorithms, particularly USA and its improved variant NUSA, achieve markedly lower relative errors for the same number of iterations while maintaining lower computational complexity. Only in exceptional cases—such as integral equations with a very small spectral radius of the kernel operator— might the biased methods exhibit comparable performance. In general, however, the USA and especially the NUSA methods are clearly preferable for applications requiring higher accuracy and stable convergence with minimal systematic bias.

Across all 1D tests, several clear patterns emerge:

• Convergence and accuracy. All stochastic methods demonstrate monotonic error reduction with increasing N and series depth i. Quasi–Monte Carlo methods, particularly MCA–MSS–2 and MCA–MSS–2–S, achieve the same accuracy as the unbiased USA with significantly fewer samples.
• Computational efficiency. Crude Monte Carlo remains the simplest but least efficient approach. Sobol-based and MSS methods provide up to two orders of magnitude improvement in accuracy at comparable runtime.
• Unbiased estimation. The USA eliminates truncation bias completely and maintains accuracy across different spatial points. Its computational cost is slightly higher than that of the MCA–MSS–2–S algorithm but yields theoretically exact expectations.
• Variance control. The symmetrization (MSS–2) and stratification (MSS–2–S) mechanisms effectively suppress variance, while the unbiased USA controls it via absorption probabilities.

In summary, for one-dimensional integral equations, quasi–Monte Carlo algorithms (especially MCA–MSS–2–S) achieve the best trade-off between accuracy and speed among the biased methods, whereas the unbiased USA and NUSA methods provides exact results with slightly higher computational effort. This justifies extending the comparison to higher-dimensional test cases, discussed in the next subsection.

Overall, these results highlight the enhanced reliability of NUSA as an unbiased stochastic solver. Its ability to achieve lower variance and smoother convergence without introducing bias demonstrates its superiority for practical applications where both accuracy and robustness are required. The modest computational overhead compared to USA makes NUSA a highly competitive method, particularly in problems involving fine discretization or strongly coupled kernels, where conventional Monte Carlo estimators typically suffer from higher variance.

3.4. Comparison of All Algorithms

Table 8 compares the Crude Monte Carlo Method (CMCM), Markov Chain Monte Carlo (MCM), three quasi-Monte Carlo variants (MCA-MSS-1, MCA-MSS-2, MCA-MSS-2-S), and the two unbiased approaches (USA and NUSA). Each algorithm used $N=1000$ trajectories or quasi-random points. All results were averaged over ten runs.

The following trends can be observed:

• The CMCM and MCM algorithms yield consistent but relatively high errors due to their purely random sampling and bias from truncation.
• The quasi-Monte Carlo methods MCA-MSS-1, MCA-MSS-2, and MCA-MSS-2-S achieve one to two orders of magnitude higher accuracy, owing to their low-discrepancy Sobol sampling. Among these, MCA-MSS-2 delivers the smallest errors due to its symmetrisation procedure.
• The classical unbiased method (USA) eliminates systematic bias but exhibits significant variance, especially near the boundary ($x_0\approx 0.1$).
• The proposed NUSA algorithm clearly provides the best overall performance: its errors are smaller than those of USA and substantially lower than those of any biased method. The adaptive absorption parameter $P_d$ effectively controls variance without compromising unbiasedness.

3.5. Comparison with Deterministic Methods

Many deterministic numerical schemes have been developed to solve Fredholm integral equations of the second kind. Among them, the Nystrom method based on the Simpson quadrature rule remains one of the most widely used due to its simplicity and effectiveness. Detailed descriptions and performance analyses of such deterministic approaches can be found in [2]. To evaluate the relative performance of stochastic algorithms against deterministic solvers, we compare the Simpson-based Nystrom method with both the classical Unbiased Stochastic Algorithm (USA) and the proposed Novel Unbiased Stochastic Algorithm (NUSA) on a set of one-dimensional test cases involving kernels and solutions of varying smoothness.

Applying the Simpson quadrature rule to the integral term of the Fredholm equation yields

$$u\left(x\right)\simeq \sum _{k=1}^n{\alpha }_{k}k(x,x_k)u\left(x_k\right)+g\left(x\right),$$

where ${\alpha }_k$ and $x_k$ denote, respectively, the weights and nodes of the quadrature. Evaluating the equation at the nodes $x_i$ results in a linear system

$$u_i=\sum _{k=1}^n{\alpha }_{k}k(x_i,x_k)u_k+g\left(x_i\right),$$

whose numerical solution approximates the continuous one through interpolation. The computational complexity of this procedure scales as $\mathcal{O}\left(n^2\right)$ due to the matrix inversion, so for a fair comparison the stochastic simulations are allocated a comparable number of samples.

The tests include both regular and discontinuous kernels and solutions:

• RKRS: Regular Kernel, Regular Solution;
• RKDS: Regular Kernel, Discontinuous Solution;
• DKRS: Discontinuous Kernel, Regular Solution;
• DKDS: Discontinuous Kernel, Discontinuous Solution.

For all experiments, the point of evaluation is $x=0.5$ with the exact value $u\left(x\right)=e^{0.5}$. The parameters $p_1=50$, $p_2=400$, and $p_3=1000$ determine the discretization level of the Simpson rule, leading to ${(2p_i+1)}^2$ stochastic simulations for comparable computational effort.

The comparison in Table 9 demonstrates the distinct performance profiles of the deterministic and stochastic algorithms. For regular kernels and smooth solutions (RKRS), the Nystrom method remains the most accurate, as expected from its reliance on high-order quadrature. However, in cases involving discontinuities in either the kernel or the solution (RKDS, DKRS, and DKDS), the performance of deterministic quadrature degrades, while the stochastic estimators remain stable. The proposed NUSA algorithm achieves up to an order of magnitude improvement in accuracy over the standard USA method, and in several discontinuous cases approaches or even surpasses the precision of the Nystrom scheme. This confirms that NUSA effectively bridges deterministic and stochastic formulations, combining the robustness of Monte Carlo sampling with accuracy levels previously attainable only by quadrature-based methods.

3.6. General USA and NUSA Algorithms

To validate the performance and generality of the unbiased algorithms, we consider a second numerical example where the kernel $k(x,x^{\prime })$ can assume both positive and negative values, and its magnitude may exceed one in absolute value:

$$u\left(x\right)={\int }_0^{1}k(x,x^{\prime })u\left(x^{\prime }\right)d{x}^{\prime }+f\left(x\right),$$

(34)

with

$$\begin{array}{cc}\hfill k(x,x^{\prime })& =\beta (2x^{\prime }-1)x^2{e}^{x^{\prime }(x-1)},\hfill \end{array}$$

(35)

$$\begin{array}{cc}\hfill f\left(x\right)& =e^x+\beta \left(2e^x-2-x{e}^x-x\right),\hfill \end{array}$$

(36)

$$\begin{array}{cc}\hfill \phi \left(x\right)& =\delta (x-x_0),\hfill \end{array}$$

(37)

where $x_0\in \{0.1,0.5,0.99\}$ is used in separate experiments. The analytical solution remains $u\left(x\right)=e^x$. The parameter $\beta$ is introduced to explore increasingly challenging configurations: as $\beta$ increases, the kernel norm approaches one, leading to numerical instability and rapidly growing variance.

Following the approach described in Section 2.6, we also compute the following weak integrals to illustrate the versatility of the unbiased framework:

$$I_1={\int }_0^{1}u\left(x\right)dx,I_2={\int }_0^{1}2xu\left(x\right)dx.$$

For $I_1$, the starting point $x_0$ is drawn uniformly from $[0,1]$. For $I_2$, two sampling strategies are considered: (i) sampling $x_0$ from the density $2x$ (which is equivalent to setting $x_0=\sqrt{U}$ for $U\sim \mathcal{U}[0,1]$); (ii) sampling $x_0$ uniformly and multiplying the obtained score by $2x_0$.

The numerical results summarized in Table 10 and Table 11 confirm the robustness and accuracy of both unbiased Monte Carlo algorithms under increasingly challenging kernel conditions. For moderate kernel strengths ($\beta =1$ and $\beta =2$), both the USA and NUSA methods yield excellent agreement with the analytical solution $u\left(x\right)=e^x$. The relative errors remain below ${10}^{-3}$, while the mean number of absorption steps i stays close to unity, reflecting the strong numerical stability of the unbiased formulation. At these parameter values, the NUSA algorithm consistently provides 30–50% smaller relative errors and noticeably reduced standard deviations compared with USA, while maintaining nearly identical computational complexity. This confirms that the adaptive variance-control feature embedded in NUSA effectively enhances estimator precision without altering the unbiased property of the solution.

As the kernel becomes more ill-conditioned at $\beta =4.15$, the variance of the stochastic estimator increases sharply, consistent with theoretical expectations for integral operators approaching the critical spectral radius. Nevertheless, the NUSA algorithm demonstrates superior numerical behavior in this near-singular regime. Both the relative errors and standard deviations remain significantly lower than those of USA, particularly for the most difficult case $x_0=0.99$, where the kernel influence is strongest. The mean number of absorption steps increases moderately, indicating that trajectories become longer as the process approaches instability, yet the convergence remains stable and unbiased.

The additional tests for the weighted functionals $I_1$ and $I_2$ further confirm the flexibility of the unbiased framework. Between the two sampling approaches considered for $I_2$, the method that starts from a uniformly distributed $x_0$ and weights the score by $2x_0$ (method (ii)) yields more accurate and stable results than direct sampling from the density $2x$ (method (i)). This observation supports the earlier conclusion that uniform initialization combined with adaptive score weighting provides a better balance between bias control and variance minimization.

Overall, these results confirm the theoretical expectations regarding stability and variance behavior and show that the NUSA algorithm systematically improves upon the standard USA, delivering higher accuracy and smoother convergence under both regular and near-critical kernel conditions.

3.6.1. Global Approximation

In this final numerical experiment, we assess the performance of the unbiased algorithms when applied to the global approximation procedure described in SubSection 2.6. Here, the objective is to approximate the solution $u\left(x\right)$ at all points of the domain by exploiting the statistical information gathered during the random walks. This experiment also demonstrates how the accuracy evolves with increasing numbers of trajectories N.

The computations were performed at two reference points, $x_1=0.5$ and $x_2=0.95$, for two different local neighborhood sizes used in the Taylor expansion, $h_1=0.1$ and $h_2=0.05$. The results for the relative error, standard deviation $\sigma$, and mean number of absorption steps i are reported in Table 12 and Table 13 for the USA and NUSA algorithms, respectively. The simulations were conducted with $N={10}^4$, ${10}^5$, ${10}^6$, and ${10}^7$ trajectories.

The results confirm that the relative error decreases monotonically as the number of trajectories N increases, in agreement with the expected Monte Carlo convergence rate. At both reference points $x_1$ and $x_2$, the NUSA algorithm consistently achieves smaller errors and standard deviations compared with USA—approximately 30–40% improvement on average—while maintaining the same mean number of absorption steps. This improvement demonstrates that the variance-reduction mechanism in NUSA remains effective even when the solution is reconstructed across the entire domain.

When the Taylor expansion is employed for local interpolation, the results remain accurate but slightly less precise than those from direct estimation. As anticipated, the parameter h controls the balance between local smoothing and variance: for larger h values ($h_1=0.1$), the averaging region becomes too broad, leading to an oversmoothing effect where the statistical information is overly integrated and the accuracy no longer improves with increasing N. Conversely, for smaller h ($h_2=0.05$), the local estimate initially suffers from a limited number of unbiased samples, but the accuracy improves significantly as N increases and the neighborhood becomes better populated.

Comparing the two algorithms, NUSA consistently produces smoother and more stable convergence curves, confirming its improved statistical efficiency. The unbiasedness of both methods is preserved throughout the computations, but NUSA’s adaptive weighting clearly mitigates stochastic variability, especially at larger N where fine-resolution averaging is required. These findings demonstrate that the proposed unbiased framework can successfully reconstruct global approximations of the solution while maintaining high numerical accuracy and robustness.

The following subsection examines the behavior of USA and NUSA in higher-dimensional domains, evaluating their scalability, numerical stability, and accuracy when applied to integral equations with multidimensional kernels. The results obtained for the global reconstruction of the one-dimensional solution confirm that both unbiased algorithms, USA and NUSA, can efficiently exploit trajectory information to recover the solution over the entire domain with high accuracy. The numerical experiments conducted also for the one-dimensional integral equations demonstrate that both the USA and NUSA algorithms deliver accurate and stable solutions over a wide range of kernel conditions, including regimes approaching the spectral instability threshold. The convergence behavior observed with increasing N and the stability of the variance reduction in NUSA indicate that the proposed approach scales favorably with problem size. Having verified the effectiveness of the proposed methods in this controlled setting, it is natural to extend the investigation to multidimensional problems, where the computational complexity and variance amplification become more pronounced. These encouraging findings motivate an extension of the analysis to multidimensional Fredholm integral equations, where the computational complexity and variance amplification are substantially more pronounced. In the following subsection, we investigate the performance of USA and NUSA in multidimensional settings, evaluating their numerical scalability, efficiency, and accuracy as the dimension and kernel complexity increase.

3.7. Multidimensional Problem

Building upon the promising one-dimensional results, we now extend the analysis of the unbiased algorithms to multidimensional Fredholm integral equations of the second kind. In higher dimensions, the numerical challenges are considerably greater due to the rapid growth of computational cost, the potential amplification of stochastic variance, and the increased difficulty of maintaining numerical stability. The goal of these experiments is to evaluate the scalability, accuracy, and robustness of both USA and NUSA algorithms when applied to multidimensional kernels with varying coupling strengths.

For this purpose, we consider test problems defined on the unit hypercube $\Omega ={[0,1]}^n$ with kernel norms scaled by a coupling parameter $\beta$. The values $\beta =0.1/{\parallel K\parallel }_{L_2}$, $\beta =0.5/{\parallel K\parallel }_{L_2}$, and $\beta =0.95/{\parallel K\parallel }_{L_2}$ are selected to represent weak, moderate, and strong interaction regimes, respectively. The solution is computed at the reference point ${\mathbf{x}}_0=(0.75,0.75,\dots ,0.75)$ for problem sizes $n=3,5,10,25,75$, and 100. For each configuration, we record the relative error, computational time, and the total number of random trajectories required for convergence.

The following tables and discussion highlight the comparative performance of the USA and NUSA algorithms across increasing dimensionality. Special attention is given to how the adaptive variance-control mechanism of NUSA influences convergence and stability in high-dimensional spaces, where unbiased stochastic methods typically face the most severe numerical difficulties.

The extension to several dimensions is considerably more demanding, since the Fredholm integral equation of the second kind becomes progressively harder to treat as n grows. Following the formulation in [6], we consider

$${\displaystyle u_n\left(\mathrm{x}\right)={\int }_0^1\dots {\int }_0^1\prod _{j=1}^n\left\{k(x^{\left(j\right)},x^{{\left(j\right)}^{\prime }})\right\}u\left({\mathrm{x}}^{\prime }\right)\prod _{j=1}^n\mathrm{d}{x}^{{\left(j\right)}^{\prime }}+f_n\left(\mathrm{x}\right),}$$

(38)

with $\mathrm{x}\equiv (x^{\left(1\right)},\dots ,x^{\left(n\right)}),{\mathrm{x}}^{\prime }\equiv (x^{{\left(1\right)}^{\prime }},\dots ,x^{{\left(n\right)}^{\prime }}),{\prod }_{j=1}^n\mathrm{d}{x}^{{\left(j\right)}^{\prime }}=\mathrm{d}{\mathrm{x}}^{\prime }.$ In this setting the ingredients are defined as

$$k(x^{\left(j\right)},x^{{\left(j\right)}^{\prime }})={\beta }_j\left(2x^{{\left(j\right)}^{\prime }}-1\right){\left\{x^{\left(j\right)}\right\}}^{2}exp\left\{x^{{\left(j\right)}^{\prime }}(x^{\left(j\right)}-1)\right\},\beta =\prod _{j=1}^n{\beta }_j,$$

(39)

$$f_n\left(x\right)={\mathrm{e}}^{{\sum }_{j=1}^n{x}^{\left(j\right)}}+{(-1)}^{(j-1)}{\beta }_n(2{\mathrm{e}}^{x^{\left(j\right)}}-2-x^{\left(j\right)}{\mathrm{e}}^{x^{\left(j\right)}}-x^{\left(j\right)}),$$

(40)

$$\phi \left(\mathrm{x}\right)=\delta (\mathrm{x}-{\mathrm{x}}_0),$$

(41)

where the reference point chosen as ${\mathrm{x}}_0=(0.75,0.75,\dots ,0.75)$ and ${\mathrm{x}}_0=(1.0,1.0,\dots ,1.0)$ for the numerical experiments.

The solution of the equation can be written as the sum of the continuous component $u\left(\mathrm{x}\right)$ and a singular part involving the Dirac delta function $\delta (\mathrm{x}-{\mathrm{x}}_0)$ centered at ${\mathrm{x}}_0$. In particular, the value of $u\left(\mathrm{x}\right)$ at ${\mathrm{x}}_0$ is obtained through the convolution of the forcing term $f\left(\mathrm{x}\right)$ with the corresponding Green function. As the dimension increases, the evaluation procedure becomes considerably more demanding, reflecting the intrinsic complexity of the problem.

Table 14 presents the detailed comparison of the classical Unbiased Stochastic Algorithm (USA) and the proposed Novel Unbiased Stochastic Algorithm (NUSA) across various dimensions and kernel strengths. The results confirm that both algorithms exhibit nearly linear growth of computational time with respect to problem dimension n, indicating that the stochastic framework is effectively free from the curse of dimensionality. This linear scaling trend is particularly evident for USA, where CPU times increase from approximately 12 s at $n=3$ to about 430 s at $n=100$, and similarly for NUSA, which remains consistently about twice as slow due to its adaptive weighting and variance-control procedures. Despite the moderate increase in runtime, NUSA achieves a substantial improvement in accuracy across all kernel regimes. For weakly coupled kernels ($\beta =0.1/{\parallel K\parallel }_{L_2}$), relative errors remain below ${10}^{-3}$ even in high dimensions, demonstrating the precision of the unbiased variance-controlled formulation. As the kernel coupling increases ($\beta =0.5$ and $0.95$), both algorithms show higher variance, but the deterioration is far less pronounced for NUSA. For example, at $\beta =0.95/{\parallel K\parallel }_{L_2}$ and $n=100$, the USA method attains a relative error of order ${10}^{-3}$, while NUSA maintains one an order of magnitude smaller, confirming its robustness under strong coupling.

The numerical results presented in Table 15 highlight the effect of dimensionality and kernel strength on the performance of the USA and NUSA algorithms. For all tested configurations, both methods reproduce the analytical solution with high precision, but the NUSA algorithm consistently achieves smaller relative errors across all problem sizes and coupling regimes. The accuracy advantage of NUSA is most evident in the moderate and strong coupling cases, where variance amplification typically degrades the performance of standard unbiased estimators. Even when the kernel norm approaches unity, NUSA maintains stable convergence and reduced stochastic dispersion, confirming the effectiveness of its adaptive variance-control mechanism. As expected, the computational time increases approximately linearly with dimension for both algorithms. The NUSA algorithm is designed to incorporate additional weighting and resampling operations, which roughly double the overall runtime compared to USA. However, this modest overhead results in substantial accuracy gains, particularly in high-dimensional cases ($n=10$) and for stronger interactions ($\beta =0.95/{\parallel K\parallel }_{L_2}$), where standard unbiased methods begin to exhibit numerical instability. The near-linear scaling of computational time with dimension demonstrates that both algorithms remain computationally feasible even for moderately high-dimensional problems.

More precisely:

• Both algorithms show the expected increase in relative error as the kernel norm approaches unity ($\beta =0.95/{\parallel K\parallel }_{L_2}$), indicating higher numerical difficulty in solving strongly coupled Fredholm equations. However, while the USA method exhibits noticeable accuracy degradation for large ${\parallel K\parallel }_{L_2}$ (up to $0.3\%$ error at $n=3$ and above $2.5\%$ for $n=75$), the NUSA algorithm maintains uniformly small errors, typically below $0.02\%$ across all tested regimes.
• The performance of the USA method remains satisfactory for low- and moderate-dimensional problems ($n\le 10$), but its accuracy declines quickly as dimensionality increases. In contrast, NUSA exhibits almost dimension-independent behaviour, maintaining sub-${10}^{-2}$ relative errors even for $n=100$, thus demonstrating excellent scalability.
• The computational time of NUSA is roughly two times larger than that of USA because of its adaptive absorption and weighted sampling mechanisms. Despite this overhead, the substantially improved accuracy results in a smaller global error per unit of CPU time, making NUSA more efficient for realistic high-dimensional applications.
• The USA algorithm becomes unstable for near-singular kernels (${\parallel K\parallel }_{L_2}\to 1$), where the variance of the estimator increases sharply. NUSA successfully mitigates this instability through controlled absorption, maintaining smooth convergence and numerical robustness for all $\beta$ values.
• For weak or moderately coupled kernels ($\beta \le 0.5$), USA achieves reasonable accuracy at minimal computational cost. However, for strongly coupled or high-dimensional systems, NUSA clearly outperforms USA, providing one to two orders of magnitude smaller errors with only modest increases in runtime.

3.8. Summary of Findings

The comparison confirms that the NUSA method offers the most balanced performance among all unbiased stochastic approaches. Its principal advantages are:

• Uniform accuracy below ${10}^{-3}$ across all tested dimensions and kernel strengths.
• Linear scaling of computational time with dimension n, preserving numerical stability.
• Significant variance reduction compared to the classical USA, especially for large ${\parallel K\parallel }_{L_2}$ values.

Across all experiments, several conclusions can be drawn:

• Accuracy hierarchy: NUSA > USA > MCA-MSS-2 > MCA-MSS-2-S > MCA-MSS-1 > MCM > CMCM.
• Variance control: NUSA achieves the lowest estimator variance through its adaptive absorption mechanism.
• Bias behaviour: All quasi-Monte Carlo methods exhibit small but non-negligible bias due to truncation, whereas USA and NUSA remain exactly unbiased.
• Computational efficiency: Although NUSA requires slightly more operations per trajectory, its reduced variance allows it to reach a target accuracy with fewer samples, making it the most efficient overall.

The numerical analysis demonstrates the following key points:

• The Crude and Markov Chain Monte Carlo methods provide consistent but relatively low-precision estimates.
• Quasi-Monte Carlo algorithms based on modified Sobol sequences (MCA-MSS-1, 2, 2-S) improve accuracy dramatically while retaining low computational cost.
• The unbiased algorithms (USA and NUSA) eliminate systematic bias completely, with NUSA offering the best balance between variance, stability, and efficiency.

The next section presents concluding remarks, summarising the theoretical and numerical insights and outlining directions for future research.

Overall, the NUSA algorithm provides a robust, scalable, and fully unbiased framework for solving multidimensional Fredholm integral equations of the second kind. Its combination of high accuracy and numerical stability makes it a strong candidate for deployment in complex high-dimensional physical and engineering models.

The multidimensional experiments confirm the theoretical expectations established in the earlier sections. Both biased and unbiased stochastic algorithms exhibit consistent convergence with increasing sample size, but their performance diverges significantly as the kernel strength and dimensionality increase. The quasi–Monte Carlo schemes (MCA–MSS–2 and MCA–MSS–2–S) remain highly efficient for smooth and weakly coupled kernels, while the unbiased approaches, particularly the proposed NUSA method, retain their accuracy and robustness even for near-singular and high-dimensional problems.

The overall results highlight the computational trade-off between the two algorithms: USA provides faster estimates with acceptable precision for weak or moderate kernels, while NUSA offers far greater stability and uniform accuracy across dimensions at roughly twice the computational cost. Importantly, both methods maintain near-linear runtime scaling, showing that the stochastic approach remains efficient even in high-dimensional integral equation problems.

The multidimensional analysis confirms that the proposed NUSA algorithm successfully combines high accuracy, unbiasedness, and computational scalability within a unified stochastic framework. Although approximately twice as computationally demanding as the classical USA method, NUSA consistently achieves one to two orders of magnitude smaller relative errors and remains stable under strong kernel coupling and increasing dimensionality. The linear scaling of CPU time with respect to dimension further demonstrates that the proposed variance-controlled stochastic strategy can be effectively applied to large-scale problems without exponential growth in computational cost. These findings validate NUSA as a robust and efficient approach for solving high-dimensional Fredholm integral equations of the second kind, providing a solid foundation for further extensions and hybrid variance-reduction schemes.

4. Discussion

The numerical experiments conducted in both one- and multi-dimensional settings provide strong evidence supporting the theoretical foundations of the proposed stochastic framework. The comparative analysis among deterministic, biased stochastic, and unbiased stochastic algorithms demonstrates that while classical methods remain highly effective for smooth low-dimensional problems, their efficiency and accuracy degrade rapidly as kernel complexity and dimensionality increase. In contrast, the Novel Unbiased Stochastic Algorithm (NUSA) maintains stable convergence, low variance, and near-constant accuracy across all tested regimes, validating the robustness of its adaptive variance-control mechanism and unbiased formulation.

The purpose of this discussion is to synthesise the theoretical and computational findings, highlight the main insights obtained from the comparative analysis, acknowledge the remaining limitations of the approach, and outline promising avenues for future research and methodological development.

The multidimensional results clearly demonstrate that the proposed NUSA algorithm unifies high accuracy, numerical stability, and scalability within a single unbiased stochastic framework. Its performance advantage over the classical USA becomes more pronounced with increasing kernel strength and dimensionality, confirming that adaptive variance control is crucial for mitigating estimator dispersion in complex integral systems. Moreover, the near-linear computational scaling observed for both algorithms highlights the suitability of stochastic approaches for problems where deterministic quadrature-based solvers suffer from exponential complexity growth. These observations motivate a broader reflection on the methodological implications, practical limitations, and prospective developments of the stochastic framework, which are discussed below.

4.1. Key Insights

The results obtained across all experiments provide compelling evidence for the accuracy, stability, and scalability of the unbiased stochastic algorithms proposed in this study. Both USA and NUSA consistently converge toward the analytical solution with decreasing variance as the number of trajectories increases, confirming the theoretical unbiasedness of the estimators. In all tested configurations—from simple one-dimensional kernels to multidimensional cases—the relative errors exhibit the expected Monte Carlo rate of convergence, while the computational cost scales linearly with the number of samples.

A key observation is the superior variance control achieved by the NUSA algorithm. Through adaptive weighting and a more efficient sampling strategy, NUSA systematically reduces stochastic dispersion without introducing bias. This improvement is particularly pronounced in challenging conditions: near the spectral stability threshold ($\beta \to 1$) and in higher-dimensional domains, where conventional unbiased estimators often suffer from unstable variance growth. The results show that NUSA achieves 30–50% lower errors than USA on average, confirming that its variance-reduction mechanism provides a substantial accuracy gain at a moderate computational cost. The ability to extend both algorithms to weak and global approximation regimes further highlights their flexibility and broad applicability to integral equations with complex structure.

4.2. Limitations

Despite these advantages, several limitations should be acknowledged. First, the computational overhead introduced by the NUSA algorithm—approximately twice that of the standard USA—may become significant in extremely high-dimensional problems or when applied to real-time simulations. While this overhead is justified by the accuracy gains, future implementations could benefit from parallelization and variance reuse strategies to reduce computational burden. Second, both methods rely on sufficient smoothness of the kernel and the solution when employing Taylor-based local averaging in the global approximation. For discontinuous or highly oscillatory kernels, the convergence rate may degrade unless suitable local adaptivity or kernel regularization is applied. Finally, although unbiasedness ensures the absence of systematic bias, achieving very low variance in near-singular configurations still requires large trajectory counts, which may limit efficiency for ill-conditioned systems.

4.3. Future Directions

Several promising research avenues arise from this work. First, integrating quasi–Monte Carlo techniques or correlated sampling schemes within the NUSA framework could further enhance convergence rates, especially for high-dimensional integrals. Second, adaptive selection of the random walk termination criteria and local averaging parameter h could improve efficiency by automatically balancing stochastic and systematic errors. A third direction involves extending NUSA to nonlinear or time-dependent integral equations, where unbiased estimation remains a critical challenge. Finally, combining the unbiased stochastic approach with deterministic preconditioners or spectral methods may yield hybrid algorithms capable of tackling large-scale problems with both accuracy and computational efficiency.

4.4. Concluding Remarks

Overall, the unbiased stochastic framework developed here, and particularly its enhanced NUSA variant, provides a robust and general methodology for solving Fredholm integral equations of the second kind. It effectively bridges the gap between traditional deterministic solvers and conventional biased Monte Carlo approaches by delivering unbiasedness, high accuracy, and scalability within a unified framework. The results confirm that unbiased algorithms can serve as a viable and efficient alternative for multidimensional integral equations, laying a solid foundation for future developments in stochastic numerical analysis and computational mathematics.

In summary, the NUSA algorithm unifies the benefits of unbiased estimation and variance-reduced sampling, offering a numerically stable and scalable solution technique for high-dimensional integral equations. Its demonstrated robustness and accuracy make it a versatile tool for stochastic numerical analysis and a foundation for future developments in probabilistic computational methods.

The discussion presented above consolidates both the theoretical and numerical perspectives of the proposed stochastic framework. Through extensive comparative analysis, it has been shown that the NUSA algorithm achieves a unique balance between accuracy, scalability, and unbiasedness, outperforming classical Monte Carlo and deterministic methods in complex and high-dimensional settings. The empirical evidence validates the theoretical foundations of the algorithm, demonstrating its capacity to maintain low variance and high precision across a broad range of kernel strengths and problem dimensions. The following section synthesises these findings into a concise set of conclusions, emphasising the methodological contributions, computational insights, and the most promising directions for future research and application.

The collective results from one-dimensional, global, and multidimensional experiments clearly demonstrate the robustness and generality of the unbiased stochastic framework. Both USA and NUSA algorithms successfully maintain unbiasedness while providing accurate approximations across a wide range of kernel strengths and dimensional settings. In higher dimensions, where deterministic methods become prohibitively expensive and biased stochastic estimators suffer from variance explosion, the proposed unbiased approaches continue to perform reliably. The NUSA algorithm, in particular, consistently delivers lower relative errors and smoother convergence profiles while preserving near-linear computational scaling with respect to dimension. Although it requires approximately twice the runtime of USA, this additional cost is offset by a marked improvement in numerical stability and precision, especially near the spectral limit of the kernel. Taken together, these findings confirm that the NUSA method represents a significant advancement in the class of unbiased Monte Carlo algorithms, providing a practical, scalable, and highly accurate solution strategy for multidimensional Fredholm integral equations.

5. Conclusions and Future Work

This paper has presented a unified theoretical and computational framework for solving Fredholm integral equations of the second kind using stochastic algorithms. This study also presented a comprehensive analysis of unbiased stochastic algorithms for solving Fredholm integral equations of the second kind, focusing on the general Unbiased Stochastic Algorithm (USA) and its enhanced variant, the Novel Unbiased Stochastic Algorithm (NUSA). Through extensive numerical experiments—including one-dimensional, global, and multidimensional test cases—the proposed methods were shown to provide accurate, stable, and scalable estimations of the solution without introducing systematic bias.

Summary of Theoretical Contributions

The work extends stochastic solution theory through the following key advances:

• Development of a probabilistic representation of the Fredholm equation based on random trajectories and absorption processes.
• Formulation of unbiased stochastic estimators (USA and NUSA) that preserve exact expectation equivalence with the analytical solution, eliminating truncation bias.
• Comparative analysis of biased algorithms—Crude Monte Carlo, Markov Chain Monte Carlo (MCM), and quasi–Monte Carlo (MCA–MSS–1, 2, and 2–S)—under a common Liouville–Neumann framework.
• Introduction of the Novel Unbiased Stochastic Algorithm (NUSA), incorporating adaptive absorption control and kernel-weight normalization to achieve substantial variance reduction while maintaining unbiasedness.

Summary of Numerical Findings

The numerical experiments validate the theoretical framework across both one-dimensional and multidimensional test cases:

• In one-dimensional problems, quasi–Monte Carlo methods achieve errors of order ${10}^{-3}$ at minimal computational cost, while the unbiased USA method attains similar accuracy when variance is properly controlled.
• The proposed NUSA algorithm consistently delivers superior precision, achieving uniform sub-${10}^{-3}$ errors across all tested domains and removing truncation bias entirely.
• In multidimensional settings ($n=3$–100), NUSA maintains stable accuracy even for strongly coupled kernels (${\parallel K\parallel }_{L_2}\approx 1$), where other algorithms exhibit pronounced variance growth.
• Although computationally more intensive than USA, NUSA’s lower variance yields the smallest overall error-to-cost ratio.

Computational and Methodological Implications

The comparative evaluation leads to several general insights:

• Low-discrepancy quasi–Monte Carlo sampling substantially accelerates convergence and should be preferred for smooth kernels and moderate dimensions.
• The unbiased USA and NUSA algorithms are indispensable for strongly coupled or high-dimensional equations, where truncation bias dominates.
• Hybrid strategies that combine quasi–Monte Carlo initialization with unbiased stochastic continuation represent a promising path toward optimal variance reduction.

Future Work

Future research will aim to broaden the applicability and theoretical foundation of the proposed framework:

• Hybrid variance-minimizing schemes: Integrating quasi–Monte Carlo initialization with NUSA dynamics to enhance convergence speed without compromising unbiasedness.
• Parallel and GPU implementations: Exploiting trajectory independence to achieve significant computational speedups on parallel hardware architectures.
• Nonlinear and time-dependent extensions: Adapting NUSA to nonlinear Fredholm and Volterra equations as well as to stochastic integro-differential problems in physics and engineering.
• Analytical variance bounds: Establishing rigorous convergence theorems and asymptotic variance estimates for unbiased stochastic estimators under general kernel conditions.

Overall, the unbiased stochastic methodology developed in this work provides a powerful foundation for future research in stochastic numerical analysis. Its demonstrated flexibility and robustness suggest broad potential for applications in scientific computing, uncertainty quantification, and complex physical modeling, particularly in domains where deterministic methods become computationally infeasible.