Partitioning Inverse Monte Carlo Iterative Algorithm for Finding the Three Smallest Eigenpairs of Generalized Eigenvalue Problem

Abstract

A new Monte Carlo approach for evaluating the generalized eigenpair of real symmetric matrices will be proposed. Algorithm for the three smallest eigenpairs based on the partitioning inverse Monte Carlo iterative (IMCI) method will be considered.

1. Introduction

It is well known that the problem of calculating the largest or smallest generalized eigenvalue problem is one of the most important problems in science and engineering [1, 2]. This problem arises naturally in many applications. Mathematically, it is a generalization of the symmetric eigenvalue problem, and it can be reduced to an equivalent symmetric eigenvalue problem. Let 𝐴,𝐵∈ℜ^𝑛×𝑛 be real symmetric matrices and the matrix 𝐵 a positive definite matrix. Consider the problem of evaluating the eigenvalues of the pencil (𝐴,𝐵), that is, the values for which𝐴𝑥=𝜆𝐵𝑥. (1.1) A generalized eigenvalue problem (1.1) is said to be symmetric positive definite (S/PD) if 𝐴 is symmetric and 𝐵 is positive definite.

2. Inverse Vector Iteration Method

Another procedure for eigenvalue prediction is to use the Rayleigh quotient given by [3](𝑥)=𝑥𝜆=𝜇^𝑇𝐴𝑥𝑥^𝑇.𝐵𝑥 (2.1) since 𝐵 is positive definite, then (2.1) is well defined.

Theorem 2.1. Suppose that 𝜆_min=𝜆₁<𝜆₂≤⋯≤𝜆_𝑛−1<𝜆_𝑛=𝜆_max are 𝑛 eigenvalues for pencil (𝐴,𝐵) and 𝑣₁,…,𝑣_𝑛, corresponding eigenvectors. Then for arbitrary initial vector 𝑥 one has [3] 𝜆₁<𝜇(𝑥)<𝜆_𝑛, (2.2) where 𝜇(𝑥) is as introduced in (2.1).

Theorem 2.2. The inverse vector iteration method for arbitrary choice vector 𝑥₀ is convergent to the smallest eigenvalue and corresponding eigenvector for pencil (𝐴,𝐵). Also, the rate of convergence depends on 𝒪(𝜆₂/𝜆₁)^𝑘, where 𝑘 is number of iterations [3].

Algorithm 1 evaluates the smallest eigenpair based on the inverse vector iteration [4].

Input:  Initial vector 𝑥0 begin  Set 𝜆1(1)=𝑥𝑇0𝐴𝑥0𝑥𝑇0𝐵𝑥0  For  𝑗=0,1,2,…  begin   Solve linear system 𝐴𝑧𝑗+1=𝐵𝑥𝑗 for 𝑧𝑗+1   Set 𝜆1(𝑗+1)=𝑧𝑇𝑗+1𝐴𝑧𝑗+1𝑧𝑇𝑗+1𝐵𝑧𝑗+1   Set 𝑥𝑗+1=𝑧𝑗+1‖𝑧𝑗+1‖   Output:  𝑥𝑗,𝜆1(𝑗)  end end

Algorithm 1 

3. Monte Carlo Method for Matrix Computations

Suppose that the matrix 𝐴∈ℜ^𝑛×𝑛 and two vectors 𝑓,ℎ∈ℜ^𝑛 are given. Consider the following Markov chain 𝑇_𝑖 with length 𝑖:𝑇_𝑖∶𝑘₀⟶𝑘₁⋯⟶𝑘_𝑖, (3.1) where for 𝑗=1,…,𝑖,𝑘_𝑗∈{1,2,…,𝑛}. The statistical nature of constructing the chain (3.1) follows as𝑝𝑘₀=𝛼=𝑝_𝛼,𝑝𝑘_𝑗=𝛽∣𝑘_𝑗−1=𝛼=𝑝_𝛼𝛽, (3.2) where 𝑝_𝛼 and 𝑝_𝛼𝛽 show the probability of starting chain at 𝛼 and transition probability from state 𝛼 to 𝛽, respectively.

In fact^𝑛_𝛼=1𝑝_𝛼=1,^𝑛_𝛽=1𝑝_𝛼𝛽=1,𝛼=1,2,…,𝑛. (3.3) Define the random variable 𝑊_𝑗 using the following recursion for𝑊₀=1,𝑊_𝑗=𝑊_𝑗−1𝑎_{𝑘𝑗−1𝑘𝑗}𝑝_{𝑘𝑗−1𝑘𝑗},𝑗=1,2,…,𝑖. (3.4) Now, define the following random variable:Θ[ℎ]=ℎ_𝑘0𝑝_𝑘0^∞_𝑗=0𝑊_𝑗𝑓_𝑘𝑗. (3.5)

