A Robust and Accurate Quasi-Monte Carlo Algorithm for Estimating Eigenvalue of Homogeneous Integral Equations

Abstract

We present an efficient numerical algorithm for computing the eigenvalue of the linear homogeneous integral equations. The proposed algorithm is based on antithetic Monte Carlo algorithm and a low-discrepancy sequence, namely, Faure sequence. To reduce the computational time we reduce the variance by using the antithetic variance reduction procedure. Numerical results show that our scheme is robust and accurate.

1. Introduction

In Monte Carlo (MC) methods the random variables are simulated by computer generated pseudorandom sequences, and the numerical solution is performed by averaging over a large number of simulations. Pseudorandom numbers mimic the realizations of independent identically distributed random variables. Due to a fundamental result in probability theory, Monte Carlo methods based on such pseudorandom numbers can converge only proportionally to the square root of the number of simulations. It was noticed that the randomness of the simulating sequence is not essential for numerical integration. Deterministic sequences which fill the space uniformly can also be used. One of the most powerful methods for improving Monte Carlo simulation is to use low-discrepancy numbers [1]. Therefore, instead of using random numbers, we can employ a deterministic sequence of numbers. Such a sequence is called a low-discrepancy sequence. The great advantage of low-discrepancy sequences is that their rate of convergence is 𝑂(𝑁⁻¹) rather than 𝑂(𝑁^−1/2), which greatly increases the competitiveness of Monte Carlo simulations. Using low-discrepancy sequences to carry out Monte Carlo algorithm is sometimes called quasi-Monte Carlo (QMC) algorithm [1]. The Faure sequence is one of the well-known quasirandom sequences used in quasi-Monte Carlo applications. In [2, 3] details of Faure sequence and its various scrambling methods are described.

The problem of using Monte Carlo and quasi-Monte Carlo methods for finding an eigenpair of matrices has been extensively studied [4]. In [5] a numerical method for finding the eigenvalue of the linear homogenous integral equation is proposed by applying antithetic Monte Carlo method to obtain a reliable scheme.

In this paper we present a new numerical method based on antithetic variance reduction method and quasi-Monte Carlo algorithm on scrambled Faure sequence for estimating the eigenvalue of the following linear homogenous integral equation: 𝑢(𝑥)=𝜆_𝐷𝑘𝑥,𝑦𝑢𝑦𝑑𝑦,𝜆≠0,𝑥∈𝐷. (1) It can be written in the form of 𝑢(𝑥)=𝜆𝐾𝑢(𝑥),𝜆≠0. (2) If 𝑢≠0, then 𝜆 and 𝑢 are called the eigenvalue and eigenfunction of the above equations, respectively. Also, 𝐾 is called the kernel of this integral equation. Throughout this paper, we suppose that the kernel 𝐾 is symmetric and positive definite; that is, 𝐾𝑥,𝑦=𝐾𝑦,𝑥,∀𝑥,𝑦∈𝐷. (3) And ⟨𝐾𝜓,𝜓⟩>0, for all 𝜓≠0.

Now, we consider the following continuous Markov Chain 𝑇_𝑘 with length, 𝑘: 𝑥₀⟶𝑥₁⟶⋯⟶𝑥_𝑘, (4) where 𝑥_𝑖∈𝐷⊆ℝ, 𝑖=1,2,…,𝑘. The statistical nature of constructing the chain (4) is 𝑝(𝑥₀) and 𝑝(𝑥,𝑦), which are the probability of starting chain from 𝑥₀∈𝐷 and transition probability from state 𝑥∈𝐷 to 𝑦∈𝐷, respectively [1].

Theorem 1. Under the above assumptions, the random variable Θ_𝑛(ℎ)=ℎ𝑥₀𝑝𝑥₀𝑊_𝑛𝜓𝑥_𝑛, (5) where 𝑊_𝑛=𝑊_𝑛−1𝐾𝑥_𝑛−1,𝑥_𝑛𝑝𝑥_𝑛−1,𝑥_𝑛,𝑊₀=1, (6) is unbiased estimator for ⟨𝐾^𝑛𝜓,ℎ⟩; that is, 𝐸Θ_𝑛(ℎ)=𝐾^𝑛𝜓,ℎ. (7)

