Introduction to Numerical Integration

• The problem of finding I = ∫_{a}^{b} f(x) g(x) dx is closely related to both differentiation rules and to formulae for expanding f(x) in [a, b]. If the expansion functions are nice, we can plug in expansion for f and do term by term integral. It divides into

• NI.1,2: f(x) assumed known at prescribed points (cf. Lagrange interpolation)

• NI.3: We can choose points where f(x) to be calculated (cf. Chebyshev expansion for interpolation)

• NI.5: Monte Carlo Integration

• NI.5: is most important as it is what you get in big problems with multidimensional integrals!

• NI.4: has remarks on use of parallel computers. A good compendium of rules will be found in M. Abramowitz and J. Stegun, “Handbook of Mathematical Functions.”

NI.1: Newton-Cotes Formulae

• (i) We write

I = ∫_{a}^{b} f(x) g(x) dx  (NI.1)

where f(x) is assumed known at x_{0} … x_{n} and g(x) is sufficiently simple so that integrals (NI.3), involving it, can be calculated. Of course, we split off g(x) ≠ 1 to make f(x) smoother and more reasonably approximated as a polynomial.

• We assume the form:

I = n ; ∑ ; r=0 w_{r}f(x_{r})  (NI.2)

where the n + 1 w_{r}’s are to be found by requiring I to be exact for all polynomials of degree ≤ n.

• This is easy-for put

I_{m} = ∫_{a}^{b} x^{m}g(x) dx  (NI.3)

which are calculable constants.

∫_{a}^{b} f(x)dx = _{i}^{∑} w_{i}f (x_{i})

Linear Equations for Coefficients

• Our condition becomes

I_{m} = n ; ∑ ; r=0 w_{r} x_{r}^{m} : 0 ≤ m ≤ n  (NI.4)

(linear nature of integral ensures that (NI.4) guarantees exactness of (NI.2) for all polynomials of degree ≤ n).

• We can write (NI.4) as:

A\underline{w} = \underline{b}  (NI.5)

\underline{w} is column vector of w_{r}’s, \underline{b} is column vector of I_{m}’s and A is (n + 1) × (n + 1) matrix of coefficients A_{mr} = (x_{r})^{m}. It is easy to show that A is nonsingular (if x_{r} distinct) and hence that (NI.5) can be solved and (NI.2) established for any g(x) in (NI.1).

• Usually (NI.5) is solved—once and for all for given g(x)—and then (NI.2) can be used quickly for any f(x) that crops up.

Examples of Newton Cotes Rules I Examples of Newton Cotes Rules II

Formulae for Examples: Trapezoidal and Simpson’s Rules

• (a) Trapezoidal Rule n = 1, g(x) = 1 x_{0} = a, x_{1} = b

∫_{a}^{b} f(x) dx ≈ \frac{h}{2}[+f(a) + f(b)] h = b − a (NI.6)

• (b) Simpson’s Rule

n = 2, g(x) = 1 x_{0} = a, x_{1} = \frac{1}{2}(a + b), x_{2} = b

∫_{a}^{b} f(x) dx = \frac{h}{3}[f(x_{0}) + 4f(x_{1}) + f(x_{2})] (NI.7)

where h = \frac{1}{2}(b − a)

Summary of Newton Cotes Rules

• Both these are equi-spaced (in x_{i}) rules with g(x) = 1 and the generalization to arbitrary n are called Newton-Cotes formulae (see Abramowitz and Stegun to find them explicitly written out).

• We can relax restrictions:

• (a) x_{i} equidistant

• (b) x_{0} = a and/or x_{n} = b

• (c) g(x) = 1

• The reader need only solve (NI.5) to find relevant formulae for any n.

• You will find in the books, discussions of so called open formulae where g(x) = 1 as above, but one omits one or both of the end-points. The remaining n or n − 1 points (x_{1} … x_{n} or x_{0} … x_{n−1} or x_{1} … x_{n−1}) are still equidistant.

Errors in Newton Cotes Formulae

• (ii) The errors in (NI.6,7) are written in the form:

Trapezoidal—(NI.6) −h^{3} f′′(ξ)/12
Simpson—(NI.7) −h^{5} f^{iv}(ξ)/90  (NI.8)

where ξ is some value of x in [a, b].

• Note that, as expected, trapezoidal formula is exact for linear functions (error ∝ f′′(ξ) vanishes for f(x) = αx + β).

• But Simpson’s rule—which as derived is exact for quadratic f(x)—is unexpectedly exact for cubics.

• In other words, its error is particularly small and it is a widely used formula in its extended form:

∫_{x_{0}}^{x_{2n}} f(x) dx = \frac{h}{3}[f_{0} + 4f_{1} + 2f_{2} + 4f_{3} + …
+ … + f_{2n}] − \frac{nh^{5}}{90}f^{(iv)}(ξ) (NI.9)

where f_{i} = f(x_{i}) and x_{i} are equidistant.

• One gets this formula by applying Simpson to intervals x_{0} to x_{2}, x_{2} to x_{4} … x_{2n−2} to x_{2n} in succession.

CPS615NI G.C.Fox, Oct 95  8  gcf@npac.syr.edu

Use of High Order Newton Cotes

• One could also generate a formula in terms of the same function values f(x_{i}) by using the Newton-Cotes formula appropriate to 2n + 1 points.

• This gains an extra derivative in error just like Simpson, and error is ch^{2n+3} f^{2n+2}(ξ) where c is some known number.

• In practice, this is not often used because

• (a) the implied interpolation of f(x) over a wide range by a high order polynomial can produce disastrous results when f(x) not polynomial (e.g., has similarities in its value \underline{or} derivatives).

• On the other hand, local interpolation by f(x) by a low order polynomial (used in (NI.9)) is much safer;

• (b) it is much easier to estimate error in iterated Simpson of equation (NI.9) as we will describe in Section NI.2.

Strategies for Manipulating Integrals

before using

Standard Numerical Integration

• (iii) Points of Technique

