# Presentation: 585279

## Monte Carlo Methods

** **

***Jerome Spanier***
***Beckman Laser Institute & Medical Clinic***
***University of California, Irvine***

**Notes:**
I thank the organizers for the invitation; it’s a pleasure and an honor to participate in these lectures honoring Ely Gelbard.  He was a good friend and wonderful colleague for many years.  This tribute to Ely is richly deserved.
I’d also like to recognize another friend and colleague at Bettis – Al Henry – who passed away not long ago.  My relationships with Al & Ely are among the most treasured, both professionally and personally.

## 0. Outline

**Introduction**
**Organizing Principles**
**Historical Perspectives (****≤~1968)**
**Landmarks**
**State of the Art**
**Challenges**

- Neutronics & Shielding External Review September 21, 2000

## 1.  Estimate  by Monte Carlo

**Monte Carlo methods estimate averages – whether naturally stochastic or not.**
**                                 **
**	**
**	**

**               ****Figure 1:**
**    ****Buffon Needle Tossing**

**by tossing a needle of length ***** l*****   on a plane ruled by parallel lines a distance ***** d*****  apart**

**   ****Buffon (1777), Laplace (1886):**

- Neutronics & Shielding External Review September 21, 2000

**Notes:**
Buffon's needle

Here we have a plane with parallel straight lines drawn on it, the distance between two neighbouring lines being d. We throw at random a needle of length l (which is less than the space d) onto the plane N times and count how often the needle crosses one of the lines. Denote this count by n. One can prove that the probability of such crossing is 
q = 2l/(pi d). 
By the law of large numbers, n/N is close to q for large N, and therefore one can take 
2lN/dn 
as an estimate for pi. 
We have no circles here at all, but still the number pi can be estimated from this random experiment!   
Fig. 1 is a picture after N=50 needle tossings with d=1 and l=0.8. Here, n=25 crossed lines, so the estimated pi = 3.2, quite good for only 50 needles. 

## 1. Estimate e by Monte Carlo

**estimate      = average number of steps made by a random walking “particle” entering a unit length slab and making steps to the right uniformly distributed on [0,1]**
**                                                 **

**or**

- Neutronics & Shielding External Review September 21, 2000

**Notes:**
This problem, although extremely simple, is much closer to the conventional Monte Carlo uses in, e.g., nuclear applications.
(Use Math 369 notes for sketch of proofs)

## 0. Outline

**Introduction**
**Organizing Principles**
**Historical Perspectives (****≤~1968)**
**Landmarks**
**State of the Art**
**Challenges**

- Neutronics & Shielding External Review September 21, 2000

## 2.  Organizing Principles

**Integrate lecture material around four primary themes**

- generating sequences

- error analysis

- error reduction

- foundations/theoretical developments

**    ****Apology/Disclaimer: I have deliberately ignored **

**Hardware/software advances**

**Overly specialized techniques**

**Growing list of important (non-nuclear) applications**

**     ****But please forgive any other important oversights**

**The slide title ****n. subject****, n = 1,2,3,4,5,6 identifies the ****subject**** of this slide according to the chronology of topic ****n**** in the outline**

- Neutronics & Shielding External Review September 21, 2000

**Notes:**
I.e., we’ll move through the different chronological periods of our review discussing for each period these four main topics, in the same order each time.
It would be presumptuous to try to cover such a broad sweep of research over so many years in just an hour.  So I’ve been necessarily selective and the choice of topics will, no doubt, reflect my own idiosyncrasies and biases.

## 0. Outline

**Introduction**
**Organizing Principles**
**Historical Perspectives (****≤~1968)**
**Landmarks**
**State of the Art**
**Challenges**

- Neutronics & Shielding External Review September 21, 2000

## 3.  The Early Days (≤~1968)

**simulation + computers  =  success**
**Los Alamos focus on weapons provided impetus for rapid initial development (Fermi, Ulam, von Neumann, Richtmyer, Metropolis); rudimentary but effective computations in support of design; general lack of published work**
**post-WWII initiatives in support of nuclear navy, reactor design, etc. at BAPL, KAPL, ORNL, BNL, ANL; “hand calculations” (1955) ****→ Marchant, Friden → ENIAC, EDVAC → IBM 650 → 704 → 7040 etc. (1967); publication mainly in government reports**

