Radiological mapping and uncertainty quantification by a fast Microcanonical Langevin Monte Carlo sampler

This paper introduces a fast Microcanonical Langevin Monte Carlo (MCLMC) sampler for radiological image reconstruction that, when accelerated by GPUs, provides both high-accuracy activity estimates and rapid uncertainty quantification, outperforming traditional methods like ML-EM in speed and reliability for nuclear emergency response.

The Gist

Imagine you are trying to find hidden radioactive "hot spots" in a large, foggy field. You have a detector that counts radiation, but it's like trying to figure out where the heat is coming from in a dark room just by listening to the crackling of a fire. You know that there is fire, and you know how loud it is, but figuring out exactly where it is and how big it is, while also knowing how sure you are about your guess, is incredibly difficult.

This is the challenge of radiological mapping.

The Old Way: The "Best Guess" with a Blindfold

Traditionally, scientists use a method called ML-EM (Maximum Likelihood Expectation-Maximization). Think of this as a student trying to solve a math problem by guessing, checking the answer, and adjusting their guess over and over.

The problem? The student doesn't know when to stop.

• If they stop too early, they haven't solved it well (under-fitting).
• If they keep going too long, they start memorizing the "noise" or random errors in the data instead of the real signal (over-fitting).
• Worst of all, the student gives you a single answer but no confidence score. They say, "The fire is here," but they don't tell you if they are 99% sure or just guessing.

The New Way: The "Microcanonical Langevin" Explorer

This paper introduces a new, super-fast tool called MCLMC (Microcanonical Langevin Monte Carlo).

Imagine instead of one student guessing, you have a swarm of 10,000 tiny, energetic explorers released into the foggy field at the same time.

1. The Exploration: These explorers run around, testing different possible maps of where the fire could be. They don't just guess; they use physics-inspired rules to move efficiently through the "fog" of possibilities.
2. The Map: After a short time, you look at where all the explorers ended up. If 9,000 of them cluster in one spot, you know with high confidence the fire is there. If they are scattered everywhere, you know you aren't sure yet.
3. The Result: You get not just a map of the fire, but a map of your uncertainty. You can see exactly where the explorers were confident and where they were confused.

Why is this a Big Deal?

The paper highlights three main superpowers of this new method:

1. Speed: The "GPU Super-Express"
Usually, running a swarm of 10,000 explorers takes hours or days on a standard computer. The authors found that by using a GPU (the powerful graphics card usually found in gaming computers), they could get the whole swarm to finish its job in about 10 seconds.

• Analogy: It's like upgrading from a bicycle courier to a high-speed maglev train. They can map a complex area in the time it takes to brew a cup of coffee.

2. No "Stop Button" Needed
With the old method, you had to guess when to stop the process to avoid over-fitting. With the MCLMC swarm, you just let them run until they settle down. Because they are sampling the whole picture of possibilities, they naturally avoid getting stuck in bad guesses. They find the true image without needing a human to say, "Okay, stop now!"

3. Smarter Guessing with "Spatial Memory"
The paper tested two types of "rules" (priors) for the explorers:

• Independent Guessing: Each explorer guesses a spot without caring about its neighbors.
• Gaussian Process Prior (GPP): The explorers talk to each other. They know that if there is a fire in one spot, the spots right next to it are likely to be hot too.
• Result: The "talking" explorers (GPP) created a much clearer, sharper map with less uncertainty, even when the data was messy or limited.

Real-World Impact

The researchers tested this on both fake data and real data from a flight over a field with hidden radioactive sources.

• The Result: The new method produced maps that looked almost exactly like the "ground truth" (the actual hidden sources).
• The Bonus: It provided a "confidence map" alongside the image.

Why Should You Care?

In a nuclear emergency (like a dirty bomb or a leak), time is everything. First responders need to know:

1. Where is the danger?
2. How bad is it?
3. How sure are we?

This new tool gives them a high-quality map and a confidence score in seconds, not hours. It allows robots or humans to make faster, safer decisions, like "Go this way, it's safe," or "Stay back, we aren't sure about that area yet."

In short: This paper turns a slow, uncertain guessing game into a fast, confident, and highly accurate exploration, helping us handle nuclear risks with much better eyes.

Technical Summary

1. Problem Statement

Radiological mapping is critical for nuclear emergency response and environmental management. It involves reconstructing the spatial and intensity distribution of radioactive sources from detector count data and contextual information (e.g., LiDAR, camera data).

• The Challenge: The reconstruction is a high-dimensional, often highly underdetermined inverse problem (many more voxels/pixels than measurements).
• Limitations of Current Methods:
  • ML-EM (Maximum Likelihood Expectation-Maximization): The standard method provides only point estimates without uncertainty quantification (UQ). It suffers from overfitting or underfitting depending on the arbitrary number of iterations chosen, lacking a well-defined stopping criterion.
  • Existing UQ Methods: Traditional Bayesian approaches (e.g., Iterative Bayesian Unfolding) are computationally prohibitive for high-dimensional problems ($O(I^2J^3)$ scaling). Other Markov Chain Monte Carlo (MCMC) samplers (like Hamiltonian Monte Carlo) are too slow for near-real-time operation. Approximations like Laplace approximation fail when the posterior is non-Gaussian (e.g., multimodal or skewed).
• Goal: Develop a method that enables near-real-time image reconstruction with robust uncertainty quantification for high-dimensional radiological maps.

2. Methodology

The authors propose using the Microcanonical Langevin Monte Carlo (MCLMC) sampler, implemented via the Blackjax library, to solve the radiological inverse problem.