• Note, when evaluating I = ∫_{a}^{b} f(x) g(x) dx, one can choose between two strategies

• (a) Change of variable. For example, dy = g(x) dx converts integral to g(y) = 1 form.

• (b) Newton-Cotes (or Gauss—see NI.3) formula with weight function g(x).

• A third strategy is to change the integration range (e.g., to take b from ∞ to finite value).

• How do we choose between strategies and how to implement them?

Criteria for Choosing Manipulation

• In each case, one must simply use criterion.

• Is f(x) or f(y) better approximated by a polynomial?

• One can use same criterion both to choose g(x) and/or subtract off some particularly nasty part (e.g., a singularity) of f(x).

• For example

I = ∫_{a}^{b} f(x) g(x) dx
= I_{1} + I_{2}
where f(x) = f_{1}(x) + f_{2}(x)
and I_{i} = ∫_{a}^{b} f_{i}(x) g(x) dx

• where f_{1}(x) may be well approximated by a polynomial and I_{1} = ∫_{a}^{b} f_{1}(x) g(x) dx could be done numerically.

• f_{2}(x) may be a nasty function to integrate numerically (e.g., x^{α} with α non-integer) but chosen so that we can do I_{2} = ∫_{a}^{b} f_{2}(x) g(x) dx analytically.

Dangers in Integrating

Non-Analytic Functions

• (iv) Further Examples of Practical Technique

• a) I = ∫_{0}^{x} x^{α} exists for α > −1. Namely, I = \frac{x^{α+1}}{α+1}, but f(x) = x^{α} is not a polynomial unless α is an integer.

• In typical integration formula

I = \frac{h}{3} [f(0) + 4f(h) + 2f(2h) + …] + 0(h^{4})

The error is 0 (h^{4}f^{(iv)}).

• Clearly this numerical procedure doesn’t work if α < 0 even though integral well defined as −1 < α < 0 ⇒ f(0) infinite, but I well defined.

• We could use an “open” formula which does not involved f(0) but even then, errors will be large. For example ∫_{0}^{h} f(x) dx has error something like h^{5} f^{(iv)}(ζ) where ζ in range 0 ≤ x ≤ h. If f(x) = x^{α}, f^{(iv)}(ζ) ∼ ζ^{α−4}, ζ ∼ h implies that error is h^{1+α} for each interval.

Why Do We Have Problems with

Non-Analytic Functions?

• What is the problem?

• f(x) = x^{α} is just not well approximated by a polynomial. Note this is true even when α is > 0, but not an integer.

• For example

f(x) = x^{2} , f^{(1)}(x) = x , f^{(2)}(x) = 1
f^{(3)...} = 0 ⇒ errors small (actually zero)

• But for α = 1/2

f(x) = \sqrt{x}
f^{(1)} = \frac{1}{2\sqrt{x}}
f^{(2)} = \frac{−1}{4x^{3/2}} ⇒ errors large

• 1/\sqrt{x} is certainly a dreadful polynomial, but \sqrt{x} is not good either as derivative is 1/\sqrt{x} and so not well behaved.

Combining Analytic and

Numerical Integration

• (b) Well, we don’t need numerical integration to perform ∫ f(x) dx with f(x) = x^{α} as we can do analytically.

• However, what about α = ±1/2 examples.

f(x) = \frac{1}{\sqrt{x}} exp (−x^{3})
or = \frac{1}{\sqrt{x}} + exp (−x^{3}) .

• There are three possible strategies which in a practical case might need all to be combined.

• (c) Strategy I for integrating non-polynomial approximatable functions.

f(x) = f_{1}(x) + f_{2}(x)

• where for example

f_{1}(x) = exp (−x^{3}) f_{2} = \frac{1}{\sqrt{x}}

• Do f_{1}(x) numerically

• Do f_{2}(x) analytically

Cunning Choice of Weight Functions

• Strategy II is used when difficult function is multiplicative, not additive.

• Suppose

I =  ∫ f(x)dx where
f(x) = f_{1}(x) f_{2}(x) so that
I =  ∫ w(x) f_{2}(x)
with f_{1}(x) viewed as a weight function w(x).

• Use f_{1}(x) as a prescribed weight function to derive special integration formulae for doing ∫ w(x) polynomials. We described these earlier.

• This strategy applies to the case f(x) = \frac{1}{\sqrt{x}} exp (−x^{3}) with w(x) = 1/\sqrt{x}.

Change of Variable Before

Numerical Integration—I

• Strategy III

• We can also use change of variables in ∫ f(x) with f(x) = \frac{1}{\sqrt{x}} exp (−x^{3}).

• Namely, if

\frac{1}{\sqrt{x}} dx = dy
y = 2\sqrt{x}
\frac{1}{\sqrt{x}} exp (−x^{3}) dx = ∫ exp (− (\frac{y}{2})^{6}) dy

• which can be done reliably numerically by approximating integrand by a polynomial, as in standard formulae.

Change of Variable Before

Numerical Integration—II

• Let us consider implications of the transformation x → y = 2\sqrt{x}.

• Integration formulae for x (which assume integrand f(x) ∼ polynomial) ⇒ equally spaced (or distribution of Gaussian formulae) grid points.

• But we have seen y is a better variable and we will rather get equally spaced points in y.

• Now dx = \sqrt{x} dy. A fixed spacing dy ⇒ spacing in x space ∝ \sqrt{x}.

• That is, x points are concentrated near the origin.

• This is intuitively reasonable as the integrand \frac{1}{\sqrt{x}} exp (−x^{3}) is rapidly varying near the origin.

How Do You Choose Correct Strategy?

• d) General Advice

• The user must “stare” at his or her function and decide on “correct” choice of strategy.

• There will be many correct choices.

• One can mix strategies and also may need different strategies in different regions.

• Consider

I = ∫_{0}^{∞} \frac{1}{\sqrt{x}} [exp (−x^{3}) + exp(−x)] dx = I_{1} + I_{2}

• with I_{1} given by

