Regularized Kaczmarz Solvers for Robust Inverse Laplace Transforms

Abstract

Inverse Laplace transforms (ILTs) are fundamental to a wide range of scientific and engineering applications—from diffusion NMR spectroscopy to medical imaging—yet their numerical inversion remains severely ill-posed, particularly in the presence of noise or sparse data. The primary objective of this study is to develop robust and efficient numerical methods that improve the stability and accuracy of ILT reconstructions under challenging conditions. In this work, we introduce a novel family of Kaczmarz-based ILT solvers that embed advanced regularization directly into the iterative projection framework. We propose three algorithmic variants—Tikhonov–Kaczmarz, total variation (TV)–Kaczmarz, and Wasserstein–Kaczmarz—each incorporating a distinct penalty to stabilize solutions and mitigate noise amplification. The Wasserstein–Kaczmarz method, in particular, leverages optimal transport theory to impose geometric priors, yielding enhanced robustness for multi-modal or highly overlapping distributions. We benchmark these methods against established ILT solvers—including CONTIN, maximum entropy (MaxEnt), TRAIn, ITAMeD, and PALMA—using synthetic single- and multi-modal diffusion distributions contaminated with 1% controlled noise. Quantitative evaluation via mean squared error (MSE), Wasserstein distance, total variation, peak signal-to-noise ratio (PSNR), and runtime demonstrates that Wasserstein–Kaczmarz attains an optimal balance of speed (0.53 s per inversion) and accuracy (MSE = $4.7\times {10}^{-8}$), while TRAIn achieves the highest fidelity (MSE = $1.5\times {10}^{-8}$) at a modest computational cost. These results elucidate the inherent trade-offs between computational efficiency and reconstruction precision and establish regularized Kaczmarz solvers as versatile, high-performance tools for ill-posed inverse problems.

1. Introduction

Inverse Laplace transforms (ILTs) are basic math transforms in many scientific and engineering fields. They facilitate the manipulation and the solution of differential equations in diverse areas such as heat conduction [1,2], electrical circuit analysis [3], mechanical vibration studies [4], and control systems design [5]. In applications ranging from magnetic resonance imaging (MRI) [6] to industrial process control, the accuracy of ILT calculations is critical to ensuring reliable and robust results.

Despite their broad utility, obtaining accurate and efficient numerical ILT solutions remains challenging. Many traditional approaches struggle with precision and computational speed, especially for complex or high-frequency functions. These difficulties are further compounded when the underlying problems are ill-posed—a common occurrence in real-world scenarios where sparse or noisy data can induce instability [7]. Such ill-posedness is frequently encountered in

Image reconstruction: Recovering detailed images from incomplete or noisy measurements (e.g., in medical or astronomical imaging) [8].

Inverse heat conduction: Estimating temperature distributions or heat fluxes from limited thermal data [9].

Seismic inversion:Inferring subsurface properties from seismic signals [10].

Machine learning: Fitting models in the presence of overfitting or high-dimensional instability [11].

Robust and efficient inverse Laplace transform (ILT) solvers are critical in a wide range of scientific and engineering fields, where they play a big role in the analysis of real and complex data. For instance, in the domain of medical imaging, such as diffusion magnetic resonance imaging (MRI), they are indispensable for generating precise tissue diffusivity maps, which are fundamental for both diagnostics and research. In the domain of polymer science, these techniques facilitate the precise determination of molecular weight distribution, thereby directly impacting the prediction of material properties. Furthermore, their role in engineering applications, such as reverse heat conduction and seismology, where data are often incomplete or noisy, is indispensable. Therefore, to mitigate these challenges, regularization techniques have been extensively developed and applied [12,13,14,15,16]. Landmark advances in computational regularization methods—such as those discussed in [17,18,19]—have provided stable frameworks for tackling ill-posed inverse problems. Yet, many of these methods require delicate parameter tuning and may still suffer from convergence issues in complex settings.

Potential strategies to address these challenges are based on the Kaczmarz method which emerged as a highly promising approach [20]. The original Kaczmarz algorithm was introduced in 1937 by Kaczmarz as an iterative method to solve the consistent systems of linear equations by sequentially projecting the current iterate onto the hyperplanes defined by each equation. Although simple and easy to implement, its performance is highly sensitive to the ordering of the equations. In its classical (cyclic) form, the convergence can be slow when successive rows are “nearly parallel” or when the system is large and ill-posed. The randomized Kaczmarz (RK) algorithm, proposed by Strohmer and Vershynin [21], was a major breakthrough. In their formulation, each row is randomly selected with a probability that increases proportionally to the square of its norm. This choice, which often leads to exponential convergence on average, helps mitigate problems caused by cyclic misordering. More recently, Gower and Richtárik [22] developed a unifying framework for randomized iterative methods in linear systems that not only recovers the RK method as a special case, but also shows that, if a positive definite ‘geometry’ matrix and sampling matrix are appropriately chosen, variants including coordinate descent and block methods can be derived. Their analysis provides a unique convergence theorem from which many known complexity results are derived and opens the door to novel algorithmic variants and preconditioning strategies. In the work of He, Dong, and Li [23], a new probability distribution is proposed for the RK algorithm, in which the probability of selecting a given row in the current iteration becomes proportional to the square of the sine of the angle between that row and the row chosen in the previous iteration. This angle-based selection is based on the observation that, when consecutive hyperplanes are nearly aligned, the benefit of an additional projection diminishes. The method aims to avoid redundant projections and improve the speed of convergence by ‘rewarding’ the selection of rows that are more ‘orthogonal’ (or, in other words, form a greater angle) relative to the previously used row. In the research area, the authors proposed acceleration strategies that, instead of projecting onto a single hyperplane, project onto the intersection of two hyperplanes. These dual-projection or intersection-based schemes have the potential to reduce the error even faster at each iteration. The convergence analysis presented in their study shows that these accelerated variants enjoy better convergence rates than the classical RK method with uniform or norm-based probabilities.

This method has gained popularity for its simplicity, speed, robustness, and convergence properties. These attributes have led to its successful application in computed tomography [24], signal processing [25], and data reconciliation [26]. Its inherent adaptability and flexibility render it an attractive candidate for the exploration of its applicability to the ILT problem.

Nevertheless, the use of the Kaczmarz method for ILT remains largely unexplored, leaving its potential for solving ILTs an open question in numerical analysis. The primary objective of this work is to adapt the Kaczmarz method specifically for ILT, thereby extending its advantages to the domain of ill-posed inverse problems. By incorporating regularization directly into the iterative scheme (whether via Tikhonov, total variation, or more advanced approaches), our goal is to enhance both the computational efficiency and the noise resilience of ILT solvers.

Several established approaches underscore the current efforts to resolve ILT challenges:

MaxEnt (maximum entropy): Utilizes non-parametric solutions that integrate prior knowledge, ensuring stability and robustness [27,28,29].

TRAIn (trust-region algorithm): Dynamically adjusts the trust-region to reconstruct ILT with superlinear convergence [30,31].

ITAMeD (iterative thresholding algorithm for multiexponential decays): Combines ${\ell }_1$-based sparsity with fast iterative methods (e.g., FISTA) to resolve discrete ILT components [32].

PALMA: Merges ${\ell }_1$ regularization and maximum entropy priors to automatically balance smoothness and sparsity [33].

While each of these methods, see Appendix A for all details of these algorithms, offers distinct advantages, they often require careful parameter tuning and may not guarantee robust convergence under challenging conditions. The Kaczmarz method, in contrast, provides an intuitive row-by-row projection strategy and can integrate seamlessly with a variety of regularization frameworks. Its iterative refinement process is similar to the Gram–Schmidt orthogonalization principle, reinforcing geometric insights that often translate into stable, noise-tolerant performance. By adapting the Kaczmarz method to ILT, we aim to exploit these strengths and develop an approach that is both computationally efficient and robust to the complexities of ill-posed Laplace inversion.