- Neutronics & Shielding External Review September 21, 2000

**Notes:**
Reference for review article on history of computing.
Perhaps discuss brief brush with ENIAC, rooms full of marching Marchants at BAPL, etc. 

## Generating Sequences 

**Desiderata:  maximal periods, minimal “structure”**

**BUT “Random numbers fall mainly in the planes” (Marsaglia, 1968) exposed undesirable higher dimensional structure and led to growing concerns over the routine use of conventional pseudorandom sequences**** **

**“****Random number generation is too important to be left                         to chance,” R.R. Coveyou, 1969        **

**Pseudorandom = computer substitute for “truly random”; e.g., linear**
**congruential (LCG) (Lehmer, 1949) **

- Neutronics & Shielding External Review September 21, 2000

**Notes:**
S.K. Park and K.W. Miller. Random number generators: good ones are hard to find. Comm. ACM, 31:1192-1201, 1988. 
Sketch of basic theory about linear congruential generators.
Quote [101] Many multiplicative linear congruential generators are descendents of the infamous RANDU ( System/360 Scientific Subroutine Package, Version III, Programmer's Manual. IBM, White Plains, New York, 1968, p. 77.). This generator was first introduced in the early 1960s; its use soon became widespread ( examples are given ) . In retrospect RANDU was a mistake. The non-prime modulus was selected to facilitate the mod operation and the multiplier was selected primarily because of the simplicity of its binary representation. Research and experience has now made it clear that RANDU represents a flawed generator with no significant redeeming features. It does not have a full period and it has some distinctly non-random characteristics. Knuth calls it ``really horrible'' [67, p. 173,]. 

## 3.  Structure of LCG’s

**Figure 3: Hyperplane Structure of Pseudorandom *****N*****-tuples**

- Neutronics & Shielding External Review September 21, 2000

**Notes:**
Illustrating the linear structure of LCG’s in higher dimensions, ala Marsaglia.

## Error Analysis/Error Reduction

**ERROR ANALYSIS:**
**        ****Minimal information was output:**
**	**
**	**

**ERROR REDUCTION:**
**  ****Surge of interest in limitations of deterministic methods for nuclear reactor design in 1950’s, ’60’s gave rise to accelerated work on Monte Carlo.**
**  ****Initial emphasis on analog Monte Carlo and “last event estimation,” but shielding design (deep penetration) problems pointed to nonanalog sampling.**
**  ****Limited use of exponential transform and roulette/splitting, survival biasing.**

**and little else**

- Neutronics & Shielding External Review September 21, 2000

**Notes:**
More about the exponential transformation on the next slide.

## Foundations/Theoretical Developments

**Relative paucity of published work**
- Kahn’s handbook (1954)
- Cashwell & Everett’s practical manual (1959)
- Hammersley & Handscomb’s book (1964)
**See excellent review by Kalos, Nakache and Celnik in *****Computing Methods in Reactor Physics,***** Gordon & Breach, 1968).**

**                              ****Early Impressions, Concerns:  **
**Shielding design issues generated a critical need for departure from analog sampling and spawned use of exponential transform (uses exponential importance function; tries to maintain spatially constant particle population with exponentially falling weights from source to detector).  But ***effective*** bias can result from undersampling certain regions, rare occurrences of high weight particles.**

- Neutronics & Shielding External Review September 21, 2000

**Notes:**
Add complete references

## Theme

***Prevailing theme of this early period:***
**Monte Carlo is a method of last resort **
**(diffusion theory still dominated reactor computations)**

- Neutronics & Shielding External Review September 21, 2000

**Notes:**
E.g., when I came to BAPL in 1955 essentially all of the numerical work relied on diffusion theory.

## 0. Outline

**Introduction**
**Organizing Principles**
**Historical Perspectives (****≤~1968)**
**Landmarks**
**State of the Art**
**Challenges**

- Neutronics & Shielding External Review September 21, 2000

## Generating Sequences

**Two important developments characterize past 30 years:**
**exhaustive, rigorous analyses of linear generators and the emergence of modern nonlinear congruential, inversive congruential, shift-register algorithms and extensions (H. Niederreiter, *****Random Number Generation & Quasi-Monte Carlo Methods,***** SIAM, 1992)**

**Figure 4: Inversive Congruential**

- Neutronics & Shielding External Review September 21, 2000

