Application and evaluation of quasi-Monte Carlo method in illumination optical systems

doi:10.1364/OE.20.009692

« Show journal navigation

Application and evaluation of quasi-Monte Carlo method in illumination optical systems

Shuhei Yoshida, Shuma Horiuchi, Zenta Ushiyama, and Manabu Yamamoto »View Author Affiliations

Optics Express, Vol. 20, Issue 9, pp. 9692-9697 (2012)
http://dx.doi.org/10.1364/OE.20.009692

View Full Text Article

Acrobat PDF (1915 KB)

Journal Search
Article Lookup

Article Tools

Citations

Alert me when this article is cited
Export Citation/Save

Article Navigation

▸ Article Outline
– Abstract
1. Introduction
2. QMC method and simulation model
3. Results
4. Conclusion
– References and links

Article Tools
Full Text
Article Info
References
Cited By
Figures
Metrics
Download PDF

Viewing Equations in the
HTML Article

Optimal Viewing Recommendations

Full Text
Article Info
References (12)
Cited By
Figures (4)
Metrics

Abstract

In this article, we evaluate a quasi-Monte Carlo (QMC) method with various low-discrepancy sequences (LDS) in illumination optical systems which are adopted in some commercial products, and clarify the method’s effectiveness quantitatively. We assumed the evaluated systems were an illumination optical system with a perfectly diffusing surface, and we compared them against the theoretical irradiance distribution. The evaluation results indicate that the QMC method delivers higher asymptotic convergence rate than the MC method does, and there is little difference between each LDS. In evaluation of simple optical systems that can be boiled down to low-dimensional numerical integration problems, the QMC method was found to be extremely effective.

1. Introduction

The purpose of this study was to evaluate the asymptotic convergence performance of the quasi-Monte Carlo (QMC) method with various low-discrepancy sequences (LDS) in illumination optical systems, and clarify the effect of the QMC method quantitatively. A low-discrepancy sequence (LDS) is a sequence with high uniformity, and has a low discrepancy. LDSs are designed to speed up the convergence of numerical integration [1

H. Niederreiter, Random Number Generation and Quasi-Monte Carlo Methods (Society for Industrial and Applied Mathematics, 1992).

–3

J. Dick and F. Pillichshammer, Digital Nets and Sequences (Cambridge University Press, 2010).

]. The Monte Carlo (MC) method that uses pseudorandom numbers has an error of the order of N^−1/2 for N samples. However, the QMC method that applies a LDS has an error of the order of (logN)^k/N in the case of k dimensional numerical integration. The QMC method that applies a LDS has been used in various fields, since its convergence property is better than that of the MC method. Major applications of the QMC method include high-order numerical integration in financial engineering [4

S. Ninomiya and S. Tezuka, “Toward real-time pricing of complex financial derivatives,” Appl. Math. Finance 3, 1–20 (1996).

] and distributed rendering in computer graphics [5

A. Keller, “Instant radiosity,” in Proceedings of the 24th annual conference on Computer graphics and interactive techniques, G. S. Owen, T. Whitted and B. Mones-Hattal ed. (ACM Press, 1997).

, 6

T. Kollig and A. Keller, “Efficient multidimensional sampling,” Comput. Graph. Forum 21, 557–563 (2002).

]; it has been proven that the convergence performance of the QMC method is better than that of the MC method in these fields. Considering these advantages of the QMC method, it can be expected that the QMC method that applies various types of LDSs is an effective method to improve the convergence performance in the evaluation of illumination optical systems. Here, we perform a detailed evaluation of the asymptotic convergence performance of the QMC method that applies four types of LDSs, and report on the results.