As a proof of concept, our investigation targets diffusion NMR (whether PFGSE—pulsed field gradient spin echo or DOSY—diffusion-ordered spectroscopy). In diffusion NMR, the measured signal attenuation due to progressively stronger magnetic field gradients results in a multi-exponential decay curve. This decay is mathematically equivalent to a Laplace transformation of the underlying distribution of diffusion coefficients. Traditional ILT approaches applied in this context are vulnerable to instability and inaccuracy. Nonetheless, the utilization of algebraic reconstruction for multi-exponential decays has previously been established in its most fundamental form. This has been achieved through the application of ILT reconstruction to molecular weight estimation in globular proteins and exopolysaccharides [34,35].

By adapting the Kaczmarz method for ILT and integrating robust regularization, our approach aims to recover complex diffusion profiles, accurately decompose overlapping components, and handle diverse sample behaviors (e.g., polydisperse systems). The enhanced numerical stability and efficiency of the proposed solver hold significant promise for advancing ILT applications in chemistry, biophysics, materials science, and chemical engineering.

This research has far-reaching implications for real-world applications of ILT. The development of an enhanced, reliable, and efficient Kaczmarz solver has the potential to substantially enhance the accuracy of problem-solving and system modelling across a broad spectrum of scientific and engineering domains. The present study has implemented several variations of the Kaczmarz algorithm, incorporating various regularization approaches, including Tikhonov, total variation, and Wasserstein regularization. The enhanced computational capabilities of such regularized Kaczmarz algorithms have the potential to transform numerous research and industrial domains.

Structure of this work. Section 2 details the Kaczmarz method for the ILT problem. In Section 3, focus on discretizing the Fredholm integral formulation and incorporating regularization. Section 4 presents a Wasserstein–Kaczmarz method from solving inverse Laplace transform. Finally, in Section 5, a comparative evaluation of the algorithm against established ILT methods in terms of accuracy, convergence speed, and noise robustness. In Section 6, discusses practical limitations, potential extensions, and the broader implications of the work for solving ill-posed inverse problems.

Objectives and Contributions

This work aims to

We adapt the Kaczmarz method to the specific requirements of ILT by discretizing the Fredholm integral formulation and integrating advanced regularization techniques.

We incorporate regularization methods (such as Tikhonov, total variation, and Wasserstein) directly into the Kaczmarz iterative process to effectively address the ill-posedness inherent in ILT problems.

We conduct an empirical evaluation of the proposed algorithm’s accuracy, convergence speed, and robustness against noise, benchmarking it against existing ILT methodologies.

We discuss the broader implications, limitations, and future directions for extending Kaczmarz-based solvers to other classes of inverse problems.

2. Kaczmarz Method for Solving Linear Systems

The Kaczmarz method is an iterative algorithm for solving the systems of linear equations of the form:

$\mathbf{K}x=\mathbf{Y},$

where $\mathbf{K}\in {\mathbb{R}}^{m\times n}$ is the system matrix, $x\in {\mathbb{R}}^n$ is the unknown solution vector, and $\mathbf{Y}\in {\mathbb{R}}^m$ is the right-hand side vector. First introduced by Stefan Kaczmarz in 1937, this algorithm has found applications in various fields such as image reconstruction, signal processing, and numerical optimization.

The core idea of the Kaczmarz method is to iteratively project the current estimate of the solution onto the hyperplane defined by one row of the Hilbert space. This projection-based approach guarantees convergence to a solution, even for inconsistent systems, by minimizing the residual error in the least-squares sense. This concept of projection is closely related to the Gram–Schmidt orthogonalization process. While the Gram–Schmidt process orthogonalizes a set of vectors, the Kaczmarz method applies successive orthogonal projections to iteratively converge to the solution. This geometric insight underscores the method’s elegance and its potential for stable, noise-tolerant performance.

This Gram–Schmidt orthogonalization, a classical method in linear algebra, systematically constructs an orthonormal basis for a vector space by projecting vectors onto subspaces spanned by previously orthogonalized vectors. Drawing parallels to Gram–Schmidt orthogonalization helps elucidate the geometric and algebraic principles underlying the Kaczmarz method.

2.1. Gram–Schmidt Orthogonalization

The Gram–Schmidt process is a method for orthogonalizing a set of linearly independent vectors in a finite-dimensional vector space. Given a set of m linearly independent vectors $v_1,v_2,\dots ,v_m\in {\mathbb{R}}^n$, the algorithm constructs an orthonormal basis $u_1,u_2,\dots ,u_m\in {\mathbb{R}}^n$ (see Algorithm 1).

Algorithm 1: Gram–Schmidt Orthogonalization

The resulting set $\{u_1,u_2,\dots ,u_m\}$ forms an orthonormal basis for the subspace spanned by the original vectors $\{v_1,v_2,\dots ,v_m\}$.

Each step in the Gram–Schmidt process involves projecting a vector $v_i$ onto the subspace spanned by the previously orthogonalized vectors. The residual vector, obtained after subtracting these projections, lies orthogonal to the subspace, ensuring orthogonality in the resulting basis.

2.2. Kaczmarz Algorithm as an Orthogonalization Process

The Kaczmarz method can be interpreted as a sequential orthogonalization process where the current estimate of the solution vector is projected onto the hyperplanes defined by the rows of A. Each row $a_i$ of A represents a hyperplane in ${\mathbb{R}}^n$, defined by the equation:

$\langle {\mathbf{K}}_i,x\rangle ={\mathbf{Y}}_i.$

The Kaczmarz method can be described step-by-step as follows (see Algorithm 2):

Algorithm 2: Kaczmarz Method for Solving Linear Systems

The algorithm successively reduces $d_i-{\mathrm{Proj}}_{a_i}\left(x_k\right)$ by projecting the current estimate onto the hyperplanes defined by the rows of A. After multiple iterations, these projections converge the estimate to a solution $x^{*}$ that minimizes the residual error in the least-squares sense:

$x^{*}=arg\underset{x}{min}{\parallel \mathbf{K}x-\mathbf{Y}\parallel }^2.$

In the Kaczmarz method, each update step, represented by ${\mathrm{Proj}}_{{\mathbf{K}}_i}\left(x_k\right)={\displaystyle \frac{\langle {\mathbf{K}}_i,x_k\rangle }{\parallel {\mathbf{K}}_i{\parallel }^2}}{\mathbf{K}}_i$, is analogous to the projection operation in Gram–Schmidt orthogonalization. However, rather than orthogonalizing vectors to construct an orthonormal basis, the method iteratively adjusts the solution vector $x_k$ to satisfy the equation of each hyperplane. By iteratively combining these projections, the Kaczmarz method ensures convergence to the optimal solution that minimizes the residual error in the least-squares sense.

3. Adaption of the Kaczmarz Method with Regularization for Solving Inverse Laplace Transforms

In this section, we consider the following Fredholm integral of the first kind as Laplace Transform and it is defined as

(1)$\mathcal{L}\left\{x\left(s\right)\right\}\left(t\right)={\int }_a^b{e}^{-ts}x\left(s\right)ds,$

where $t\ge 0$ and $s\in [a,b]$. We view $\left\{x\right(s\left)\right\}$ as the unknown function we wish to recover (i.e., the “inverse Laplace transform” of $f\left(t\right)$), and the kernel $k(t,s)=e^{-ts}$ is the classical Laplace transform kernel. To solve (1) via the Kaczmarz method, we first discretize both the integral domain in s and the function values in t. The goal is to recast the integral equation into a linear system of the form

$\mathbf{K}x=\mathbf{Y},$

where $\mathbf{K}$ approximates the continuous operator $\mathcal{L}\left\{x\left(s\right)\right\}\left(t\right)={\int }_a^b{e}^{-ts}x\left(s\right)ds$, $\mathbf{Y}$ corresponds to sampled values of $f\left(t\right)$, and x is the discretized representation of $x\left(s\right)$.

Step 1:. Partition the Interval $[a,b]$.