**Notes:**
Inversive Congruential Generators 
Sometimes the Parallel Hyperplanes phenomenon inherent in LCGs may cause adverse effects to certain simulation applications because the space between the hyperplanes will never be hit by any point of the generator, and the simulation result may be very sensitive to this kind of regularities. Inversive Congruential Generators (ICG) are designed to overcome this difficulty. They are variants of LCGs introduced by Eichenauer and Lehn in 1986. Although the extra inversion step eliminates Parallel Hyperplanes (see Figure 3), it also changes the intrinsic structures and correlation behaviors of LCGs [5]. ICGs are promising candidates for parallelization, because unlike LCGs, ICGs do not have long-range autocorrelation problems. However, with current state of the art, ICGs are substantially (8X) slower than LCGs due to the inversion process. ICGs have the same period as their LCG counterparts.

## Generating Sequences

**     ****the parallel development of quasirandom, or low discrepancy, sequences (finally a way around the assumption of randomness of nonrandom quantities): ignore randomness; concentrate on *****uniformity over *****[0,1]****s****, s = “dimension” of problem**
**      ****the latter constitutes a paradigm shift whose foundations lie in number theory (theory of uniform distribution mod 1) (see Niederreiter book cited above)**

**Figure 5: Halton Pairs**

- Neutronics & Shielding External Review September 21, 2000

**Notes:**
The fundamental idea behind quasirandom (also called low discrepancy) sequences is very different.  All pretense at random behavior is abandoned in favor of greater uniformity in higher dimensions.  The underlying theory is number theoretic.

## 4.  Pseudorandom vs. Quasirandom

**Comparison of 2000 pairs: **

**                                ****Figure 6:**
**Pseudorandom vs. Low Discrepancy Sequences**

- Neutronics & Shielding External Review September 21, 2000

**Notes:**
Side by side comparison accentuates the differences between pseudo- and quasi-random pairs.

## Radically Different Error Analyses

**(Koksma, 1942; Hlawka, 1961) for integrals and similar results have been established for the transport equation (Chelson, 1976; Maize, 1980; Spanier-Li, 1998), which is shown in a rather natural way to be an infinite dimensional extension of the finite dimensional quadrature problem.**

** ****For pseudorandom implementations,**
**  ****increasing focus on confidence interval estimation**

** ****where *****L******N,****** U******N*****  are determined by *****N, m******N****** ,s******N****** a ,b*****, **
** ****has  led to use of higher moments (VOV), multiple statistical tests, as in MCNP.**

** ****When low discrepancy sequences are used as generators, *****deterministic *****error bounds are required:**

- Neutronics & Shielding External Review September 21, 2000

