- Email: [email protected]
On the magnetisation of the ground states in two dimensional Ising spin glasses
Computer Physics Communications 49 (1988) 417—421 North-Holland, Amsterdam
417
ON THE MAGNETISATION OF THE GROUND STATES IN TWO DIMENSIONAL ISING SPIN GLASSES * Francisco BARAHONA Department of Combinatorics and Optimization, University of Waterloo, Waterloo, Ontario, Canada N2L 3G1
and Adolfo CASARI Depto. de Matematica.s, Universidad de Chile, Casilla 5272, Santiago 3, Chile Received 30 August 1987
We study the configurations of minimum energy among those with zero magnetisation in Ising spin glasses. We formulate that problem as a graph partitioning problem and we give an algorithm for it. We determine the concentration of negative interactions such that those configurations give a ground state.
1. Introduction We study two dimensional Ising spin glasses on the square lattice with periodic boundary conditions. The exchange interactions J and J have the respective probabilities p and 1 p. Several authors have observed the existence of a critical concentration p~for which ferromagnetism disappears. Kirkpatrick [1] used the Monte Carlo method, that only gives an approximation of a ground state, to study 80 X 80 grids. He claimed that 0.15
—
*
Research partially supported by the Volkswagenwerk Stiftung, and by Departaniento de Investigación y Bibliotecas, Universidad de Chile.
in a different way. Our aim is to find for each p the configuration of minimum energy among those with magnetisation equal to zero. We compare this energy with the ground state energy. As we shall see in section 2 there is a value ~ such that both energies coincide for p ~ We shall also see that ~ is different from p~.Our algorithm can handle 9 X 9 (81 spins) grids. This can also be considered an approximation. Therefore we do not pretend to have computed an exact value of We want to present this new approach to the study of the magnetisation of the ground states, and to show what can be done with the present state of the art of integer programming. There are only two cases of spin glasses problems known to be polynomially solvable. One of them is the 2-dimensional case which reduces to a Chinese Postman Problem, cf. refs. [2,3]. The other case corresponds to the Random Field Ising ferromagnets which reduces to a max flow mm cut problem, this follows from the work of Picard and Ratliff [4] as pointed out in Barahona [5]. The problem we study in this paper is NP-complete, hence the existence of a polynomial algorithm is very unlikely. Our algonthm uses cutting planes m a first phase and then Branch and Bound. The recent
0010-4655/88/$03.50 © Elsevier Science Publishers B.V. (North-Holland Physics Publishing Division)
~.
~.
.
418
F. Barahona, A. Casari
/
Two dimensional Ising spin glasses
n
success of cutting planes methods in combinatorial optimization may suggest that one could handle problems of large size. Padberg and Rinaldi [6] have solved to optimality one Travelling Salesman Problem with 2392 cities. Barahona and Macciom [7] have found true ground states of 3-dimensional spin glasses in 5 X 5 x 5 grids (125 spins). Barahona, Grötschel, Junger and Reinelt [8] have found true ground states of 2-dimensional spin glasses with an exterior magnetic field in 40 x 40 grids (1600 spins). In the last two cases the spin glass problem reduces to a max cut problem in a graph. In all these examples, the algorithms are based on a good knowledge of the facial structure of the Travelling Salesman polytope, cf. GrOtschel and Padberg [9},and of the cut polytope, cf. Barahona and Mahjoub [10]. In this paper we require the magnetisation to be zero, then our problem reduces to a discrete quadratic programming problem with one constraint, it also reduces to a “graph partitioning problem”. Unfortunately, the facial structure of the graph partitioning polytope is not yet well understood. We could handle problems with 81 nodes, and it seems that some years will pass before one can solve problems that are substantially bigger. For discrete quadratic problems with one constraint, Gab, Hammer and Simeone [11] have reported results with problems up to 70 variables, Also Hartwig, Daske and Kobe [12] have reported results with up to 60 variables. For graph partitioning problems, Christophides and Brooker [13] report results with graphs up to 46 nodes, and Barnes [14] mention examples with 20 nodes. Therefore, we have slightly increased the size of the graph partitioning problems that can be solved at present. Given a square grid (graph) G = (V, E) and a function w: E { —J, J }, the energy of a spin configuration is given by the Hamiltonian
subject to
s,
~
=
0,
s, E (—1, 1), i = 1
n
This is a quadratic discrete optimization problem with one constraint. In what follows we shall reduce it to a combinatorial optimization problem m a graph. Since H= —~(w~~:s1=s~} + >{w1~:s~=—sj}, and ~
(w~:s1= s~)+ =
=
~
c,
(w1~:s1=
~Sj}
(j,j)EE
H= 2L(w,~:s,= —s~} C. —
Thus minimizing H is equivalent to minimize ~
(w~:s,= Our problem reduces to the following.
Equipartition of graphs: Given a graph G = (V, E), and a weight function w: E (—J, J) find a note set Uc V, such that UI =[IVI/2], and —~
~
j E U,
jE
V\ U
}
is minimum. Note that U = (i I s~= 1 }, V\ U = { i I = —1) and I U I = [I V 1/2] means E~ s = 0. This is an NP-hard problem even for planar graphs, hence we do not expect to find a polynomial time algorithm. This also suggests that we shall not be able to solve problems of large size. In section 2 we present an integer programming formulation of this problem and the numerical results.
—,
H=
~
—
w11s~s1.
2. Integer programming fonnulation and numerical results
(i,j)EE
Our problem can be stated as minimize
H
=
—
~
Given a graph G = (V, E) and a weight function w: E R, we associate to each node I E V a variable y, and to each edge ~j~ E a variable x~. Our problem can be formulated as follows: —~
F. Barahona~A. Ca.sari
/
Two dimensional Ising spin glasses
minimize ~
+
x,, ~ 2,
Before doing that we shall introduce some graph theoretic terminology. G = (V, a stable is in a set of Given nodes aS graph ~ V, such that E) there is no set edge E with both endnodes in S. If S is a stable set each
j’~_j’~—
x,1 ~ 0,
edge has at least one endnode in V\S.
ii
UE
subject to (a) y, + (b)
419
(c) —y, +y3
—
x,, ~ 0,
(2.1)
Lemma 2.3. Let S be a stable set of G, the
(d) —y,—y~+x,~
integrality of the variables (y1) for IE V\S implies the integrality of all the other variables.
:1
Proof. Let ij be an edge of E. If i and j belong to V\ S it is clear that x~,will be integer valued.
x,,
E
(0, 1
} for each edge ~/e E,
(0, 1) for each node i E V,
Let us assume that i ~S and j ~ V\S. If 0 inequalities 2.1(b) and (d) imply x,1 =~• ~f = 1 inequalities 2.1(a) and (c) imply x.j = 1 y,. Hence given an assignment of integer values to the variables { y~)for ~/~ V\ S our problem reduces to minimize ~ =
where n = I V I and k = [n/2]. The system of inequalities (2.1) was suggested to us by GrOtschel [15]. Let G’ be the graph obtained from G by adding a universal node (a node joined to all the other nodes of G), the inequalities (2.1) corresponds to the triangle inequalities of the cut polytope of G’. The triangle inequalities have been introduced by Barahona and Mahjoub [10]. This formulation represents a partition of V into U and V\U, y1=l if iEU and y1=O if i V\ U. For an edge u the variable ~ takes the value 1 if and only if i and j do not belong to the same set of the partition. We shall use integer programming techniques to solve this problem. As a first step we shall prove that we do not need to require all the variables to be integer valued, the integrality of some of them implies the integrality of the others. As a second step we shall introduce some new inequalities that cut fractional vectors but leave all the integer solutions, those inequalities are called cutting planes in the integer programming terminology. Remark 2.2. The integrality of the variables { y,) implies the integrality of the variables { x,1). Thus we can drop the constraints x,1 E (0, 1) for if E E. Hence the problem has I V I + E variables where I V I of them are required to be integer valued. We shall prove a Lemma that allows us to decrease again the number of integer variables.
—
subject to
~
iES ~
0 where k’
=
k—
=
k’
~ ~,~E 5, V\S
y1, and a,
=
E( w,1 x,~=
It is well known that this kind of linear program has an integer valued optimal solution. 0 Hence if S is a maximum stable set, out problem has I V I + E I variables but only I V\ S I are required to be integer valued. Good formulations are essential to solve integer programming problems efficiently. Integer programniing algorithms require a lower bound on the value of the minimum, and the efficiency of the algorithm is very dependent on the sharpness of the bound. For the integer programming problem .
mimmize subject to
cx Ax ~ b x {0, 1 }“, a lower bound is determined by solving the linear program minimize subject to
cx Ax ~ b 0 ~ x ~ 1.
F. Barahona, A. Casari / Two dimensional Ising spin glasses
420
Hence, given two systems of inequalities Ax ~ b and Dx ~ f such that
)
x I Ax ~ b, x E (0, 1) = {x I Dx ~f, x E (0, 1}”), and (x I Ax ~ b, 0 ~ x ~ 1) c
{ x I Dx ~
.~
0~
~ 1)
we say that the first system gives a better (or a tighter) formulation. The lower bound given by the first system is sharper than the bound given by the second system. In order to tighten our formulation we shall introduce three new sets of inequalities. Let G = (V, E) be a graph. Given U ç V the set C= (41EE liE U, jE V\U) is called a cut. A cycle is a set of edges i1i2, i2i3,. ..,i11i1, i1i~ we say that the cycle has length 1. A subiree is a set of edges that induces a connected subgraph and does not contain any cycle, Any integer vector x that satisfies the inequa]ities of (2.1) represents a cut of G. Hence x must satisfy >
1
for each subtree T of cardinality k + 1. These inequalities are called the “subtree inequalities”, For a cycle C of length greater than k + 1 the inequality ~ x~~ 2 Ii
C
must also be satisfied. These are called the “long cycle inequalities”, Since a cut has an even intersection with a cycle, the vector x must satisfy x, zjF
x~~
—
IFI
1,
yC\F
for each cycle C and F ç C with I F I odd. These are called the “cycle inequalities”. Now we are ready to describe the method of solution. First we apply the heuristic of Kemighan and Lin [16] to obtain an upper bound of the optimum value. Then we solve the linear programming problem obtained from (2.1) by dropping
the integrality requirement for the variables. This gives a lower bound for the optimum value. In order to decrease the gap between those two bounds we apply the cutting plane phase. The cutting plane phase consists of two steps that are applied iteratively: Step 1. Given a fractional vector x find a list of subtree inequalities, long cycle inequalities and cycle inequalities violated by x. Step 2. Add this set of inequalities to the current linear program and solve it again, this gives a new lower bound and a new vector x that is used in Step 1. We do not have a polynomial algorithm to find subtree and long cycle inequalities, which are NPcomplete problems. For the subtree inequalities, we look for a minimum weighted spanning tree, and remove edges until it has the required size. For the long cycle inequalities, we first look for a minimum weighted spanning tree, and then add one edge to it, if the cycle thus obtained has the required size we keep it, otherwise we try again with a new edge. These two procedures are just heuristics. For the third class of inequalities, the cycle inequalities, we have a polynomial algorithm, it appears in Barahona and Mahjoub [10]. It finds an inequality by solving a sequence of shortest path problems. When we cannot find any violated inequality in Step 1 we apply the Branch and Bound phase. Branch and Bound is the classical way of solving integer programming problems, cf. Papadimitriou and Steiglitz [17]. We used this algorithm in 9 x 9 grids, i.e. graphs with 81 nodes and 162 edges. The concentration of negative interactions was between 0.05 and 0.24. We used MPSX/MIP on an IBM 3031 computer. In a typical case the cutting plane phase took 26 mm of CPU time, this lead to a reduction of 50% of the grap between the lower and the upper bound. In this same case the branch and bound phase took 11 mm of CPU time. Table 1 shows our numerical results, they consist of averages of three samples. In fig. 1 we plot
F. Barahona A. Casari
/
Two dimensional Ising spin glasses
421
coincide for p It seems that 0.22 ~ ~ 0.25. It also seems that ~ is different from the critical
Table 1
~
p
Energy per spin
0.05 0.10 0.12 0.14 0.16 0.18 0.20 0.22 0.24
—1.51 —1.41 —1.41 —1.43 1.36 —1.39 —1.39 —1.41 —1.41
concentration p,~= 0.145 obtained by Bieche et al. [2]. This is reasonable, because for computing p~ they ask the magnetisation to be near zero, and for ~ we ask it to be exactly zero. Computing ~ with accuracy for bigger grids is a challenging open problem, but new developments in combinatorial optimization need to be made first. We conjecture that ~ is greater than p,~for grids of any size.
E
References -1.3 [11 S. Kirkpatrick, Phys. Rev. B16 (1977) 4630. [2] I. Bieche, R. Maynard, R. Rammal and J.P. Uhry, J. Phys. A 13 (1980) 2553. [3] F. Barahona, R. Maynard, R. Rainmal and J.P. Uhry, J. Phys. A 15 (1982) 673. [4] J.C. Picard and H.D. Ratliff, Networks 5 (1974) 357. [5] F. Barahona, J. Phys. A 18 (1985) L673. [6] M. Padberg and G. Rinaldi, Operations Res. Lett. 6
0
-1.4
0
-1.5
0 0
0
0
00
0
0 0 -1.6
0
0
-1.7 -1.8
(1987) 1. Barahona and E. Maccioni, J. Phys. A 15 (1982) L611. [81 F. Barahona, M. Grötschel, M. Junger and G. Reinelt, Report CORR 86-17, University of Waterloo (1986) to
[71F.
-1.9 I
0.1
I
I
0.2
0.3
P
Fig. 1.
with circles the energy per spin of the ground state obtained by Bieche et al. [2], we do not have their numerical results, therefore we have reproduced fig. 7 of their paper. They compute the groundstate energy for p in the interval 0.10 ~ p ~ 0.20 and for p = 0.30. We plot with squares the energy per spin obtained after imposing the magnetisalion to be zero. This graphic seems to indicate that there is a concentration ~ such that both energies
appear in Operations Research. [9] M. Grötschel and M. Padberg, Math. Programming 16 (1979) 265, 281. [10] F. Barahona and A.R. Mahjoub, Mathematical Programming 36 (1986) 157. [111 0. Gab, P.L. Hammer and B. Simeone, Math. Programming Study 12 (1980) 132. [12] A. Hartwig, F. Daske mun. 32 (1984) 133. and S. Kobe, Comput. Phys. Corn[13] N. Christophides and P. Brooker, SIAM J. Appl. Math. 30 (1976) 55. [14] E. Barnes, SIAM J. Alg. Disc. Methods 3 (1982) 541. [15] M. Grdtschel, private communication (1984). [16] 291. B.W. Kernighan and S. Lin, Bell System Tech. J. 49(1970) [17] C.H. Papadimitriou and K. Steigjitz, Combinatorial Optimization (Prentice Hall, Englewood Cliffs, 1982).
417
ON THE MAGNETISATION OF THE GROUND STATES IN TWO DIMENSIONAL ISING SPIN GLASSES * Francisco BARAHONA Department of Combinatorics and Optimization, University of Waterloo, Waterloo, Ontario, Canada N2L 3G1
and Adolfo CASARI Depto. de Matematica.s, Universidad de Chile, Casilla 5272, Santiago 3, Chile Received 30 August 1987
We study the configurations of minimum energy among those with zero magnetisation in Ising spin glasses. We formulate that problem as a graph partitioning problem and we give an algorithm for it. We determine the concentration of negative interactions such that those configurations give a ground state.
1. Introduction We study two dimensional Ising spin glasses on the square lattice with periodic boundary conditions. The exchange interactions J and J have the respective probabilities p and 1 p. Several authors have observed the existence of a critical concentration p~for which ferromagnetism disappears. Kirkpatrick [1] used the Monte Carlo method, that only gives an approximation of a ground state, to study 80 X 80 grids. He claimed that 0.15
—
*
Research partially supported by the Volkswagenwerk Stiftung, and by Departaniento de Investigación y Bibliotecas, Universidad de Chile.
in a different way. Our aim is to find for each p the configuration of minimum energy among those with magnetisation equal to zero. We compare this energy with the ground state energy. As we shall see in section 2 there is a value ~ such that both energies coincide for p ~ We shall also see that ~ is different from p~.Our algorithm can handle 9 X 9 (81 spins) grids. This can also be considered an approximation. Therefore we do not pretend to have computed an exact value of We want to present this new approach to the study of the magnetisation of the ground states, and to show what can be done with the present state of the art of integer programming. There are only two cases of spin glasses problems known to be polynomially solvable. One of them is the 2-dimensional case which reduces to a Chinese Postman Problem, cf. refs. [2,3]. The other case corresponds to the Random Field Ising ferromagnets which reduces to a max flow mm cut problem, this follows from the work of Picard and Ratliff [4] as pointed out in Barahona [5]. The problem we study in this paper is NP-complete, hence the existence of a polynomial algorithm is very unlikely. Our algonthm uses cutting planes m a first phase and then Branch and Bound. The recent
0010-4655/88/$03.50 © Elsevier Science Publishers B.V. (North-Holland Physics Publishing Division)
~.
~.
.
418
F. Barahona, A. Casari
/
Two dimensional Ising spin glasses
n
success of cutting planes methods in combinatorial optimization may suggest that one could handle problems of large size. Padberg and Rinaldi [6] have solved to optimality one Travelling Salesman Problem with 2392 cities. Barahona and Macciom [7] have found true ground states of 3-dimensional spin glasses in 5 X 5 x 5 grids (125 spins). Barahona, Grötschel, Junger and Reinelt [8] have found true ground states of 2-dimensional spin glasses with an exterior magnetic field in 40 x 40 grids (1600 spins). In the last two cases the spin glass problem reduces to a max cut problem in a graph. In all these examples, the algorithms are based on a good knowledge of the facial structure of the Travelling Salesman polytope, cf. GrOtschel and Padberg [9},and of the cut polytope, cf. Barahona and Mahjoub [10]. In this paper we require the magnetisation to be zero, then our problem reduces to a discrete quadratic programming problem with one constraint, it also reduces to a “graph partitioning problem”. Unfortunately, the facial structure of the graph partitioning polytope is not yet well understood. We could handle problems with 81 nodes, and it seems that some years will pass before one can solve problems that are substantially bigger. For discrete quadratic problems with one constraint, Gab, Hammer and Simeone [11] have reported results with problems up to 70 variables, Also Hartwig, Daske and Kobe [12] have reported results with up to 60 variables. For graph partitioning problems, Christophides and Brooker [13] report results with graphs up to 46 nodes, and Barnes [14] mention examples with 20 nodes. Therefore, we have slightly increased the size of the graph partitioning problems that can be solved at present. Given a square grid (graph) G = (V, E) and a function w: E { —J, J }, the energy of a spin configuration is given by the Hamiltonian
subject to
s,
~
=
0,
s, E (—1, 1), i = 1
n
This is a quadratic discrete optimization problem with one constraint. In what follows we shall reduce it to a combinatorial optimization problem m a graph. Since H= —~(w~~:s1=s~} + >{w1~:s~=—sj}, and ~
(w~:s1= s~)+ =
=
~
c,
(w1~:s1=
~Sj}
(j,j)EE
H= 2L(w,~:s,= —s~} C. —
Thus minimizing H is equivalent to minimize ~
(w~:s,= Our problem reduces to the following.
Equipartition of graphs: Given a graph G = (V, E), and a weight function w: E (—J, J) find a note set Uc V, such that UI =[IVI/2], and —~
~
j E U,
jE
V\ U
}
is minimum. Note that U = (i I s~= 1 }, V\ U = { i I = —1) and I U I = [I V 1/2] means E~ s = 0. This is an NP-hard problem even for planar graphs, hence we do not expect to find a polynomial time algorithm. This also suggests that we shall not be able to solve problems of large size. In section 2 we present an integer programming formulation of this problem and the numerical results.
—,
H=
~
—
w11s~s1.
2. Integer programming fonnulation and numerical results
(i,j)EE
Our problem can be stated as minimize
H
=
—
~
Given a graph G = (V, E) and a weight function w: E R, we associate to each node I E V a variable y, and to each edge ~j~ E a variable x~. Our problem can be formulated as follows: —~
F. Barahona~A. Ca.sari
/
Two dimensional Ising spin glasses
minimize ~
+
x,, ~ 2,
Before doing that we shall introduce some graph theoretic terminology. G = (V, a stable is in a set of Given nodes aS graph ~ V, such that E) there is no set edge E with both endnodes in S. If S is a stable set each
j’~_j’~—
x,1 ~ 0,
edge has at least one endnode in V\S.
ii
UE
subject to (a) y, + (b)
419
(c) —y, +y3
—
x,, ~ 0,
(2.1)
Lemma 2.3. Let S be a stable set of G, the
(d) —y,—y~+x,~
integrality of the variables (y1) for IE V\S implies the integrality of all the other variables.
:1
Proof. Let ij be an edge of E. If i and j belong to V\ S it is clear that x~,will be integer valued.
x,,
E
(0, 1
} for each edge ~/e E,
(0, 1) for each node i E V,
Let us assume that i ~S and j ~ V\S. If 0 inequalities 2.1(b) and (d) imply x,1 =~• ~f = 1 inequalities 2.1(a) and (c) imply x.j = 1 y,. Hence given an assignment of integer values to the variables { y~)for ~/~ V\ S our problem reduces to minimize ~ =
where n = I V I and k = [n/2]. The system of inequalities (2.1) was suggested to us by GrOtschel [15]. Let G’ be the graph obtained from G by adding a universal node (a node joined to all the other nodes of G), the inequalities (2.1) corresponds to the triangle inequalities of the cut polytope of G’. The triangle inequalities have been introduced by Barahona and Mahjoub [10]. This formulation represents a partition of V into U and V\U, y1=l if iEU and y1=O if i V\ U. For an edge u the variable ~ takes the value 1 if and only if i and j do not belong to the same set of the partition. We shall use integer programming techniques to solve this problem. As a first step we shall prove that we do not need to require all the variables to be integer valued, the integrality of some of them implies the integrality of the others. As a second step we shall introduce some new inequalities that cut fractional vectors but leave all the integer solutions, those inequalities are called cutting planes in the integer programming terminology. Remark 2.2. The integrality of the variables { y,) implies the integrality of the variables { x,1). Thus we can drop the constraints x,1 E (0, 1) for if E E. Hence the problem has I V I + E variables where I V I of them are required to be integer valued. We shall prove a Lemma that allows us to decrease again the number of integer variables.
—
subject to
~
iES ~
0 where k’
=
k—
=
k’
~ ~,~E 5, V\S
y1, and a,
=
E( w,1 x,~=
It is well known that this kind of linear program has an integer valued optimal solution. 0 Hence if S is a maximum stable set, out problem has I V I + E I variables but only I V\ S I are required to be integer valued. Good formulations are essential to solve integer programming problems efficiently. Integer programniing algorithms require a lower bound on the value of the minimum, and the efficiency of the algorithm is very dependent on the sharpness of the bound. For the integer programming problem .
mimmize subject to
cx Ax ~ b x {0, 1 }“, a lower bound is determined by solving the linear program minimize subject to
cx Ax ~ b 0 ~ x ~ 1.
F. Barahona, A. Casari / Two dimensional Ising spin glasses
420
Hence, given two systems of inequalities Ax ~ b and Dx ~ f such that
)
x I Ax ~ b, x E (0, 1) = {x I Dx ~f, x E (0, 1}”), and (x I Ax ~ b, 0 ~ x ~ 1) c
{ x I Dx ~
.~
0~
~ 1)
we say that the first system gives a better (or a tighter) formulation. The lower bound given by the first system is sharper than the bound given by the second system. In order to tighten our formulation we shall introduce three new sets of inequalities. Let G = (V, E) be a graph. Given U ç V the set C= (41EE liE U, jE V\U) is called a cut. A cycle is a set of edges i1i2, i2i3,. ..,i11i1, i1i~ we say that the cycle has length 1. A subiree is a set of edges that induces a connected subgraph and does not contain any cycle, Any integer vector x that satisfies the inequa]ities of (2.1) represents a cut of G. Hence x must satisfy >
1
for each subtree T of cardinality k + 1. These inequalities are called the “subtree inequalities”, For a cycle C of length greater than k + 1 the inequality ~ x~~ 2 Ii
C
must also be satisfied. These are called the “long cycle inequalities”, Since a cut has an even intersection with a cycle, the vector x must satisfy x, zjF
x~~
—
IFI
1,
yC\F
for each cycle C and F ç C with I F I odd. These are called the “cycle inequalities”. Now we are ready to describe the method of solution. First we apply the heuristic of Kemighan and Lin [16] to obtain an upper bound of the optimum value. Then we solve the linear programming problem obtained from (2.1) by dropping
the integrality requirement for the variables. This gives a lower bound for the optimum value. In order to decrease the gap between those two bounds we apply the cutting plane phase. The cutting plane phase consists of two steps that are applied iteratively: Step 1. Given a fractional vector x find a list of subtree inequalities, long cycle inequalities and cycle inequalities violated by x. Step 2. Add this set of inequalities to the current linear program and solve it again, this gives a new lower bound and a new vector x that is used in Step 1. We do not have a polynomial algorithm to find subtree and long cycle inequalities, which are NPcomplete problems. For the subtree inequalities, we look for a minimum weighted spanning tree, and remove edges until it has the required size. For the long cycle inequalities, we first look for a minimum weighted spanning tree, and then add one edge to it, if the cycle thus obtained has the required size we keep it, otherwise we try again with a new edge. These two procedures are just heuristics. For the third class of inequalities, the cycle inequalities, we have a polynomial algorithm, it appears in Barahona and Mahjoub [10]. It finds an inequality by solving a sequence of shortest path problems. When we cannot find any violated inequality in Step 1 we apply the Branch and Bound phase. Branch and Bound is the classical way of solving integer programming problems, cf. Papadimitriou and Steiglitz [17]. We used this algorithm in 9 x 9 grids, i.e. graphs with 81 nodes and 162 edges. The concentration of negative interactions was between 0.05 and 0.24. We used MPSX/MIP on an IBM 3031 computer. In a typical case the cutting plane phase took 26 mm of CPU time, this lead to a reduction of 50% of the grap between the lower and the upper bound. In this same case the branch and bound phase took 11 mm of CPU time. Table 1 shows our numerical results, they consist of averages of three samples. In fig. 1 we plot
F. Barahona A. Casari
/
Two dimensional Ising spin glasses
421
coincide for p It seems that 0.22 ~ ~ 0.25. It also seems that ~ is different from the critical
Table 1
~
p
Energy per spin
0.05 0.10 0.12 0.14 0.16 0.18 0.20 0.22 0.24
—1.51 —1.41 —1.41 —1.43 1.36 —1.39 —1.39 —1.41 —1.41
concentration p,~= 0.145 obtained by Bieche et al. [2]. This is reasonable, because for computing p~ they ask the magnetisation to be near zero, and for ~ we ask it to be exactly zero. Computing ~ with accuracy for bigger grids is a challenging open problem, but new developments in combinatorial optimization need to be made first. We conjecture that ~ is greater than p,~for grids of any size.
E
References -1.3 [11 S. Kirkpatrick, Phys. Rev. B16 (1977) 4630. [2] I. Bieche, R. Maynard, R. Rammal and J.P. Uhry, J. Phys. A 13 (1980) 2553. [3] F. Barahona, R. Maynard, R. Rainmal and J.P. Uhry, J. Phys. A 15 (1982) 673. [4] J.C. Picard and H.D. Ratliff, Networks 5 (1974) 357. [5] F. Barahona, J. Phys. A 18 (1985) L673. [6] M. Padberg and G. Rinaldi, Operations Res. Lett. 6
0
-1.4
0
-1.5
0 0
0
0
00
0
0 0 -1.6
0
0
-1.7 -1.8
(1987) 1. Barahona and E. Maccioni, J. Phys. A 15 (1982) L611. [81 F. Barahona, M. Grötschel, M. Junger and G. Reinelt, Report CORR 86-17, University of Waterloo (1986) to
[71F.
-1.9 I
0.1
I
I
0.2
0.3
P
Fig. 1.
with circles the energy per spin of the ground state obtained by Bieche et al. [2], we do not have their numerical results, therefore we have reproduced fig. 7 of their paper. They compute the groundstate energy for p in the interval 0.10 ~ p ~ 0.20 and for p = 0.30. We plot with squares the energy per spin obtained after imposing the magnetisalion to be zero. This graphic seems to indicate that there is a concentration ~ such that both energies
appear in Operations Research. [9] M. Grötschel and M. Padberg, Math. Programming 16 (1979) 265, 281. [10] F. Barahona and A.R. Mahjoub, Mathematical Programming 36 (1986) 157. [111 0. Gab, P.L. Hammer and B. Simeone, Math. Programming Study 12 (1980) 132. [12] A. Hartwig, F. Daske mun. 32 (1984) 133. and S. Kobe, Comput. Phys. Corn[13] N. Christophides and P. Brooker, SIAM J. Appl. Math. 30 (1976) 55. [14] E. Barnes, SIAM J. Alg. Disc. Methods 3 (1982) 541. [15] M. Grdtschel, private communication (1984). [16] 291. B.W. Kernighan and S. Lin, Bell System Tech. J. 49(1970) [17] C.H. Papadimitriou and K. Steigjitz, Combinatorial Optimization (Prentice Hall, Englewood Cliffs, 1982).