Divide the interval $[a,b]$ into n subintervals by points $s_1,s_2,\dots ,s_n$. Often, a uniform spacing is chosen, i.e., $\Delta s={\displaystyle \frac{b-a}{n}}$, but non-uniform partitions may be used if certain regions of $[a,b]$ require finer resolution.

Step 2:. Approximate the Integral.

Replace the integral by a numerical quadrature:

${\int }_a^b{e}^{-ts}x\left(s\right)ds\approx \sum _{j=1}^n{w}_j{e}^{-t{s}_j}x\left(s_j\right),$

where $w_j$ is the quadrature weight for the j-th subinterval. Common choices for $w_j$ include the rectangle rule ($w_j=\Delta s$) or more accurate schemes such as trapezoidal or Gaussian quadrature.

Step 3:. Sample $\mathcal{L}\left\{x\right(s\left)\right\}\left(t\right)$ at Discrete Points.

Select m distinct values of t, say $t_1,t_2,\dots ,t_m$. For each $t_i$, approximate $\mathcal{L}\left\{x\left(s\right)\right\}\left(t_i\right)$ by the given data or by sampling from the continuous function $\mathcal{L}\left\{x\right(s\left)\right\}\left(t\right)$.

Step 4:. Assemble the Linear System.

Define the matrix $\mathbf{K}\in {\mathbb{R}}^{m\times n}$ by

${\mathbf{K}}_{i,j}=w_j{e}^{-t_i{s}_j},i=1,\dots ,m,j=1,\dots ,n.$

Next, form the vector $\mathbf{Y}\in {\mathbb{R}}^m$ by

${\mathbf{Y}}_i=\mathcal{L}\left\{x\left(s\right)\right\}\left(t_i\right),i=1,\dots ,m.$

Finally, let $x\in {\mathbb{R}}^n$ represent the unknown discrete values $x_j=x\left(s_j\right)$. Thus, the approximate equation

${\int }_a^b{e}^{-t_{i}s}x\left(s\right)ds\approx \mathcal{L}\left\{x\left(s\right)\right\}\left(t_i\right)$

becomes the matrix problem

$\mathbf{K}x=\mathbf{Y},$

where $\mathbf{K}\in {\mathbb{R}}^{m\times n}$ is the system matrix and $\mathbf{Y}\in {\mathbb{R}}^m$ is the right-hand side. By sequentially projecting the current estimate $x_k$ onto the hyperplanes defined by each row of $\mathbf{K}$, the Kaczmarz method converges toward a solution that satisfies $\mathbf{K}x\approx \mathbf{Y}$. In the presence of noise or ill-posedness, however, regularization may be needed to stabilize and improve the quality of the solution. Below, we consider two common forms of regularization:

Tikhonov regularization, which penalizes large norms of the solution.

Total variation (TV) regularization, which promotes piecewise-smooth solutions by penalizing large discrete gradients.

3.1. Tikhonov Regularization

Tikhonov regularization (sometimes called ${\ell }_2$-regularization) introduces a parameter $\lambda >0$ to counteract noise amplification and large solution components. To efficiently incorporate Tikhonov regularization into our Kaczmarz-type iterations, we employ an incremental proximal-gradient (forward–backward) splitting scheme first popularized in signal-processing [36,37]. Concretely, at each step k, we perform

$d_i={\displaystyle \frac{Y_i{K}_i-\lambda \left(x^k\odot K_i\right)}{\parallel K_i{\parallel }_2^2+\lambda }}$

followed by a relaxation $x^{k+{\displaystyle \frac{1}{2}}}=x^k+\omega d_i$, $\omega$ is relaxation parameter that used to be 1 and the global ${\ell }_2$-proximal step $x^{k+1}={\mathrm{prox}}_{{\ell }_2,\lambda }\left(x^{k+{\displaystyle \frac{1}{2}}}\right)$. This per-row approximation converges to the unique minimizer, and is closely related to generalized forward–backward splitting [38]. Together, these steps implement an incremental proximal-gradient method whose fixed point satisfies

$(K^{T}K+\lambda I)x=K^{T}Y$

the normal equations of the Tikhonov-regularized problem. Hence, it is fully consistent with classical Tikhonov regularization.

3.2. Total Variation Regularization

Total Variation (TV) regularization is well suited for problems where the solution x is expected to be piecewise constant or smooth. Instead of penalizing the size of x directly, we penalize its discrete gradient. In a one-dimensional setting (i.e., each component $x_i$ for $i=1,\dots ,n$), a simple finite-difference approximation of the TV penalty involves terms of the form

$\alpha \left(2x^j-x^{j-1}-x^{j+1}\right),$

where $\alpha >0$ controls the strength of the regularization. Thus, the “distance and direction” step can be modified to

$d_i={\displaystyle \frac{{\mathbf{Y}}_i}{\parallel {\mathbf{K}}_i{\parallel }^2}}{\mathbf{K}}_i-\alpha \left[\mathrm{TV}‐\mathrm{penalty}\mathrm{on}{x}_k\right].$

In practice, one updates each component $x_i$ with information from its neighbors, promoting locally constant or smoothly varying solutions. Boundary conditions (e.g., for $x_1$ and $x_n$) may need special handling.

3.3. Algorithm: Kaczmarz Method with Regularization

This subsection presents the Kaczmarz method with regularization algorithm to solve the ILTs. In many applications (e.g., inverse problems or image reconstruction), the solution may also be constrained to be non-negative. The pseudocode below illustrates how to incorporate these features into each Kaczmarz iteration, and we subsequently discuss practical tips for parameter selection, boundary conditions, and stopping criteria (see Algorithm 3).

Algorithm 3: Kaczmarz method for solving linear systems (with optional regularization)

In practice, it is essential to choose the regularization parameters carefully. Larger values of $\lambda$ or $\alpha$ impose stronger regularization, leading to more stabilization but also potentially slower convergence and a biased estimate. Boundary conditions matter especially for the TV case, where one must decide how to handle edges (e.g., via forward/backward differences or reflection). If physical considerations dictate non-negativity (such as concentration profiles in inversion problems), projecting each iterate onto the positive orthant provides a simple and effective enforcement of $x_k\ge 0$. In terms of stopping criteria, one can monitor $\parallel x_{k+1}-x_k\parallel$, the residual norm $\parallel \mathbf{K}{x}_{k+1}-\mathbf{Y}\parallel$, or both. In summary, the Kaczmarz method furnishes a flexible row-by-row updating scheme; by directly embedding Tikhonov or total variation regularization into its “distance and direction” step, it effectively mitigates noise amplification while supporting structurally relevant constraints such as smoothness, sparsity, or non-negativity.

3.4. Discussion and Practical Considerations

Choice of parameters. Both Tikhonov ($\lambda$) and TV ($\alpha$) regularization parameters significantly affect the trade-off between fitting the data and enforcing smoothness/small norms. These typically require some experimentation or cross-validation strategy to select optimal values.

Computational cost. Each inner iteration of the Kaczmarz method only accesses one row of A at a time, which can be advantageous for large problems. However, TV regularization can introduce additional computation because one must update each component $x_i$ with information from its neighbors.

Convergence behavior. The Kaczmarz method often converges quickly in practice, especially if the rows of A are processed in a randomized or optimal order. Regularization may slow down convergence slightly but typically leads to more stable and physically meaningful solutions in noisy or ill-conditioned scenarios.

By integrating either Tikhonov or TV regularization into the classical Kaczmarz method, we obtain a family of algorithms capable of handling ill-posed and large-scale problems, making them suitable for applications in signal/image processing, tomography, and inverse problems where stability and interpretability are crucial.

4. Wasserstein–Kaczmarz Method for Solving Inverse Laplace Transforms

ILTs for numerical solving can become particularly difficult when the parameter s contains limited information, the signal-to-noise ratio is low, or the sampling in the t-domain is sparse. Such scenarios exacerbate the ill-posed nature of the inverse problem, rendering standard regularization approaches ($L_1$ or $L_2$ norms) insufficient. These traditional techniques often fail to capture the underlying distribution accurately, resulting in oversmoothed or excessively noisy reconstructions under sparse or noisy data conditions.