Proof. See [5].

Theorem 2. Under the above conditions, one has 𝜆=lim_𝑛→∞𝐸Θ_𝑛(ℎ)𝐸Θ_𝑛+1(ℎ)=lim_𝑛→∞𝐾^𝑛𝜓,ℎ𝐾^𝑛+1.𝜓,ℎ (8)

Proof. See [5].

2. Monte Carlo and Quasi-Monte Carlo Methods

The Monte Carlo method was invented by Stanislaw Ulam, a Polish born mathematician who worked for John von Neumann on the United States’ Manhattan Project during World War II, in 1946 [6]. Monte Carlo methods are based on the simulation of stochastic processes whose expected values are equal to computationally interesting quantities. MC methods offer simplicity of construction and are often designed to mirror some process whose behavior is only understood in a statistical sense. However, there are wide classes of problems where MC methods are the only known computational method of solution. Despite the universality of MC methods, a serious drawback is their slow convergence, which is based on the 𝑂(𝑁^−1/2) behavior of the size of statistical sampling errors. One generic approach to improve the convergence of MC methods has been the use of highly uniform random numbers in place of the usual pseudorandom numbers. While pseudorandom numbers are constructed to mimic the behavior of truly random numbers, these highly uniform numbers, called quasirandom numbers (or low-discrepancy sequences), are constructed to be as evenly distributed as mathematically possible. The use of quasirandom numbers in MC leads to quasi-Monte Carlo methods [7, 8]. Quasi-Monte Carlo methods are a variant of ordinary Monte Carlo methods that employ highly uniform quasirandom numbers in place of Monte Carlo’s pseudorandom numbers. QMC methods can be viewed as deterministic versions of Monte Carlo methods [7]. Determinism enters in two ways: by working with deterministic points rather than random samples and by the availability of deterministic error bounds instead of probabilistic MC error bounds. It could be argued that most practical implementations of MC methods are, in fact, quasi-Monte Carlo methods since the purportedly random samples that are used in Monte Carlo calculation are often generated in the computer by the deterministic algorithm. In QMC methods deterministic nodes are selected in such a way that the error bound is as small as possible. The very nature of the QMC methods with its completely deterministic procedures implies that we get deterministic and thus guaranteed error bounds [7]. In principle, it is therefore always possible to determine in advance an integration rule that yields a given accuracy. For example, for MC it is necessary to increase the number of simulations100 times to reduce the error by a factor of 10, whereas the QMC requires less (in general much less) than 100 times, and only 10 times in optimal cases. Monte Carlo methods offer statistical error estimates; however, while quasi-Monte Carlo has a faster convergence rate than normal Monte Carlo, one cannot obtain error estimates from quasi-Monte Carlo sample values by any practical way. A recently proposed method, called randomized quasi-Monte Carlo method, takes advantage of MC and QMC methods. Randomness can be brought to bear on quasirandom sequences through scrambling and other related randomization techniques in randomized QMC methods, which provide an elegant approach to obtain error estimates for quasi-Monte Carlo based on treating each scrambled sequence as a different and independent random sample. The core of randomized quasi-Monte Carlo method is to find an effective and fast algorithm to scramble (randomize) quasirandom sequences. Besides obtaining error estimates for quasi-Monte Carlo, scrambling techniques provide a natural way to parallelize quasirandom sequences. This scheme is especially suitable for distributed or grid computing. By scrambling a quasirandom sequence we can produce a family of related quasirandom sequences.

3. Quasirandom or Low-Discrepancy Sequences

