Algebraic multi-grid method for disordered systems

Nuclear Physics B (Proc. Suppl.) 9 (1989) 521-524 North-Holland, Amsterdam

521

Algebraic Multi-Grid Method for Disordered Systems Robert G. Edwards [email protected]

J o n a t h a n Goodman [email protected]

Alan D. Sokal [email protected]~fU.EDU

Departments of Physics and Mathematics New York University

The computation of the fermion (quark) propagator in a random gauge field is a key step in Monte Carlo studies of hadron masses in quantum chromodynamics [1] and in nearly allalgorithms for dynamical fermions [2]. For the latticesizes used in present-day hadron-mass studies, each propagator computation takes several hours of supercomputer time [3]. It is thus essentialto seek new and more efficientalgorithms for this problem. From a mathematical point of view, the computation of the fermion propagator is a deterministic problem of numerical linear algebra: the solution of a large, sparse system of linear equations arising from the discretizationof an elliptic differentialoperator (the Euclidean Dirac operator in a background gauge field). This problem is thus part of a wen-establlshed field of numerical analysis. But the equations studied by numerical analysts usually have coefficients which are smooth on the characteristiclength scale of the problem, e.g.

physics as a simplifiedmodel of electricalconduction in composite materials [5,6]. Consider a d-dimensional lattice with L -4-1 sites in the vertical direction (free boundaries) and L sitesin each of the d - 1 horizontal directions (periodic boundary conditions). Let each nearest-neighbor bond be a 1-ohm resistor with probability p, or an infinite resistor with probability 1 - p (independently on each bond). N o w attach conducting plates to the top and bottom faces, connect a 1-volt battery, and measure the current flow (i.e.the conductance GL). Then z =

defines an "effective conductivity" of the random medium; we would like to know how Z depends on p. Clearly ~ vanishes whenever p is less than the percolation threshold Pc; and it can be proven that ~ is nonvanishing for p > pc [7]. It is thus reasonable to expect that

S ~ (p-pc) t

--A~o+(l+sinzl)~O = f

or have at worst jump discontinuities,e.g. -A~o + (1 + sgn¢l)~p = f . In QCD, by contrast, the coefficients are the (random) values of the gauge field, which in general are VerY rough on the characteristic length scale m~ 1 or m~,1. This is because quantum fields are very singular objects: in the continuum limit they are distributions in the sense of Schwartz [4]. We emphasize that this would be the case even if gauge-fixing were employed: the point is that the curvature of the gauge field over one fermion correlation length is far from small. On the other hand, ff no gauge-fixing is imposed, then the gauge field will vary rapidly from one link to the next even ff the curvature is zero. The fermion-propagator problem is thus a problem on the frontiers of numerical analysis: a linear partial differential equation with wildly discontinuous coefficients. As a first step, it would be valuable to have a "warm-up problem" in which the key issue of disordered coefficients can be dealt with separately from the technical complications caused by gauge invariance and spin. One such problem is that of computing currents in a random resistor network near percolation threshold, which is of interest in solld-state

0920-5632/89/$03.50 Q Elsevier Science Publishers B.V. (North-tIolland Physics Publishing Division)

~ m
aspJ, pc

for some criticalexponent t which we would like to compute. One obvious approach to this problem is direct Monte Carlo simulation combined with finite-size scaling. The Monte Carlo aspect of this problem is trivial (independent bond percolation); the hard part is computing the conductance for a given resistor configuration. This boil~ down to solving Kirchhoff's equations for the current flow, or equivalently a Laplace-like equation for the node voltages -- a large, sparse, and highly ill-conditioned system of linear equations. Direct methods of solution, such as sophisticated variants of Gauss elimination, have been tried with modest success [6,8],but the computational labor grows as L 3d-a when p = 1 (albeit more slowly when p ~ pc) [9]. Iterative methods of solution require a labor of only order L d per iteration (when p = 1), but typically suffer from critical slowing-down, so that 'the number of iterations required grows as a power of L. In fact, the criticalslowing-down is much more severe in disordered problems (p ~ Pc) than in regular problems (p = 1).

In recent years m u c h efforthas been devoted to devising improved iterative algorithms with reduced criticalslowingdown (see Table 1). In particular, Fourier preconditioning [19,20] has been shown to reduce appreciably the critical