Wasserstein regularization, inspired by optimal transport theory, is a robust method for aligning the reconstructed $x\left(s\right)$ with a reference distribution $x_{\mathrm{ref}}\left(s\right)$. This method is particularly effective for preserving the geometric structure of $x\left(s\right)$, even under high noise conditions. The p-Wasserstein distance between two distributions f and g is defined as

$W_p(f,g)={\left(\underset{\gamma \in \Pi (f,g)}{inf}\int {|x-y|}^{p}d\gamma (x,y)\right)}^{1/p},$

where $\Pi (f,g)$ represents all transport plans $\gamma (x,y)$ with marginals f and g.

The Wasserstein-2 distance ($W_2$) is typically used. It provides a robust alternative by taking advantage of the geometry of the solution space. By aligning the reconstructed solution with a prior reference $x_{\mathrm{ref}}$, it compensates for missing information, preserves physical interpretability, and effectively handles sparse datasets. Moreover, the Wasserstein technique is particularly well-suited for sparse datasets because of its ability to incorporate geometric priors and maintain robustness against noise and missing data. The objective function is defined as

$\underset{x}{min}{\parallel \mathbf{K}x-\mathbf{Y}\parallel }^2+\lambda W_2^2(x,x_{\mathrm{ref}}),$

where

$W_2^2(x,x_{\mathrm{ref}})$: The Wasserstein-2 distance aligns the reconstructed x with the reference distribution $x_{\mathrm{ref}}$.

$x_{\mathrm{ref}}$: A prior distribution reflecting expected physical properties or sparsity patterns.