**Notes:**
There is no longer a probabilistic/statistical basis for analyzing the limiting behavior of quasirandom sequences so new error analysis methods have to be developed.  These were established for finite dimensional integrals by Koksma (1942) and Hlawka (1961 and this theory was extended to integral equations by Paul Chelson (1976) and Earl Maize (1980) in their dissertations.  See also Spanier-Li (1998).

## 4.  Comparison of Convergence Rates

**                                                                      ****Figure 7:**
**       ****Comparison of Pseudorandom, Quasirandom (*****s***** =1 & 2) Convergence Rates**

**convergence rate (Halton, 1960; Faure, 1981, 1982; Niederreiter, 1987) but  error bounds**

- Neutronics & Shielding External Review September 21, 2000

**Notes:**
Here we see that the slope of the pseudorandom line is –1/2, the limiting slope (for s = 1) of the quasirandom line is –1 and the higher dimensional quasirandom sequences lie (mostly) between, depending on the relative sizes of s and N.   Note the scale of the horizontal axis; most practical MC simulations are crowded into the left part of the graph. 

## 4.  General Error Reduction

**GENERAL PURPOSE METHODS:**
**Maturation of theoretical importance sampling:  Goertzel (1949), Kahn (1954), Goertzel & Kalos (1958), Kalos (1963), Coveyou (1967), Spanier & Gelbard (1969), use of “intelligent” random numbers Booth (1986).**
**Most recent development:  adaptive importance sampling and related techniques for achieving  *****geometric convergence*****.**
**      ****The idea behind adaptive methods is:  Each random walk simulated adds information.  How can this information be used to improve the estimates steadily?**
**      ****Alter sampling/weighting strategy after a block of  *****N*****  histories ( = “stage” in a multistage algorithm) has been generated; expect ***** n*****th  stage estimates to improve on those from stage****  *****n-1.***

**The past three decades have produced a remarkable growth in both general and special purpose methods for reducing error in Monte Carlo and quasi-Monte Carlo calculations.**

- Neutronics & Shielding External Review September 21, 2000

**Notes:**
Booth’s 1986 paper in NSE discusses the differences between use of intelligent random numbers and conventional biasing/weighting done in the random walk space.  He also discusses some of the potential advantages of performing the learning in the random number space.

## 4.  Exponential Error Reduction Rates

**I.e., goal is to achieve**

***n***** = stage number**

- Neutronics & Shielding External Review September 21, 2000

**Notes:**
This initiative was stimulated by an interest, at LANL, in greatly increasing the efficiency of MCNP.

## 4.  Context for Adaptive Methods

**                          ****Figure 8:**
**Illustration of Geometric Convergence**

**(****Kong, 1999; Spanier, Lai & Spanier, Kong and Spanier, Hayakawa & Spanier, MCQMC98)**

- Neutronics & Shielding External Review September 21, 2000

**Notes:**
Add full references.  Note dramatic effect of geometric convergence.

## 4.  Geometric Error Reduction

**Figure  9:  Geometric Convergence for 1D Model Problem**

- Neutronics & Shielding External Review September 21, 2000

**Notes:**
We’ve developed three distinct algorithms that achieve geometric convergence for global transport problems.  One of these is based on the use of correlated sampling, a second on the use of (unbiased) importance sampling, and a third, generalizing the second, on a class of (generally) biased importance sampling estimators. Time does not permit delving into them further here, but each depends upon learning in some form, of course.  Fig. 9 compares the behavior of these three methods for a simple 1-D model transport problem.

## 4.  Practical Error Reduction

**CONTROL OF WEIGHTS**
**One can use analytically correlated estimates of variance to select nearly minimal variance parameter values from a parametrized family of sampling schemes (Spanier, 1971).**
**A very effective technique incorporated into MCNP is the use of weight windows, which combines Russian roulette and splitting in order to control both the density of particles throughout phase space and their weights and avoids catastrophic contributions to variance in that way.**
**For quasi-Monte Carlo methods, there has been a parallel development of *****variation***** reduction methods that are adapted from v*****ariance***** reduction procedures established for conventional pseudorandom Monte Carlo.**

**and**

- Neutronics & Shielding External Review September 21, 2000

**Notes:**
Roughly speaking, the idea in my 1971 paper is to define a low dimensional subspace of the (infinite dimensional) space of all importance sampling transformations (e.g., use of exponential transformations for deep penetration problems) and then choose the best among these using correlated estimates of the variance (considered as a function of these parameters).
References to MCNP for weight windows.
Parallel between variance reduction for pseudorandom MC and variation reduction for quasi-MC.

## 4.  Next Event Estimation

**Another effective method is to use information-rich **
**estimators (improved scoring techniques)  ********  getting more**
** ****information out of each random walk *****without altering the walk   probabilities*****.**

**This idea was applied very productively recently and led to creation of  a special  version of MCNP for estimators useful in well-logging applications (Mosher et. al., 2002).**

- Neutronics & Shielding External Review September 21, 2000

**Notes:**
Add full references

## 4.  Use of Superposition 

**The first of two other methods that have been successfully used to convert one transport problem to an equivalent problem easier to solve makes use of the linearity of the transport operator**
**                                                          **

- Neutronics & Shielding External Review September 21, 2000

## 4.  Illustration of Superposition 

**Figure 10: 2 Region Illustration**

- Neutronics & Shielding External Review September 21, 2000

**Notes:**
The definition of QII makes the 2-region problem source constant and isotropic, hence the sum of the fluxes is also constant and isotropic (it’s easy to see that this provides a valid solution to the equation for the sum; by uniqueness, this is the only solution.

## 4. Use of Reciprocity

**Simulation of the adjoint transport equation has been shown to be very **
**effective for estimating reaction rates in small volumes (even at a point!)**

**        ****The second method relies on the reciprocity principle for **
**        ****the transport operator (Maynard, 1960):**

- Neutronics & Shielding External Review September 21, 2000

## 4.  Illustration of Reciprocity 

**      ****Superposition and reciprocity can lead to very different simulations and results, each with its own advantages and disadvantages**

**      ****The smaller is region III, the**** greater is the advantage of reciprocity.**

**Figure 11: 3 Region Illustration**

- Neutronics & Shielding External Review September 21, 2000

**Notes:**
The small red region could identify a local “hot spot” of an operating reactor (cell geometry).

## Other Error Reduction Methods

**Perturbation techniques have been valuable in estimating small differences accurately by Monte Carlo methods.  **
**Methods have been devised for obtaining highly correlated estimates of small differences, and this idea can also be used to provide estimates of rates of change (sensitivity coefficients).  **
**When these are coupled with optimization algorithms, an effective Monte Carlo inverse problem-solving algorithm can result (Hayakawa et. al., 2001).**

- Neutronics & Shielding External Review September 21, 2000

## Illustration of Inverse Problems

**      ****Figure  12: Illustration of Inverse Problem**

- Neutronics & Shielding External Review September 21, 2000

**Notes:**
Fig. 12 illustrates use of pMC (and dMC) plus optimization to solve a rather general class of inverse problems.

## Foundations/Theoretical Developments

**Prompted by need for improved understanding of puzzling and counterintuitive results:**
**1. erratic behavior of sample means**
**2. possibility of theoretically infinite variance importance   sampling strategies that seem reasonable otherwise).  **
**Role of “model problem” calculations in debugging and enriching understanding of transport processes**.

- Neutronics & Shielding External Review September 21, 2000

## 4.  Other Theoretical Advances

**Theoretical connections of transport problems, matrix problems and finite dimensional quadrature problems (Spanier & Gelbard, 1969, Spanier, 1995, Spanier & Li, 1998).**

**Derivation of coupled integral equations for expected values of higher moments of estimators (Amster & Djomehri (1976), Lux & Koblinger (1991).**

- Neutronics & Shielding External Review September 21, 2000

## 0. Outline

**Introduction**
**Organizing Principles**
**Historical Perspectives (****≤~1968)**
**Landmarks**
**State of the Art**
**Challenges**

- Neutronics & Shielding External Review September 21, 2000

## Generating Sequences 

**Pseudorandom sequences still dominate traditional applications but low discrepancy and hybrid mixtures of pseudorandom and low discrepancy sequences are gaining popularity (especially for financial modeling, computer graphics).**

- Neutronics & Shielding External Review September 21, 2000

## Comparison of Generating Sequences

**Figure 13:  Mixed, Quasi- and Pseudo-random Sequences**

- Neutronics & Shielding External Review September 21, 2000

**Notes:**
Fig. 13 compares plots of pairs of mixed, quasi- and pseudo- random sequences.  The differences here are hard to detect but they can produce substantially different results in higher dimensions.

## Comparison of Convergence Rates

**                             ****Figure 14: **
**    ****Comparison of Mixed, Quasi- and Pseudo Convergence for Different Dimensions**

- Neutronics & Shielding External Review September 21, 2000

**Notes:**
Fig. 14 provides dramatic evidence of the advantage of using mixed (some quasi- components, some pseudo- components) sequences for high dimensional problems.  The upper left graph is for a 40_D problem treated either using pseudorandom numbers in all 40 dimensions, or quasirandom numbers in all 40, or half in each (with those chosen for quasirandom chosen judiciously).  The lower left graph is mixed (100, 40) and the right graph is mixed (200, 40).

## Theory

**Equivalence between analytic model (transport equation) and either a probabilistic model  (for pseudorandom simulation) or a deterministic model (quasirandom simulation) is beneficial. **
**Utility of coupled integral equations for expected values of higher moments of estimators (Amster & Djomehri, 1976, Lux and Koblinger, 1991) is well understood.  **
**Applying these models and derivations to achieve optimal performance is still as much art as it is science.**

- Neutronics & Shielding External Review September 21, 2000

## 6.  Challenges for the Future

**Optimum use of pseudorandom and quasirandom sequences for transport problems and very high dimensional integrals**

** ****Complexity analysis of Monte Carlo; optimality with respect to    **
** ****complexity theory in various settings**

**Monte Carlo methods for inverse problems**

** ****Optimization of adaptive, zero variance algorithms**

**IT SHOULD BE A GREAT CENTURY!**

- Neutronics & Shielding External Review September 21, 2000