# A Cost-Effective Smoothed Multigrid with Modified Neighborhood-Based Aggregation for Markov Chains

Authors: Zhao-Li Shen; Ting-Zhu Huang; Bruno Carpentieri; Chun Wen
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

Smoothed aggregation multigrid method is considered for computing stationary distributions of Markov chains. A judgement which determines whether to implement the whole aggregation procedure is proposed. Through this strategy, a large amount of time in the aggregation procedure is saved without affecting the convergence behavior. Besides this, we explain the shortage and irrationality of the Neighborhood-Based aggregation which is commonly used in multigrid methods. Then a modified version is presented to remedy and improve it. Numerical experiments on some typical Markov chain problems are reported to illustrate the performance of these methods.


---

## Body

## 1. Introduction

Markov chains have a large number of applications for scientific research and technological optimization in many areas, including queuing systems and statistical automata networks, web page ranking [1, 2] and gene ranking [3], risk management, and customer relationship analysis. Markov chains also have an important application for the noise and signal analysis, processing of voice, image, and screen, analog of communications network and transport phenomena, economic calculation methods, and so forth. For the needs of these practical applications, there arise different types of Markov chain models, such as the hidden Markov chain [4], multivariate Markov chain, and high-order Markov chain [5]. For Markov models, computing the stationary distribution is an important issue.The transition matrix of a homogeneous irreducible Markov chain withn states can be denoted as B. B∈Rn×n is a stochastic matrix and commonly defined as a column-stochastic matrix; that is, for every 1≤i, j≤n,(1)0≤bi,j≤1,eTB=eT,where column vector e is equal to one.Then the stationary distribution problem is mathematically given by(2)Bx=x,x1=1,x≥0.In this study, we only consider the situation whereB is also an irreducible matrix, which means, in its directed graph, there exists a directed path from any vertex i to any other vertex j. If B is irreducible, according to the Perron-Frobenius theorem [6], there exists a unique solution x to the linear system (2).This linear problem (2) is often reformulated as(3)Ax=0,x1=1,x≥0with(4)A=I-B,where I is an identity matrix. As Sterck et al. showed in their work [7], this formulation has two specific advantages. First, the solution x is not only the right eigenvector of A corresponding to eigenvalue 0, but also the right singular vector of A corresponding to singular value 0. The second advantage of formulation (3) is that the M-matrix structure of A is amenable to additive algebraic multigrid solvers designed for nonsingular linear problems [8, 9]. Hence, our interest is to solve the linear system (3), corresponding to a homogeneous irreducible Markov chain.In many applications, the size of the Markov chain is very large, leading to the fact that direct solvers are impractical. Hence, the iterative procedures are widely used to calculate the stationary probability of Markov chains. However, traditional one-level iterative methods, for example, the Power method, (weighted) Jacobi method, and (weighted) Gauss-Seidel method, converge very slowly for calculating the solutionx in the situations where the subdominant eigenvalue λ2 of B satisfies λ2~1 [10, 11]. The Markov chain in this case is a slowly mixing type.Recently, multigrid iterative methods aiming to accelerate the convergence rate for these types of Markov chains have received much attention. The multigrid methods for Markov chains have been developed and applied very successfully since the earliest work of Takahashi’s two-level iterative aggregation/disaggregation method for Markov chains [12]. Direct extension of two levels to multilevels was first explored in [13]. Then, adaptive algebraic multigrid methods whose aggregates are formed algebraically based on the strength of connection in the problem matrix column-scaled by the current iterate were developed in [14, 15], and for Markov chains [16]. This method shows good potential. Following this idea, there have arisen some techniques including smoothing the interpolation and restriction operators [17], using acceleration on all recursive levels [18] or on the top level [2, 19], the on-the-fly interactive strategy [20], and so forth.The numerical experiments in [18] show that the aggregation strategy has a certain effect on the performance of these adaptive algebraic multigrid methods. Employing reasonable and efficient aggregation strategy can improve the performance and applicability of these multigrid methods. However, in some cases, we found that the Neighborhood-Based aggregation used commonly by the above-mentioned methods generates aggregates with many states of very weak connection between each other. These unreasonable aggregates often lead to a poor convergence rate of the algebraic multigrid methods.To remedy these problems and improve the applicability of Neighborhood-Based aggregation, this paper propose a modified version which changes some rules.Besides the effect on the convergence rate, the computational cost of aggregation procedure accounts for a large proportion of the whole cost of these adaptive algebraic multigrid methods [20]. We study the aggregation method which is based on the connection strength and propose a judge mechanism to capture in advance the situations where the aggregates will be the same as previous. Then there is no need to implement the whole aggregation method. This strategy keeps the convergence rate the same with computing time reduced.The remainder of this paper is organized as follows. In Section2, the smoothed aggregation multigrid (SAM) method for the Markov chains from [17] is reviewed. In Section 3, we analyse and demonstrate the drawback of the Neighborhood-Based aggregation method and present the modified version. The judge mechanism and the cost-effective SAM are presented in Section 4. Numerical tests are presented in Section 5. Finally, conclusions are reported in Section 6.


## 2. Smoothed Aggregation Multigrid (SAM) Method for Markov Chains

First, we briefly review the framework of SAM method for Markov chains from [17], where the hierarchy of coarse grids is developed based on the structure of the scaled problem matrix Adiag⁡(xi), instead of relying on the geometry of the original problem. Here, we follow the presentation of [17].Multilevel methods aim to reduce the algebraically smooth error left by the classical iterative methods. So classical iterative methods like Power method, (weighted) Jacobi method, and (weighted) Gauss-Seidel method are often used as relaxation techniques for multilevel methods.In this study, the weighted Jacobi iterative method is used as relaxation for (3) and it is given by(5)x⟵I-ωD-1Ax,where D is the diagonal matrix of A.System (3) can be rewritten into the following form:(6)Adiag⁡xiei=0,where xi is the approximate got by last multilevel cycle and current relaxations and ei is the multiplicative error vector and should be all ones when convergence is achieved.Then the aggregation matrixQ∈Rn×m needs to be constructed, where qij=1 if the node i belongs to aggregate j and qij=0 otherwise. Each column of Q represents an aggregate, and each row has only one component equal to 1 while others are 0. For example, if the nodes are ordered according to the aggregates they belong to, the matrix Q has the form like(7)Q=100⋯100⋯100⋯010⋯010⋯010⋯001⋯001⋯001⋯⋮⋮⋮⋱.The aggregation algorithm for constructing Q is discussed in Section 3.WhenQ is constructed, the coarse-level equation of (6) is given by(8)QTAdiag⁡xiQec=0,where ec is the coarse-level multiplicative error corresponding to ei.Then, the restriction operatorR and the prolongation operator P are defined by R=QT and P=diag⁡xiQ. Equation (8) is rewritten as(9)RAPec=0,and the coarse-level operator Ac can be defined as(10)Ac=RAP.Ifec is known, then an improved coarse-level approximation can be obtained by(11)xc=diag⁡PT1ec.With (11), (9) is rewritten as the coarse-level probability equation:(12)Acdiag⁡PT1-1xc=0.Finally, (12) is used to accelerate the classical one-level iterative methods like Power method, (weighted) Jacobi method, or (weighted) Gauss-Seidel method for (3). We take the two-level case, for example.A two-level aggregation/disaggregation method is applied after some relaxations on (3) at the fine level. Through constructing the aggregation Q, the coarse-level equation (12) is formatted, and we solve (12) to get the coarse-level approximation xc. Then, a coarse-level correction(13)xi+1=Pec=Pdiag⁡PT1-1xcis implemented followed by some relaxations again to obtain the improved approximation at fine level.Two-level iterative aggregation/disaggregation can be naturally extended to multiple levels by recursively applying the two-level method to solve the coarse-level probability equation (12). There are many types of recursive cycles such as V cycle which solves the coarse-level probability equation (12) once at each intermediate level of a cycle and W cycle which solves it twice. Here we only describe V cycle. It is illustrated by Figure 1.Figure 1
Multilevel V cycle: the finest-level operations are represented at the top, the coarsest-level operations are at the bottom, and intermediate-level operations are in between. The black dots represent relaxation operations on corresponding levels and the open dots represent direct solvers.Performance of this aggregation multigrid method for Markov chains is often unsatisfactory. Inspired by smoothed aggregation (SA) multigrid methods for linear systems, smoothing the interpolation and restriction operators,P and R, is proposed to improve the convergence behavior of the aggregation multigrid method [17].All columns of the interpolation operator,P, are smoothed by (weighted) Jacobi method with weight ω:(14)Ps=I-ωD-1AP.For the restriction operator, R, all rows are smoothed as (15)Rs=RI-ωAD-1.The smoothed operators Ps and Rs still hold the properties: (16)1cTRs=1T,Ps1c=x,forxi=x.When Ps and Rs are generated, the smoothed coarse-level operator is naturally defined as (17)Acs=RsAPs.Then the coarse-level probability equation is given by(18)Acsdiag⁡PsT1-1xc=0.Note that because the smoothed coarse-level operatorAcs may lose the properties of an irreducible singular M-matrix, a lumping technique is necessary to modify Acs to keep these properties. This technique is implemented before the coarse-level probability equation; refer to [17] for details. This smoothed aggregation multigrid (SAM) method for Markov chain can be formulated as Algorithm 1.Algorithm 1:SAM method for Markov chain (V cycle), x← SAM(A,x,v1,v2).
Input:
A: problem matrix,
x: iteration solution,
v1, v2: pre and post relaxation steps.
Output:
x: iteration solution.
(1) if not on the coarsest level then
(2)  x  ← Relax(A,x) v1 times.
(3)  Construct Q.
(4)  R←QT and P←diag⁡(x)Q.
(5)  Rs  ← smooth R, Ps← smooth P.
(6)  Ac  ← lumping RsAPs, xc←Ps′1.
(7)  xc=SAMAcdiag⁡(xc)-1,xc,v1,v2
(8)  x←Psdiag⁡PsT1-1xc.
(9)  x  ← Relax(A,x) v2 times.
(10) else
(11)  x← the solution of Ax=0, x1=1, x>0 by a direct solver.
(12) end if
(13) return  x;It is necessary to apply a direct solver for the coarsest level and the direct solver used in Algorithm1 is based on the following theorem.Theorem 1 (Theorem  4.16 in [21]).
IfA is an irreducible singular M-matrix, then each of its principal submatrices other than A itself is a nonsingular M-matrix.LetAL denote the coarsest-level operator; then we use the direct method presented in the coarsest-level Algorithm 2 to solve the coarsest-level equation ALxL=0.Algorithm 2:Coarsest-level algorithm.
Input:
A: problem matrix,
Output:
x: the solution with x1=1.
(1) Compute NL≔size(AL,1);
(2) Determine ALp≔AL(1:NL-1,1:NL-1);
(3) Determine bLp≔-AL(1:NL-1,NL);
(4) Compute xLp(1:NL-1)≔ALp∖bLp; let xLp(NL)=1;
(5) Set the coarsest-level solution x=xLp/xLp1.
(6) return  x;


