Automated variance reduction

[vitmog]

Hello all,
This issue is related to the following topic mentioned in "EGSnrc development project ideas #790":

Deep learning for Monte Carlo of radiation transport

• Variance reduction parameter selection

I am currently working on the design of a portable open source C++ library for adaptive variance reduction techniques extending some efficient techniques already implemented in a shielding Monte Carlo Fortran code that I maintain. The library aims to provide a concise API for collecting information about Monte Carlo sampling by sending current particle tracking and scoring data (aka "learning") and returning the optimal parameters of variance reduction techniques on request. I would like to investigate the possibility and relevance of linking it to the EGSnrc family codes and adjust the method conception at an early stage if needed.

While I am familiar with neutral particle transport Monte Carlo codes organization principles and their nuclear reactor and radiation shielding applications, I am less familiar with the accelerator and medical physics areas of EGSnrc applicability. Therefore, it would be highly appreciated if the principal developers could clarify several points:

1. What are the most relevant codes from the EGSnrc family where variance reduction is desirable?
   EGSnrc is a quite comprehensive system, so I plan to dive into a specific code for investigation of possible linking issues. From my study of the documentation, it seems that it could be BEAMnrc implementing the split/roulette and implicit capture techniques, but I can be wrong.
2. Is calling external C/Fortran routines from the Mortran code supported?
   I am not experienced with Mortran, but since Mortran is just a preprocessor it seems to be possible to call routines of a specially designed wrapper, as a part of EGSnrc, from specific tracking and scoring points.
3. Where can I find real-world problem examples for variance reduction efficiency investigation and debugging?

Also, it would be helpful to know thoughts of community members about potential utility of automated variance reduction in EGSnrc. For example, radiation shielding problems of nuclear power facilities and spent nuclear fuel and radioactive waste transportation, as a rule, cannot be solved in the analog Monte Carlo mode due to very high attenuation. Therefore, these calculations are always performed using variance reduction techniques which provide several orders of magnitude computational gain. What is the status of variance reduction in treatment planning, accelerator physics, and other problems for which EGSnrc is used?

Thanks for your time.

Best regards,
Vitaly Mogulian

[ftessier]

Thank you for your interest in this topic!

On your specific questions

1. Relevant applications for variance reduction (VR):

   The most relevant EGSnrc apps for VR research are probably BEAMnrc and egs_chamber. The broad classes of VR methods are indeed particle splitting and roulette, with notable examples such as Directional Bremsstrahlung Splitting and Range-Based Roulette, which can yield efficiency improvements of several orders of magnitude, as you said. EGSnrc apps also include more specialized schemes, such as Temporary Phase Space Storage and Cross-Section Enhancement.

   Note that EGSnrc does not have dedicated VR methods for high-attenuation scenarios, which remains a limitation for shielding problems. For exploring adaptive VR concepts, I recommend starting with a minimal egs++ application built from the egs_app stub. This allows full control and clarity while testing your proof-of-principle implementations.

2. Calling external routines:

   Modern EGSnrc applications (such as egs_chamber) are written in C++ using the egs++ class library. These applications interface with the Mortran-based physics engine but handle geometry, sources, and VR primarily on the C++ side.

   While it is technically possible to call external C or Fortran routines from the Mortran backend, I don't recommend it. Development in C++ via the egs++ framework provides much better flexibility and integration. The long-term direction for EGSnrc development is towards re-factoring in modern languages (such as C++), so at any rate following the egs_chamber model is the most future-proof approach.

3. Examples for VR efficiency and debugging:

   For real-world examples, refer to the EGSnrc manual sections on VR, and to the following key papers, which illustrate how VR techniques are implemented and benchmarked within EGSnrc:

   • Directional Bremsstrahlung Splitting (Med. Phys. 31, 2004)

   • egs_chamber: A C++ class library for ion chamber simulations (Med. Phys. 35, 2008)

On the potential utility of adaptive or automated VR

In shielding problems, adaptive VR would be a significant win! In medical physics and industrial radiation processing applications, VR methods tend to fall in generic classes, for example electron beam bremsstrahlung generation, or dose scoring in small volumes within large phantoms.