Theorem 3.1. Consider the following system: 𝐴𝑥=𝑏. (3.6)Let the nonsingular matrix 𝑀∈ℜ^𝑛,such that 𝑀𝐴=𝐼−𝐿, then the system (3.6) can be presented in the following form: 𝑥=𝐿𝑥+𝑓, (3.7) where 𝑓=𝑀𝑏. Then under condition Max_𝑖∑^𝑛_𝑗=1|𝑙_𝑖𝑗|<1,one has [5] 𝐸{Θ[ℎ]}=⟨ℎ,𝑥⟩. (3.8)

Suppose that 𝑥^(𝑖) is the 𝑖th iterative solution of the following recursion relation with 𝑥⁽⁰⁾=𝑓. If we set the random variableΘ_𝑖[ℎ]=ℎ_𝑘0𝑝_𝑘0^𝑖_𝑗=0𝑊_𝑗𝑓_𝑘𝑗, (3.9) then𝐸Θ_𝑖[ℎ]=ℎ,𝑥^(𝑖+1). (3.10)

By simulating 𝑁 random paths with length 𝑖𝑇_𝑖^(𝑠)∶𝑘₀^(𝑠)⟶𝑘₁^(𝑠)⟶⋯⟶𝑘_𝑖^(𝑠),𝑠=1,2,…,𝑁, (3.11) we can findΘ_𝑖^(𝑠)(ℎ)=ℎ_𝑘0𝑝_𝑘^(𝑠)₀^𝑖_𝑗=0𝑊_𝑗^(𝑠)𝑓_𝑘𝑗,𝑠=1,…,𝑁. (3.12) The Monte Carlo estimation can be evaluated byΘ_𝑖=1𝑁^𝑁_𝑠=1Θ_𝑖^(𝑠)(ℎ) (3.13) which is an approximation of ⟨ℎ,𝑥^(𝑖+1)⟩.

From all possible permissible densities, we apply the following:𝑝_𝛼=ℎ_𝛼∑^𝑛_𝛼=1ℎ_𝛼,𝑝_𝛼𝛽=𝑎_𝛼𝛽∑^𝑛_𝛽=1𝑎_𝛼𝛽,𝛼=1,2,…,𝑛. (3.14) The choice of the initial density vector and the transition probability matrix leads to an almost Optimal Monte Carlo (MAO) algorithm.

Theorem 3.2. Using the above choice 𝑝={𝑝_𝛼}^𝑛_𝛼=1 and 𝑃={𝑝_𝛼𝛽}^𝑛_𝛼,𝛽=1 the variance of the unbiased estimator for obtaining the inverse matrix is minimized [4].

There is a global algorithm that evaluates the solution of system (3.6) for every matrix 𝐴. The complexity of algorithm is 𝑂(𝑛²𝑙𝑁), where 𝑙 and 𝑁 are the average length of Markov chian and the number of simulated paths, respectively [2].

4. Inverse Monte Carlo Iterative Algorithm (IMCI)

Inverse Monte Carlo iterative algorithm can be applied when 𝐴 is a nonsingular matrix. In this method, we calculate the following functional in each steps:𝐴𝑓_𝑗,ℎ_𝑗𝐵𝑓_j,ℎ_𝑗=𝐵𝑓_𝑗−1,ℎ_𝑗𝐵𝑓_𝑗,ℎ_𝑗. (4.1) It is more efficient that we first evaluate the inverse matrix using the Monte Carlo algorithm [1, 2, 4]. The algorithm can be realized as in Algorithm 2.

input:  𝐴∈𝑅𝑛×𝑛,𝑓0∈𝑅𝑛 begin  Starting from initial vector 𝑓0  For  𝑗=1,2,…  begin   Using global algorithm, calculate the sequence of Monte Carlo   iterations by solving the following system           𝐴𝑓𝑗=𝐵𝑓𝑗−1   Set        𝜆(𝑗)=⟨𝐴𝑓𝑗,ℎ𝑗⟩⟨𝐵𝑓𝑗,ℎ𝑗⟩=⟨𝐵𝑓𝑗−1,ℎ𝑗⟩⟨𝐵𝑓𝑗,ℎ𝑗⟩   Output:   Smallest eigenvector 𝜆1(𝑗), and corresponding eigenvector 𝑓𝑗.  end end

Algorithm 2 

5. Partitioning IMCI