## 3. Neighborhood-Based Aggregation and Our Modified Version

How to construct the aggregation matrixQ is a crucial issue in algebraic multigrid method. Here, we consider the aggregation methods which determine aggregates based on strength of connection in matrices. There exist some corresponding aggregation methods like distance-one aggregation, distance-two aggregation [16, 17], pairwise aggregation, double pairwise aggregation [22], Neighborhood-Based aggregation [23], bottom-up aggregation [24], and some other types of aggregation. Serving as the aggregation method in [7, 18, 19], Neighborhood-Based aggregation has some advantages over others. It usually generates well-balanced aggregates of approximately equal size and reduces the number of unknowns quickly in coarsening process. Coarse-level stencil sizes tend to be rather uniform and do not grow too quickly [18].Here, the aggregation procedure is based on matrixA-=Adiag⁡(xi) and a node i is said to be strongly connected to node j in the graph of A- if(19)-a-ij≥θmaxk≠i⁡-a-ikor-a-ji≥θmaxk≠j⁡-a-jk,where θ is a strength of connection parameter. The strong neighborhood Ni of node i is the set of all points which are strongly connected to i in the graph of A- including node i itself. The Neighborhood-Based aggregation can be represented as Algorithm 3.Algorithm 3:Neighborhood-Based aggregation, {Qj}j=1J←Neighborhood-Based(A-,θ).
Input:
A-: scaled problem matrix,
θ: tolerance.
Output:
Q: aggregation matrix.
1st process: generate the initial aggregates
(1) Set R←{1,2,…,n} and J←0.
(2) for  i∈{1,2,…,n}  do
(3)  construct strong neighborhoods Ni based on (19).
(4)  if  Ni∩R=Ni  then
(5)   J←J+1.
(6)   QJ←Ni
(7)   R←R∖Ni.
(8)  end if
(9) end for
2nd process: connect the remaining nodes to their most connected aggregates
(10) while  R≠∅  do
(11) pick i∈R and set j=argmaxk=1,…,J⁡card(Ni∩Qk).
(12) Qj←Qj∪{i} and R←R∖{i}.
(13) end while
3rd process: construct the aggregation matrix  Q∈Rn×J
(14) for  i∈{1,2,…,n}  do
(15) if  i∈Qj, j=1,2,…,J  then
(16)  Qij=1.
(17) else
(18)  Qij=0.
(19) end if
(20) end for
(21) return  Q;Though the Neighborhood-Based aggregation has many advantages, it is somewhat unreasonable in some detail.In Algorithm3, it picks i∈{1,2,…,n} in a natural order and defines a new aggregate denoting the (J+1)th aggregate as the strong neighborhood of node i, if such strong neighborhood is contained in the left set: R∖⋃j=1JQj.As for a strong neighborhoodNi of some node i, from (19) we can see that if the maxk≠i{-a-ik} is relatively small, then the strength of connection from other nodes in Ni to node i will be weak. Even worse, the nodes in Ni except i may be weakly connected with many others universally; that is to say, -a-jk is relatively small for many j,k∈Ni.And if such nodei is in the beginning of {1,2,…,n}, it would be aggregated earlier than most other nodes, and then it is very likely to construct an aggregate using such a neighborhood Ni described above.This aggregate constructed throughNi is unreasonable. First, it contains many nodes which are neither strongly connected with i nor strongly connected with each other. Secondly, because it is aggregated earlier, it may separate some nodes which are truly strongly connected.Here, we take a toy example to illustrate.Example 2.
LetA- be an 8×8 matrix:(20)0.07-0.01-0.01-0.01-0.01-0.01-0.01-0.01-0.010.94-0.3-0.3-0.3-0.01-0.01-0.01-0.01-0.30.94-0.3-0.3-0.01-0.01-0.01-0.01-0.3-0.30.94-0.3-0.01-0.01-0.01-0.01-0.3-0.3-0.30.94-0.01-0.01-0.01-0.01-0.01-0.01-0.01-0.010.95-0.45-0.45-0.01-0.01-0.01-0.01-0.01-0.450.95-0.45-0.01-0.01-0.01-0.01-0.01-0.45-0.450.95.We can write this matrix as its block form:(21)a-11OOT2–5OOT6–8,where T2–5 denotes the square block containing elements in the columns and rows from 2nd to 5th of A-, T6–8 is defined correspondingly, and O means the elements of it are relatively very small.This block form (21) can directly drive the reasonable aggregates for matrix A-. With the restriction that each aggregate must have at least two nodes, the reasonable aggregates should be constructed as follows: 2nd to 5th nodes form an aggregate; the 6th to 8th nodes form another aggregate; the first node is added to either of these two aggregates and where it is added makes little difference.However, let us check which aggregates of this matrix will be constructed by Neighborhood-Based algorithm. We define the strong connection matrixS as follows:Fori≠j,(22)sij=1,if-a-ij≥θmaxk≠i⁡-a-ikor-a-ji≥θmaxk≠j⁡-a-jk,0,otherwise.Fori=j,(23)sii=1.So,Sij=1 means that node j is strongly connected with node i and is included in Ni. S is a symmetric matrix and the ith row or ith column represents the strong neighborhood Ni. In this example, S is(24)S=1111111111111000111110001111100011111000100001111000011110000111.In Neighborhood-Based Algorithm3, nodes are picked in a natural order, so the 1st node is first picked and its strong neighborhood N1 is used as 1st aggregate. It is seen from S in (24) that all the other nodes are strongly connected to node 1. Adding 1st node itself, the strong neighborhood N1 contains all the nodes of A-. The 1st aggregate is constructed containing all the nodes and the algorithm ends. It is too aggressive and unreasonable.We consider that picking nodes in a natural order to construct aggregates is unreasonable. Here we consider picking nodes in another order.We define a measurementsize_nei(i) of node i as(25)size_neii≔lengthNi,which means the number of nodes in the strong neighborhood Ni of i.We can suppose that the more the nodes the strong neighborhoodNi contains, the more likely the probabilities to be distributed evenly and be relatively smaller. So, we preferentially pick the nodes with smaller size_nei value to construct aggregates. We apply this strategy by reordering R by size_nei value to make the reordered set increasing in size_nei. This strategy can avoid generating the unreasonable aggregate in the 1st process: generate the initial aggregates, as the example shows.However the Neighborhood-Based aggregation also has some defects in the 2nd process that may lead to bad aggregates.As seen from the 2nd process of Algorithm3, it picks the remaining nodes of 1st process and adds them to the “most connected” aggregates. We conclude that the reason one node i is left is because its strong neighborhood Ni has some common nodes with at least one of the existing aggregates and the “most connected” aggregate is defined as the aggregate which has the most common nodes with Ni.However, when the algorithm considers which aggregatei should be added to, it ignores Ni itself. Consider a situation that Ni has a small number of common nodes with the existing aggregates compared to size of Ni. If we still add i to the “most connected” aggregate as Algorithm 3 does, it may cause some bad aggregates containing many nodes which should be separated because of the weak connection between each other.We use Example2 to illustrate again.We change the order of picking nodes to avoid the problem caused by the 1st process discussed above. First the 6th node should be picked to construct aggregate in 1st process of Algorithm3; because N6 contains the 1st node and nodes from 6th to 8th, the 1st aggregate is constructed as {1,6,7,8}. Then we pick nodes from 2nd to 5th. It can be seen from S in (24) that the strong neighborhoods of them are the same and contain the nodes from 1st to 5th. These strong neighborhoods have a common node with the 1st aggregate: 1st node, so the nodes from 2nd to 5th are all remaining nodes and should be added to the “most connected” aggregate. According to 2nd process of Algorithm 3, each node should be added to the aggregate which has most common nodes with its strong neighborhood. Thus, the result of the algorithm is also an aggregate with all the nodes, still too aggressive and unreasonable.Here we propose an improved scheme to the 2nd process in Algorithm3. Let RNi denote the remaining strong neighborhood of the remaining note i as(26)RNi=R∩Ni,which means removing the common nodes with the existing aggregates from Ni. And we consider RNi as a candidate aggregate. Then we find the aggregates with most common nodes with Ni as previous. If the aggregate founded is in the existing aggregates, we add the node i to the corresponding aggregate. Otherwise, we use the remaining strong neighborhood RNi  which added the node i as a new aggregate, the (J+1)th aggregate (J is the number of existing aggregates).The whole algorithm is represented in Algorithm4.Algorithm 4:Modified Neighborhood-Based aggregation, Q←MN-BA(A-,θ).
Input:
A-: scaled problem matrix,
θ: tolerance.
Output:
Q: aggregation matrix.
1st process: generate the initial aggregates
(1) Set R←{1,2,…,n} and J←0.
(2) for  i∈{1,2,…,n}  do
(3)  compute size_nei(i) based on (25).
(4) end for
(5) Reorder R by increasing size_nei values:
(6) R← new reordered R.
(7) for  i∈R  do
(8)  construct strong neighborhoods Ni based on (19)
(9)  if  Ni∩R=Ni  then
(10)  J←J+1.
(11)  QJ←Ni.
(12)  R←R∖Ni.
(13)  end if
(14) end for
2nd process: connect the remaining nodes to their most connected aggregates
(15) while  R≠∅  do
(16)  Pick i∈R and set RNi according to (26)
(17)  QJ+1←RNi∖{i}
(18)  Set j=argmaxk=1,…,J+1⁡card(Ni∩Qk).
(19)  if  j==J+1  then
(20)   J=J+1.
(21)    Qj←Qj∪{i} and R←R∖RNi.
(22)  else
(23)   Qj←Qj∪{i} and R←R∖{i}.
(24)   QJ+1={}
(25)  end if
(26) end while
3rd process: construct the aggregation matrix  Q∈Rn×J
(27) for  i∈{1,2,…,n}  do
(28)  if  i∈Qj,j=1,2,…,J  then
(29)   Qij=1.
(30)  else
(31)   Qij=0.
(32)  end if
(33) end for
(34) return  Q;Note that the main extra work of our modified Neighborhood-Based Algorithm4 compared to the old version consists of reordering the set R which can be solved in O(log2⁡n) work in the first process, and the loss of the parallelization of the second process.