I_{1} = ∫_{0}^{x∼1} \frac{1}{\sqrt{x}} [exp (−x^{3}) + exp(−x)] dx.

• Use change of variables y = 2\sqrt{x} \underline{or} weight w(x) = 1/\sqrt{x} formula.

Hybrid Strategy Example (continued)

• The rest of the integrand

I_{2} = ∫_{1}^{∞} f(x) dx f(x) = \frac{1}{\sqrt{x}} [exp (−x^{3}) + exp(−x)] again.

• Here, the problem is “infinite range” of integral.

• So we let either

y = exp(−x/2) x = (−2 log y)
dy = −\frac{1}{2} exp(−x/2) dx
I_{2} = ∫_{0}^{1/\sqrt{e}} \frac{2}{\sqrt{(−2 log y)}} (y + \frac{1}{y} exp ((2 log y)^{3})) dy

• or

y = \frac{1}{x}
I_{2} = ∫_{0}^{1} (\frac{1}{y})^{3/2} [exp (−\frac{1}{y^{3}}) + exp (−\frac{1}{y})] dy

Note these integrals are well behaved.

CPS615NI G.C.Fox, Oct 95  19  gcf@npac.syr.edu

NI.2: Romberg or

Adaptive Integration

• A very convenient way of evaluating integrals is to repeatedly apply Simpson’s rule—halving interval each time.

a  b
x  x  x Step 1
•  x  •  x  •  2
• x • x • x • x •  3
→ x

• The point is that, as diagram indicates, one need only evaluate half as many (x in diagram) functions as you expect because you use old values (• in figure).

• As the error is proportional to nh^{5} and n doubles and h halves each go, error should go down by ≈ 1/16 each iteration and convergence is fast. So, in practice, one halves h until change in integral between steps is less than some assigned value.

NI.3: Gaussian Integration

• i) We are evaluating the same formula

I = ∫_{a}^{b} g(x) f(x) dx  (NI.1)

in terms of the same ansatz

I = n ; ∑ ; r=0 w_{r} f(x_{r})  (NI.2)

• Now we choose \underline{both} w_{r} \underline{and} x_{r} (a grand total of 2n + 2 parameters) to make (NI.2) exact for polynomials of degree ≤ 2n + 1. (2n + 2 equations analogous to (NI.3)).

• With Newton-Cotes, one only chooses n + 1 parameters w_{r} viewing x_{r} as fixed.

• Thus, Newton-Cotes formulae are exact for polynomials of degree ≤ n.

Some Mathematics Underlying

Gaussian Integration

• This can be ignored!

• Let P_{n}(x) be the polynomials orthogonal with respect to weight function g(x)

∫_{a}^{b} g(x) P_{n}(x) P_{m}(x) dx = h_{n} δ_{nm}  (NI.10)

h_{n} is a known irrelevant constant.

• Let x_{i}^{(n)} be the i th zero of P_{n}(x).

• Define K_{n}(x, y) by the formula

K_{n}(x, y) = n ; ∑ ; k=0 P_{k}(x) \frac{P_{k}(y)}{h_{k}}  (NI.11)

• If Π_{n}(x) is any polynomial of degree n, then one can easily show:

∫_{a}^{b} g(x) K_{n}(x, y) Π_{n}(x) dx = Π_{n}(y)  (NI.12)

• It is not difficult to see that (NI.2) is satisfied by

x_{i} = x_{i}^{(n+1)} , w_{i} = \frac{1}{K_{n+1} (x_{i}^{(n+1)}; x_{i}^{(n+1)})} (NI.13)

CPS615NI G.C.Fox, Oct 95  22  gcf@npac.syr.edu

More Mathematics Deriving Gaussian

Integration Formulae

• To convince yourself, choose a particular complete set of polynomials of degree 2n + 1.

• a) Firstly, P_{n+1}(x) multiplied into 1, x, x^{2}, … x^{n}

• For these, I is zero because x_{i} = x_{i}^{(n+1)} are at zeros of P_{n+1}(x).

• b) Secondly,

K_{n+1} (x_{i}^{(n+1)}, x) (x − x_{0}^{(n+1)}) (x − x_{1}^{(n+1)})

… (x − x_{i−1}^{(n+1)}) (x − x_{i+1}^{(n+1)}) … (x − x_{n}^{(n+1)})

Note: Omit (x − x_{i}^{(n+1)}) factor

• By construction, the only term that contributes to I for this polynomial is that with r = i and this gives formula for w_{i}.

• This is quick derivation of Gaussian Integration. Key point is that x_{i} are zeroes of the orthogonal polynomial P_{n+1}(x).

Basic Gaussian Integration

• For g(x) = 1, we get normal Gauss formulae with:

x_{i} = i^{th} zero of P_{n+1}(x)
w_{i} = \frac{2}{[1 − x_{i}^{2}]} [P′_{n+1}(x_{i})]^{−2}  (NI.14)

P_{n+1}(x) is Legendre Polynomial of degree n + 1.

• Clearly, Gaussian integration is best when f(x) is very hard to evaluate, and one must minimize value of n.

• However, it is very hard to estimate error in calculation as the simple iterative scheme described in Section NI.2 for Simpson’s rule cannot work.

• Not only is hard to evaluate weights and x_{i}’s as n iterates → ∞, but you have to recalculate all n functions each time. Note: There is no overlap in the x_{i} values as we double grid points (n → 2n + 1).

When Should You Use Gaussian?

• (iv) Comparing Gauss with fixed grid point method, Gauss “saves” a factor of two by doubling the polynomial order for fixed number of function evaluations. However, this advantage is lost when compared to Romberg as latter also gains a factor of two as reuses points. Romberg is further much easier to code.

• (v) So use Romberg if you have a few integrals which are needed to prescribed accuracy. Use Gaussian Integration when you have a lot of integrals with similar f(x) and the right value of n can be determined before hand.

Three Choices of x_{i} in Numerical Integration

Express I = n ; ∑ ; r=0 w_{r}f (x_{r})

