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

Authors: Behrouz Fathi Vajargah; Farshid Mehrdoust
Journal: Advances in Numerical Analysis
 (2011)

---

## 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.


---

## Body

## 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 A,B∈ℜn×n be real symmetric matrices and the matrix B a positive definite matrix. Consider the problem of evaluating the eigenvalues of the pencil (A,B), that is, the values for which(1.1)Ax=λBx.
A generalized eigenvalue problem (1.1) is said to be symmetric positive definite (S/PD) if A is symmetric and B is positive definite.


## 2. Inverse Vector Iteration Method

Another procedure for eigenvalue prediction is to use the Rayleigh quotient given by [3](2.1)λ=μ(x)=xTAxxTBx.
since B is positive definite, then (2.1) is well defined.Theorem 2.1.
Suppose thatλmin=λ1<λ2≤⋯≤λn-1<λn=λmax are n eigenvalues for pencil (A,B) and v1,…,vn, corresponding eigenvectors. Then for arbitrary initial vector x one has [3]
(2.2)λ1<μ(x)<λn,
where μ(x) is as introduced in (2.1).Theorem 2.2.
The inverse vector iteration method for arbitrary choice vectorx0 is convergent to the smallest eigenvalue and corresponding eigenvector for pencil (A,B). Also, the rate of convergence depends on 𝒪(λ2/λ1)k, where k is number of iterations [3].Algorithm1 evaluates the smallest eigenpair based on the inverse vector iteration [4].Algorithm 1
Input:Initial vectorx0beginSetλ1(1)=x0TAx0x0TBx0For  j=0,1,2,…beginSolve linear systemAzj+1=Bxj for zj+1Setλ1(j+1)=zj+1TAzj+1zj+1TBzj+1Setxj+1=zj+1∥zj+1∥Output:  xj,λ1(j)endend


## 3. Monte Carlo Method for Matrix Computations

Suppose that the matrixA∈ℜn×n and two vectors f,h∈ℜn are given. Consider the following Markov chain Ti with length i:(3.1)Ti:k0⟶k1⋯⟶ki,
where for j=1,…,i,kj∈{1,2,…,n}. The statistical nature of constructing the chain (3.1) follows as(3.2)p(k0=α)=pα,p(kj=β∣kj-1=α)=pαβ,
where pα and pαβ show the probability of starting chain at α and transition probability from state α to β, respectively.In fact(3.3)∑α=1npα=1,∑β=1npαβ=1,α=1,2,…,n.
Define the random variable Wj using the following recursion for(3.4)W0=1,Wj=Wj-1akj-1kjpkj-1kj,j=1,2,…,i.
Now, define the following random variable:(3.5)Θ[h]=hk0pk0∑j=0∞Wjfkj.Theorem 3.1.
Consider the following system:(3.6)Ax=b.Let the nonsingular matrix M∈ℜn,such that MA=I-L, then the system (3.6) can be presented in the following form:
(3.7)x=Lx+f,
where f=Mb. Then under condition Maxi∑j=1n|lij|<1,one has [5]
(3.8)E{Θ[h]}=〈h,x〉.Suppose thatx(i) is the ith iterative solution of the following recursion relation with x(0)=f. If we set the random variable(3.9)Θi[h]=hk0pk0∑j=0iWjfkj,
then(3.10)E{Θi[h]}=〈h,x(i+1)〉.By simulatingN random paths with length i(3.11)Ti(s):k0(s)⟶k1(s)⟶⋯⟶ki(s),s=1,2,…,N,
we can find(3.12)Θi(s)(h)=hk0pk0(s)∑j=0iWj(s)fkj,s=1,…,N.
The Monte Carlo estimation can be evaluated by(3.13)Θi=1N∑s=1NΘi(s)(h)
which is an approximation of 〈h,x(i+1)〉.From all possible permissible densities, we apply the following:(3.14)pα=|hα|∑α=1n|hα|,pαβ=|aαβ|∑β=1n|aαβ|,α=1,2,…,n.
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 choicep={pα}α=1n and P={pαβ}α,β=1n the variance of the unbiased estimator for obtaining the inverse matrix is minimized [4].There is aglobal algorithm that evaluates the solution of system (3.6) for every matrix A. The complexity of algorithm is O(n2lN), where l and N 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 whenA is a nonsingular matrix. In this method, we calculate the following functional in each steps:(4.1)〈Afj,hj〉〈Bfj,hj〉=〈Bfj-1,hj〉〈Bfj,hj〉.
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.Algorithm 2
input:A∈Rn×n,f0∈RnbeginStarting from initial vectorf0For  j=1,2,…beginUsingglobal algorithm, calculate the sequence of Monte Carloiterations by solving the following systemAfj=Bfj-1Setλ(j)=〈Afj,hj〉〈Bfj,hj〉=〈Bfj-1,hj〉〈Bfj,hj〉Output:Smallest eigenvectorλ1(j), and corresponding eigenvector fj.endend