## 4. Cost-Effective SAM for Markov Chains

The smoothed aggregation multigrid for Markov chains described in Algorithm1 is one type of adaptive algebraic multigrid algorithm. It constructs aggregation matrix Q based on the strength of connection in the scaled problem matrix A-=Adiag⁡xi, without any need for advance knowledge of the topology of the Markov chain. The aggregation process is fully adaptive, with aggregates constructed adaptively in each iteration and at all recursive levels. This type of multigrid method is powerful and scalable for many Markov chain problems compared to the classical AMG which use the operators without any update. But they are relatively more expensive because they construct the whole multigrid operators in every cycle and unfortunately the computational cost of constructing operators is accounted for a large proportion of the whole cost in a cycle. In [24], it is noted that, empirically, more than 50% of the computation time of each cycle is spent on coarse matrix construction. And through our observation, the most costly procedure is the aggregation procedure which constructs the matrix Q.In the numerical experiments, there are some cases where the aggregation matrixQ at a certain level of current cycle is the same as the corresponding one in the previous cycle, so implementing the aggregation method to update Q is totally a waste. The Neighborhood-Based aggregation and our modified version depend on the strong neighborhood Ni of each node i and the ordering of picking nodes. And because the strong neighborhood Ni just depends on the strong connection matrix S in (22) and so does the reordering in our modified version, it can conclude that, with the same aggregation matrix S, the results after Neighborhood-Based aggregation or our modified version would be the same.Actually, in Neighborhood-Based aggregation, it is necessary to construct strong neighborhoodNi of each node i and that is equal to constructing the strong connection matrix S because Ni=S(i,:) according to the definition of Ni and S. So there is no extra work to construct S and there is only a need to add a judgment step. If the current S is the same as the corresponding one in previous cycle, there is no need to implement the remaining part of Neighborhood-Based aggregation and the corresponding aggregation matrix Q can still be used, leading to much time being saved in this procedure. Meanwhile, the storage will be increased; for that we also need to store the strong connection matrix S and the aggregation matrix Q at each level of previous cycle. The aggregation matrix Q is very sparse with only n nonzeros; n is the size of current level. The strong connection matrix S is symmetric, only the upper (lower) triangular part and the diagonal elements should be stored. The whole algorithm is presented in Algorithm 5.Algorithm 5:Cost-effective SAM for Markov chains (CESAM), x←CESAM(A,x,v1,v2).
Input:
A: problem matrix,
x: iteration solution,
v1, v2: pre and post relaxation steps.
Output:
x: iteration solution.
(1) if not on the coarsest level then
(2)  x← Relax(A,x) v1 times.
(3)  Construct S according to (22).
(4)  if in the first cycle then
(5)   Spre←S.
(6)   Construct Q based on S, and Qpre←Q.
(7)  else if  S==Spre  then
(8)   Q←Qpre.
(9)  else
(10)  Spre←S.
(11)  Construct Q based on S, and Qpre←Q.
(12) end if
(13) R←QT and P←diag⁡(x)Q.
(14) Rs← smooth R, Ps← smooth P.
(15) Ac← lumping RsAPs, xc←Ps′1.
(16) xc=CESAMAcdiag⁡xc-1,xc,v1,v2
(17) x←Psdiag⁡PsT1-1xc
(18) x← Relax(A,x) v2 times.
(19) else
(20)x← the solution of Ax=0, x1=1, x>0 by a direct solver.
(21) end if
(22) return  x;Note that the aggregation methods based on the connection strength including Neighborhood-Based aggregation, our modified Neighborhood-Based aggregation, and distance-one (distance-two) aggregation are all suitable for Algorithm5. And this cost-effective SAM (Algorithm 5) is mathematically equal to SAM (Algorithm 1).


## 5. Numerical Experiments

In this section, we report numerical results obtained using a Matlab R2014a implementation on 64-bit Windows-7 with a core i5-3470 processor and 8 GB RAM memory.One purpose of our numerical tests is to examine whether our modified Neighborhood-Based aggregation can generate aggregates in a more reasonable way and bring accuracy improvement compared to the original version and to see if it is applicable for some case where the original Neighborhood-Based aggregation performs badly. Another purpose is to test the efficiency of our cost-effective SAM (referred as CESAM) for Markov chains. We also compared CESAM method with two standard methods, Power method and Gmres(50) method, to demonstrate the efficiency of it. Four Markov chain problems studied in [7] including 2 structured problems and 2 random walk on unstructured planar graph, two Markov chain examples from queueing models, and a typical Markov chain constructed by ourselves like Example 2 in Section 3 have been employed in our experiments. In our experiments, we use V cycle only. For convenience, we let CESAM_O and SAM_O denote the applications of CESAM and SAM with original Neighborhood-Based aggregation, respectively, and CESAM_M and SAM_M denote the applications of CESAM and SAM with our modified Neighborhood-Based aggregation, respectively.Without loss of generality, special sets of the parameters are employed and they are taken from [7]. As mentioned early, the weighted Jacobi method is used as the pre- and postsmoothing in SAM (Algorithm 1) and CESAM (Algorithm 5). Let v1=v2=1 and set the relaxation parameter ω=0.7 in the experiments. Note that, even though the values of parameters work well for all tests considered here, they are likely to be problem-dependent. The coarsest-level solution is implemented by the coarsest-level Algorithm 2 when the size is below 50; the strength of connection parameter is chosen as θ=0.25. The initial guess is generated by random sampling with a uniform (0,1) distribution and normalized to one in the one norm for the first four Markov chains and is all equal to 1/n (n is the dimension of the Markov chain) in the last Markov chain. Iterations are terminated when res=Axk1/Ax01≤10-12 with xk the current approximate solution and x0 the initial guess solution. The max number of iterations is set as 200 for SAM and CESAM methods, 106 for Power method, and 5000 for Gmres(50) method.Numerical results are reported in the tables, where “n” is the problem size, “it” denotes the iteration counts, and “CPUt” denotes the total computing time in seconds. “Cop” is the operator complexity of the last cycle, which is defined as the sum of the number of nonzero entries in all operator Al on all levels divided by the number of nonzero entries in the fine-level operator A1; that is, Cop=∑l=1L‍nnz(Al)/nnz(A1) with nnz(Al) being the number of nonzero entries in Al.


### 5.1. Structured Problems

#### 5.1.1. Uniform 2D Lattice

The first test problem is a Markov chain on 2D lattice with uniform weights. MatrixA is essentially a scaled graph-Laplacian on a 2D uniform quadrilateral lattice with 5-point stencil (see Figure 2). It is a typically structured problem. Tables 1 and 2 give the numerical results.Table 1
Numerical results for 2D lattice with uniform weights.
n
SAM_O
SAM_M
CESAM_O
CESAM_M
C
op
it
CP
U
t
C
op
it
CP
U
t
it
CP
U
t
it
CP
U
t
4096
1.34
32
1.532
1.39
31
1.643
32
0.554
31
0.642
16384
1.33
36
6.338
1.38
31
6.121
36
1.751
31
2.094
65536
1.34
37
26.477
1.38
31
24.882
37
8.992
31
8.837Table 2
Numerical results for 2D lattice with uniform weights.
n
CESAM_M
Power
Gmres(50)
it
CP
U
t
it
CP
U
t
it
CP
U
t
4096
31
0.642
>106
>254.7
734
2.264
16384
31
2.094
>106
>33.74
2087
17.44
65536
31
8.837
>106
>254.7
>5000
>108.6Figure 2
Graph for a 2D lattice with uniform weights.Table1 shows that, for this simple and highly structured problem, there is little difference in iterations between the Neighborhood-Based aggregation and our modified version. The operator complexity “Cop” of SAM_M is a little larger than that of SAM_O for the fact that our modified Neighborhood-Based aggregation is less aggressive than the original version. This larger operator complexity and the addition work of our modified Neighborhood-Based aggregation lead to a larger computing time per cycle. For our CESAM, it has the same “Cop” (which is not shown in the table) and iterations with SAM but converges in a much smaller time. The CESAM can save the calls of much work in aggregation procedure and this is problem-dependent. For example, for the size of 4096, the CESAM_O needs 35 times of implementing the remaining work of aggregation procedure while SAM_O needs 96 times. In addition, these savings most happen in finer level, resulting in a decrease of 63.8% proportion in computing time by CESAM_O when compared to SAM_O. It is similar for other sizes and for CESAM_M compared to SAM_M. Seen from Table 2, when compared to Power and Gmres(50) method, the CESAM_M shows better performance in both iterations and computing time.