At present, the MC method is popularly used for the evaluation of illumination optical systems, and it is also incorporated in commercial software programs. The QMC method is also adopted in some commercial products such as ZEMAX® and LightTools®. In these products, the Sobol’ sequence [7

I. M. Sobol', “On the distribution of points in a cube and the approximate evaluation of integrals,” USSR Comput. Math. Math. Phys. 7, 86–112 (1967).

] is used as a LDS. In addition to the Sobol’ sequence, we evaluate Richtmyer [8

R. D. Richtmyer, “The evaluation of definite integrals, and a quasi-Monte-Carlo method based on the properties of algebraic numbers,” LA-1342, Los Alamos Scientific Laboratories, (1951).

], Halton [9

J. H. Halton, “On the efficiency of certain quasi-random sequences of points in evaluating multi-dimensional integrals,” Numer. Math. 2, 84–90 (1960).

, 10

J. M. Hammersley, “Monte Carlo methods for solving multivariable problems,” Ann. N. Y. Acad. Sci. 86, 844–874 (1960).

], Faure [11

H. Faure, “On the star-discrepancy of generalized Hammersley sequences in two dimensions,” Monatsh. Math. 101, 291–300 (1986).

] sequences. We believe the evaluation provides a better understanding of the effect of the QMC method in the field of illumination optics. Generally, to evaluate illumination optical systems on a computer, the radiant intensity distribution on the illuminated surface is calculated by focusing incoherent incident light onto the system from multiple light sources. The rays produced from the light sources are generated using the statistical sampling. In cases where the light sources have a certain area, the coordinates on the light source plane at which the rays are formed are determined using sampling points; simultaneously, the direction in which the rays advance is also determined. The distribution of sampling points is set such that it reflects the physical characteristics of the light sources, and the elements of the ray coordinates and direction are determined independently from one another. Here, we can expect the convergence to be improved by adopting LDSs as a sampling method. In this study, we assumed an illumination optical system with a perfectly diffusing surface as a simple optical system. We clarify the convergence performance of the QMC method is better than MC method, quantitatively.

2. QMC method and simulation model

This section describes the QMC method and the simulation model for illumination optical systems.

2.1 QMC method with LDS

First, the discrepancy and LDS are defined. Consider x_n = (x_n⁽¹⁾, x_n⁽²⁾, …, x_n⁽^k⁾), n = 0, 1, …, N − 1 as the set of points in the k-dimensional unit cube [0, 1]^k and A(E; N) as the number of points in an interval E out of x_n. When coordinates in the k-dimensional unit cube are t = (t₁, t₂, …, t_k), the discrepancy D_N⁽^k⁾ with the maximum norm value (L_∞ norm) is defined as

D_{N}^{(k)} = \sup_{t \in {[0, 1]}^{k}} | \frac{A ([0, t); N)}{N} - λ_{k} ([0, t)) | .

(1)

Here, sup means least upper bound of set, λ_k is the Lebesgue measure of the k^th dimension. The meaning of Eq. (1) is a discrepancy between the volume of the k-dimensional cuboid and its evaluation by numerical integration with the statistical sampling. The LDS is the infinite sequence satisfying the following formula for all values of N > 1

D_{N}^{(k)} \leq C_{k} {\frac{(\log N)}{N}}^{k} .

(2)

Here, C_k is a constant that depends only on dimension k. The types of LDSs adopted in this study are Richtmyer, Halton, Faure and Sobol’ sequence.

Figure 1 shows distribution of two dimensional LDSs with 1000 sampling points. Here, to compare the QMC sampling with the MC sampling, distribution of pseudorandom numbers using a mersenne twister (MT) [12

M. Matsumoto and T. Nishimura, “Mersenne twister: a 623-dimensionally equidistributed uniform pseudo-random number generator,” ACM Trans. Model. Comput. Simul. 8, 3–30 (1998).

] is shown in Fig. 1(a). The distribution represents the uniformity of LDSs is higher than that of pseudorandom numbers.

Fig. 1 Distribution of two dimensional LDSs and pseudorandom numbers. (a) MT. (b) Richtmyer. (c) Halton. (d) Faure. (e) Sobol’.

The period of MT we used is p = 2¹⁹⁹²⁷ − 1. By contrast, there is no period in the LDS. This is the feature of the LDS in comparison with pseudorandom numbers.

2.2 Simulation model

A perfectly diffusing surface is a plane with a stable radiance in all directions for viewing a light source. The radiant intensity in the direction that forms the angle θ with the normal of the perfectly diffusing surface is expressed as follows, according to Lambert’s cosine law:

I_{θ} = I_{n} \cos θ,

(3)

Here, I_n is the radiant intensity in the direction vertical to the light source plane. In this case, the irradiance E_θ is

E_{θ} = I_{θ} \frac{\cos θ}{d^{2}} = \frac{I_{n} \cos^{2} θ}{d^{2}},

(4)

Here, d is the distance between the light source and a position on the irradiance plane. Next, the method for reproducing a perfectly diffusing surface on a computer, using a uniformly distributed point sequence is described. We express the light emission angle as (ϕ, θ) in the spherical coordinate system. Here, the ranges of ϕ and θ are set up as 0 ≤ ϕ < 2π and 0 ≤ θ ≤ π/2, respectively. The cumulative distribution function G(θ) for θ is

G (θ) = \frac{1}{2} (1 - \cos 2 θ),

(5)

Therefore, the ray emission direction θ can be determined as follows, using an inverse transform method:

θ = G^{- 1} (x) = \frac{1}{2} arc \cos (1 - 2 x),

(6)

Here, x is a point sequence uniformly distributed in [0, 1]. Since the distribution of ϕ is uniform in [0, 2π], it is determined as ϕ = 2πx.

The optical system in a simulation is shown in Fig. 2 . In a simulation, we generated the rays based on the method described in the previous section from the coordinates (x, y, z) = (0, 0, 0). A position 1.0 mm away from the light source was set up as the irradiance plane, and precision was evaluated by comparing the irradiance distribution with the theoretical values Eq. (4). In this case, the irradiance distribution was determined depending on the number of rays in each area by dividing the 2.0 mm × 2.0 mm area centered on the irradiance plane into 512 × 512 portions.

Fig. 2 Optical system in the simulation.

3. Results

This section provides the results of a performance evaluation of the QMC method. We compared the irradiance distributions obtained using the MC and QMC methods as qualitative evaluation methods, and calculated the asymptotic convergence rate for each LDS as quantitative evaluation methods.

The irradiance distribution when the number of rays is N = 10⁷ is shown in Fig. 3 . It is evident from these results that the QMC method that applies a LDS shows smaller fluctuations in irradiance distribution than the MC method does.

Fig. 3 Radiant intensity distributions obtained through simulations. (a) MT. (b) Richtmyer. (c) Halton. (d) Faure. (e) Sobol’.

Next, we evaluated the precision of each LDS by comparing the irradiance distribution calculated by the QMC method with the theoretical values. The root-mean-square error (RMSE), defined by the following equation, was used to evaluate the precision

RMSE = \sqrt{\frac{{‖ E - \hat{E} ‖}^{2}}{N}} = \sqrt{\frac{\sum_{i} \sum_{j} {[E (i, j) - \hat{E} (i, j)]}^{2}}{N}} .

(7)

Here, E and

\hat{E}

indicate the irradiance distributions obtained by theory and simulation, respectively. The evaluation results are shown in Fig. 4 .

Fig. 4 Root-mean-square error evaluation results.

These results show that the QMC method delivers higher precision than the MC method does. Although slight vibrations are observed for the Richtmyer and Faure sequences, the precision is better than that of MT for all N values in all cases. Next, to evaluate the asymptotic convergence rate of the QMC method, we fitted the graph shown in the figure, which indicates the relationship between N and the RMSE, using the function cN⁻^a. The results are shown in Table 1 .

Table 1 Asymptotic Convergence Rate

Sequence	c	a
MT	298.3	0.500
Richtmyer	2017	0.670
Halton	2327	0.691
Faure	2189	0.680
Sobol’	2212	0.687

A larger a value indicates better asymptotic convergence. MT gives a = 0.5, which exactly matches the theoretical value. The asymptotic convergence of LDS shows an a value higher than that of MT by more than 0.17, indicating extremely high asymptotic convergence rate. Since the evaluation of the irradiance distribution is a low-dimensional numerical integration problem in simple illumination optical systems, as with perfectly diffusing surfaces, the QMC method delivers high precision, matching the theory. In addition, there is little difference between each LDS, but it might be said that the Halton and Sobol’ sequence are some better than other LDSs. Based on these results, the QMC method is expected to deliver good results in simple illumination optical systems, and it may safely be said that there is little reason to employ multiple LDS for simulation software.