as ∫_{a}^{b} f(x)dx = _{i}^{∑} w_{i}f (x_{i})

NI.4: Parallel Computing

and Integration

• If have n + 1 point integration (n = 15 above) divide into N_{proc} (N_{proc} = 4 above) equal parts.

• M = (n + 1)/N_{proc} (= 4 above) integration points in each processor. Each processor does an internal sum.

I_{l} = F_{l} ; ∑ ; i=B_{l} w_{i} f(x_{i})  l = 0 … N_{proc−1}
index B_{l} = Ml
F_{l} = B_{l} + M − 1

• Then total integral is I = I_{0} + I_{1} … + I_{N_{proc}−1}

• There is ONLY communication in one place—a global COMBINE (reduction).

Issues for Parallel Computing

• i) If N_{proc} does not exactly divide n + 1, then one (or more) processors have less work than the others—minor load imbalance and slightly tedious programming.

• ii) Adaptive algorithms (e.g., Romberg) could lead to load imbalance if one introduces new integration points selectively (so that one gets more integration points in one processor than another).

• Parallel random numbers were discussed in statistics module.

Two-Dimensional Integrals

• iii) For most one-dimensional integrals would not need parallel computing as rarely is this a big calculation.

• However, in multidimensional integrals, one can need more computing power.

∫_{a}^{b} dx ∫_{c}^{d} dy f(x, y)

• In first term, divide the x points between processors, as illustrated, for one-dimensional case. In second term, do internally to each processor.

• iv) We will return to multidimensional case when we discuss Monte Carlo methods, which are a more general and useful approach.

CPS615NI G.C.Fox, Oct 95  29  gcf@npac.syr.edu

NI.5: Monte Carlo Integration

• (i) “All” integration methods can be written as

I = ∫_{a}^{b} g(x) f(x) dx = n ; ∑ ; i=0 w_{i}f (x_{i})

• where the n + 1 points x_{i} and n + 1 weights w_{i} are determined in sundry ways to optimize various properties, as described earlier.

• If we increase n, the dependence of error is

Error ∝  \frac{1}{n^{2}}  iterated Trapezoidal
\frac{1}{n^{4}}  iterated Simpson
\frac{1}{n^{2m}}  iterated m-order Gauss

• with iterated groups of k points.

Basic Formulation of Integral

as a Mean (Expectation Value)

• Taking g(x) = 1 for convenience, we write our basic integral as

I = (b − a) ∫_{−∞}^{+∞} f(x) p(x) dx
p(x) = \frac{1}{(b − a)}  a ≤ x ≤ b
= 0 outside [a, b]

• where p(x) has carefully been constructed to be a properly normalized probability distribution for a random variable x which is uniformly distributed in the range a ≤ x ≤ b.

Basic Monte Carlo Approach

to Integration

• Now consider the Monte Carlo integration method, which writes I = ⟨f̃⟩.

⟨f̃⟩ = ∫_{a}^{b} \frac{f̃(x) dx}{b − a}  f̃ = (b − a) f  (NI.15)

I is mean of f̃(x) where x is uniformly distributed random variable between a and b.

• Defining

σ^{2} = ∫_{a}^{b} \frac{[f̃(x) − I ]^{2} dx}{b − a}  (NI.16)

• We know from central limit theorem that

\frac{1}{n} n ; ∑ ; i=1 f̃(x_{i}) → I  (NI.17)

with error \frac{σ}{\sqrt{n}}.

• (NI.15) is of the same form as our previous formulae: all the w_{i} are equal, while the positions x_{i} are random. However, it seems a bad idea as error → 0 much slower than even iterated trapezoidal.

Why Monte Carlo Methods Are Best

in Multidimensional Integrals

• But wait! Consider an integral in k dimensions – for other methods, we must make a grid in each direction with n^{1/k} points for each coordinate.

• Thus, error for Trapezoidal is ∝ 1/n^{2/k} which has a slower n dependence as k increases.

• However, Monte Carlo method still keeps its 1/\sqrt{n} behavior and for k ≥ 5, converges faster than Trapezoidal method.

• The other methods appear better than Monte Carlo until larger k.

• However, in practice, this is not so because there the implicit condition n^{1/k} >> 1 before the error formulae even apply.

• For instance, n^{1/k} = 10, k = 5 gives large (10^{5}) values of n while n^{1/k} = 1000 gives ridiculous values for n (10^{15}) with k = 5.

• This example shows that the problem of multidimensional integrals is not a trivial generalization of that in one dimension.

Best Multidimensional

Integration Formulae

• You can pose general problem of best method in k dimensions: the answer is not known in general and is probably not either Monte Carlo \underline{or} Gaussian for each coordinate.

• Once k is bigger than around 3, Monte Carlo is usually best pragmatic choice. The distribution of points in multidimensional integrals is contrasted in figures MC.1 and MC.2.

Distribution of Points in Two-dimensional Integral

Done by Newton-Cotes Style Formulae

Figure MC.1:

Two-dimensional integral done by Newton-Cotes style formulae

Distribution of Points in Two-dimensional Integral

Done by Monte Carlo

Figure MC.2:

Two-dimensional integral done by Monte Carlo

Use of Bounding Boxes

to Calculate π — I

• (ii) Complicated Domains

• One advantage of the Monte Carlo method is that it can very easily take account of complicated integration domains.

• We illustrate this first with a simple example—calculation of π.

• Suppose we wished to calculate the integral over the inside of the circle

I = ∫_{x^{2}=y^{2}≤1} g(x, y) dx dy.

• One way of doing this is to put

\begin{aligned}g′(x, y) & = g(x, y) & x^{2} + y^{2} ≤ 1 \\ & = 0 & x^{2} + y^{2} > 1 \end{aligned}

to Calculate π — II

• And calculate

I = ∫_{−1}^{1} dx ∫_{−1}^{1} dy g′(x, y)

which is an elementary Monte Carlo integral with x and y generated independently between −1 and 1.

• Actually, the brute force method (generate x uniformly between −1 and 1, y uniformly between ±\sqrt{1 − x^{2}}) is more efficient—in the language of (NI.16) this second method has a smaller variance.