#### 5.1.2. Tandem Queuing Network

The next test problem is an open tandem queuing network from [10]. Tandem queuing network exists in some practical applications, for example, in wireless networks [25] and blood screening procedures [26]. Following the settings and description in [17], we consider the situation that there are two finite queues with single servers, customers arrive according to a Poisson distribution with rate μ, and the service time distribution at the two single-server stations is Poisson with rates μ1 and μ2. In our numerical simulation, we set the number of customers in the queues as N=63,127,255. We set μ=10, μ1=11, and μ2=10. The states of this system can be represented by tuples (n1,n2), with n1 denoting the number of customers waiting in the first queue and n2 in the second queue. Then the total number of states is given by (N+1)2. It is also a structured problem and illustrated in Figure 3. The numerical results for this problem are presented in Tables 3 and 4.Table 3
Numerical results for tandem queuing network.
n
SAM_O
SAM_M
CESAM_O
CESAM_M
C
op
it
CP
U
t
C
op
it
CP
U
t
it
CP
U
t
it
CP
U
t
4096
1.24
34
1.459
1.27
34
1.576
34
0.450
34
0.442
16384
1.24
33
5.007
1.27
32
5.494
33
1.188
32
1.226
65536
1.24
43
26.51
1.25
33
22.69
43
6.529
33
5.190Table 4
Numerical results for tandem queuing network.
n
CESAM_M
Power
Gmres(50)
it
CP
U
t
it
CP
U
t
it
CP
U
t
4096
34
0.442
>106
>9.469
1652
4.631
16384
32
1.226
>106
>28.44
>5000
>40.53
65536
33
5.190
>106
>193.0
>5000
>86.82Figure 3
Graph for tandem queuing network.Seen from Table3, like the previous example, “Cop” of SAM with our modified Neighborhood-Based aggregation is still larger than that of SAM with Neighborhood-Based aggregation. This is mainly caused by the difference in the 2nd process of these two aggregation methods because these two test problems are very regular which means the strong neighborhood of each node is almost with the same size. For this problem, the number of iterations is the same for these methods when the size is 16384 or 65536. For the same reason explained above, our SAM_M costs more time per cycle than SAM_O and this leads to a larger computing time of sizes 4096 and 16384 where the iterations are almost the same. For the size of 65536, the number of iterations is reduced by a larger proportion and this offsets the larger computing time per cycle, resulting in a smaller computing time. Similar to the previous example, CESAM converges in a much smaller time than SAM, with 69.2%–76.3% proportion of time saved by CESAM_O compared to SAM_O and 72.0%–77.7% proportion of time saved by CESAM_M compared to SAM_M. Note that because the CESAM reduces the calls of much work in aggregation procedure, the higher cost of modified Neighborhood-Based aggregation has less impact. If the saving in calls is almost the same, the CESAM_M with less iterations will be more efficient than CESAM_O; this can be verified by “CPUt” of these two methods. The CESAM_M method still shows better performance than that of Power method and Gmres(50) method, which can be seen in Table 4.


#### 5.1.3. M/H2/1 Queues

The next test problem is M/H2/1 queues from [10]. M/H2/1 queues are used in mobile communications [27]. M/H2/1 queue can be illustrated in Figure 4. Arrivals are generated according to a Poisson distribution at rate λ, while the service is represented by a two-phase hyperexponential distribution. With probability α a customer entering service receives service at rate μ1, while with probability 1-α this customer receives service at rate μ2. The numerical results for this problem are presented in Tables 5 and 6.Table 5
Numerical results for M/H2/1 queues.
n
SAM_O
SAM_M
CESAM_O
CESAM_M
C
op
it
CP
U
t
C
op
it
CP
U
t
it
CP
U
t
it
CP
U
t
4001
1.43
32
1.706
1.25
28
1.507
32
0.438
28
0.374
20481
1.44
33
8.042
1.25
27
6.691
33
1.825
27
1.191Table 6
Numerical results for M/H2/1 queues.
n
CESAM_M
Power
Gmres(50)
it
CP
U
t
it
CP
U
t
it
CP
U
t
4001
28
0.374
>106
>9.906
>5000
>16.73
20481
27
1.191
>106
>39.36
>5000
>52.60Figure 4
Graph for M/H2/1 queues.Table5 shows that the iterations needed for SAM_M are less than those needed for SAM_O. This means the modified Neighbour-Based aggregation can generate more suitable aggregates than the original version. Less iterations also result in a smaller computing time of SAM_M and CESAM_M compared to SAM_O and CESAM_O, respectively. For the efficiency of the cost-effective strategy, we can see that CESAM converges in a much smaller time than SAM, with 74.3%–75.2% proportion of time saved by CESAM_O compared to SAM_O and 75.2%–82.2% proportion of time saved by CESAM_M compared to SAM_M. Finally, as Table 6 shows, the CESAM_M method still obtains much better performance for this problem than that of Power and Gmres(50) methods.


#### 5.1.4. H2/E3/1 Queues

The final structured test problem is H2/E3/1 queues from [10]. H2/E3/1 queues are used in network engineering including computer communications networks [28]. The arrival process of the H2/E3/1 queue is a two-phase hyperexponential distribution while the service process is an Erlang-3 distribution. This queueing system can be shown graphically in Figure 5. For more details, see [10]. The numerical results for this problem are presented in Tables 7 and 8.Table 7
Numerical results for H2/E3/1 queues.
n
SAM_O
SAM_M
CESAM_O
CESAM_M
C
op
it
CP
U
t
C
op
it
CP
U
t
it
CP
U
t
it
CP
U
t
6146
1.95
29
2.382
1.80
24
1.999
29
0.775
24
0.794
24578
1.95
39
12.10
1.81
23
7.368
39
3.386
23
2.915
120002
2.00
61
101.01
1.83
52
93.07
61
30.16
52
21.01Table 8
Numerical results for H2/E3/1 queues.
n
CESAM_M
Power
Gmres(50)
it
CP
U
t
it
CP
U
t
it
CP
U
t
6146
24
0.794
>106
>13.55
>5000
>18.77
24578
23
2.915
>106
>46.98
>5000
>64.08
120002
52
21.01
>106
>469.2
>5000
>201.3Figure 5
Graph for H2/E3/1 queues.As seen from Table7, the iterations needed for SAM_M are still less than those needed for SAM_O, which results in a smaller computing time of SAM_M and CESAM_M compared to SAM_O and CESAM_O, respectively. For the efficiency of the cost-effective strategy, we can see that CESAM converges in a much smaller time than SAM, with 67.5%–72.0% proportion of time saved by CESAM_O compared to SAM_O and 60.3%–60.4% proportion of time saved by CESAM_M compared to SAM_M. Finally, as Table 8 shows, the CESAM_M method still obtains much better performance for this problem than that of Power and Gmres(50) methods. For these four unstructured test problems, the SAM_M shows higher numerical scalability than SAM_O since the iterations needed almost keep the same when the sizes of the problems increase.


### 5.2. Unstructured Problems

#### 5.2.1. Random Planar Graph (Undirected)

Random walks on graphs play important roles in several fields, including information retrieval and statistical physics. For example, many ranking problems, for example, PageRank [2] and GeneRank [3] corresponding to directed and undirected graphs, respectively, can be modeled by random walks on graphs which represent relations between the items. First, we consider an unstructured planar (undirected) graph and calculate the stationary probability distribution of the random walk on the graph. We randomly distribute n nodes in (0,1)2. Then a planar graph connecting these nodes is constructed by using Delaunay triangulation and each node that shares an edge within the triangulation is connected by bidirectional links. The probability of transition from node i to node j is given by the reciprocal of the number of outward links from node i. See Figure 6(a) for a small version of this example.Figure 6
Graphs for small versions ofrandom planar graph (a) andrandom planar graph,nonsymmetric (b). Black dots represent nodes, and light gray arrows represent bidirectional links. For (b), black arrows represent unidirectional links and triangles with “+” inside have a single link that was made unidirectional. For easier visualization, the graphs shown here have a more regular distribution of points than the actual points used to build the Markov chains.
(a)
(b)Table9 shows that, for this unstructured problem, much iterations are saved using our modified Neighborhood-Based aggregation and this saving is more remarkable than that in previous examples of structured problems. Also note that the difference between the operator complexity “Cop” of SAM_O and SAM_M is larger than that in previous examples. We think the reason is that, for unstructured problems, modified Neighborhood-Based aggregation can generate more reasonable aggregates with less strong connections that are separated and less weak connections that are gathered, which leads to a much faster convergence rate and usually much more aggregates with smaller size. For this problem, 23.3%–35.3% proportion of computing time is reduced by our modified Neighborhood-Based aggregation compared to Neighborhood-Based aggregation within SAM. Similar to previous examples, there is a considerable saving in computing time by our CESAM, and this saving is more remarkable in the case of using modified Neighborhood-Based aggregation. For the size of 65536, the computing time of SAM_M is 69.7% as large as that of SAM_O, but the computing time of CESAM_M is 58.4% as large as that of CESAM_O. As shown in Table 10, when compared with Power and Gmres(50) methods, the CESAM_M still gets the best performance.Table 9
Numerical results for random planar graph (undirected).
n
SAM_O
SAM_M
CESAM_O
CESAM_M
C
op
it
CP
U
t
C
op
it
CP
U
t
it
CP
U
t
it
CP
U
t
4096
1.25
59
2.506
1.66
30
1.923
59
1.034
30
0.748
16384
1.25
81
13.08
1.72
31
8.464
81
5.190
31
3.030
65536
1.25
92
64.68
1.74
31
45.05
92
27.94
31
16.31Table 10
Numerical results for random planar graph (undirected).
n
CESAM_M
Power
Gmres(50)
it
CP
U
t
it
CP
U
t
it
CP
U
t
4096
30
0.748
14434
2.586
430
1.309
16384
31
3.030
50070
36.59
995
8.634
65536
31
16.31
>106
>545.8
2831
47.70