4. Conclusion

We have evaluated the QMC method that applies various LDSs in the field of illumination optics, and clarified the effect of the method quantitatively. Evaluation results show the asymptotic convergence rate of the QMC method is extremely higher than the MC method in evaluation of simple optical systems that can be described by low-dimensional numerical integration problems, and there is little difference between each LDS.

We believe this study provides a better understanding of the effect of the QMC method in the field of illumination optics.

References and links

1.	H. Niederreiter, Random Number Generation and Quasi-Monte Carlo Methods (Society for Industrial and Applied Mathematics, 1992).
2.	M. Drmota and F. Robert, Tichy, Sequences, Discrepancies and Applications (Springer, 1997).
3.	J. Dick and F. Pillichshammer, Digital Nets and Sequences (Cambridge University Press, 2010).
4.	S. Ninomiya and S. Tezuka, “Toward real-time pricing of complex financial derivatives,” Appl. Math. Finance 3, 1–20 (1996).
5.	A. Keller, “Instant radiosity,” in Proceedings of the 24th annual conference on Computer graphics and interactive techniques, G. S. Owen, T. Whitted and B. Mones-Hattal ed. (ACM Press, 1997).
6.	T. Kollig and A. Keller, “Efficient multidimensional sampling,” Comput. Graph. Forum 21, 557–563 (2002).
7.	I. M. Sobol', “On the distribution of points in a cube and the approximate evaluation of integrals,” USSR Comput. Math. Math. Phys. 7, 86–112 (1967).
8.	R. D. Richtmyer, “The evaluation of definite integrals, and a quasi-Monte-Carlo method based on the properties of algebraic numbers,” LA-1342, Los Alamos Scientific Laboratories, (1951).
9.	J. H. Halton, “On the efficiency of certain quasi-random sequences of points in evaluating multi-dimensional integrals,” Numer. Math. 2, 84–90 (1960).
10.	J. M. Hammersley, “Monte Carlo methods for solving multivariable problems,” Ann. N. Y. Acad. Sci. 86, 844–874 (1960).
11.	H. Faure, “On the star-discrepancy of generalized Hammersley sequences in two dimensions,” Monatsh. Math. 101, 291–300 (1986).
12.	M. Matsumoto and T. Nishimura, “Mersenne twister: a 623-dimensionally equidistributed uniform pseudo-random number generator,” ACM Trans. Model. Comput. Simul. 8, 3–30 (1998).

OCIS Codes
(080.1753) Geometric optics : Computation methods
(220.2945) Optical design and fabrication : Illumination design

ToC Category:
Optical Design and Fabrication

History
Original Manuscript: March 1, 2012
Revised Manuscript: April 5, 2012
Manuscript Accepted: April 9, 2012
Published: April 12, 2012

Citation
Shuhei Yoshida, Shuma Horiuchi, Zenta Ushiyama, and Manabu Yamamoto, "Application and evaluation of quasi-Monte Carlo method in illumination optical systems," Opt. Express 20, 9692-9697 (2012)
http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-20-9-9692

Sort: Author | Year | Journal | Reset

References