Our experience is that the currently implemented VR techniques have efficiency profiles that exhibit broad, smooth peaks around their optimal parameters. Because of this, we typically tune VR parameters by running a few short simulations on a logarithmic scale to identify the optimal parameter space. An adaptive VR system would need to outperform this manual tuning process to yield practical advantages.

That said, your proposed research direction is valuable; it would provide insight and perhaps advances in optimizing Monte Carlo radiation transport simulations.

[vitmog]

@ftessier Thank you very much for such a comprehensive explanation! It would have taken me a lot of time to obtain just a small piece of the information that you provided here. I am really happy that this topic is relevant!

I will prepare the response a bit later, after diving into the details.

[vitmog]

Thank you again for the useful information and insights!

On specific questions

1. Relevant applications for variance reduction (VR):

The egs_chamber app and its target problems seem quite interesting to me and the code is written in C++. I will focus on the egs++ applications. Thank you!

2. Calling external routines:

If I realize it right, the the_stack handling is organized by the LIFO (Last In, First Out) straightforward principle that is desirable. Also, all the_stack handling is on the C++ side, so there is no internal Mortran manipulation with the_stack records. Consequently, any modification of the Mortran side of the code will not be required anyway, so it will be further considered as a black box.

Perhaps it will even be not needed to modify the egs++ side, because:

1. Both EGS_SimpleApplication and EGS_AdvancedApplication classes provide the entry point augsab() which allows to collect particle tracking data and scoring, and crucial history events are implemented right in the overridden methods too;
2. Probably, it is possible to pull information about completed track branching events between two neighboring augsab() callings from the relative change of the the_stack array if it is not possible to do from augsab();
3. At least a part of the VR techniques is implemented on the top level part of the code, such as the egs_chamber app.

I will likely use the tutor4pp example to learn more about handling the derived particles storage the_stack when I will manage to compile run the apps.

3. Examples for VR efficiency and debugging:

Thank you! My representation of possible VR practical applications to particle transport became wider than before. I am also glad that there is no conflict with my library concept.

On the potential utility of adaptive or automated VR

Problems like the specific small volume-averaged dose calculation with the cascades of particles of diverse types seem to me to be very interesting, because they provide permanent demand on the VR techniques with charged particles unlike neutral ones.

On the outperforming standard schemes

If it is enough just to hit the right order of magnitude with the manual parameter selection in practice, then the optimization of exactly that parameter cannot outperform it significantly, but can be used for the control of the implementation. However, I suppose that just increasing the number of the variance reduction entities with individual parameters, whose selection is not a problem when the scheme becomes adaptive, should already provide the desired gain. Also, the library concept will provide implementation of almost fully-automated schemes with large number of such parameters.

UPD. 21.01.26. Actually, the Mortran side manipulates the_stack, but I still hope that it is possible to recover enough information about branching in augsab().
UPD. 22.01.26. According to PIRS-701, the LIFO stack can be re-sorted on the Mortran side using a special user-defined macro that makes the linking more problematic. It will be considered as an incompatible feature for a while.
Finally, I see a possibility to prepare a minimal proof-of-principle example as an user app, either new one or an existent app modification. Of course, it is still possible that I may meet some unexpected troubles at the implementation stage. In the case of success, designing of the automated VR schemes noted above would obviously require strong support, if not guidance, of the physicists/maintainers who have deep intuitive representation about the target physical systems, the understanding of code internals and user habits.

[ftessier]

Looking further ahead, I can imagine eliminating VR parameters entirely via a scheme I'll call retrospective splitting. The insight is that the simulation already knows when scoring events occur, so it should be able to concentrate efforts in the relevant regions of phase space. In terms of integrating the Boltzmann transport equation, the simulation would automatically allocate more samples around tracks that contribute to the score. Visually: when a particle history reaches the scoring region, the simulation would "rewind" and enrich the upstream cascade to extract more information from that productive trajectory.

The challenge is to do this without introducing bias. For example, resampling the last interaction multiple times leads to oversampling with no obvious way to adjust the weights. The solution is to split particles upstream in the cascade. Since splitting a particle into $n$ copies (each with weight $1/n$) is always valid at any point in the simulation, one can retroactively apply this operation to ancestors of a scoring particle! The point is that you don't predict which particles will contribute to scoring; you react to scoring events by enriching the trajectory that led there. This is reminiscent of Monte Carlo adjoint-based importance sampling, but carried out on the fly per track, so to speak, without requiring a separate adjoint calculation.