#### 5.2.2. Random Planar Graph (Directed)

Then, we consider directed unstructured problems. The unstructured planar graphs from the previous example are used to form a similar problem with nonsymmetric sparsity structure. Based on the graphs described in the previous example, we select a subset of triangles from the triangulation such that no two triangles in the set are neighbored. This is done by randomly selecting a triangle, marking it with “+” and marking all of its three neighbors with “−”. Then repeat this process for the next unmarked triangle until all triangles are marked. Then we randomly delete one of the six directed arcs that connect the three nodes in each “+” triangle to destroy the symmetry. This process ensures that the resulting Markov chain is still irreducible. The probability of transition from nodei to node j is given by the reciprocal of the number of outward links from node i. See Figure 6(b) for a small version of this example with the “+” triangles marked.Seen from Table11, for this nonsymmetric unstructured case, the behavior of these methods is very similar to that in previous example. Also much iterations are saved using our modified Neighborhood-Based aggregation. And the difference between the operator complexity “Cop” of SAM_O and SAM_M is still large. For this problem, 28.4%–44.2% proportion of computing time is reduced by our modified Neighborhood-Based aggregation compared to Neighborhood-Based aggregation within SAM. When CESAM is used, computing time is reduced by 51.3%–65.0% and this proportion increases with the problem size. As shown in Table 12, when the problem size is as small as 1024, the CESAM_M gets the worst performance compared with Power and Gmres(50) methods. When the size comes to 4096 and 16384, the CESAM_M outperforms the Power and Gmres(50) methods and this advantage increases with the problem size since CESAM_M is much more scalable.Table 11
Numerical results for random planar graph (directed).
n
SAM_O
SAM_M
CESAM_O
CESAM_M
C
op
it
CP
U
t
C
op
it
CP
U
t
it
CP
U
t
it
CP
U
t
1024
1.20
69
0.947
1.61
27
0.528
69
0.461
27
0.246
4096
1.24
69
2.917
1.66
29
1.872
69
1.184
29
0.742
16384
1.24
72
11.50
1.70
31
8.232
72
4.549
31
2.883Table 12
Numerical results for random planar graph (directed).
n
CESAM_M
Power
Gmres(50)
it
CP
U
t
it
CP
U
t
it
CP
U
t
1024
27
0.246
4375
0.242
173
0.139
4096
29
0.742
17030
3.077
377
1.286
16384
31
2.883
63051
45.93
1146
10.72


### 5.3. Extreme Example

The last test problem is a Markov chain matrix constructed by ourselves. It is similar to Example2 in Section 3 but we expand the size to 1024. It is a typical NCD (nearly completely decomposable) Markov chain which consists of groups of nodes that are strongly connected with each other and very weakly connected to the nodes in other groups. Here we suppose that there exists a node which is very weakly connected to all the other nodes, and the transition probabilities between it and other nodes are distributed evenly. For convenience, we constructed this matrix as(27)
0.9
1
e
-
4
1
e
-
4
0.1
1023
A
O
1
0.1
1023
O
2
B
,where(28)A=3.9e-3⋯9.9599e-1⋮⋱⋮9.9599e-1⋯3.9e-3,with size of 600×600,(29)B=1.9e-3⋯9.9422e-1⋮⋱⋮9.9422e-1⋯1.9e-3,with size of 423×423,(30)O1=4300e-3⋯4300e-3⋮⋱⋮4300e-3⋯4300e-3,with size of 600×423,(31)O2=6423e-3⋯6423e-3⋮⋱⋮6423e-3⋯6423e-3,with size of 423×600.Although this situation seems extreme, it represents the similar cases which will possibly happen in the matrices of Markov chains locally. So we hope to evaluate the rationality and the universality of the original Neighborhood-Based aggregation and our modified version through this simple test matrix. Result of numerical experiments is presented in Table5. Since we focus on the convergence rate affected by original Neighborhood-Based aggregation and modified Neighborhood-Based aggregation and there is no difference between CESAM and SAM, only the CESAM with these two aggregation methods are implemented for this example.From Table13, we can observe that the computing time and the iteration counts are reduced tremendously when our modified Neighborhood-Based aggregation is used. The performance of the original Neighborhood-Based aggregation over this problem is unacceptable which means it is not applicable for this problem, but our modified Neighborhood-Based aggregation can figure it out efficiently. The “coarsest-size” in Table 5 shows that the original Neighborhood-Based aggregation generates 1 aggregate containing all the nodes for this test matrix and our modified Neighborhood-Based aggregation generates 2 aggregates. Seen from (27), this test matrix should be divided into 2 aggregates as our modified Neighborhood-Based aggregation did. Because many nodes with weak connections are gathered together by Neighborhood-Based aggregation, the restriction and prolongation operators are very rough, leading to a poor performance as shown in Table 5. In application examples, especially in unstructured problems, the situations like this may exist locally and our modified Neighborhood-Based aggregation will be more suitable.Table 13
Numerical results of this extreme example.
n
CESAM_O
CESAM_M
C
op
Coarsest-size
it
CP
U
t
C
op
Coarsest-size
it
CP
U
t
1024
1.00
1
>200
>50.93
1.00
2
6
1.66


## 5.1. Structured Problems

### 5.1.1. Uniform 2D Lattice

The first test problem is a Markov chain on 2D lattice with uniform weights. MatrixA is essentially a scaled graph-Laplacian on a 2D uniform quadrilateral lattice with 5-point stencil (see Figure 2). It is a typically structured problem. Tables 1 and 2 give the numerical results.Table 1
Numerical results for 2D lattice with uniform weights.
n
SAM_O
SAM_M
CESAM_O
CESAM_M
C
op
it
CP
U
t
C
op
it
CP
U
t
it
CP
U
t
it
CP
U
t
4096
1.34
32
1.532
1.39
31
1.643
32
0.554
31
0.642
16384
1.33
36
6.338
1.38
31
6.121
36
1.751
31
2.094
65536
1.34
37
26.477
1.38
31
24.882
37
8.992
31
8.837Table 2
Numerical results for 2D lattice with uniform weights.
n
CESAM_M
Power
Gmres(50)
it
CP
U
t
it
CP
U
t
it
CP
U
t
4096
31
0.642
>106
>254.7
734
2.264
16384
31
2.094
>106
>33.74
2087
17.44
65536
31
8.837
>106
>254.7
>5000
>108.6Figure 2
Graph for a 2D lattice with uniform weights.Table1 shows that, for this simple and highly structured problem, there is little difference in iterations between the Neighborhood-Based aggregation and our modified version. The operator complexity “Cop” of SAM_M is a little larger than that of SAM_O for the fact that our modified Neighborhood-Based aggregation is less aggressive than the original version. This larger operator complexity and the addition work of our modified Neighborhood-Based aggregation lead to a larger computing time per cycle. For our CESAM, it has the same “Cop” (which is not shown in the table) and iterations with SAM but converges in a much smaller time. The CESAM can save the calls of much work in aggregation procedure and this is problem-dependent. For example, for the size of 4096, the CESAM_O needs 35 times of implementing the remaining work of aggregation procedure while SAM_O needs 96 times. In addition, these savings most happen in finer level, resulting in a decrease of 63.8% proportion in computing time by CESAM_O when compared to SAM_O. It is similar for other sizes and for CESAM_M compared to SAM_M. Seen from Table 2, when compared to Power and Gmres(50) method, the CESAM_M shows better performance in both iterations and computing time.


### 5.1.2. Tandem Queuing Network