H. Niederreiter, Random Number Generation and Quasi-Monte Carlo Methods (Society for Industrial and Applied Mathematics, 1992).
M. Drmota and F. Robert, Tichy, Sequences, Discrepancies and Applications (Springer, 1997).
J. Dick and F. Pillichshammer, Digital Nets and Sequences (Cambridge University Press, 2010).
S. Ninomiya and S. Tezuka, “Toward real-time pricing of complex financial derivatives,” Appl. Math. Finance3, 1–20 (1996).
A. Keller, “Instant radiosity,” in Proceedings of the 24th annual conference on Computer graphics and interactive techniques, G. S. Owen, T. Whitted and B. Mones-Hattal ed. (ACM Press, 1997).
T. Kollig and A. Keller, “Efficient multidimensional sampling,” Comput. Graph. Forum21, 557–563 (2002).
I. M. Sobol', “On the distribution of points in a cube and the approximate evaluation of integrals,” USSR Comput. Math. Math. Phys.7, 86–112 (1967).
R. D. Richtmyer, “The evaluation of definite integrals, and a quasi-Monte-Carlo method based on the properties of algebraic numbers,” LA-1342, Los Alamos Scientific Laboratories, (1951).
J. H. Halton, “On the efficiency of certain quasi-random sequences of points in evaluating multi-dimensional integrals,” Numer. Math.2, 84–90 (1960).
J. M. Hammersley, “Monte Carlo methods for solving multivariable problems,” Ann. N. Y. Acad. Sci.86, 844–874 (1960).
H. Faure, “On the star-discrepancy of generalized Hammersley sequences in two dimensions,” Monatsh. Math.101, 291–300 (1986).
M. Matsumoto and T. Nishimura, “Mersenne twister: a 623-dimensionally equidistributed uniform pseudo-random number generator,” ACM Trans. Model. Comput. Simul.8, 3–30 (1998).

Faure, H.
Halton, J. H.
Hammersley, J. M.
Keller, A.
Kollig, T.
Matsumoto, M.
Ninomiya, S.
Nishimura, T.
Sobol', I. M.
Tezuka, S.

ACM Trans. Model. Comput. Simul.

M. Matsumoto and T. Nishimura, “Mersenne twister: a 623-dimensionally equidistributed uniform pseudo-random number generator,” ACM Trans. Model. Comput. Simul.8, 3–30 (1998).

Ann. N. Y. Acad. Sci.

J. M. Hammersley, “Monte Carlo methods for solving multivariable problems,” Ann. N. Y. Acad. Sci.86, 844–874 (1960).

Appl. Math. Finance

S. Ninomiya and S. Tezuka, “Toward real-time pricing of complex financial derivatives,” Appl. Math. Finance3, 1–20 (1996).

Comput. Graph. Forum

T. Kollig and A. Keller, “Efficient multidimensional sampling,” Comput. Graph. Forum21, 557–563 (2002).

Monatsh. Math.

H. Faure, “On the star-discrepancy of generalized Hammersley sequences in two dimensions,” Monatsh. Math.101, 291–300 (1986).

Numer. Math.

J. H. Halton, “On the efficiency of certain quasi-random sequences of points in evaluating multi-dimensional integrals,” Numer. Math.2, 84–90 (1960).

USSR Comput. Math. Math. Phys.

I. M. Sobol', “On the distribution of points in a cube and the approximate evaluation of integrals,” USSR Comput. Math. Math. Phys.7, 86–112 (1967).

Other

R. D. Richtmyer, “The evaluation of definite integrals, and a quasi-Monte-Carlo method based on the properties of algebraic numbers,” LA-1342, Los Alamos Scientific Laboratories, (1951).
A. Keller, “Instant radiosity,” in Proceedings of the 24th annual conference on Computer graphics and interactive techniques, G. S. Owen, T. Whitted and B. Mones-Hattal ed. (ACM Press, 1997).
H. Niederreiter, Random Number Generation and Quasi-Monte Carlo Methods (Society for Industrial and Applied Mathematics, 1992).
M. Drmota and F. Robert, Tichy, Sequences, Discrepancies and Applications (Springer, 1997).
J. Dick and F. Pillichshammer, Digital Nets and Sequences (Cambridge University Press, 2010).

2002, Kollig, Comput. Graph. Forum
1998, Matsumoto, ACM Trans. Model. Comput. Simul.
1996, Ninomiya, Appl. Math. Finance
1986, Faure, Monatsh. Math.
1967, Sobol', USSR Comput. Math. Math. Phys.
1960, Halton, Numer. Math.
1960, Hammersley, Ann. N. Y. Acad. Sci.

Cited By

Alert me when this paper is cited

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.

Click here to see a list of articles that cite this paper