Quasirandom numbers are constructed to minimize a measure of their deviation from uniformity called discrepancy. A quasirandom sequence is normally generated in the unit 𝑠-dimensional hypercube, 𝐼^𝑠=[0,1)^𝑠, and attempts to fill the hypercube as uniformly as possible. The original construction of quasirandom sequences was related to the van der Corput sequence, which is a one-dimension quasirandom sequence based on digital inversion. This digital inversion method is a central idea behind the construction of many current quasirandom sequences in arbitrary bases and dimensions. Following that, Halton [9] generalized the van der Corput sequence to 𝑠 dimensions, and Sobol [10] introduced the sequences that now bear his name. As significant generalization of these methods was proposed by Faure [6]. Niederreiter [7] generalized the existing construction of the Sobol and Faure sequences to arbitrary bases. These are now called Niederreiter sequences. Tezuka [11] further generalized Niederreiter sequences by using the polynomial arithmetic analogue of Halton sequences. The construction of quasirandom sequences is based on minimizing their discrepancy. Quasirandom sequences aim to have the fraction of their points within any subset 𝐽(𝜈)=[0,𝜈₁)×[0,𝜈₂)×⋯×[0,𝜈_𝑠) as close as possible to the subset’s volume fraction. Based on stardiscrepancy, we expose that the following definition of a low-discrepancy sequence in [0,1)^𝑠 is as follows.

Definition 3. For any 𝑁>1 and sequence {𝑥_𝑖}, let {𝑥_𝑖}_{1≤𝑖≤𝑁} denote the first 𝑁 points of the sequence {𝑥_𝑖}. If one has 𝐷^∗_𝑁≤𝐶_𝑠log𝑁^𝑠𝑁, (9) where the constant 𝑐_𝑠 depends only on the dimension 𝑠, then the sequence {𝑥_𝑖} is called a low-discrepancy sequence.

There are some variations for random and quasirandom sequences. Tezuka formulated the generalized Faure sequence, based on Halton sequence but using polynomials for reordering. According to Traub, generalized Faure method converges significantly faster than Monte Carlo and simple Faure method. In [2] we represent two methods that mix Faure sequence with the best scrambling Faure sequence and called them MFaure and M2Faure sequences. Figures 1, 2, and 3 represent Faure, MFaure, and M2Faure sequences for dimension 40. Figure 1 illustrates a disadvantage of the Faure sequence: the Faure construction leads to a sequence that has correlations between its individual coordinates. This leads, among others, to bad two-dimensional projections and also has its consequences when the sequence is used for numerical integration. We must stress that its performance in low dimensions is better. In Figures 2 and 3, we see the advantages of the MFaure and M2Faure sequences [3]. In this paper, we employ M2Faure sequence in the proposed algorithm.

Figure 1 

Projection of the first 1024 points from a 40-dimensional Faure sequence in base 41.

Figure 2 

Projection of the first 1024 points from a 40-dimensional MFaure sequence in base 41.

Figure 3 

Projection of the first 1024 points from a 40-dimensional M2Faure sequence in base 41.

4. Variance Reduction Techniques

As Monte Carlo simulation is slow to converge and there are various methods for increasing the speed of convergence, this essentially comes down to simulate the estimator in such a way that their variance is reduced. One simple method is antithetic sampling. With antithetic sampling one draws samples in pairs. If 𝑈 is uniformly distributed, then so is 1−𝑈. For every random, draw two paths that are therefore simulated: one with 𝑈 and one with 1−𝑈. This ensures that the mean of the drawn paths is zero, and the symmetry of the normal distribution is achieved by construction [12].

Assume that we want to compute 𝐸[𝑓(𝑋)] with 𝑋 a random variable uniformly distributed on [0,1]. The Naive Monte Carlo (NMC) is as follows: 𝑓(𝑋)=1𝑁^𝑁_𝑖=1𝑓𝑋_𝑖, (10) where 𝑋_𝑖, 𝑖=1,…,𝑁, are independent copies of 𝑋. In the antithetic variates (AV) method we use the numbers 1−𝑋₁,…,1−𝑋_𝑁 and introduce the antithetic Monte Carlo estimator 𝑓_AV(𝑋)=12𝑁^𝑁_𝑖=1𝑓𝑋_𝑖+𝑓1−𝑋_𝑖. (11) We note that the AV estimator is unbiased [13]. Suppose that 𝜎²=var(𝑓(𝑋)). Then the variance of AV estimator is given by 𝑓_AV(𝑋)=𝜎²+12𝑁𝑓(𝑋)().2𝑁cov,𝑓1−𝑋 (12) If 𝑓(𝑋) and 𝑓(1−𝑋) are negatively correlated, then we have a reduction of the variance compared to Naive Monte Carlo estimator based on 2𝑁 random numbers.