R.G. Edwards et al./ Algebraic multi-grid method

522

Algorithm

r (p = Pc cluster) ~" (P = Pc backbone)

r (p = 1)

Jacohi, Gauss-Seidel Optimal S O R Conjugate Gradient (CG) Fourier-Preconditioned Jacobi Fourier-Preconditioned C G Algebraic Multi-Grid ( A M G )

L ~z'sz [8,10,11] L ~1'44 [8,10,11]* L ~1"44 [8,10,11]*

L 2 [13] L 1 [13] L 1 [14] L° L° L ° [lS]

L ~'1"°

L°

[16] [17]

L~'s L ~l"s L ~t'3 L ~1"4 L ~°'7 L°

[8,10,11,12] [8,10,11,12]* [8,10,11,12]* [15] [15]* [17]

Table 1: Critical slowing-down in selected iterative algorithms for the two-dimensional random resistor problem. Asterisk indicates t h a t critical exponent is inferred from cited article together with standard SOR/CG theory [13,14].

slowing-down in the random-resistor problem [15,16], but severe critical slowing-down still rer-a~us. This is because the preconditioning operator - - a fractional power of the Laplacian - - mim~cs the ensemble-averaged current flow in the lattice, but the current flow in any particular realization of the random resistor network is affected strongly by the local and global topology of interconnections. This reasoning suggests that an improved strategy should take account of the topology of the particular resistor network at hand. The multi-grid method [18,21]is known to be an extraordinarily effective approach for solving large linear systems arising from the discretizationof ellipticpartial differential equations (PDEs). Usually the coarse grids and interpolation operators are defined geometrically, e.g. cubical (2 × 2) blocks with piecewise-constant or piecewise-linearinterpolation. This approach, which is suitable for P D E s with smooth coefficients,will clearly not be appropriate in disordered systems such as the random-resistor network: just because two sites are close geometrically does not mean that they are close in the topology of the resistor network and hence in voltage. Rather, a successfulmulti-grid algorithm will have to define its coarse grids and interpolation operators in accordance with the connection structure of the particular resistor network. Algebraic multi-grid ( A M G ) [22] is just such a strategy: the A M G codes are are "black-box" solvers which, given an arbitrary matrix (of a certain type), attempt to choose appropriate coarse grids and interpolation operators. Of course, there is no guarantee as to how well the algorithm will perform on any given problem. Heretofore A M G has been applied to P D E s with strong jump discontinuities,but not to problems as singular as the random-resistor network. W e have tried a standard A M G code ( A M G I R 4 [23])on the two-dimensional random-resistor problem, and have found (somewhat to our surprise) that it succeeds in eliminating entirely the criticalslowing-down. Our computer experiments are as follows [17]: Using an L × (L + 1) latticeas explained above, we set the electricpotential to be ~ = 0 on the bottom row, ~ = 1 on the top row, and compute ~ on the interior sites using Kirchhoff's and Ohm's laws. T w o prel{miuary reductions are applied before solving the resulting system of linear equations. First, we

L 100 200 400

Cluster 0.348 (0.002) 0.403 (0.001)

Backbone 0.264 (0.002) 0.319 (0.001) 0.362 (0.001)

Table 2: Mean measured convergence factors for algebraic multigrid (AMG) algorithm on two-dimensional random resistor problem at percolation threshold, computed on either connected cluster or current-carrying backbone. Standard error is shown in parentheses.

compute the "connected cluster" using the Fischer-GallerHoshen-Kopelman algorithm [24]. A site is in the connected cluster if there is a path of bonds from it to the top and to the bottom. Next we compute the "blconnected cluster" or "current-carrying backbone". A bond is not part of the backbone if for topological reasons it could not possibly carry current, i.e.if it is part of a "dangling end". Tarjan's depth-first-searchalgorithm [251 finds the backbone in time proportional to the number of sitesin the connected cluster. The C P U time per iteration of the A M G I R 4 code is also proportional to the number of unknowns. W e present results for L = i00,200,400 at p = Pc = 1/2. W e started with an initial guess of ~ -- 0 in the interior of the lattice and iterated until the error ][e][ [~,i~e, (current loss)2] 1/2 < 10 -12. This corresponds to a reduction in ]lell by about 13 orders of maguitude. Such high accuracy may be unnecessary for physical applications but it gives a better understanding of the numerical method. The measured convergence .factor is then defined to be the factor by which the error is reduced in the last multigrid iteration. The true (asymptotic) convergence factor is slightly higher titan this, but the practically relevant (average) convergence factor is significantly lower. The mean (over realizations) values of the measured convergence factor are given in Table 2. Thus, for the 400 × 400 lattice we can obtain the solution to six-decimal-place accuracy in about 10 iterations of the A M G method, while it would take roughly 10 million iterations using the Gauss-Seidel method! Although the convergence factors increase slightly with L, there is no evidence that they are approaching 1 as L --~ co. W e conclude that criticalslowing-down is completely (or

R.G. Edwards et M. / Algebraic multi-grid method

almost completely) absent. Of course, this extraordinary performance does not come for free: setting up the coarse grids and interpolation operators takes a significant CPU time, and both memory usage and CPU time per iteration are high due to the overhead associated with linked lists. 1 Nevertheless, we feel that at this stage it is most urgent to concentrate on the physics of critical slowing-down. Later, after a good algorithm has been identified, one can work on making the implementation more efficient, taking advantage of the structure of a particular problem. For lack of space, we can give here only a very rough sketch of how the AMG codes choose the coarse grids and interpolation operators; for more details, see [22]. The existing AMG algorithms are designed primarily for diagonally dominant symmetric M-matrices, i.e. real symmetric matrices that are negative off-diagonal and have non-negative row sums. Such matrices are essentially"Laplacian-like": a typical example is the SchxSdinger operator - A + V on an arbitrary fiuitegraph with an arbitrary potential V > 0. The first step in the A M G setup phase is to divide the latticepoints into coarse-gridpoints (C-points)and allothers (F-points), so that every F-point is "strongly connected" (by the given matrix A) to at least one C-point, using as few Cpoints as possible consistent with this constraint. There are heuristic algorithms to do this [22]. The next step is to define the interpolation weights {w~k}ieFuc, k~C: here wlh is the strength of interpolation from a coarse-grid point k to a fine-grid point i. If i is a Cpoint, we make the obvious choice wlk : 61k. On the other hand, if i is an F-point, the idea is that wik should be roughly proportional to the "strength of connection" between i and k in the original matrix A. The simplest implementation of this idea is interpolation along direct connections: Let Ni (the "neighborhood" of i) be the set of all sites j ~ i such that aij ~ O; and let Ci be some chosen subset of Ni N C. Then we define --aik wlk - ~ a~j j~v~ for k E Ci, and zero otherwise. Finally, we define the coarsegrid matrix by the variational (Galerkin) formula: a ~ "'~ = y~w~a~jv~jl %J