• However, use of bounding boxes is an important technique for it is often impossible to find ranges analytically while it is easy to decide if a given (x, y) lies within range of integration.

Use of Bounding Boxes

in Complicated Geometries — I

• Example: Integrate g(x, y) over region h(x, y) ≥ 0. h(x, y) = x^{3} + 3xy^{2} + 23y^{3}

or something complicated like this.

• It is indeed \underline{easy} to decide if h(x, y) ≥ 0 or < 0 for a given value (x, y).

• However, it is hard to solve equation x^{3} + 3xy^{2} + 23y^{3} = 0 (the boundary of domain).

• Note you cannot use Simpson/Gauss/Trapezoidal to integrate

∫_{−1}^{+1} ∫_{−1}^{+1} dx dy F (x, y)

\begin{aligned}where F (x, y) & = f(x, y) & x^{2} + y^{2} ≤ 1 \\ & = 0 & x^{2} + y^{2} > 1 \end{aligned}

in Complicated Geometries — II

• F (x, y) is discontinuous at x = ±\sqrt{1 − y^{2}}. This violates assumptions in analytic methods which require F (x, y) be approximated by a linear (trapezoidal), quadratic (Simpson) or higher order polynomial.

• This assumption is violated by a function F (x, y) which is discontinuous. Even if f(x, y) happens to be zero at x^{2} + y^{2} = 1, it would rarely have required zero derivatives (of higher order for Gauss methods).

• On the other hand, Monte Carlo method only requires finite standard deviation. If has NO requirement on continuity of values; derivatives need not even exist.

• Thus, only Monte Carlo methods can use this “bounding box” approach to complicated integration domains. In practice, the latter are very common.

IMPORTANCE Sampling

Basic Theory

• (iii) Of course, it is obvious from (NI.16) that Monte Carlo integration gets better as standard deviation σ gets smaller.

• Using definition of σ^{2} (or remembering its name “variance”), it follows that method is best when f(x) is slowly varying over domain.

• For instance, if f(x) ≡ constant, σ = 0 and method is exact.

• Importance sampling (“sampling where integrand important”) is a technique for reducing variation of f(x).

• Thus if

I = ∫ f(x) dx
change variables  y = q(x) , dy = q′(x) dx
put \frac{f(x)}{q′(x)} = F (y)
then  I = ∫ F (y) dy.

• We now apply Monte Carlo method generating y uniformly.

Choice of Importance Sampling

Weight Function — I

• We try to choose q′(x) so f(x)/q′(x) has a smaller standard deviation than f(x).

• The best choice is

q′(x) = f(x)

• which implies

• 1. F (y) is constant with zero standard deviation • 2. Probability of getting a particular x value is

∝ q′(x) = f(x). That is, one “chooses” x values to congregate where function f is large.

• This is, of course, not realistic as to find q′(x) = f(x) requires solving

q(x) = ∫_{a}^{x} f(x) dx

(i.e., solving indefinite integral and that’s harder than definite integral which is q(b)).

Choice of Importance Sampling

Weight Function — II

• In general, one chooses q′(x) so that • a) q′(x) is “approximately” f(x)

• b) ∫ q′(x) can be performed easily.

• Example: f(x) = e^{−x} (⊗ a mess) where mess slowly varying.

• Choose

q′(x) = e^{−x}
q(x) = 1 − e^{−x}

• Then F (y) = “a mess”.

Monte Carlo Approach

to Discrete Integrals

• (iv) Discrete Integrals

• One can apply Monte Carlo methods to sums as well as continuous integrals. For example, suppose we have to calculate

∫ dx_{a} dx_{b} dx_{c} … 2 ; ∑ ; l=1 f_{l} (x_{a}, x_{b}, x_{c}…)

• Then generate using normal Monte Carlo approach with random numbers to generate x_{a}, x_{b}, x_{c}.

• Now, generate another random number r with 0 ≤ r ≤ 1 uniformly and then code

if (r < 0.5) l = 1
else l = 2

• and obvious modifications for other discrete summations.

Why Use Monte Carlo

for Summations?

• This can have important consequences. Suppose our integral involves

S = 1000 ; ∑ ; l=1 f_{l} (x_{a}, x_{b}, x_{c})

• Then if we direct calculate the sum, it has computational cost, which is 1000 times that of calculating f_{l}.

• Suppose we generate r and set l_{r} = 1 + INT(1000 × r) (INT(y) gives you integer part of floating point number y).

• And then let S = 1000 × f_{l_{r}} (x_{a}, x_{b}, x_{c}) where x_{a}, x_{b}, x_{c} chosen by continuous Monte Carlo procedure.

Example of using Monte Carlo

for Summations

• Then computational complexity of whole problem is reduced by a factor of 1000! The price, of course, is that error is increased by variance (over l) of f_{l}.

• This variance is added in to variance of anything else in integrand including variance of f_{l} over x_{a}, x_{b}, x_{c} which is already present.

• Thus, we can increase total number of points n (where x_{a}, x_{b}, x_{c}, l are chosen separately n times) by 1000 and reduce total error by \sqrt{1000} ∼ 32 if we generate l randomly rather than using full calculation of sum.

• This strategy is often a big win!

• See following discussion for a related comment.

The Wrong Way to Perform

Multiple Monte Carlo

• (v) When one does multiple integrals ∫ dx ∫ dy f(x, y) by Monte Carlo and chooses a set of n points \underline{x}_{i} = (x_{i}, y_{i}) where for \underline{each} i, one chooses separately x_{i}, y_{i} randomly. This gives \underline{x}_{i} distributed randomly over the plane.

• It would be \underline{very} wrong to view as in Newton Cotes as an integral over x of F (x) = ∫ dy f(x, y) where we think of doing x and y integrals by Monte Carlo.

• This incorrect approach gives a set of points x_{i}, y_{i} where one first chooses x_{1} randomly, then fixing this value of x_{1}, chooses n_{y} ∼ \sqrt{n} random values of y; then change to a new random x and so on.

