CS294-13: Special Topics Lecture #4Advanced Computer GraphicsUC Berkeley Monday, 1 4 Sep. 2009Monte Carlo IntegrationLecture #4: Monday, 1 4 Sep. 2009Lecturer: Ravi RamamoorthiScribe: Fu-Chung Huang1 Introduction and Quadrature MethodsIn rendering we have a problem to determine the intensity of a certain 3D point. This isdone by gathering lighting from all directions to that point, or mathematically speaking,integrating the incoming lights in the range of a unit hemisphere around that point,which could be hard. Here we introduce a very simple but effective method called MonteCarlo method, which utilizes the power of randomness to compute the expected valuefor the intensity o f that point.1.1 Integration in 1D with Quadrature MethodsBefore diving into complex problem domain, we first illustrate integratio n over certain1D function f(x). The problem is defined to find the integrated value I of function f(x)over some range x ∈ [a, b]:I =Zbaf(x)dxThe easiest way is to find its integral function F (x) =Rx−∞f(u)du, and plug in the rangeI = F (b) − F (a). However, finding an analytic form f or that integral function is thebiggest challenge, or even worse when the function is not continuous and the there is nosuch analytical form.One way to go around this is to equally divide the range [a, b] into many smallintervals, says n intervals with width h = (b − a)/n, which we call step size. For eachinterval, we evaluate the function value, and make the value represents that interval. Theother way to think about it is we are calculating the area covered by the function f(x).Now we can add up the contribution from each interval, and say the summation shouldapproximate the real value I. This is called the Quadrature Integration, and va riousrules specify how to evaluate the ar ea under f(x) to represent that interval.Rectangle rule uses one(or two) evaluated value(s) to calculate the area, Tra pezoidalrule uses two, a nd Simpson’s rule uses three, as illustrated in Fig. 1. Overall the methodhas a general form looks like:ˆI =nXi=1wif(xi) (1)2 CS294-13: Lecture #4Figure 1: Trapezoidal rule and Simpson’s rule for quadrature method The trape-zoidal method uses linear approximation, and Simpson’s method assumes a parabolic arc(polynomial of degree 2) to calculate the area.Since various methods use different way to approximate area, there are still uncovered/over-covered errors. In general the error decreases as we decrease the step size(by increasing n).Providing that f(x) has at least two continuous derivatives within [a, b] and |f′′(x)| ≤ M,for the rectangle rule the error is given by:ˆIR− I ≤(b − a)24h2Mand by substituting the step size(width) h with (b − a)/n, we have the error:ˆIR− I ≤(b − a)324n2M = O(n−2)Trapezoidal rule has similar error bound since it also uses linear function to approximatethe covered area, a nd its error is given byˆIT− I ≤(b − a)12h2M =(b − a)312n2M = O(n−2)Simpson’s rule uses a more sophisticated quadratic polynomial function(assuming f(x)has at least 4-th derivatives) to calculate the covered area, and thus has better errorbound:ˆIS− I ≤(b − a)180h4M =(b − a)5180n4M = O(n−4)With 1D function the above methods can converge rapidly(or error reduces rapidlyas n increases), but they do not have such a good behavior when it comes to higherdimension.CS294-13: Lecture #4 31.2 Integration in higher dimensionsAs in 1D the definition Eq. 1, we can extend the approximation using tensor productrule to s-dimensional function, which is redefined:ˆI =nXi1=1nXi2=1...nXis=1w1w2...wisf(x1, x2, ..., xis) (2)We again simply use quadrature rules to evaluate everything with n samples in eachdimension, and requiring nssamples in total. In the case of 1D function, the error b oundis O(n−r)(where r = 2 or 4). However with s-dimensional space using N = nssamples,the error bound is O(N−r) = O(n−r/s), which degrades rapidly because of the curse ofdimensionality. This means that even if we increase the number of sampled region byslicing the region more finely, the error doesn’t go away as fast as we increase samples.Furthermore, if the function f(x) presents discontinuities, the error bound is at bestO(n−1) for 1D case, then with higher dimensional space is at best O(n−1/s). There is animportant result which limits the convergence rate(as we have already shown) for anydeterministic quadrature rule, called Bakhalov’s theorem (which we are not going to talkabout it here).Since the error bound(or convergence rate) is not so good as quadrature method goingto higher dimension, we need something better, and it should also deals with disconti-nuities. Monte Carlo method has the properties that error bound is always O(n−1/2)regardless of dimensionality, and it is also immune t o discontinuities.2 Monte Carlo Met hodMonte Carlo method in plain words is simply randomness. Note there is another classof ra ndom algorithm called Las Vegas, which always leads to correct result, fo r examplequicksort picks random pivot. Monte Carlo method does not provide 100% correctness,but in general the expected results will be correct. Before talking how to use MonteCarlo method to integrate function, we first review some probability concepts that areuseful as building block.2.1 Probability Reviews2.1.1 CDF and PDFCumulative Distribution Function (CDF) P (x) for random variable X describes the prob-ability that r andom variable X is less than or equal to some value x. Let’s take the diceas an example with possible outcome X = {1, 2, 3, 4, 5, 6}. The CD F P(2) for rolling thedice is 1/3, since the possibility of the outcomes smaller or equal to 2 out of 6 possiblefaces are 1/3; similarly P (3) = 1/2 and P (4) = 2/ 3.4 CS294-13: Lecture #4Probability Density Function (PDF) p( x) is the po ssibility of each outcome for randomvariable X, defined as dP (x)/dx. In the dice example, the PDF p(x) for the occurrenceof each face is 1/6.The possibility of cumulated outcome over some range x ∈ [a, b] is defined asP (x ∈ [a, b]) =Zbap(x)dx = P (b) −P (a)2.1.2 Expected Value and VarianceGiven the definitions for CDF and PDF , we want to ask: what is the averaged outcomein general or in the long run? The expected value is to answer such question and definedby:Ep(f(x)) =ZΩf(x)p(x)dx (3)The expected value Ep() for some function f(x) is computed via drawing random samplex with some probability distribution p(x). For discrete cases, the expected value
View Full Document
Unlocking...