This procedure (choice of C-points, w and Ac°ar'e) is then repeated to define the sequence of coarse grids. For the A M G iterationto work well, it is necessary that the "interpolatory connections" aij (j E Ci) comprise a reasonable fraction of i's total connection strength, for each i E F; but one wants to choose C~ as small as possible consistent with this constraint,in order to keep the coarse-grid a T h e t o t a l t i m e needed to solve the L = 200 l i n e a r equations on

the cluster (~esp. backbone) was roughly 5.3 (~esp. 1.5) minutes per configuration on a Sun 3/160 workstation with floating-polnt accelerator and 16 Mbytes of memory. Memory usage is approximately 60 floatingpoint words per cluster (or backbone) site•

523

matrix as sparse as possible. It turns out that the best overall performance is obtained by using interpolationsmore sophisticated than the one explained above, such as interpolation along second-order connections (i.e.along matrix dements of A 2 as well as of A). Such an interpolation is used in the AMG1R4 code. Unfortunately, the foregoing theory does not extend immediately to complez matrices such as the lattice Laplace or Dirac operator in a background gauge field. Indeed, the failure of the symmetric M-matrix condition is intimately related to the curvature of the gauge field. We are currently working on developing gauge-covariant extensions of AMG ideas. W e remark, finally,that a stochastic generalization [26] of the A M G algorithm can be used for updating the scalar sector in the standard model of dynamically triangulated random surfaces (= discretizedbosonic Polyakov string) [27]. This work was supported in part by N S F grants PHY8413569 and DMS-8705599. One of the authors (J.G.) was partially supported by a Sloan Foundation Fellowship and by D A R P A .

References [1] I. Montvay, Rev. Mod. Phys. 59, 263 (1987). [2] D. Weingarten, review talk at this conference. [3] S.R. Sharpe, J. Stat. Phys. 43, 1129 (1986), see p. 1133. [4] N. Bohx and L. Rosenfeld, Dansk Vid. Selsk. Mat. Fys. Medd. 12, no. 8 (1933); A.S. Wightman, Ann. Inst. Henri Poincar~ A1, 403 (1964); Z. Wizimirski, Bull. Acad. Polon. Sci. Set. Math. Astr. Phys. 14, 91 (1966); P. Colella and O.E. Lanford, in ConstructiveQuantum Field Theory (Lecture Notes in Physics #25), ed. G. Velo and A. Wightman (Springer, Berlin, 1973). [5] S. Kirkpatrick, Rev. Mod. Phys. 45, 574 (1973). [6] S. Kirkpatrick, in Ill-Condensed Matter (Les Houches 1978), ed. R. Balian et al. (North-Holland, Amsterdam,

1979). [7] J.T. Chayesand L. Chayes,Commun.Math. Phys. 1{}5, 133 (1986). [8] D.J. Frank and C.J. Lobb, Phys. Rev. B37, 302 (1988). [9] R.J. Lipton, D.J. Rose and R.E. Tarjan, SIAM J. Numet. Anal. 16, 346 (1979). [10] J.-M. Normand, H.J. Herrmann and M. Hajjar, J. Stat. Phys. 52,441 (1988). [11] D.C. Hong, S. Havlin, H.J. Herrmann and H.E. Stanley, Phys. Rev. B30, 4083 (1984). [12] H.J. Herr~ann~ D.C. Hong and H.E. Stanley, J. Phys. A17, L261 (1984); T. Nagatard, J. Phys. A19, Ll165 (1986).

524

R.G. Edwards et aJ. ~Algebraic multi-grid method

[13] R.S. Varga, Matri~ Iterative Analysis (Prentice-Hall, Englewood Cliffs NJ, 1962).

[21] W.L. Briggs, A Multigrid Tutorial (SIAM, Philadelphia,

[14] L.A. Hageman and D.M. Young, Applied Iterative Methods (Academic Press, New York, 1981).

[22] J. Ruge and K. Stfiben, in Multigrid Methods, ed. S.F. McCormick (SIAM, Philadelphia, 1987).

[15] G.G. Batrouni, A. Hansen and M. Nelkin, Phys. Rev. Left. 571 1336 (1986).

[23] R. Hempel, private commlmication via J. Ruge (cecrces @csugold.bitnet).

[16] G.G. Batrouni and A. Hansenl J. Stat. Phys. 52, 747 (1988).

[24] B.A. Galler and M.J. Fischer, Commun. ACM 7, 301, 506 (1964); J. Hoshen and R. Kopelman, Phys. Rev. B14, 3438 (1976); R.G. Edwards and A.D. SokalI in preparation.

[17] R.G. Edwards, J. GoodmAn and A.D. Sokal, Phys. Rev. Lett. 61, 1333 (1988).

1987).

[18] W. Hackbusch I M~Iti-Grid Methods and Applications (Springer, Berlin, 1985).

[25] A.V. Alto, J.E. Hopcroft and J.D. Ullman, The Design and Analysis of Computer Algorithms (Addison-Wesley, Reading, Mass., 1974), pp. 176-187.

[19] P. Concus, G.H. Golub and D.P. O'Leary, in Sparse Matriz Computations Ied. J.R. Bunch and D.J. Rose (Academic Press, New York, 1976).

[26] J. Goodman and A.D. Sokal, Phys. Rev. Left. 56, 1015 (1986); and in preparation.

[20] G.G. Batrouni et al., Phys. Rev. D321 2736 (1985).

[27] A.A. Migdal, review talk at this conference.