• This scenario gives use to a set (x_{i}, y_{i}) where

x_{1} = x_{2} … = x_{n_{y}}
x_{n_{y}+1} = x_{2n_{y}} etc.

• In this case error is proportional to 1/\sqrt{n_{x}} ∼ 1/n^{1/4} and NOT 1/\sqrt{n} where all x_{i} and all y_{i} are chosen independently.

Stock Market Example of

Multiple Monte Carlos — I

• You may think of this as obvious, but let us study its consequences in a generalized example. We wish to make money on the stock market using a Monte Carlo model of their behavior on certain variables I = ∫ f (x_{1}, x_{2} …) dx_{1} . . ..

• We need an accuracy of say 1%. Then we need a number of events n_{f} where σ_{f}/n_{f}^{1/2} ∼ 0.1⟨f⟩ where desired integral is up to a constant ⟨f⟩ = μ_{f}. Then

n_{f} = 10000 (\frac{σ_{f}}{μ_{f}} )^{2}

• Now suppose your friend had another set of variables whose Monte Carlo J = ∫ g (y_{1}, y_{2},…) dy_{1}, dy_{2} required

n_{g} = 10000 (\frac{σ_{g}}{μ_{g}} )^{2}

Stock Market Example of

Multiple Monte Carlos — II

• Now we realize that we need to combine the two approaches and integrate over both x_{i} and y_{j}!

• Then naively this is hopeless as it requires

n = n_{f} · n_{g} = 10^{8} (\frac{σ_{f} σ_{g}}{μ_{f} μ_{g}} )^{2} points.

• But this is regarding joint Monte Carlo independently, but rather we should think about problem as

K = ∫ f(x_{1} …) g(y_{1} …) dx_{1} … dy_{1} …

where fresh x_{i}, y_{j} are chosen each point.

Number of Points Needed

in Joint Monte Carlo

• Then assuming special case written above where integrand has product form f g and f and g are independent, we need

n = n_{k} = 10^{4} (\frac{σ_{fg}}{μ_{f} μ_{g}} )^{2} points

σ_{fg}^{2} = σ_{f}^{2} σ_{g}^{2} + μ_{f}^{2} σ_{g}^{2} + μ_{g}^{2} σ_{f}^{2}

as can easily be proven.

\frac{σ_{fg}^{2}}{μ_{f}^{2} μ_{g}^{2}} ∼ \frac{σ_{f}^{2}}{μ_{f}^{2}} + \frac{σ_{g}^{2}}{μ_{g}^{2}}.

• We see that we need NOT n_{k} = n_{f} · n_{g} points, but rather n_{k} = n_{f} + n_{g} points!

• The exact calculation requires one to relax the unrealistic assumption that f and g are independent, but the point is general and the estimate of n_{k} is reasonably accurate!

Accept/Reject Method for Generating

General Probability Distributions

• (vi) Accept/Reject Method

• This is an important practical method for generating random numbers according to particular probability distributions, which cannot be integrated so that ∫_{x_{1}}^{x} p(x) can be used as described earlier.

• We wish to generate x according to a given probability distribution p(x) in range x_{1} ≤ x ≤ x_{2}.

• Put

p_{max} = max_{x} p(x)

Let w_{max} be (any number) ≥ p_{max}.

• Generate x uniformly in [x_{1}x_{2}]. (*) Generate r uniformly in [01].

If r > \frac{p(x)}{w_{max}} reject x and start again at (*)
r ≤ \frac{p(x)}{w_{max}} accept x

Estimate of Maximum in

Accept/Reject Method

• The resultant accepted x is indeed distributed according to p(x).

• Often, it is hard to find p_{max}—remember we might apply method to p(x_{1}, x_{2},… x_{r}) with a very complicated formula for p. One can estimate a w_{max} sometimes.

• The method works for any w_{max} ≥ p_{max} but if w_{max} is chosen as twice the real p_{max}, then one unnecessarily rejects half of the possible events.

• The method works, but is inefficient.

• Actually, one can use w_{max} < p_{max} if either

• weight events by p(x)/w_{max} if p(x) > w_{max}

• just smile and accept error if doesn’t happen often

Introduction to Metropolis Method

• (vii) Metropolis Method

• This very important method was introduced by N. Metropolis in Journal of Chemical Physics, 21, 1087 (1953).

• It is used in nearly all large scale statistical physics simulations.

• We are again doing “importance” sampling for ∫ p(\underline{x}) f(\underline{x}) d\underline{x} where we wish to choose \underline{x} according to p(\underline{x}).

• In the simplest physics applications p(x) = exp (−(Energy = E(\underline{x}))/kT) and is (due to exponential and large energy E(\underline{x})) typically tiny!

• Direct application of Monte Carlo without importance sampling or use of “accept/reject” method will “never” find interesting values of \underline{x} with p(\underline{x}) ∼ p_{max}.

The Metropolis Procedure

• The method starts with a point \underline{x}_{i} which can choose at random or with some prejudice. We now wish to generate the next point \underline{x}_{j}.

• Choose \underline{x}_{t} uniformly in a neighborhood of \underline{x}_{i}, i.e., −Δ\underline{x} + \underline{x}_{i} ≤ \underline{x}_{t} ≤ Δx + \underline{x}_{i}.

• Then examine p (\underline{x}_{i}) and p (\underline{x}_{t})

• if p (\underline{x}_{t}) ≥ p (\underline{x}_{i}), accept \underline{x}_{t} and put \underline{x}_{i+1} = \underline{x}_{t} and

• if p (\underline{x}_{t}) < p (\underline{x}_{i}), generate a new random number r uniformly in [0, 1].

• if r < p (\underline{x}_{t}))/p (x_{i}) < 1, set \underline{x}_{i+1} = \underline{x}_{t}

• if r ≥ p (\underline{x}_{t}))/p (x_{i}) < 1, set \underline{x}_{i+1} = \underline{x}_{i}