• Forward Model:

  • The problem is modeled as a Poisson process where measured counts $x$ follow a Poisson distribution with mean $\lambda$.
  • $\lambda = Vw + bt$, where $w$ is the vector of voxel intensities, $V$ is the system matrix (geometric efficiency), $b$ is the background rate, and $t$ is dwell time.
  • The negative log-likelihood (Poisson loss) is minimized in traditional methods but serves as the likelihood term in the Bayesian framework.

• Bayesian Framework:

  • The goal is to sample from the joint posterior distribution $p(w, b | x) \propto p(x|\lambda) \cdot p(w, b)$.
  • Priors: Two priors were tested:
    1. Truncated Gaussian Prior: Independent Gaussian distributions for each voxel (truncated at zero).
    2. Gaussian Process Prior (GPP): Incorporates spatial correlations between voxels, automatically tuning hyperparameters via empirical Bayes.
  • Link Function: Since intensities must be non-negative, a link function (element-wise exponential) maps latent samples $\phi$ to the physical variables $w$ and $b$.

• The MCLMC Sampler:

  • Unlike traditional MCMC, MCLMC uses physics-inspired dynamics (microcanonical ensemble) to sample high-dimensional distributions efficiently.
  • It features automatic hyperparameter tuning and extremely fast "burn-in" (convergence).
  • Posterior Reduction: To generate 2D maps, the authors compare reducing the posterior to a point estimate via:
    • Marginal Mean
    • Marginal Mode
    • Joint Mode (approximated by the sample with the highest log-posterior or via L-BFGS optimization).
  • Uncertainty Estimation: The 68% Highest Density Interval (HDI) of the marginal distribution is used as the uncertainty metric.

• Hardware Acceleration: The method leverages GPU parallelization (NVIDIA CUDA) to run multiple MCMC chains simultaneously, drastically reducing computation time.

3. Key Contributions

1. Application of MCLMC to Radiology: First application of the Microcanonical Langevin Monte Carlo sampler to distributed radiological source imaging and UQ.
2. Speed and Efficiency: Demonstrated that MCLMC converges orders of magnitude faster than state-of-the-art samplers like Hamiltonian Monte Carlo (HMC) and No-U-Turn Sampler (NUTS) in high dimensions.
3. Robustness to Overfitting: Unlike ML-EM, MCLMC does not require a user-defined stopping criterion; with sufficient data, it naturally converges to the true image without overfitting.
4. Superiority of GPP: Showed that Gaussian Process Priors significantly outperform independent Gaussian priors by leveraging spatial correlations, leading to sharper images and lower uncertainty.
5. Real-World Validation: Validated the method on both synthetic data and real-world aerial mapping data from the Johns Hopkins University Applied Physics Laboratory (JHU APL).

4. Results

Synthetic Data Experiments

• Convergence: MCLMC achieved good convergence ($\hat{R} < 1.1$, ESS > 200) with $\sim 10^4$ samples.
• Speed Comparison:
  • For $J=480$ voxels, MCLMC took 11.8 seconds to achieve a mean Effective Sample Size (ESS) of 112.
  • HMC required 5105 seconds to achieve a similar ESS (0.5 ES/s vs. 9.5 ES/s for MCLMC).
  • To match an ESS of 240, MCLMC took 12.9 seconds, while HMC took 474 seconds.
• Image Quality:
  • MCLMC with GPP achieved PSNR = 29.05 and SSIM = 0.95 with $T=100$ min measurement time.
  • This outperformed the Truncated Gaussian prior (PSNR = 25.42) and even surpassed the Gaussian prior with a much longer measurement time ($T=1000$ min).
• Uncertainty: The uncertainty maps (68% HDI) were sharper and more confined to source areas when using GPP.

Real-World Data (JHU APL Aerial Campaign)

• Dataset: Aerial mapping of a grid of Cs-137 sources ($J \approx 30,000$ pixels).
• Performance:
  • MCLMC with GPP yielded PSNR = 16.38 and SSIM = 0.66, significantly better than the Truncated Gaussian prior (PSNR = 15.38, SSIM = 0.35).
  • The reconstructed total intensity was within 14% of the ground truth.
• Scalability:
  • Despite increasing dimensionality from $J=480$ (synthetic) to $J \approx 30,000$ (real), computation time did not increase dramatically when using GPU parallelization.
  • GPU vs. CPU: Running 10 chains in parallel on an NVIDIA RTX 4090 reduced the time for 10,000 samples to ~3.7 seconds (second run), compared to ~16.6 seconds on CPU multi-threading.

5. Significance

• Near-Real-Time Decision Making: The ability to generate accurate radiation maps with full uncertainty quantification in ~10 seconds enables rapid decision-making in nuclear emergencies (e.g., path planning for responders, defining safe zones).
• Operational Feasibility: The method overcomes the "curse of dimensionality" that has historically prevented Bayesian UQ in large-scale imaging.
• Future Applications: The authors suggest this technology can drive autonomous path planning (moving toward high-uncertainty areas) and define data sufficiency criteria (stopping measurements when uncertainty drops below a threshold).
• Limitations: The method currently requires storing the full system matrix $V$ in memory, which can be challenging for very large scenes on consumer hardware, though running on smaller regions or lower fidelity can mitigate this.

In conclusion, the paper establishes MCLMC as a transformative tool for radiological imaging, offering a rare combination of high accuracy, robust uncertainty quantification, and computational speed suitable for real-world emergency response.