Let the matrix 𝐴 be partitioned into four blocks 𝐴₁,𝐴₂,𝐴₃, and 𝐴₄, where 𝐴₁ and 𝐴₄ are square matrices of order 𝑝 and 𝑞 such that 𝑝+𝑞=𝑛:𝐴𝐴=¹𝐴₂𝐴₃𝐴₄. (5.1) By assumption that all the indicated matrix inversions are realized, it is easy to verify that𝐴⁻¹=𝐵𝐿𝑀𝑁, (5.2) where𝑁=𝐴₄−𝐴₃𝐴₁⁻¹𝐴₂⁻¹,𝑀=−𝑁𝐴₃𝐴₁⁻¹,𝐿=−𝐴₁⁻¹𝐴₂𝑁,𝐾=𝐴₁⁻¹−𝐴₁⁻¹𝐴₂𝑀. (5.3) Thus inverting a matrix of order 𝑛 comes down to inverting four matrices, of which two have order 𝑝 and two have order 𝑞, and several matrix multiplications. Therefore the basic Monte Carlo for solving 𝐴𝑓_𝑗=𝐵𝑓_𝑗−1 will be called as the dimension of matrix 𝐴andequals to threshold. This action causes the convergence acceleration. Now, we can use the following recursion algorithm to obtain the inverse of matrix 𝐴 (see Algorithm 3).

Partitioning inverse  (𝑆,𝑛) begin:   𝑛=rank(𝑆);𝑝=𝑛/2   𝐴=𝑆[1∶𝑝,1∶𝑝];𝐵=𝑆[1∶𝑝,𝑝+1∶𝑛]   𝐶=𝑆[𝑝+1∶𝑛,1∶𝑝];𝐷=𝑆[𝑝+1∶𝑛,𝑝+1∶𝑛]   𝑚= size (𝐴)  if 𝑚≤ threshold   𝐴𝐴  =  Monte Carlo procedure  (𝐴)  else begin:    𝐴𝐴 =  Partitioning  inverse  (𝐴,𝑚)    𝑁 =  Partitioning  inverse  (𝐷−𝐶∗𝐴𝐴∗𝐵)    𝑀=−𝑁∗𝐶∗𝐴𝐴;𝐿=−𝐴𝐴∗𝐵∗𝑁    𝐾=𝐴𝐴−𝐴𝐴∗𝐵∗𝑀    𝑆𝑆[1∶𝑝,1∶𝑝]=𝐾;𝑆𝑆[1∶𝑝,𝑝+1∶𝑛]=𝐿    𝑆𝑆[𝑝+1∶𝑛,1∶𝑝]=𝑀;𝑆𝑆[𝑝+1∶𝑛,𝑝+1∶𝑛]=𝑁  end end

Algorithm 3 

6. Finding More Than One Generalized Eigenvalues

Assume that an eigenvalue 𝜆₁ and its corresponding eigenvector 𝑣₁ have been computed using the partitioning IMCI algorithm. In the first step of the above algorithm, we deflate the matrix 𝐴 to the matrix 𝐵. Then, we repeat again the first step of the algorithm to obtain the dominant eigenvalue of 𝐵 which is the second dominant eigenvalue of 𝐴. Let 𝑝 values of eigenvalues of pencil (𝐴,𝐵) be computed. Suppose that 𝑉_𝑝 is a matrix such that the columns of 𝑉_𝑝 are 𝑝 vector of eigenvector of pencil (𝐴,𝐵), that is,𝑉_𝑝=𝜈₁,…,𝜈_𝑝, (6.1) where 𝜈_𝑖 is eigenvector corresponding eigevalue 𝜆_𝑖,𝑖=1,2,…,𝑝.

Now, let𝐴_𝑝,𝐵=𝐴+𝐵𝑉_𝑝Λ𝐵𝑉_𝑝^𝑇,𝐵, (6.2) whereΛ=diag𝛿_𝑖−𝜆_𝑖,𝛿_𝑖>𝜆_𝑝,𝑖=1,…,𝑝. (6.3) Hence, if we find the 𝑝th smallest eigenpair of pencil (𝐴,𝐵), then we can evalute 𝜆_𝑝+1, that is, (𝑝+1)th smallest eigenvalue of pencil (𝐴,𝐵).

7. Numerical Results

In this section, the experimental results for obtaining the three smallest eigenpairs outlined in Tables 1, 2, and 3. The numerical tests are performed on Intel(R) (Core(TM)2 CPU, 1.83 GHz) personal machine.

Table 1 

Number of chains 𝑁=80.

Table 2 

The solution when the number of chains increases.

Table 3 

Total computational time for general and partitioning methods.

8. Conclusion and Future Study

We have seen that Monte Carlo algorithms can be used for finding more than one eigenpair of generalized eigenvalue problems. We analyze the computational complexity, speedup, and efficiency of the algorithm in the case of dealing with sparse matrices. Finally, a new method for computing eigenpairs as the partitioned method is presented. In Figure 1 the comparison of computational times between general Monte Carlo algorithm and partitioning algorithm is shown. The scatter diagram in Figure 2, shows that there is a linear relationship between matrix dimension (equivalently, the number of matrix elements) and total computational time for partitioning IMCI.

Figure 1 

Computational times for general and partition methods.

Figure 2 

Regression function 𝑦=0.0000658008𝑥+0.1048607485.