5. The Proposed Algorithm (AVQMC Algorithm)

(1)Choose a value of Markov chain length 𝑛 and a value of number of Markov chain 𝑁.(2)Create a continuous quasirandom walk (Markov chain) 𝑥₀→𝑥₁→⋯→𝑥_𝑛, where 𝑥_𝑖 for 𝑖=0,…,𝑛 are elements of M2Faure sequence.(3)For each quasirandom walk 𝑠=1,…,𝑁  simulate the following random variable: Θ_𝑛^(𝑠,1)(ℎ)=ℎ𝑥₀^(𝑠)𝑝𝑥₀^(𝑠)𝑊_𝑛^(𝑠)𝜓𝑥_𝑛^(𝑠). (13)(3.1)Create𝑦₀→𝑦₁→⋯→𝑦_𝑛 such that 𝑦_𝑖=1−𝑥_𝑖 for 𝑖=0,…,𝑛.(3.2)Set the following random variable: Θ_𝑛^(𝑠,2)(ℎ)=ℎ𝑦₀^(𝑠)𝑝𝑦₀^(𝑠)𝑊_𝑛^(𝑠)𝜓𝑦_𝑛^(𝑠). (14)(4)Obtain the quasi-Monte Carlo estimation Θ_𝑛=12𝑁^𝑁_𝑠=1Θ_𝑛^(𝑠,1)(ℎ)+^𝑁_𝑠=1Θ_𝑛^(𝑠,2)(ℎ). (15)(5)PrintΘ_𝑛≈𝜆.(6)End of algorithm.

6. Simulation Results

In this section we compare three approaches of Monte Carlo methods, namely, Naive Monte Carlo (NMC), antithetic variates Monte Carlo (AVMC) [13], and the proposed algorithm (AVQMC) methods. We run algorithms for evaluating the eigenvalue of the following examples. We compare the relative error of AVQMC and AVMC algorithms and the Naive Monte Carlo (NMC) algorithm for each example. Figures 4 and 5 illustrate the relative error against number of sample sizes with 𝑁=1000 for Examples 4 and 5, respectively. All figures show that, when the sample size is small, AVQMC and AVMC algorithms are better than NMC algorithm. Also, when the sample size is large, AVQMC is better than all, and another advantage is that relative errors of method with Faure sequence are deterministic. Also, we obtain the estimation of the eigenvalue, with standard deviation (STD) estimation for different values of the number of Markov chains (𝑁). For each example the Monte Carlo simulation result is summarized in a figure. The numerical results in Figures 6 and 7 confirm that the AVQMC algorithm is a reliable stochastic algorithm for computing the eigenvalue of the homogeneous linear integral equations. Also, we outline the standard deviation (STD) of the unbiased Monte Carlo estimator, the approximated eigenvalue (𝜆), and the relative errors of three methods for Examples 4 and 5 in Table 1.

Table 1 

Relatives errors, STD, and 𝜆.

Figure 4 

Relative error against the number of sample points with 𝑁=1000 for Example 4.

Figure 5 

Relative error against the number of sample points with 𝑁=1000 for Example 5.

Figure 6 

Comparison between NMC, AVMC, and AVQMC methods regarding Example 4.

Figure 7 

Comparison between NMC, AVMC, and AVQMC methods regarding Example 5.

Example 4. Consider equation (2) with kernel 𝐾(𝑥,𝑦)=𝑒^𝑥−𝑦, 𝐷=[0,1]. We can easily show that the exact eigenvalue is 1.

Example 5. Consider equation (2) with kernel 𝐾(𝑥,𝑦)=sin(𝑥)sin(𝑦), 𝐷=[0,𝜋/2]. We can easily show that the exact eigenvalue is 4/𝜋=1.2732.

7. Conclusion

In this paper we have proposed an efficient algorithm and analyzed the performance of the algorithm for computing the eigenvalue of the linear homogenous integral equations. We studied the accuracy and robustness of the proposed algorithm. The numerical results in Table 1 show that AVQMC algorithm has better improvement, in both error and convergence, than AVMC and NMC algorithms.