The next test problem is an open tandem queuing network from [10]. Tandem queuing network exists in some practical applications, for example, in wireless networks [25] and blood screening procedures [26]. Following the settings and description in [17], we consider the situation that there are two finite queues with single servers, customers arrive according to a Poisson distribution with rate μ, and the service time distribution at the two single-server stations is Poisson with rates μ1 and μ2. In our numerical simulation, we set the number of customers in the queues as N=63,127,255. We set μ=10, μ1=11, and μ2=10. The states of this system can be represented by tuples (n1,n2), with n1 denoting the number of customers waiting in the first queue and n2 in the second queue. Then the total number of states is given by (N+1)2. It is also a structured problem and illustrated in Figure 3. The numerical results for this problem are presented in Tables 3 and 4.Table 3
Numerical results for tandem queuing network.
n
SAM_O
SAM_M
CESAM_O
CESAM_M
C
op
it
CP
U
t
C
op
it
CP
U
t
it
CP
U
t
it
CP
U
t
4096
1.24
34
1.459
1.27
34
1.576
34
0.450
34
0.442
16384
1.24
33
5.007
1.27
32
5.494
33
1.188
32
1.226
65536
1.24
43
26.51
1.25
33
22.69
43
6.529
33
5.190Table 4
Numerical results for tandem queuing network.
n
CESAM_M
Power
Gmres(50)
it
CP
U
t
it
CP
U
t
it
CP
U
t
4096
34
0.442
>106
>9.469
1652
4.631
16384
32
1.226
>106
>28.44
>5000
>40.53
65536
33
5.190
>106
>193.0
>5000
>86.82Figure 3
Graph for tandem queuing network.Seen from Table3, like the previous example, “Cop” of SAM with our modified Neighborhood-Based aggregation is still larger than that of SAM with Neighborhood-Based aggregation. This is mainly caused by the difference in the 2nd process of these two aggregation methods because these two test problems are very regular which means the strong neighborhood of each node is almost with the same size. For this problem, the number of iterations is the same for these methods when the size is 16384 or 65536. For the same reason explained above, our SAM_M costs more time per cycle than SAM_O and this leads to a larger computing time of sizes 4096 and 16384 where the iterations are almost the same. For the size of 65536, the number of iterations is reduced by a larger proportion and this offsets the larger computing time per cycle, resulting in a smaller computing time. Similar to the previous example, CESAM converges in a much smaller time than SAM, with 69.2%–76.3% proportion of time saved by CESAM_O compared to SAM_O and 72.0%–77.7% proportion of time saved by CESAM_M compared to SAM_M. Note that because the CESAM reduces the calls of much work in aggregation procedure, the higher cost of modified Neighborhood-Based aggregation has less impact. If the saving in calls is almost the same, the CESAM_M with less iterations will be more efficient than CESAM_O; this can be verified by “CPUt” of these two methods. The CESAM_M method still shows better performance than that of Power method and Gmres(50) method, which can be seen in Table 4.


### 5.1.3. M/H2/1 Queues

The next test problem is M/H2/1 queues from [10]. M/H2/1 queues are used in mobile communications [27]. M/H2/1 queue can be illustrated in Figure 4. Arrivals are generated according to a Poisson distribution at rate λ, while the service is represented by a two-phase hyperexponential distribution. With probability α a customer entering service receives service at rate μ1, while with probability 1-α this customer receives service at rate μ2. The numerical results for this problem are presented in Tables 5 and 6.Table 5
Numerical results for M/H2/1 queues.
n
SAM_O
SAM_M
CESAM_O
CESAM_M
C
op
it
CP
U
t
C
op
it
CP
U
t
it
CP
U
t
it
CP
U
t
4001
1.43
32
1.706
1.25
28
1.507
32
0.438
28
0.374
20481
1.44
33
8.042
1.25
27
6.691
33
1.825
27
1.191Table 6
Numerical results for M/H2/1 queues.
n
CESAM_M
Power
Gmres(50)
it
CP
U
t
it
CP
U
t
it
CP
U
t
4001
28
0.374
>106
>9.906
>5000
>16.73
20481
27
1.191
>106
>39.36
>5000
>52.60Figure 4
Graph for M/H2/1 queues.Table5 shows that the iterations needed for SAM_M are less than those needed for SAM_O. This means the modified Neighbour-Based aggregation can generate more suitable aggregates than the original version. Less iterations also result in a smaller computing time of SAM_M and CESAM_M compared to SAM_O and CESAM_O, respectively. For the efficiency of the cost-effective strategy, we can see that CESAM converges in a much smaller time than SAM, with 74.3%–75.2% proportion of time saved by CESAM_O compared to SAM_O and 75.2%–82.2% proportion of time saved by CESAM_M compared to SAM_M. Finally, as Table 6 shows, the CESAM_M method still obtains much better performance for this problem than that of Power and Gmres(50) methods.


### 5.1.4. H2/E3/1 Queues

The final structured test problem is H2/E3/1 queues from [10]. H2/E3/1 queues are used in network engineering including computer communications networks [28]. The arrival process of the H2/E3/1 queue is a two-phase hyperexponential distribution while the service process is an Erlang-3 distribution. This queueing system can be shown graphically in Figure 5. For more details, see [10]. The numerical results for this problem are presented in Tables 7 and 8.Table 7
Numerical results for H2/E3/1 queues.
n
SAM_O
SAM_M
CESAM_O
CESAM_M
C
op
it
CP
U
t
C
op
it
CP
U
t
it
CP
U
t
it
CP
U
t
6146
1.95
29
2.382
1.80
24
1.999
29
0.775
24
0.794
24578
1.95
39
12.10
1.81
23
7.368
39
3.386
23
2.915
120002
2.00
61
101.01
1.83
52
93.07
61
30.16
52
21.01Table 8
Numerical results for H2/E3/1 queues.
n
CESAM_M
Power
Gmres(50)
it
CP
U
t
it
CP
U
t
it
CP
U
t
6146
24
0.794
>106
>13.55
>5000
>18.77
24578
23
2.915
>106
>46.98
>5000
>64.08
120002
52
21.01
>106
>469.2
>5000
>201.3Figure 5
Graph for H2/E3/1 queues.As seen from Table7, the iterations needed for SAM_M are still less than those needed for SAM_O, which results in a smaller computing time of SAM_M and CESAM_M compared to SAM_O and CESAM_O, respectively. For the efficiency of the cost-effective strategy, we can see that CESAM converges in a much smaller time than SAM, with 67.5%–72.0% proportion of time saved by CESAM_O compared to SAM_O and 60.3%–60.4% proportion of time saved by CESAM_M compared to SAM_M. Finally, as Table 8 shows, the CESAM_M method still obtains much better performance for this problem than that of Power and Gmres(50) methods. For these four unstructured test problems, the SAM_M shows higher numerical scalability than SAM_O since the iterations needed almost keep the same when the sizes of the problems increase.


## 5.1.1. Uniform 2D Lattice

The first test problem is a Markov chain on 2D lattice with uniform weights. MatrixA is essentially a scaled graph-Laplacian on a 2D uniform quadrilateral lattice with 5-point stencil (see Figure 2). It is a typically structured problem. Tables 1 and 2 give the numerical results.Table 1
Numerical results for 2D lattice with uniform weights.
n
SAM_O
SAM_M
CESAM_O
CESAM_M
C
op
it
CP
U
t
C
op
it
CP
U
t
it
CP
U
t
it
CP
U
t
4096
1.34
32
1.532
1.39
31
1.643
32
0.554
31
0.642
16384
1.33
36
6.338
1.38
31
6.121
36
1.751
31
2.094
65536
1.34
37
26.477
1.38
31
24.882
37
8.992
31
8.837Table 2
Numerical results for 2D lattice with uniform weights.
n
CESAM_M
Power
Gmres(50)
it
CP
U
t
it
CP
U
t
it
CP
U
t
4096
31
0.642
>106
>254.7
734
2.264
16384
31
2.094
>106
>33.74
2087
17.44
65536
31
8.837
>106
>254.7
>5000
>108.6Figure 2
Graph for a 2D lattice with uniform weights.Table1 shows that, for this simple and highly structured problem, there is little difference in iterations between the Neighborhood-Based aggregation and our modified version. The operator complexity “Cop” of SAM_M is a little larger than that of SAM_O for the fact that our modified Neighborhood-Based aggregation is less aggressive than the original version. This larger operator complexity and the addition work of our modified Neighborhood-Based aggregation lead to a larger computing time per cycle. For our CESAM, it has the same “Cop” (which is not shown in the table) and iterations with SAM but converges in a much smaller time. The CESAM can save the calls of much work in aggregation procedure and this is problem-dependent. For example, for the size of 4096, the CESAM_O needs 35 times of implementing the remaining work of aggregation procedure while SAM_O needs 96 times. In addition, these savings most happen in finer level, resulting in a decrease of 63.8% proportion in computing time by CESAM_O when compared to SAM_O. It is similar for other sizes and for CESAM_M compared to SAM_M. Seen from Table 2, when compared to Power and Gmres(50) method, the CESAM_M shows better performance in both iterations and computing time.


## 5.1.2. Tandem Queuing Network

The next test problem is an open tandem queuing network from [10]. Tandem queuing network exists in some practical applications, for example, in wireless networks [25] and blood screening procedures [26]. Following the settings and description in [17], we consider the situation that there are two finite queues with single servers, customers arrive according to a Poisson distribution with rate μ, and the service time distribution at the two single-server stations is Poisson with rates μ1 and μ2. In our numerical simulation, we set the number of customers in the queues as N=63,127,255. We set μ=10, μ1=11, and μ2=10. The states of this system can be represented by tuples (n1,n2), with n1 denoting the number of customers waiting in the first queue and n2 in the second queue. Then the total number of states is given by (N+1)2. It is also a structured problem and illustrated in Figure 3. The numerical results for this problem are presented in Tables 3 and 4.Table 3
Numerical results for tandem queuing network.
n
SAM_O
SAM_M
CESAM_O
CESAM_M
C
op
it
CP
U
t
C
op
it
CP
U
t
it
CP
U
t
it
CP
U
t
4096
1.24
34
1.459
1.27
34
1.576
34
0.450
34
0.442
16384
1.24
33
5.007
1.27
32
5.494
33
1.188
32
1.226
65536
1.24
43
26.51
1.25
33
22.69
43
6.529
33
5.190Table 4
Numerical results for tandem queuing network.
n
CESAM_M
Power
Gmres(50)
it
CP
U
t
it
CP
U
t
it
CP
U
t
4096
34
0.442
>106
>9.469
1652
4.631
16384
32
1.226
>106
>28.44
>5000
>40.53
65536
33
5.190
>106
>193.0
>5000
>86.82Figure 3
Graph for tandem queuing network.Seen from Table3, like the previous example, “Cop” of SAM with our modified Neighborhood-Based aggregation is still larger than that of SAM with Neighborhood-Based aggregation. This is mainly caused by the difference in the 2nd process of these two aggregation methods because these two test problems are very regular which means the strong neighborhood of each node is almost with the same size. For this problem, the number of iterations is the same for these methods when the size is 16384 or 65536. For the same reason explained above, our SAM_M costs more time per cycle than SAM_O and this leads to a larger computing time of sizes 4096 and 16384 where the iterations are almost the same. For the size of 65536, the number of iterations is reduced by a larger proportion and this offsets the larger computing time per cycle, resulting in a smaller computing time. Similar to the previous example, CESAM converges in a much smaller time than SAM, with 69.2%–76.3% proportion of time saved by CESAM_O compared to SAM_O and 72.0%–77.7% proportion of time saved by CESAM_M compared to SAM_M. Note that because the CESAM reduces the calls of much work in aggregation procedure, the higher cost of modified Neighborhood-Based aggregation has less impact. If the saving in calls is almost the same, the CESAM_M with less iterations will be more efficient than CESAM_O; this can be verified by “CPUt” of these two methods. The CESAM_M method still shows better performance than that of Power method and Gmres(50) method, which can be seen in Table 4.


