Open Access
Add to BibSonomyAbstract
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.
©2012 Optical Society of America
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–3]. 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] and distributed rendering in computer graphics [5, 6]; 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] is used as a LDS. In addition to the Sobol’ sequence, we evaluate Richtmyer [8], Halton [9, 10], Faure [11] 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 xn = (xn(1), xn(2), …, xn(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 xn. When coordinates in the k-dimensional unit cube are t = (t1, t2, …, tk), the discrepancy DN(k) with the maximum norm value (L∞ norm) is defined as
Here, sup means least upper bound of set, λk is the Lebesgue measure of the kth 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 > 1Here, Ck 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] 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 = 219927 − 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:
Here, In is the radiant intensity in the direction vertical to the light source plane. In this case, the irradiance Eθ isHere, 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 θ isTherefore, the ray emission direction θ can be determined as follows, using an inverse transform method: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.
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 = 107 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
Here, E and indicate the irradiance distributions obtained by theory and simulation, respectively. The evaluation results are shown in Fig. 4 .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
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).