• \underline{x}_{i+1} is next random point (note it may be identical to \underline{x}_{i}). Eventually, this produces a set of points distributed “according to p(x)”.

Why Metropolis Method Works

• Reasons:

• 1. The neighborhoods Δ\underline{x} − \underline{x}_{i} ≤ \underline{x}_{t} ≤ \underline{x}_{i} + Δ\underline{x} eventually cover the whole space.

• 2. Take any two points \underline{x}_{i}, \underline{x}_{j}. Suppose p (\underline{x}_{i}) > p (\underline{x}_{j}), i, j label states or points in system.

• The Metropolis Method is such that probabilities of transition are P r(i → j), which is p (x_{j}) /p (x_{i}) and P r(j → i) is 1. That is,

\frac{P r(i → j)}{P r(j → i)} = \frac{p (x_{j})}{p (x_{i})}

• 3. In equilibrium,

N_{i} P r(i → j) = N_{j} P r(j → i).

• 4. Therefore,

\frac{N_{i}}{N_{j}} ∝ \frac{p (x_{i})}{p (x_{j})}

• 5. Here N_{i} and N_{j} are number (density) of points generated at \underline{x}_{i} and \underline{x}_{j}.

Monte Carlo Examples

Example 1: An Experimental

Physics Application

• The SSC near Dallas was to have been built at a cost of > ∼$10B. Our experiment has looked for a Higgs Particle and after a one-year run, it has seen no candidates?

• What does this mean for the underlying knowledge of physics?

• Theory:

• Lots of Particles (including Higgs).

A High Energy Experiment Scenario

• Each particle is seen in a huge detector with some probability.

• Some particles miss apparatus. Others hit chambers/bend into magnets.

• All “multiple scatter”.

• Many produce secondary “showers” (i.e., final detection is a “mess”).

• Further theory is • a) An approximation;

• b) Hard to calculate.

• However it is “just” an integral.

CPS615NI G.C.Fox, Oct 95  57  gcf@npac.syr.edu

An Experimental Physics

Monte Carlo

• Chance of seeing a Higgs

= _{∫∗} (Cross-Section(= Probability) to

Produce Higgs)

⊗ _{∫∗∗} (Particles Interacting with
Experimental Detector) (NI.18)

• * represents all relevant theoretical parameters • ** represents all possible things that happen to

particles as they traverse real world for ∼ 50 nanoseconds (i.e., for around 50 feet as speed of light is a nanosecond per foot!).

• This is a “Double” Monte Carlo—each Monte Carlo is a multiple integral.

• Note that “Second Monte Carlo” has lots of “non-differentiable functions” as particles are lost/detected with probability 0,1 respectively at sharp geometric boundaries.

Double Monte Carlo’s Again — I

• A digression on double integrals, which we already discussed in NI.4 in general terms

I = ∫ f(x) dx ∫ g(x, y) dy

• Each integral ∫ could be several dimensions, but write as one. If

∫ f(x) dx = I_{1} needs N_{1} points
∫ g(x, y) dy = I_{2} needs N_{2} points

• Then for analytic method (such as Romberg), the calculation of I would need (approximately) N_{1} · N_{2} points. Computer time is “multiplicative”.

• What about Monte Carlo! Fractional error

\frac{δI_{i}}{I_{i}} = (\frac{σ_{i}}{I_{i}}) /\sqrt{N_{i}}

Double Monte Carlo’s Again — II

• We can estimate this error assuming f and g are independent.

σ_{fg}^{2} = ∫ (f^{2}(x) g^{2}(y) − ⟨f⟩^{2} ⟨g⟩^{2}) dx dy
= ∫ f^{2}(x) dx ∫ g^{2}(y) dy − ⟨f⟩^{2} ⟨g⟩^{2}
= (σ_{f}^{2} + ⟨f⟩^{2}) (σ_{g}^{2} + ⟨g⟩^{2}) − ⟨f⟩^{2} ⟨g^{2}⟩
\frac{σ_{fg}^{2}}{⟨f⟩^{2} ⟨g⟩^{2}} = \frac{σ_{f}^{2}}{⟨f⟩^{2}} + \frac{σ_{g}^{2}}{⟨g⟩^{2}}
(\frac{δI}{I} )^{2} = (\frac{δI_{1}}{I_{1}})^{2} + (\frac{δI_{2}}{I_{2}})^{2}

• NOTE: We add, not multiply errors.

• Suppose we use N_{1} = N_{2} = N events.

• Therefore, use of 2N events gets the same fractional error in I as in I_{1}, I_{2}, if these each use N events and have σ_{f}/⟨f⟩ = σ_{g}/⟨g⟩.

• So do NOT think of (NI.18) as a Monte Carlo over theory and for each theory (model) a Monte Carlo over experiment.

• Rather, (NI.18) is a single Monte Carlo where “random x” plus “random y” form a single random point in multidimensional space.

CPS615NI G.C.Fox, Oct 95  60  gcf@npac.syr.edu

A Monte Carlo Event

• We call the resultant point “an event” in analogy to what happens in real world. Monte Carlo Integrals are “just” averages over events. Call

• “f” = Cross-section to produce Higgs

• “g” = Particles interacting with experimental detector.

• In the real world, each event has “unit weight” (God does perfect importance sampling).

• In the Monte Carlo,

• g is typically 0 or 1, as it is in real world,

• f is a widely varying function reflecting theoretical expectation (model).

• Calculating g is very time consuming as need to track many (100’s at SSC) particles through apparatus.

• Also, one typically has several different g’s to calculate for each f.

• These different g’s are the responses of different parts of apparatus.

Uniform Weight Events

• This seems a terrible “waste” if weight f(\underline{x}) (particular event) small compared to average, i.e.,

\frac{f (\underline{x})}{⟨f⟩} small.

• So we use “accept-reject” method, accepting event \underline{x} with probability

\frac{f (\underline{x})}{max\underline{x} f (\underline{x}) |_{guess/correct}}

• In practice, often hard to find and one guesses the maximum—often incorrectly (i.e., less than true maximum).