Implementing this would require storing the full cascade history (positions, directions, energies, interactions, perhaps even rng states!) to enable "rewinding" to ancestor particles. That's a lot of bookkeeping, but otherwise would not significantly complicate current variance estimators if carried out per history. I wouldn't want to shoehorn this into the existing EGSnrc codebase! 😨 But I trust one could demonstrate the idea in a minimal egs++ application with a small purpose-built state machine.

Of course, this doesn't eliminate the optimization question. You'd still have parameters like recursion depth and splitting ratios to choose. But my hunch is that these would be more problem-agnostic than current VR parameters; there might emerge a reasonable universal heuristic for "how far back to split," for example. And since the simulation can track its own efficiency in real time, could it potentially tune these on the fly, walking toward an optimum as it runs?

Just thinking out loud! 😃

[vitmog]

If I got it right, the described idea of retrospective splitting is potentially valid technically, but exactly this splitting scheme would provide a biased estimator, because it will systematically underestimate target functionals due to reducing higher already sampled scoring depending on its values. That is, we divide the score by $n$ and, instead of the remaining part $(n-1)/n$, we use something different. If $n$ was permanently fixed, it would be OK, but the proportionality to the scoring value will affect the means. Thus, the choice of the splitting parameter must be made before the sampling. Alternatively, it could be implemented using something like such a test branch(-es) which must be always excluded from the statistics.

A scheme, quite close to what you described, can be implemented using diverse applications of Markov Chain Monte Carlo (MCMC) with the Metropolis-Hastings (MH) algorithm for sampling arbitrary distributions, closer to the optimal ones, using dependent samples. This family of methods allows us to do the desirable -- to mutate (resample) an existing sampled history depending on its scoring (retrospective splitting is an example of such a mutation). It is not difficult in final algorithmic implementation, but to substantively describe it in a few words is not easy (actually harder). A famous example of a MCMC application to particle transport is the Metropolis Light Transport (MLT) approach, but now there are lots of works in the physically based rendering. The idea mentioned by you with keeping random variables is used too as replacing of phase-space coordinates with, so called, Primary Sample Space coordinates which are the random variables used for history sampling. This is a very large topic!

I also made some experimental implementations of promising schemes (unlike MLT) 1.5-2 years ago for neutral particle transport to study it and test the fundamental possibility of obtaining the efficiency gain. I really see diverse benefits of such schemes application with neutral particles, but I am currently not sure about the charged ones.

I absolutely agree with all that you wrote about the bulkiness of working with entire existing history, especially with the general nuclear physics applications rather than the nuclear power engineering ones. For this reason, my idea is based on implementation of a scheme which just returns the estimated parameters of traditional VR techniques and does not require changing of algorithms that provide its portability, as expected. It can be used in the single-stage (on-line) mode, but also allows restarting from previously generated data. I hope that preparation of the concept description in the form of a preprint will not take too much time, but I am also doing this in parallel with coding and studying.

Regarding the optimal splitting parameters, their definition was formulated in a publication of Ogibin in 1967 and proved by Medvedev not later 2011. Roughly speaking, the optimal splitting rate equals to square root of the ratio of the expected variance of further scoring at this point to the one at the previous point. You can find the full definition in Section 2.3.5 (Equation (32)) of my working paper. I am sorry for its (and this) language and know about many other issues which I have planned to fix in the further release. Note that the optimal parameter formula is recurrent and does not depend on the chain size, i.e. all is much simpler than it might seem!

In addition, the currently proposed method can be considered as a generalization of the adaptive splitting used in the calculations of the working paper, but it is not described there.

[vitmog]

I would also to clarify the idea about MCMC schemes application to the particles transport which is mentioned above. I meant that the MCMC schemes for charged particles can be efficient, but I am not sure that the efficiency will be higher than simply optimized splitting scheme like that I am working on. In contrast, MCMC can provide cardinal improvement with neutral particles, because possible mutations strategies for them seems more flexible. However, this topic requires development anyway.