## 5.1.3. M/H2/1 Queues

The next test problem is M/H2/1 queues from [10]. M/H2/1 queues are used in mobile communications [27]. M/H2/1 queue can be illustrated in Figure 4. Arrivals are generated according to a Poisson distribution at rate λ, while the service is represented by a two-phase hyperexponential distribution. With probability α a customer entering service receives service at rate μ1, while with probability 1-α this customer receives service at rate μ2. The numerical results for this problem are presented in Tables 5 and 6.Table 5
Numerical results for M/H2/1 queues.
n
SAM_O
SAM_M
CESAM_O
CESAM_M
C
op
it
CP
U
t
C
op
it
CP
U
t
it
CP
U
t
it
CP
U
t
4001
1.43
32
1.706
1.25
28
1.507
32
0.438
28
0.374
20481
1.44
33
8.042
1.25
27
6.691
33
1.825
27
1.191Table 6
Numerical results for M/H2/1 queues.
n
CESAM_M
Power
Gmres(50)
it
CP
U
t
it
CP
U
t
it
CP
U
t
4001
28
0.374
>106
>9.906
>5000
>16.73
20481
27
1.191
>106
>39.36
>5000
>52.60Figure 4
Graph for M/H2/1 queues.Table5 shows that the iterations needed for SAM_M are less than those needed for SAM_O. This means the modified Neighbour-Based aggregation can generate more suitable aggregates than the original version. Less iterations also result in a smaller computing time of SAM_M and CESAM_M compared to SAM_O and CESAM_O, respectively. For the efficiency of the cost-effective strategy, we can see that CESAM converges in a much smaller time than SAM, with 74.3%–75.2% proportion of time saved by CESAM_O compared to SAM_O and 75.2%–82.2% proportion of time saved by CESAM_M compared to SAM_M. Finally, as Table 6 shows, the CESAM_M method still obtains much better performance for this problem than that of Power and Gmres(50) methods.


## 5.1.4. H2/E3/1 Queues

The final structured test problem is H2/E3/1 queues from [10]. H2/E3/1 queues are used in network engineering including computer communications networks [28]. The arrival process of the H2/E3/1 queue is a two-phase hyperexponential distribution while the service process is an Erlang-3 distribution. This queueing system can be shown graphically in Figure 5. For more details, see [10]. The numerical results for this problem are presented in Tables 7 and 8.Table 7
Numerical results for H2/E3/1 queues.
n
SAM_O
SAM_M
CESAM_O
CESAM_M
C
op
it
CP
U
t
C
op
it
CP
U
t
it
CP
U
t
it
CP
U
t
6146
1.95
29
2.382
1.80
24
1.999
29
0.775
24
0.794
24578
1.95
39
12.10
1.81
23
7.368
39
3.386
23
2.915
120002
2.00
61
101.01
1.83
52
93.07
61
30.16
52
21.01Table 8
Numerical results for H2/E3/1 queues.
n
CESAM_M
Power
Gmres(50)
it
CP
U
t
it
CP
U
t
it
CP
U
t
6146
24
0.794
>106
>13.55
>5000
>18.77
24578
23
2.915
>106
>46.98
>5000
>64.08
120002
52
21.01
>106
>469.2
>5000
>201.3Figure 5
Graph for H2/E3/1 queues.As seen from Table7, the iterations needed for SAM_M are still less than those needed for SAM_O, which results in a smaller computing time of SAM_M and CESAM_M compared to SAM_O and CESAM_O, respectively. For the efficiency of the cost-effective strategy, we can see that CESAM converges in a much smaller time than SAM, with 67.5%–72.0% proportion of time saved by CESAM_O compared to SAM_O and 60.3%–60.4% proportion of time saved by CESAM_M compared to SAM_M. Finally, as Table 8 shows, the CESAM_M method still obtains much better performance for this problem than that of Power and Gmres(50) methods. For these four unstructured test problems, the SAM_M shows higher numerical scalability than SAM_O since the iterations needed almost keep the same when the sizes of the problems increase.


## 5.2. Unstructured Problems

### 5.2.1. Random Planar Graph (Undirected)

Random walks on graphs play important roles in several fields, including information retrieval and statistical physics. For example, many ranking problems, for example, PageRank [2] and GeneRank [3] corresponding to directed and undirected graphs, respectively, can be modeled by random walks on graphs which represent relations between the items. First, we consider an unstructured planar (undirected) graph and calculate the stationary probability distribution of the random walk on the graph. We randomly distribute n nodes in (0,1)2. Then a planar graph connecting these nodes is constructed by using Delaunay triangulation and each node that shares an edge within the triangulation is connected by bidirectional links. The probability of transition from node i to node j is given by the reciprocal of the number of outward links from node i. See Figure 6(a) for a small version of this example.Figure 6
Graphs for small versions ofrandom planar graph (a) andrandom planar graph,nonsymmetric (b). Black dots represent nodes, and light gray arrows represent bidirectional links. For (b), black arrows represent unidirectional links and triangles with “+” inside have a single link that was made unidirectional. For easier visualization, the graphs shown here have a more regular distribution of points than the actual points used to build the Markov chains.
(a)
(b)Table9 shows that, for this unstructured problem, much iterations are saved using our modified Neighborhood-Based aggregation and this saving is more remarkable than that in previous examples of structured problems. Also note that the difference between the operator complexity “Cop” of SAM_O and SAM_M is larger than that in previous examples. We think the reason is that, for unstructured problems, modified Neighborhood-Based aggregation can generate more reasonable aggregates with less strong connections that are separated and less weak connections that are gathered, which leads to a much faster convergence rate and usually much more aggregates with smaller size. For this problem, 23.3%–35.3% proportion of computing time is reduced by our modified Neighborhood-Based aggregation compared to Neighborhood-Based aggregation within SAM. Similar to previous examples, there is a considerable saving in computing time by our CESAM, and this saving is more remarkable in the case of using modified Neighborhood-Based aggregation. For the size of 65536, the computing time of SAM_M is 69.7% as large as that of SAM_O, but the computing time of CESAM_M is 58.4% as large as that of CESAM_O. As shown in Table 10, when compared with Power and Gmres(50) methods, the CESAM_M still gets the best performance.Table 9
Numerical results for random planar graph (undirected).
n
SAM_O
SAM_M
CESAM_O
CESAM_M
C
op
it
CP
U
t
C
op
it
CP
U
t
it
CP
U
t
it
CP
U
t
4096
1.25
59
2.506
1.66
30
1.923
59
1.034
30
0.748
16384
1.25
81
13.08
1.72
31
8.464
81
5.190
31
3.030
65536
1.25
92
64.68
1.74
31
45.05
92
27.94
31
16.31Table 10
Numerical results for random planar graph (undirected).
n
CESAM_M
Power
Gmres(50)
it
CP
U
t
it
CP
U
t
it
CP
U
t
4096
30
0.748
14434
2.586
430
1.309
16384
31
3.030
50070
36.59
995
8.634
65536
31
16.31
>106
>545.8
2831
47.70


### 5.2.2. Random Planar Graph (Directed)

Then, we consider directed unstructured problems. The unstructured planar graphs from the previous example are used to form a similar problem with nonsymmetric sparsity structure. Based on the graphs described in the previous example, we select a subset of triangles from the triangulation such that no two triangles in the set are neighbored. This is done by randomly selecting a triangle, marking it with “+” and marking all of its three neighbors with “−”. Then repeat this process for the next unmarked triangle until all triangles are marked. Then we randomly delete one of the six directed arcs that connect the three nodes in each “+” triangle to destroy the symmetry. This process ensures that the resulting Markov chain is still irreducible. The probability of transition from nodei to node j is given by the reciprocal of the number of outward links from node i. See Figure 6(b) for a small version of this example with the “+” triangles marked.Seen from Table11, for this nonsymmetric unstructured case, the behavior of these methods is very similar to that in previous example. Also much iterations are saved using our modified Neighborhood-Based aggregation. And the difference between the operator complexity “Cop” of SAM_O and SAM_M is still large. For this problem, 28.4%–44.2% proportion of computing time is reduced by our modified Neighborhood-Based aggregation compared to Neighborhood-Based aggregation within SAM. When CESAM is used, computing time is reduced by 51.3%–65.0% and this proportion increases with the problem size. As shown in Table 12, when the problem size is as small as 1024, the CESAM_M gets the worst performance compared with Power and Gmres(50) methods. When the size comes to 4096 and 16384, the CESAM_M outperforms the Power and Gmres(50) methods and this advantage increases with the problem size since CESAM_M is much more scalable.Table 11
Numerical results for random planar graph (directed).
n
SAM_O
SAM_M
CESAM_O
CESAM_M
C
op
it
CP
U
t
C
op
it
CP
U
t
it
CP
U
t
it
CP
U
t
1024
1.20
69
0.947
1.61
27
0.528
69
0.461
27
0.246
4096
1.24
69
2.917
1.66
29
1.872
69
1.184
29
0.742
16384
1.24
72
11.50
1.70
31
8.232
72
4.549
31
2.883Table 12
Numerical results for random planar graph (directed).
n
CESAM_M
Power
Gmres(50)
it
CP
U
t
it
CP
U
t
it
CP
U
t
1024
27
0.246
4375
0.242
173
0.139
4096
29
0.742
17030
3.077
377
1.286
16384
31
2.883
63051
45.93
1146
10.72