What Happens If You Miscalculate

Maximum Weight in

Accept/Reject Approach

• It is possible to cope with an f (\underline{x}) > guessed max\underline{x} f(x) in two ways.

• 1. Weight event with factor w (\underline{x}) in integrand

w (\underline{x}) = \frac{f (\underline{x})}{max\underline{x} f (\underline{x}) |_{guess}}

• 2. Accept

\frac{f (\underline{x})}{max\underline{x} f (\underline{x}) |_{guess}}

copies of \underline{x}.

• This approach (2.) is possible even if not an integer using probability to accept > 1 particles!

Example 2: Parallel Computing

for “Event” Monte Carlos

• The easiest way to use parallel computing for Monte Carlo is to generate separate events on separate processors.

• This just needs independent random number sequences on each node (we discussed in statistics module).

• These “independent” events are naturally implemented in HPF with an INDEPENDENT DO loop.

• This is often called the “farm” concept.

• A key hardware feature is that little communication is needed between processors.

• One needs to accumulate the central limit theorem sums such as ∑_{events x} “f(x) g(y)”.

• This “reduction” included in both MPI (message passing) and HPF.

• Note “loosely-coupled” systems (networks of workstations) are suitable as little communication.

Example 3: Lattice Monte Carlo

Theoretical Physics

• A typical “theoretical” calculation in condensed matter or high energy physics involves integrals of form

I = ∫ exp [\frac{−H}{kT} ] Θ ds_{1},… ds_{n}

where s_{i} is spin (gluon field) at site i.

• Typically, the sites are in a regular lattice in two, three, or four dimensions.

• This integral would involve solve 10^{4} → 10^{6} dimensions and so must be evaluated by Monte Carlo.

• Θ is not a very rapidly varying function (typically a polynomial in the s_{i}).

• Also, we want to do integral for several different choices for Θ.

• The Hamiltonian H is proportional to size of system and so the weighting function exp[−H/kT] is very small and one must be very careful to choose one’s integration points appropriately.

Choice of Points in Lattice Monte Carlo

• Assume you have a particular choice s_{1}^{(0)} … s_{n}^{(0)} with H^{(0)} = H (s_{1}^{(0)} … s_{n}^{(0)}).

• Remember exp[−H/kT] is a tiny number (as H is energy of full system H/kT is approximately the number of objects in system) and more importantly, varies rapidly.

• Thus, we cannot choose points “uniformly”.

• Rather, we assume point \underline{s}^{(1)} “near” \underline{s}^{(0)}. So that exp[−H/kT] only changes a little. This can be achieved by only changing one of the spins at a time, i.e., we move s_{k}^{(0)} → s_{k}^{(1)} and s_{i}^{(0)}i ≠ k unchanged.

• Typically, this produces a change ΔH ∼ kT as kT is energy per site (up to numerical constants).

• Remember classical result of energy of 1/2kT per degree of freedom.

• So change in exp[−H/kT] is roughly exp[I kT/kT] (i.e., modest).

• Thus, we will not move off “peak” of integrand.

Pictorial View of

Lattice Monte Carlo Integrands

• Essential point is that integral looks like multidimensional version of figure below, with a single peak becoming a ridge in many dimensions.

• Our strategy ensures that we move along ridge keeping integrand large.

• The naive strategy of uniform s_{i} will always choose s_{i} in the valley with “zero” integrand.

Metropolis and Heat Bath Methods

• One usually considers one of two methods to change s_{k}.

• (a) Metropolis — generate s_{k}^{(1)} uniformly about s_{k}^{(0)} and accept or reject change based on

r = exp [− [\frac{H (… s_{k}^{(1)} …)}{kT} − \frac{H (… s_{k}^{(0)} …)}{kT}]]

• if r > 1 (energy decreased), accept s_{k}^{(1)}

• if r < 1 (energy increased), accept s_{k}^{(1)} with probability r.

• (b) Heat Bath — sometimes one can use “exact” method to generate a variable according to a given probability distribution.

• That is, if distribution function

F_{k}(x) = ∫ ds_{k} exp (−\frac{H (… s_{k}…)}{kT})

can be calculated and inverted.

• Here, we view H as a function of s_{k} with other spins fixed.

Calculation of Observables

• In each case, (a) or (b), spins are generated according to distribution exp[−H/kT] and so values of observables are calculated as

⟨Θ⟩ = \frac{1}{#spin_{{configurations \atop =N}}} N ; ∑ ; l=1 Θ (s_{1}^{(l)} … s_{m}^{(l)})

• If we had not been careful about generating the spin values, but say chosen them uniformly, then we must use

⟨Θ⟩ = \frac{N ; ∑ ; l=1 exp [−H[s^{(l)}]/kT ] Θ (s^{(l)} …)}{N ; ∑ ; l=1 exp [−H[s^{(l)}]/kT ]}

• This is correct, but numerically hopeless. The wide variation in exponential (remember H proportional to size) would give a huge standard deviation in ⟨Θ⟩.

Example 4: Parallel Computing

for Lattice Theory

• One can in some cases do Lattice Theory as in “event” case with separate lattice’s on each node of farm.

• Parallelism comes from the different events (here states of system) gotten with different random numbers on each node.

• However, this is not usually possible as system is too big to fit on a node.

• Rather, we take (two, three or four dimensional) lattice and domain decompose (i.e., get data parallelism from space on which “spins” defined).

• This leads to parallel issues similar to partial differential equations.

• The decomposition is illustrated in the following two figures.

A 12^{3} Problem Lattice

Decomposed Onto a 64-node Machine

Arranged as a 4^{3} Machine Lattice

A 12^{3} Problem Lattice Decomposed onto a 64-node Machine Arranged as a 4^{3} Machine Lattice. Points

labeled X (“red”) or • (“black”) can be updated at the same time.

The 16 Time and Eight Internal Gluon Degrees of Freedom

The 16 time and eight internal gluon degrees of freedom stored at each point shown in previous figure.