## 5. Partitioning IMCI

Let the matrixA be partitioned into four blocks A1,A2,A3, and A4, where A1 and A4 are square matrices of order p and q such that p+q=n:(5.1)A=(A1A2A3A4).
By assumption that all the indicated matrix inversions are realized, it is easy to verify that(5.2)A-1=(BLMN),
where(5.3)N=(A4-A3A1-1A2)-1,M=-NA3A1-1,L=-A1-1A2N,K=A1-1-A1-1A2M.
Thus inverting a matrix of order n comes down to inverting four matrices, of which two have order p and two have order q, and several matrix multiplications. Therefore the basic Monte Carlo for solving Afj=Bfj-1 will be called as the dimension of matrix Aandequals to threshold. This action causes the convergence acceleration. Now, we can use the following recursion algorithm to obtain the inverse of matrix A (see Algorithm 3).Algorithm 3
Partitioning inverse  (S,n)begin:n=rank(S);p=n/2A=S[1:p,1:p];B=S[1:p,p+1:n]C=S[p+1:n,1:p];D=S[p+1:n,p+1:n]m= size (A)ifm≤ thresholdAA  =  Monte Carlo procedure  (A)elsebegin:AA =  Partitioning  inverse  (A,m)N =  Partitioning  inverse  (D-C*AA*B)M=-N*C*AA;L=-AA*B*NK=AA-AA*B*MSS[1:p,1:p]=K;SS[1:p,p+1:n]=LSS[p+1:n,1:p]=M;SS[p+1:n,p+1:n]=Nendend


## 6. Finding More Than One Generalized Eigenvalues

Assume that an eigenvalueλ1 and its corresponding eigenvector v1 have been computed using the partitioning IMCI algorithm. In the first step of the above algorithm, we deflate the matrix A to the matrix B. Then, we repeat again the first step of the algorithm to obtain the dominant eigenvalue of B which is the second dominant eigenvalue of A. Let p values of eigenvalues of pencil (A,B) be computed. Suppose that Vp is a matrix such that the columns of Vp are p vector of eigenvector of pencil (A,B), that is,(6.1)Vp=[ν1,…,νp],
where νi is eigenvector corresponding eigevalue λi,i=1,2,…,p.Now, let(6.2)(Ap,B)=(A+(BVp)Λ(BVp)T,B),
where(6.3)Λ=diag⁡{δi-λi},δi>λp,i=1,…,p.
Hence, if we find the pth smallest eigenpair of pencil (A,B), then we can evalute λp+1, that is, (p+1)th smallest eigenvalue of pencil (A,B).


## 7. Numerical Results

In this section, the experimental results for obtaining the three smallest eigenpairs outlined in Tables1, 2, and 3. The numerical tests are performed on Intel(R) (Core(TM)2 CPU, 1.83 GHz) personal machine.Table 1
Number of chainsN=80.
DimensionEigenvalues errorEigenvectors errorλ1(1)λ1(2)λ1(3)v1(1)v1(2)v1(3)64×641.26×10-74.5×10-71.26×10-68.36×10-4.002813.008154128×1281.18×10-74.46×10-71.33×10-66.33×10-4.002534.008074256×2566.48×10-84.37×10-71.51×10-63.16×10-4.002151.008070512×5127.91×10-82.87×10-71.72×10-61.83×10-4.001284.0077591024×10247.04×10-81.99×10-71.79×10-61.06×10-4.000712.00769Table 2
The solution when the number of chains increases.
Number of chainsCalculated eignvaluesNExactλ1(1)=.74529395Exactλ1(2)=.85537131Exactλ1(3)=.9103274020.74737713.85513704.9100823840.74529979.85537157.9103430480.74529407.85537175.9102873Table 3
Total computational time for general and partitioning methods.
DimensionTime (Sec.)General methodPartitioning method20×20.21.2040×40.40.3260×60.71.5480×801.25.60100×1001.83.75120×1202.541.18140×1403.671.57160×1604.422.04180×1805.582.54200×2007.203.18


## 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 Figure1 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 functiony=0.0000658008x+0.1048607485.



---

*Source: 826376-2011-04-11.xml*