A proximal update is performed in the Wasserstein-2 metric (this equivalence, introduced by Jordan, Kinderlehrer and Otto, inspired the so-called JKO scheme to approximate these diffusion processes via an implicit discretization of the gradient flow in the Wasserstein metric [39,40]. The Wasserstein-proximal subproblem at iterate $x_k$ is

$x_{k+1}=arg\underset{x}{min}{\displaystyle \frac{1}{\tau }}{W}_2^2(x,x_{\mathrm{ref}})+{∥x-x_{\mathrm{k}}∥}_2^2$

In scenarios where the reference distribution $x_{\mathrm{ref}}$ is unknown, Wasserstein regularization can be adapted to infer or estimate $x_{\mathrm{ref}}$. The following approaches can be applied:

$x_{\mathrm{ref}}={\displaystyle \frac{1}{N}},$

where N is the number of diffusion coefficients. This ensures no prior bias in the reconstruction. Infer $x_{\mathrm{ref}}$ dynamically during reconstruction. Starting with $x_{\mathrm{ref}}={\displaystyle \frac{1}{N}}$, update $x_{\mathrm{ref}}$ at each iteration based on the reconstructed x:

$x_{\mathrm{ref}}^{k+1}=\alpha x^{k+1}+(1-\alpha )x_{\mathrm{ref}}^k,$

where $\alpha$ controls how quickly the reference adapts. Thus, this convex combination is precisely the proximal map of the squared-Wasserstein term, i.e., the first JKO-time-step of the Wasserstein gradient flow [41], and so faithfully regularizes toward $x_{\mathrm{ref}}$ in the $W_2$ metric. This Wasserstein–Kaczmarz method is presented in the Algorithm 4.

Algorithm 4: Wasserstein–Kaczmarz method for solving inverse Laplace transforms

This Wasserstein–Kaczmarz algorithm effectively recovers the sought function $x\left(s\right)$ in the presence of noise, sparse data, or limited prior information. The Wasserstein term aligns x with $x_{\mathrm{ref}}$ to preserve geometrical or distributional characteristics even under challenging conditions, ensuring a more robust and physically interpretable solution compared to conventional methods.

5. Results

In this section, we present the performance results of the evaluated algorithms on synthetic datasets designed to assess their reconstruction accuracy and robustness under noise conditions. The methods compared include classical Kaczmarz, Tikhonov-regularized Kaczmarz, TV-regularized Kaczmarz, Wasserstein–Kaczmarz, CONTIN, maximum entropy (MaxEnt), TRAIn, ITAMeD, and PALMA.

5.1. Synthetic Data Generation and Evaluation Metrics

Synthetic datasets were generated by simulating Laplace-transformed signals from analytically prescribed “ground-truth” distributions—both single- and multi-component—employing Gaussian and exponential functional forms. In order to emulate experimental measurement variability, a controlled noise level of 1% was introduced. In the context of practical diffusion NMR experiments, the absolute value of the diffusion coefficient, denoted by D, and the reconstruction’s noise sensitivity are considered to be essential performance metrics. However, the explicit variation of these parameters is not a prerequisite for the mathematical demonstrations presented herein. The performance of each reconstruction method was quantitatively evaluated using the following metrics:

Total computation time (seconds): Measures the total runtime required for each algorithm to reach convergence.

Mean squared error (MSE): Represents the average squared difference between reconstructed and true distributions, quantifying reconstruction accuracy.

$\mathrm{MSE}={\displaystyle \frac{1}{N}}\sum _{i=1}^N{(x_{\mathrm{recon},i}-x_{\mathrm{true},i})}^2.$

Wasserstein distance: Measures the distributional difference between reconstructed and true distributions, indicating how closely the overall shapes match.

Total variation (TV): Quantifies the smoothness of the reconstructed distribution; lower TV indicates smoother solutions, while higher TV indicates sharper, potentially noisier solutions.

$\mathrm{TV}\left(x\right)=\sum _{i=1}^{N-1}|x_{i+1}-x_i|.$

Peak signal-to-noise ratio (PSNR): Provides a logarithmic measure of the reconstruction quality, with higher values corresponding to better accuracy.

$\mathrm{PSNR}=10{log}_{10}\left({\displaystyle \frac{max{\left(x_{\mathrm{true}}\right)}^2}{\mathrm{MSE}}}\right).$

All algorithms were uniformly initialized, and regularization parameters were optimized through the L-curve method.

5.2. Single-Component Data Reconstruction

Table 1 summarizes the reconstruction results for single-component synthetic data at a 1% noise level, displaying total computation time, MSE, Wasserstein distance, total variation, and PSNR.

Figure 1 illustrates the underlying true distribution along with both the true signal and its noisy counterpart. This visualization provides a qualitative baseline that helps in assessing the effectiveness of the reconstruction algorithms in preserving the original signal structure while mitigating noise.

As illustrated in Figure 2, the performance of various reconstruction algorithms on single-component synthetic data with 1% noise is demonstrated. The TRAIn algorithm attains the lowest mean squared error (MSE) of 1.5 × ${10}^{-8}$ and the smallest Wasserstein distance of 5.4 × ${10}^{-8}$ thereby signifying its superior reconstruction accuracy. However, it should be noted that this is accompanied by a higher computational cost (19.36 s). In contrast, the classical Kaczmarz method, despite its high processing speed (0.62 s), exhibits comparatively elevated errors with an MSE of 1.6 × ${10}^{-7}$ and a greater total variation. The Tikhonov–Kaczmarz and TV–Kaczmarz methods provide a balanced trade-off between accuracy and efficiency, with Tikhonov–Kaczmarz showing a particularly low MSE (2.4 × ${10}^{-8}$) and minimal total variation (3.7 × ${10}^{-7}$) at a moderate runtime (1.26 s). Additionally, MaxEnt delivers a high peak signal-to-noise ratio (35.3 dB), whereas PALMA, despite its extended computational time (261.88 s), exhibits the poorest accuracy, as evidenced by its high MSE and Wasserstein values. The results emphasize the inherent trade-offs between computational efficiency and reconstruction fidelity, guiding the selection of an appropriate algorithm based on specific application needs. The results in Table 1 indicate that the TRAIn and Wasserstein–Kaczmarz methods achieve the highest reconstruction accuracy by recording the lowest MSE and Wasserstein distance values, with competitive PSNR figures. Conversely, while the classical Kaczmarz method is the fastest in terms of computational time, its performance is compromised by higher errors. Tikhonov and TV regularizations achieve a balance between accuracy and computational cost, thus rendering them viable options when moderate reconstruction quality is deemed acceptable alongside efficiency.

5.3. Multi-Component Distribution Reconstruction

The best algorithms evaluated in the previous section were then subjected to further evaluation through the utilization of multi-modal synthetic distributions, characterized by closely overlapping diffusion peaks. This approach emulates realistic and analytically challenging scenarios. In such conditions, characterized by stringent constraints, Wasserstein–Kaczmarz and TRAIn demonstrated remarkable stability and precision, adeptly disentangling and retrieving individual distribution components in the face of substantial noise contamination. From a pragmatic standpoint, the minimum resolvable separation between adjacent diffusion coefficients ($\Delta D_{min}$) is a pivotal parameter for experimental applicability. Nevertheless, it is not necessary to explicitly investigate $\Delta D_{min}$ within the purely mathematical framework presented herein.

Table 2 presents a quantitative comparison of the reconstruction algorithms’ performance on multi-component synthetic data contaminated with 1% noise. Metrics reported include mean squared error (MSE), Wasserstein distance, total variation (TV), Peak signal-to-noise ratio (PSNR), and computational time.

Figure 3 provides a visual representation of the true distribution, alongside the clean and noisy signals. This qualitative visualization underscores the severity of noise interference and highlights the necessity for robust reconstruction methods.

Figure 4 visually demonstrates the algorithms’ reconstruction fidelity, particularly highlighting the robustness and accuracy of Wasserstein–Kaczmarz. Table 2 reveals clear differences among the evaluated methods. TRAIn achieved the lowest MSE and Wasserstein distance, indicating superior accuracy in reconstructing the original distributions despite having a higher total variation. Wasserstein-Kaczmarz also displayed excellent results with slightly higher MSE but significantly reduced computational time, showcasing its strength as an efficient and accurate alternative. MaxEnt provided rapid reconstruction but compromised somewhat on Wasserstein distance accuracy. Lastly, Tikhonov-Kaczmarz, although highly effective at noise mitigation, lagged in computational speed. Based on these observations, the algorithms are ranked as follows:

TRAIn (Best accuracy, balanced efficiency).

Wasserstein–Kaczmarz (excellent accuracy and fastest reconstruction).

MaxEnt (fast computation, moderate accuracy).

Tikhonov–Kaczmarz (good noise reduction, slower computational time).

5.4. Practical Considerations

Optimal parameter tuning was indispensable to each algorithm’s efficacy. Regularization parameters (e.g., $\lambda$ and $\alpha$ for Tikhonov and total-variation methods) were selected via L-curve analysis [42,43]. Furthermore, in practical diffusion NMR measurements, both the accurate recovery of absolute diffusion coefficient values and the reconstruction’s sensitivity to measurement noise critically affect robustness, while the minimum resolvable separation between adjacent diffusion coefficients ($\Delta D_{min}$) represents a key benchmark for peak discrimination. The incorporation of non-negativity constraints further enhanced the physical plausibility of the reconstructed distributions. However, within the purely mathematical framework presented herein, the explicit investigation of noise sensitivity and $\Delta D_{min}$ was deemed unnecessary.

Figure 4

Reconstruction of multi-component synthetic data (1% noise level). Panel (a) Tikhonov–Kaczmarz shows effective noise suppression with minimal error; Panel (b) Wasserstein–Kaczmarz demonstrates superior accuracy and computational efficiency; Panel (c) maximum entropy (MaxEnt) provides good signal preservation with acceptable noise management; Panel (d) TRAIn highlights excellent reconstruction fidelity balanced with moderate processing time.

[Figure omitted. See PDF]

6. Discussion

This study highlights the significant potential of adapting the Kaczmarz method for solving inverse Laplace transform (ILT) problems, particularly when enhanced with robust regularization techniques. The evaluation across multiple synthetic scenarios demonstrates the substantial advantages of these adaptations in terms of accuracy, stability, and computational efficiency compared to traditional methods.

The proposed Wasserstein–Kaczmarz and other algorithms already developed emerged as particularly effective under challenging conditions, such as sparse, noisy, or multi-modal distributions. Wasserstein–Kaczmarz demonstrated superior computational efficiency and strong reconstruction accuracy, making it especially valuable for time-sensitive applications or computationally constrained scenarios. Meanwhile, TRAIn provided exceptional accuracy, albeit with a slightly increased computational load, reinforcing its suitability for applications requiring high-fidelity reconstruction.

While established algorithms such as MaxEnt, ITAMeD, and PALMA each possess distinct strengths—ranging from rapid computation to specific robustness features—they also exhibit limitations, particularly in complex or highly noisy conditions. For instance, MaxEnt, despite its rapid execution, showed a compromise in accuracy, particularly evident through increased Wasserstein distances. Similarly, Tikhonov–Kaczmarz, despite its strong performance in noise suppression, was limited by slower computational speeds.

From a practical laboratory standpoint, several additional considerations must guide implementation. First, the accurate determination of the absolute diffusion coefficient D is indispensable for meaningful physical interpretation, necessitating careful calibration and reference measurements. Second, the minimum resolvable separation between adjacent diffusion peaks ($\Delta D_{min}$) must be empirically quantified to ensure the reliable discrimination of closely spaced species. Finally, a more exhaustive investigation of signal-to-noise ratio (SNR) effects—in particular, how SNR degradation influences both peak recovery and parameter estimation—will be essential to translate these algorithms from simulation to experiment. Although such detailed experimental studies fall outside the strict mathematical framework presented here, they constitute the next steps for validating and optimizing Kaczmarz-based ILT solvers in real-world analytical contexts.

Beyond the specific findings of this study, the broader implications suggest significant potential for extending Kaczmarz-based solvers to a wide variety of other inverse problems, including tomography, seismic inversion, and medical imaging. Despite the clear benefits, limitations such as parameter tuning complexity and computational scalability need further exploration. Future research directions should focus on developing adaptive parameter selection methods and investigating hybrid regularization strategies that could further enhance performance and robustness across diverse classes of inverse problems.

7. Conclusions

The adapted Kaczmarz method, enhanced with advanced regularization strategies such as Wasserstein and trust-region techniques, represents a significant advancement in addressing the complexities associated with inverse Laplace transform problems. The results underscore the Wasserstein–Kaczmarz algorithm’s optimal balance between computational efficiency and reconstruction accuracy, positioning it as a strong candidate for real-world applications demanding rapid and accurate solutions.

The TRAIn algorithm also stands out by delivering unmatched accuracy under challenging conditions, demonstrating its value for applications prioritizing precise reconstructions despite higher computational demands.

The methodological advances presented herein not only establish new performance benchmarks for ILT solvers but also pave the way for the further investigation of iterative regularization techniques in addressing ill-posed inverse problems across a wide range of scientific and engineering applications. To ensure reproducibility and to support ongoing development, the full implementation of all algorithms described in this study (software version v1.0.0) is available at the project’s GitHub repository: https://github.com/fmarrabal/kaczmarz-ilt-solver (accessed on 25 June 2025).

Author Contributions

Conceptualization, M.G.-L., I.F. and F.M.A.-C.; methodology, M.G.-L., E.V. and V.V.; software, M.G.-L. and V.V.; validation, M.G.-L., E.V. and V.V.; formal analysis, M.G.-L.; investigation, M.G.-L., E.V. and V.V.; resources, F.M.A.-C.; data curation, M.G.-L.; writing—original draft preparation, M.G.-L.; writing—review and editing, M.G.-L., E.V., V.V., I.F. and F.M.A.-C.; visualization, M.G.-L.; supervision, I.F. and F.M.A.-C.; project administration, F.M.A.-C.; funding acquisition, I.F. and F.M.A.-C. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

The data presented in this study are openly available in https://github.com/fmarrabal/kaczmarz-ilt-solver (accessed on 25 June 2025).

Conflicts of Interest

The authors declare no conflicts of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1 Visual comparison of the true distribution, the true signal, and the noisy signal.

Figure 2 Comparison of reconstruction algorithms on single-component synthetic data (1% noise level). Panel (a) shows the Wasserstein–Kaczmarz method, which achieves robust noise mitigation with minimal reconstruction error as quantified by low MSE and Wasserstein distance. Panel (b) displays the TRAIn algorithm, noted for its superior accuracy and balanced performance between error reduction and computational efficiency. Panel (c) illustrates the maximum entropy (MaxEnt) approach, offering a compromise that preserves signal details while effectively suppressing noise. Panel (d) presents the Tikhonov–Kaczmarz method, which delivers efficient regularization with an optimal trade-off between reconstruction fidelity and processing time.

Figure 3 Comparison of true and noisy signals with the underlying true distribution.

Table 1

Algorithm performance comparison on single-component synthetic data (1% noise level).

Algorithm	MSE	Wasserstein	TV	PSNR (dB)	Time (s)
Classical Kaczmarz	1.6 × ${10}^{-7}$	5.6 × ${10}^{-5}$	3.5 × ${10}^{-2}$	33.57	0.42
Tikhonov–Kaczmarz	2.4 × ${10}^{-8}$	5.2 × ${10}^{-5}$	3.7 × ${10}^{-7}$	33.94	1.26
TV–Kaczmarz	4.1 × ${10}^{-7}$	1.3 × ${10}^{-4}$	2.8 × ${10}^{-2}$	29.5	3.06
Wasserstein–Kaczmarz	4.7 × ${10}^{-8}$	6.1 × ${10}^{-5}$	6.9 × ${10}^{-3}$	21.54	0.53
CONTIN	1.6 × ${10}^{-7}$	9.0 × ${10}^{-5}$	3.2 × ${10}^{-2}$	33.55	2.13
MaxEnt	1.1 × ${10}^{-7}$	5.7 × ${10}^{-5}$	4.2 × ${10}^{-2}$	35.3	1.46
TRAIn	1.5 × ${10}^{-8}$	5.4 × ${10}^{-8}$	8.1 × ${10}^{-3}$	29.4	19.36
ITAMeD	2.2 × ${10}^{-8}$	6.2 × ${10}^{-5}$	1.9 × ${10}^{-3}$	26.7	4.35
PALMA	1.2 × ${10}^{-6}$	3.2 × ${10}^{-4}$	3.3 × ${10}^{-2}$	24.94	261.88

Table 2

Algorithm performance comparison on multi-component synthetic data (1% noise level).

Algorithm	MSE	Wasserstein	TV	PSNR (dB)	Time (s)
Tikhonov–Kaczmarz	5.4 × ${10}^{-7}$	3.5 × ${10}^{-4}$	1.2 × ${10}^{-2}$	19.19	108.22
Wasserstein–Kaczmarz	3.6 × ${10}^{-7}$	2.6 × ${10}^{-4}$	1.4 × ${10}^{-2}$	20.98	6.67
MaxEnt	4.1 × ${10}^{-7}$	23.1 × ${10}^{-4}$	1.6 × ${10}^{-2}$	20.42	2.08
TRAIn	2.1 × ${10}^{-7}$	2.3 × ${10}^{-4}$	7.9 × ${10}^{-2}$	19.49	22.35

Appendix A. Regularization Optimization for Ill-Posed Problems

Ill-posed problems are a class of mathematical problems that fail to meet one or more of the criteria for well-posedness established by Hadamard [7]. According to Hadamard, a problem is well-posed if it satisfies the three following conditions:

Existence: A solution to the problem exists.

Uniqueness: The solution is unique.

Stability: The solution’s behavior changes continuously with changes in the initial data.

If any of these criteria are not satisfied, the problem is termed ill-posed. Ill-posed problems commonly arise in inverse problems, integral equations, and partial differential equations, particularly in scenarios involving physical systems where data is noisy, incomplete, or ambiguous. Such problems are inherently unstable, meaning that small perturbations in the input data can lead to disproportionately large deviations in the solution.

These challenges necessitate additional techniques to stabilize the problem and obtain meaningful solutions. In the absence of stabilization, solutions to ill-posed problems may be non-unique, physically implausible, or overly sensitive to noise in the data.

The general form of an ill-posed problem can be expressed mathematically as(A1)$$Kx=f,K:X\to Y$$ where K is a compact operator mapping elements from the vector space X to Y, $x\in X$ represents the unknown solution, and $f\in Y$ is the data. Compact operators, by their nature, amplify noise, exacerbating the instability of the problem.

In practical scenarios, the observed data f is often corrupted by noise, which is represented as(A2)$$f_{\delta }=f+\delta ,\parallel f-f_{\delta }\parallel <\delta$$ where $\delta$ denotes the noise level. The presence of noise makes it difficult to directly solve Equation (A1) because small errors in f can lead to large errors in x. This inherent instability underscores the need for regularization techniques to constrain the solution space and ensure stability.

Given the perturbed data $f_{\delta }$, the solution set can be expressed as(A3)$$Q_{\delta }:=\{x\in X:\parallel Kx-f_{\delta }\parallel \le \delta \}$$ This formulation defines a set of possible solutions that satisfy the noise-constrained condition. However, the dependence of the solution $x_{\delta }$ on the noisy data $f_{\delta }$ introduces significant challenges. The inequality $\parallel Kx-f_{\delta }\parallel \le \delta$ does not inherently ensure either stability or uniqueness of the solution. Small variations in $f_{\delta }$ can produce large variations in $x_{\delta }$, leading to instability. Moreover, multiple values of x may satisfy the inequality, making the solution non-unique.

These issues highlight the core difficulties of solving ill-posed problems directly. Without further constraints or regularization, the resulting solutions are often impractical, unreliable, or physically meaningless. Regularization methods address these deficiencies by imposing additional structure or prior knowledge about the solution, which will be discussed in subsequent sections.

To address the instability inherent in ill-posed problems, regularization methods impose additional constraints or introduce prior knowledge about the solution. This process helps ensure that the solution is stable, unique, and physically meaningful. Regularization techniques are discussed in detail in the following sections.

To address these challenges, additional information is introduced in the form of regularization. Regularization imposes constraints or preferences on the solution, effectively narrowing the solution space to ensure stability and feasibility. This additional information can reflect physical, mathematical, or empirical considerations. For instance, in many applications, solutions are expected to be smooth or adhere to certain positivity constraints.

One widely accepted approach to regularization is Tikhonov regularization [19]. Tikhonov regularization stabilizes the solution by introducing a stabilizing functional, modifying the optimization problem to(A4)$$\underset{x\in X}{min}\parallel Kx-f_{\delta }{\parallel }^2+\lambda {\parallel Lx\parallel }^2$$ where L is a linear operator that incorporates prior knowledge about the solution (e.g., smoothness or sparsity), and $\lambda >0$ is the regularization parameter that balances the trade-off between fidelity to the data and the stabilizing functional.

The choice of L and $\lambda$ significantly impacts the solution. Common choices for L include the identity operator (standard Tikhonov regularization) or higher-order derivatives to enforce smoothness. Various methods exist to determine the optimal $\lambda$, such as the L-curve criterion [42], cross-validation, or generalized cross-validation [44].

The operator K in Equation (A1) is often a continuous compact operator. For many problems, K is expressed as a Fredholm integral equation of the first kind:(A5)$$Kx:={\int }_a^{b}k(t,s)x\left(s\right)ds=f\left(t\right),K:X\to Y,t\in (c,d)$$ Here, $k(t,s)$ is the kernel function, which defines the relationship between the unknown function $x\left(s\right)$ and the data $f\left(t\right)$. The kernel function $k(t,s)$ itself lies in a Hilbert space, providing a structured framework for its analysis.

A Hilbert space is a complete vector space equipped with an inner product that induces a norm. The inner product, denoted by $\langle u,vs.\rangle$, is a bilinear form that satisfies the following properties for all $u,v,w\in H$ and $\alpha \in \mathbb{R}$ or $\mathbb{C}$:

Linearity: $\langle \alpha u+v,w\rangle =\alpha \langle u,w\rangle +\langle v,w\rangle$.

Symmetry: $\langle u,vs.\rangle =\overline{\langle v,u\rangle }$.

Positive definiteness: $\langle u,u\rangle \ge 0$, and $\langle u,u\rangle =0$ if and only if $u=0$.

A norm is induced by the inner product as $\parallel u\parallel =\sqrt{\langle u,u\rangle }$. Completeness means that every Cauchy sequence in the Hilbert space converges to an element within the space. Examples of Hilbert spaces include ${\ell }^2$ (the space of square-summable sequences) and $L^2(a,b)$ (the space of square-integrable functions on $[a,b]$).

Compact operators such as K are characterized by their ability to map bounded sets in X to relatively compact sets in Y. This property exacerbates the ill-posedness of the problem, as small perturbations in $f\left(t\right)$ can lead to disproportionately large changes in $x\left(s\right)$. Compactness implies that the image under K of any bounded sequence in X has a subsequence that converges in Y. The integral form of K highlights its dependence on the kernel function $k(t,s)$, which encapsulates the physics or underlying principles of the problem at hand.

Discretization of the Fredholm integral operator is commonly performed to solve the problem numerically. This process leads to the matrix formulation:(A6)$$\mathbf{Y}=\mathbf{K}\mathbf{x}$$ where $\mathbf{Y}$ represents the discrete approximation of f, $\mathbf{x}$ is the discrete solution, and $\mathbf{K}$ is the discretized operator derived from $k(t,s)$. However, discretization introduces numerical challenges such as instability due to the ill-conditioning of $\mathbf{K}$. Ill-conditioning means that small errors in $\mathbf{f}$ can result in large errors in $\mathbf{x}$, making the system sensitive to perturbations. These challenges necessitate the use of regularization techniques to mitigate the effects of noise and stabilize the solution.

The Fredholm integral equation of the first kind provides a mathematical framework for many inverse problems. The Hilbert space setting and compact operator properties play a central role in the analysis and solution of such equations, while regularization techniques are crucial for addressing their inherent instability.

Herein, we delve into several widely used algorithms developed to solve such problems. These include CONTIN [45], Maximum Entropy [29], ITAMeD [32], TRAIn [31], and PALMA [33], among others. Each algorithm employs unique approaches to tackle the challenges posed by ill-posedness, noise, and instability in the context of Fredholm integral equations. We will explore their mathematical foundations, implementation strategies, and applications across various scientific fields.

Appendix A.1. The CONTIN Algorithm: A Constrained Regularization Approach

The CONTIN algorithm, developed by S.W. Provencher, is a powerful and flexible method for solving noisy linear operator equations, including Fredholm integral equations of the first kind. It is designed to handle the inherent ill-posedness of such problems by incorporating constraints and regularization techniques. CONTIN addresses equations of the form:(A7)$$f\left(t\right)={\int }_a^{b}k(t,s)x\left(s\right)ds+{ϵ}_t,t=1,\dots ,N,$$ where $f\left(t\right)$ are the observed data points, $k(t,s)$ represents the kernel function, $x\left(s\right)$ is the unknown solution, and ${ϵ}_k$ accounts for noise. This formulation encompasses a wide range of experimental contexts, such as photon correlation spectroscopy and Laplace transforms. To solve Equation (A7), CONTIN discretizes it into a system of linear algebraic equations:(A8)$$\mathbf{Y}=\mathbf{K}\mathbf{x}+ϵ,$$ where $\mathbf{Y}$ is the vector of observed data, $\mathbf{K}$ is the discretized kernel matrix, $\mathbf{x}$ is the solution vector, and $ϵ$ represents noise. The ill-conditioning of $\mathbf{K}$ necessitates regularization to ensure stable and meaningful solutions.

CONTIN minimizes the following functional to determine the optimal solution:(A9)$$V\left(\alpha \right)={\parallel \mathbf{W}(\mathbf{Y}-\mathbf{K}\mathbf{x})\parallel }^2+{\alpha }^2{\parallel \mathbf{R}\mathbf{x}\parallel }^2,$$ where $\mathbf{W}$ is a weighting matrix based on the covariance of the noise, $\mathbf{R}$ represents the regularization operator, and $\alpha$ is the regularization parameter. The second term penalizes non-physical or oscillatory solutions, promoting stability and smoothness.

The functional $V\left(\alpha \right)$ comprises two main components:

Data fidelity term: ${\parallel \mathbf{W}(\mathbf{Y}-\mathbf{K}\mathbf{x})\parallel }^2$, which quantifies the deviation of the model solution from the observed data, weighted by the noise covariance matrix $\mathbf{W}$.

Regularization term: ${\alpha }^2{\parallel \mathbf{R}\mathbf{x}\parallel }^2$, where $\mathbf{R}$ encodes prior knowledge about the solution, such as smoothness or sparsity. The parameter $\alpha$ balances these two terms.

The minimization of $V\left(\alpha \right)$ is performed using numerical optimization techniques. The gradient of $V\left(\alpha \right)$ with respect to $\mathbf{x}$ is given by(A10)$$\nabla V\left(\alpha \right)=-2{\mathbf{K}}^T{\mathbf{W}}^T(\mathbf{Y}-\mathbf{K}\mathbf{x})+2{\alpha }^2{\mathbf{R}}^T\mathbf{R}\mathbf{x}.$$ Iterative solvers, such as conjugate gradient or quasi-Newton methods, are employed to find $\mathbf{x}$, which minimizes $V\left(\alpha \right)$ under the imposed constraints.

The regularization parameter $\alpha$ plays a critical role in balancing the data fidelity and smoothness of the solution. CONTIN iteratively adjusts $\alpha$ using the following algorithm (Algorithm A1):

Algorithm A1: CONTIN algorithm for adaptive regularization

[Image omitted. Please see PDF.]

This algorithm allows for the incorporation of equality and inequality constraints, which reflect prior knowledge about the solution. For instance, non-negativity constraints can be imposed to ensure that the solution $s\left(\lambda \right)$ remains physically meaningful:(A11)$$\begin{array}{cc}\hfill \mathbf{D}\mathbf{x}& \ge \mathbf{d},\hfill \end{array}$$(A12)$$\begin{array}{cc}\hfill \mathbf{E}\mathbf{x}& =\mathbf{e},\hfill \end{array}$$ where $\mathbf{D}$ and $\mathbf{E}$ are user-defined matrices, and $\mathbf{d}$ and $\mathbf{e}$ are vectors specifying the constraints.

Appendix A.2. The Maximum Entropy Algorithm: A Regularization Method for Ill-Posed Problems

The maximum entropy (MaxEnt) algorithm is a widely used method for solving ill-posed problems. MaxEnt relies on the principle of entropy maximization to find the most probable solution consistent with the observed data, incorporating prior knowledge while avoiding overfitting to noise. We assume an objective of the form,$${\chi }^2\left(x\right)-\alpha S\left(x\right),$$ or equivalently$$J\left(x\right)={\displaystyle \frac{1}{2}}\sum _{j=1}^p{\displaystyle \frac{{\left(Y_j-{\sum }_{i=1}^n{K}_{ji}{x}_i\right)}^2}{{\sigma }_j^2}}-\alpha \sum _{i=1}^n{x}_{i}ln{x}_i,$$ where $\alpha >0$ is a Lagrange multiplier (or regularization parameter) controlling the balance between data fidelity ${\chi }^2$ and entropy $S\left(x\right)$.

For the entropy term$$S\left(x\right)=-\sum _i{x}_{i}ln\left(x_i\right),$$ its gradient with respect to $x_i$ is $-1-ln\left(x_i\right)$. For the ${\chi }^2$ term, we carefully accumulate contributions from the residual $\mathbf{Y}-\mathbf{K}x$, scaled by ${\sigma }_j^{-2}$.

After exponentiating the negative gradient, solutions are often re-scaled to respect constraints such as ${\sum }_i{x}_i=F$ (e.g., $F=1$ for a probability distribution). This step is crucial for maintaining physically or probabilistically meaningful solutions.

To solve this ill-posed problem, MaxEnt maximizes the entropy functional S:(A13)$$S=-\sum _{i=1}^n{x}_{i}log\left({\displaystyle \frac{x_i}{F}}\right),F=\sum _{i=1}^n{x}_i,$$ where $x_i$ is the discretized value of $x\left(s\right)$ and F ensures normalization. The entropy maximization is subject to constraints that enforce consistency with the observed data. The optimization problem is formulated as(A14)$$Q=S-\lambda {\chi }^2,$$ where $\lambda$ is a Lagrange multiplier, and ${\chi }^2$ quantifies the fit to the data:(A15)$${\chi }^2=\sum _{i=1}^p{\displaystyle \frac{{\left({\mathbf{Y}}_j-{\sum }_{i=1}^n{\mathbf{K}}_{ij}{x}_i\right)}^2}{{\sigma }_j^2}},$$ where ${\mathbf{Y}}_i$ is the observed data, ${\mathbf{K}}_{ij}$ represents the integral kernel, $x_j$ is the solution vector, and ${\sigma }_i$ is the standard deviation of the noise. The optimization is performed iteratively (see Algorithm A2):

Algorithm A2: MaxEnt algorithm

[Image omitted. Please see PDF.]

Appendix A.3. Trust-Region Algorithm: A Regularization Approach for Inversion Problems

The trust-region algorithm for the inversion (TRAIn) is a robust and iterative regularization technique. Traditional iterative methods often struggle to balance convergence with robustness, resulting in solutions that may either converge too slowly or become unstable in the presence of measurement errors.

To address these issues, this algorithm has been developed as a powerful iterative regularization method. TRAIn minimizes the standard least-squares objective$$\Phi \left(x\right)={\parallel \mathbf{K}x-\mathbf{Y}\parallel }_2^2,$$ while imposing two key constraints: non-negativity of the solution ($\mathbf{x}\ge 0$) and a bound on the change between consecutive iterations ($\parallel x^{(k+1)}-x^{\left(k\right)}{\parallel }_2\le r^{\left(k\right)}$), where $r^{\left(k\right)}$ denotes the trust-region radius at iteration k. These constraints ensure that each iteration remains within a reliably stable region, effectively controlling the step size and preventing divergence.

At the heart of TRAIn lies the solution of a trust-region subproblem at each iteration. A trial step $\mathbf{z}$ is computed by approximately minimizing$$\underset{\mathbf{z}}{min}\Phi (x^{\left(k\right)}+\mathbf{z})\mathrm{subject}\mathrm{to}{\parallel \mathbf{z}\parallel }_2\le r^{\left(k\right)},$$ using the truncated conjugate gradient method. The quality of this trial step is then quantified by the ratio$$\rho ={\displaystyle \frac{\Delta {\Phi }_{\mathrm{actual}}}{\Delta {\Phi }_{\mathrm{predicted}}}},$$ which compares the actual reduction in the objective to that predicted by the model. Based on the value of $\rho$, the trust-region radius is adjusted adaptively:

If $\rho >{\mu }_2$, the radius is increased ($r^{(k+1)}={\nu }_2{r}^{\left(k\right)}$).

If ${\mu }_1\le \rho \le {\mu }_2$, the radius remains unchanged.

If $\rho <{\mu }_1$, the radius is decreased ($r^{(k+1)}={\nu }_1{r}^{\left(k\right)}$).

Typical parameter values are ${\mu }_1=0.25$, ${\mu }_2=0.75$, ${\nu }_1=0.5$, and ${\nu }_2=2$. The iterative process is terminated when the residual norm satisfies$$\parallel \mathbf{K}{x}^{\left(k\right)}{-\mathbf{Y}\parallel }_2\le \tau ϵ,$$ with $ϵ$ representing the noise level and $\tau$ is a constant slightly greater than 1 (e.g., $\tau =1.02$).

The following algorithm provides a detailed exposition of TRAIn, describing the adaptive strategies that guarantee both stability and convergence even in the presence of significant noise (see Algorithm A3).

Algorithm A3: Trust-region algorithm for the inversion (TRAIn)

[Image omitted. Please see PDF.]

Appendix A.4. Iterative Thresholding Algorithm for Inversion Problems

The iterative thresholding algorithm for inversion problems is a robust and efficient regularization method. By employing sparsity-promoting regularization and leveraging the fast iterative shrinkage-thresholding algorithm (FISTA). This provides high-resolution reconstruction of the inverse Laplace transform. This task is notoriously ill-posed, as it is sensitive to noise and requires regularization. The method formulates the problem as an ${\ell }_1$-regularized least-squares minimization:(A16)$$\underset{x}{min}{\displaystyle \frac{1}{2}}{\parallel \mathbf{K}x-\mathbf{Y}\parallel }_2^2+\tau {\parallel x\parallel }_1,$$ where

$\mathbf{K}$ is the kernel compact operator mapping elements from vector space X to Y;

x is the discretized unknown solution;

$\mathbf{Y}$ is the discrete approximation of f;

${\parallel x\parallel }_1$ promotes sparsity in x;

$\tau$ is the regularization parameter.

FISTA is employed to efficiently solve the optimization problem. The algorithm iteratively updates the solution x as follows (see Algorithm A4):

Algorithm A4: Iterative update algorithm with proximal operator

[Image omitted. Please see PDF.]

Appendix A.5. PALMA Algorithm: Hybrid Regularization for Inversion Problems

The PALMA algorithm, standing for “Proximal Algorithm for ${\ell }_1$ combined with MAxEnt prior,” is a hybrid regularization method designed to address the challenges of inversion problems. This algorithm integrates the principles of sparsity and entropy maximization to provide robust and efficient solutions for the inverse Laplace transform problem. The PALMA algorithm solves the following constrained optimization problem:(A17)$$\underset{x\in {\mathbb{R}}^N}{min}\Psi \left(x\right)\mathrm{subject}\mathrm{to}{\parallel \mathbf{K}x-\mathbf{Y}\parallel }_2\le \eta ,$$ where:

$\Psi \left(x\right)$ is the hybrid regularization function combining entropy and sparsity priors;

$\eta$ is the tolerance related to the experimental noise level.

The hybrid regularization function is defined as(A18)$$\Psi \left(x\right)=\lambda \mathrm{ent}(x,a)+{(1-\lambda )\parallel x\parallel }_1,$$ where

$\mathrm{ent}(\mathbf{X},a)$ is the negative Shannon entropy with a flat prior a:(A19)$$\mathrm{ent}(x,a)=\sum _{i=1}^N{\displaystyle \frac{x_i}{a}}log\left({\displaystyle \frac{x_i}{a}}\right),$$

${\parallel x\parallel }_1={\sum }_{i=1}^N\left|x_i\right|$ is the ${\ell }_1$ norm promoting sparsity;

$\lambda \in [0,1]$ balances the trade-off between entropy and sparsity.

The PALMA algorithm employs a proximal optimization approach to solve the hybrid regularization problem. The key steps are included in the algorihtm as follows (Algorithm A5):

Algorithm A5: Proximal and projection algorithm

[Image omitted. Please see PDF.]

The PALMA algorithm requires the careful selection of parameters. The entropy prior a is typically chosen as the area under the expected spectrum, providing a meaningful baseline for the algorithm. The noise tolerance $\eta$ is estimated based on the noise level in the experimental data, ensuring the algorithm’s robustness against measurement inaccuracies. Lastly, the weight $\lambda$ balances the influence of entropy and sparsity in the regularization process, with suggested values ranging between 0.01 and 0.5 to achieve an optimal trade-off between these competing factors.