## 5.2.1. Random Planar Graph (Undirected)

Random walks on graphs play important roles in several fields, including information retrieval and statistical physics. For example, many ranking problems, for example, PageRank [2] and GeneRank [3] corresponding to directed and undirected graphs, respectively, can be modeled by random walks on graphs which represent relations between the items. First, we consider an unstructured planar (undirected) graph and calculate the stationary probability distribution of the random walk on the graph. We randomly distribute n nodes in (0,1)2. Then a planar graph connecting these nodes is constructed by using Delaunay triangulation and each node that shares an edge within the triangulation is connected by bidirectional links. The probability of transition from node i to node j is given by the reciprocal of the number of outward links from node i. See Figure 6(a) for a small version of this example.Figure 6
Graphs for small versions ofrandom planar graph (a) andrandom planar graph,nonsymmetric (b). Black dots represent nodes, and light gray arrows represent bidirectional links. For (b), black arrows represent unidirectional links and triangles with “+” inside have a single link that was made unidirectional. For easier visualization, the graphs shown here have a more regular distribution of points than the actual points used to build the Markov chains.
(a)
(b)Table9 shows that, for this unstructured problem, much iterations are saved using our modified Neighborhood-Based aggregation and this saving is more remarkable than that in previous examples of structured problems. Also note that the difference between the operator complexity “Cop” of SAM_O and SAM_M is larger than that in previous examples. We think the reason is that, for unstructured problems, modified Neighborhood-Based aggregation can generate more reasonable aggregates with less strong connections that are separated and less weak connections that are gathered, which leads to a much faster convergence rate and usually much more aggregates with smaller size. For this problem, 23.3%–35.3% proportion of computing time is reduced by our modified Neighborhood-Based aggregation compared to Neighborhood-Based aggregation within SAM. Similar to previous examples, there is a considerable saving in computing time by our CESAM, and this saving is more remarkable in the case of using modified Neighborhood-Based aggregation. For the size of 65536, the computing time of SAM_M is 69.7% as large as that of SAM_O, but the computing time of CESAM_M is 58.4% as large as that of CESAM_O. As shown in Table 10, when compared with Power and Gmres(50) methods, the CESAM_M still gets the best performance.Table 9
Numerical results for random planar graph (undirected).
n
SAM_O
SAM_M
CESAM_O
CESAM_M
C
op
it
CP
U
t
C
op
it
CP
U
t
it
CP
U
t
it
CP
U
t
4096
1.25
59
2.506
1.66
30
1.923
59
1.034
30
0.748
16384
1.25
81
13.08
1.72
31
8.464
81
5.190
31
3.030
65536
1.25
92
64.68
1.74
31
45.05
92
27.94
31
16.31Table 10
Numerical results for random planar graph (undirected).
n
CESAM_M
Power
Gmres(50)
it
CP
U
t
it
CP
U
t
it
CP
U
t
4096
30
0.748
14434
2.586
430
1.309
16384
31
3.030
50070
36.59
995
8.634
65536
31
16.31
>106
>545.8
2831
47.70


## 5.2.2. Random Planar Graph (Directed)

Then, we consider directed unstructured problems. The unstructured planar graphs from the previous example are used to form a similar problem with nonsymmetric sparsity structure. Based on the graphs described in the previous example, we select a subset of triangles from the triangulation such that no two triangles in the set are neighbored. This is done by randomly selecting a triangle, marking it with “+” and marking all of its three neighbors with “−”. Then repeat this process for the next unmarked triangle until all triangles are marked. Then we randomly delete one of the six directed arcs that connect the three nodes in each “+” triangle to destroy the symmetry. This process ensures that the resulting Markov chain is still irreducible. The probability of transition from nodei to node j is given by the reciprocal of the number of outward links from node i. See Figure 6(b) for a small version of this example with the “+” triangles marked.Seen from Table11, for this nonsymmetric unstructured case, the behavior of these methods is very similar to that in previous example. Also much iterations are saved using our modified Neighborhood-Based aggregation. And the difference between the operator complexity “Cop” of SAM_O and SAM_M is still large. For this problem, 28.4%–44.2% proportion of computing time is reduced by our modified Neighborhood-Based aggregation compared to Neighborhood-Based aggregation within SAM. When CESAM is used, computing time is reduced by 51.3%–65.0% and this proportion increases with the problem size. As shown in Table 12, when the problem size is as small as 1024, the CESAM_M gets the worst performance compared with Power and Gmres(50) methods. When the size comes to 4096 and 16384, the CESAM_M outperforms the Power and Gmres(50) methods and this advantage increases with the problem size since CESAM_M is much more scalable.Table 11
Numerical results for random planar graph (directed).
n
SAM_O
SAM_M
CESAM_O
CESAM_M
C
op
it
CP
U
t
C
op
it
CP
U
t
it
CP
U
t
it
CP
U
t
1024
1.20
69
0.947
1.61
27
0.528
69
0.461
27
0.246
4096
1.24
69
2.917
1.66
29
1.872
69
1.184
29
0.742
16384
1.24
72
11.50
1.70
31
8.232
72
4.549
31
2.883Table 12
Numerical results for random planar graph (directed).
n
CESAM_M
Power
Gmres(50)
it
CP
U
t
it
CP
U
t
it
CP
U
t
1024
27
0.246
4375
0.242
173
0.139
4096
29
0.742
17030
3.077
377
1.286
16384
31
2.883
63051
45.93
1146
10.72


## 5.3. Extreme Example

The last test problem is a Markov chain matrix constructed by ourselves. It is similar to Example2 in Section 3 but we expand the size to 1024. It is a typical NCD (nearly completely decomposable) Markov chain which consists of groups of nodes that are strongly connected with each other and very weakly connected to the nodes in other groups. Here we suppose that there exists a node which is very weakly connected to all the other nodes, and the transition probabilities between it and other nodes are distributed evenly. For convenience, we constructed this matrix as(27)
0.9
1
e
-
4
1
e
-
4
0.1
1023
A
O
1
0.1
1023
O
2
B
,where(28)A=3.9e-3⋯9.9599e-1⋮⋱⋮9.9599e-1⋯3.9e-3,with size of 600×600,(29)B=1.9e-3⋯9.9422e-1⋮⋱⋮9.9422e-1⋯1.9e-3,with size of 423×423,(30)O1=4300e-3⋯4300e-3⋮⋱⋮4300e-3⋯4300e-3,with size of 600×423,(31)O2=6423e-3⋯6423e-3⋮⋱⋮6423e-3⋯6423e-3,with size of 423×600.Although this situation seems extreme, it represents the similar cases which will possibly happen in the matrices of Markov chains locally. So we hope to evaluate the rationality and the universality of the original Neighborhood-Based aggregation and our modified version through this simple test matrix. Result of numerical experiments is presented in Table5. Since we focus on the convergence rate affected by original Neighborhood-Based aggregation and modified Neighborhood-Based aggregation and there is no difference between CESAM and SAM, only the CESAM with these two aggregation methods are implemented for this example.From Table13, we can observe that the computing time and the iteration counts are reduced tremendously when our modified Neighborhood-Based aggregation is used. The performance of the original Neighborhood-Based aggregation over this problem is unacceptable which means it is not applicable for this problem, but our modified Neighborhood-Based aggregation can figure it out efficiently. The “coarsest-size” in Table 5 shows that the original Neighborhood-Based aggregation generates 1 aggregate containing all the nodes for this test matrix and our modified Neighborhood-Based aggregation generates 2 aggregates. Seen from (27), this test matrix should be divided into 2 aggregates as our modified Neighborhood-Based aggregation did. Because many nodes with weak connections are gathered together by Neighborhood-Based aggregation, the restriction and prolongation operators are very rough, leading to a poor performance as shown in Table 5. In application examples, especially in unstructured problems, the situations like this may exist locally and our modified Neighborhood-Based aggregation will be more suitable.Table 13
Numerical results of this extreme example.
n
CESAM_O
CESAM_M
C
op
Coarsest-size
it
CP
U
t
C
op
Coarsest-size
it
CP
U
t
1024
1.00
1
>200
>50.93
1.00
2
6
1.66


## 6. Conclusions

In this paper, we have studied the Neighborhood-Based aggregation and explained the shortage and the irrationality of it. Then a modified Neighborhood-Based aggregation is presented. We also propose a CESAM method which reduces the calls of much work in aggregation procedure with the same convergence rate of SAM. The numerical results confirm that the modified Neighborhood-Based aggregation method is much more efficient than Neighborhood-Based aggregation for unstructured problems and with stronger applicability, and a significant saving in computing time is achieved by CESAM method. The CESAM_M method also shows good potential when compared with other methods for Markov chain problems.



---

*Source: 816935-2015-10-29.xml*