One Network to Approximate Them All: Amortized
Variational Inference of Ising Ground States

Sebastian Sanokowski1

Wilhelm Berghammer1

Sepp Hochreiter1,21Johannes Kofler1

Sebastian Lehner1,2

ELLIS Unit Linz and LIT AI Lab, Institute for Machine Learning, Johannes Kepler University Linz, Austria
2Institute of Advanced Research in Artificial Intelligence (IARAI), Vienna, Austria
sanokowski@ml.jku.at

Abstract
For a wide range of combinatorial optimization problems, finding the optimal
solutions is equivalent to finding the ground states of corresponding Ising Hamiltonians. Recent work shows that these ground states are found more efficiently by
variational approaches using autoregressive models than by traditional methods.
In contrast to previous works, where for every problem instance a new model has
to be trained, we aim at a single model that approximates the ground states fora whole family of Hamiltonians. We demonstrate that autoregregressive neuralnetworks can be trained to achieve this goal and are able to generalize across a classof problems. We iteratively approximate the ground state based on a representationof the Hamiltonian that is provided by a graph neural network. Our experiments
show that solving a large number of related problem instances by a single model
can be considerably more efficient than solving them individually.

1

Introduction

Ising models are of broad interest since they exhibit a remarkably wide range of interesting phenomena
in statistical mechanics and provide a universal model of classical spin physics [1; 2]. Besides their
popularity in physics, they are frequently applied in disciplines like e.g. molecular biology, operations
research, and neuroscience.
The energy of an Ising system with N spins is given by the following Hamiltonian:
H=−

X
i<j

Jij σi σj −NX

hi σi ,

(1)i=1where σi , σj ∈ {−1, 1} are the spins, i.e. the system’s degrees of freedom which are locatedat lattice sites indexed by i, j ∈ {1, . . . , N }. The interactions among these spins are given by
the couplings Jij ∈ R in the first term. The second term represents couplings to external fields
hi ∈ R. For the problems considered in this work, this term vanishes, because external fields
can be absorbed into the coupling term. A specific instance of an Ising model is thus in generaldefined by the couplings and the external fields. The problem of finding the spin configuration withthe lowest energy, i.e the ground state, is known to be equivalent to numerous NP-hard problems
including combinatorial optimization (CO) problems, like Max-Cut and the Traveling Salesperson
problem [3]. Consequently, while the problem of finding or approximating the ground states ofIsing Hamiltonians is of general interest, under the commonly held P ̸= NP assumption, it is
excluded that solutions to all problem instances can be found efficiently. Nevertheless, this does
Machine Learning and the Physical Sciences workshop, NeurIPS 2022.

not imply that it is impossible that an optimal solution or an approximation might be efficiently
computable for many practically occurring problem instances. Traditionally, methods like simulated
annealing (SA) [4] are used to approximate ground states. More recently, variational methods using
autoregressive neural networks were shown to yield better approximations of the system’s Boltzmann
distribution (or of the wave function in the corresponding quantum setting) of Ising-type problems
and CO problems [5; 6; 7; 8; 9; 10; 11]. An attractive feature of these methods is that they are
general, i.e. they do not rely on a combination with other algorithms or problem-specific heuristics
and do not require ground truth data, which is typically problematic [12; 13]. The work of [14]employs meta-learning methods to similar variational problems to obtain more efficient optimization
procedures for individual problem instances. However, also here each problem instance requires
dedicated fine-tuning. In time and computationally constrained scenarios, however, these approaches
may be infeasible since they require a potentially lengthy variational optimization procedure for eachnew problem instance. We, therefore, aim at a model that approximates solutions to a family of CO
problems with a fixed set of model parameters.
Since many CO problems naturally admit a representation as graphs, numerous Graph Neural Network(GNN) [15] applications to CO problems have been demonstrated. As pointed out in a recent review
article [16], most of these applications rely either on the availability of ground truth solutions or on
the combination with traditional algorithms or handcrafted cost functions. Notable exceptions are
the GNN-based approaches in [17; 18]. Similar to the variational approaches discussed above, theserequire an individual optimization procedure for new problem instances as well. This limitation does
not apply to the GNN-based approaches in [19; 20]. However, here the learned distribution overstates are assumed to factorize with respect to individual state variables and thus neglect correlations
among them. As pointed out in [5], these mean-field approximations are in general inferior comparedto the autoregressive neural networks used in aforementioned variational approaches and in this work.
The contribution of the present work is threefold. First, we extend the variational approaches byintroducing the additional goal of generalization over families of problem instances. This aspectis analogous to amortized variational inference where the solution to a variational optimization
problem is directly approximated by the output of a learned function and not the result of a dedicatedoptimization procedure [21]. We will refer to our method as Amortized Variational Annealing (AVA).
Second, we investigate how the size of the models affect the achievable solution quality and how themodel generalises to different problem sizes. Finally, we define a new measure of problem instancedifficulty which appears to be a useful indicator for the hardness of CO problem instances.

2

Problem Setting and Methods

We consider the problem of approximating the ground state σ0 ∈ {−1, +1}N of Ising Hamiltonians
with N spins as defined in (Eq. 1). We take a variational approach to this problem and optimize thevariational parameters θ of a probability distribution pθ (σ) over the states σ ∈ {−1, +1}N such thata high probability mass is assigned to σ0 . Once a sufficiently good approximation θ̂ of the optimalvariational parameters is found, the ground state is likely to be obtained by sampling states frompθ̂ (σ) and choosing the sample with the lowest energy. Instead of optimizing pθ (σ) for an individual
Hamiltonian Hτ we are interested in constructing a parameterization of pθ (σ|Hτ ) that adapts to
the Hamiltonian of interest without requiring an optimization of θ. That is, we aim at finding a
distribution that is conditioned on Hamiltonians which are sampled from a distribution pH (Hτ ).
In order to prevent our algorithm from being stuck in local minima, we use the Variational ClassicalAnnealing (VCA) approach, where an entropy term is added to the minimization objective andannealed from high to low temperature T (A.5.3). The resulting task is to minimize the following
variational free energy with respect to θ:


min E
Hτ (σ) + T log pθ (σ|Hτ ) .
(2)EθHτ ∼pH σ∼pθ (σ|Hτ )Our model AVA (see Fig. 1) is built up of three components (A.4). The first part consists of a GNN
that maps the graph of an Ising Hamiltonian to a node-wise embedding. The embedding of the Ising
Hamiltonian is obtained by mapping its couplings on the edges of a graph. We then use a GNN withseveral message passing steps [22] to obtain an embedding ei (Hτ ) for each node. Then, an LSTM
[23] is used which iteratively processes these embeddings along with the previously generated spinvalue σi−1 . Finally, an MLP maps the hidden state of the LSTM to the parameter of a Bernoulli
distribution for the state at node i. A spin configuration is generated by sampling these distributions
2

Figure 1: An illustration of our architecture. We first run a forward pass through a GNN in order toobtain a Hamiltonian embedding ei (Hτ ) for each node. Then, we iterate in an autoregressive mannerover each node in the graph, where the LSTM receives the GNN embedding of the current node iand the previously generated spin value i − 1 as an input. Here, i denotes the index of the currently
visited site and < i denotes indices of already visited nodes.

Figure 2: Left: ∆Erel over number of annealing steps for AVA with different types of Hamiltonian
embeddings. The network is evaluated on 100 SK-16 test Hamiltonians. Right: ∆Erel over thenumber of annealing steps for different system sizes in the training and test set. The model is trainedon N spins and evaluated on Ñ spins. ∆Erel was evaluated on 300 test Hamiltonians and 50 sampled
spin configurations. Error bars in both plots indicate the standard deviation over three independent
runs.
autoregressively. The corresponding variational distribution of our autoregressive model can therefore
be written asNYpθ (σ|Hτ ) =pθ (σi |σ<i , e≤i (Hτ )).
(3)
i=1

3

Experiments

In the following, we report results for two classes of Hamiltonians. One class are the SherringtonKirkpatrick spin-glass Hamiltonians (A.1.1) with N spins, which we denote as SK-N , and the otherclass of Hamiltonians represents the Max-Cut problem (A.1.2).

Expressivity of the Embedding Network: To investigate how the model performance depends onthe expressivity of the Hamiltonian embedding, we plot in Fig. 2 (left) the relative energy difference∆Erel (A.1.3) over the number of annealing steps for different numbers of GNN layers on SK-16.
The model tends to perform better, the more GNN layers are used. To characterize the inductive biasof the GNN, we compare the GNN embedding with an embedding generated by a comparable MLP(A.4.3). The results indicate that the GNN embeddings allow for more efficient annealing.
Generalisation to other System Sizes: Another advantage of the GNNs compared to the MLPs
is that once they were trained on some system size they can in principle be applied to systems ofany size. To test the generalization capability to other system sizes, we train the GNN-based model
on SK-16. Then, we make a prediction for different system sizes and evaluate ∆Erel using 300
Hamiltonians per system size. The results in Fig. 2 (right) show that the prediction on smaller systems
3

becomes slightly better and the prediction on larger systems becomes slightly worse compared to theperformance on SK-16 itself. The model that is trained on SK-16 and applied on SK-21 (red star) is
only slightly worse than the model that is directly trained on SK-21 (purple triangle, shown only for
Nannealing = 105 ).
Amortisation: Since AVA learns a conditional distribution over a distribution of Hamiltonians
while in Variational Neural Annealing (VNA) [8] each model is only trained on one specificHamiltonian, the task of VNA is considerably easier. Thus, it is not unexpected that for a given
amount of annealing steps, VNA performs better than AVA (see direct comparison of decay rates
in A.3). On the other hand, an important advantage of our method is that our model is able toapproximate ground states of i.i.d. Hamiltonians at inference time, whereas in VNA a new model hast
to be trained from scratch for every new Hamiltonian. Therefore, given a certain target ∆Erel we estimate in Fig. 3 (left) at which number of Hamiltonians it is faster to use AVA compared to VNA (A.3).

Hardness Measure: Since not all problem instances from a certain CO problem class are equally
hard, [24] propose a simplicity criterion which indicates whether a given CO problem instance canbe solved by simply diagonalizing the corresponding Ising Hamiltonian Hτ and discretizing the stateτvmin∈ RN that minimizes the energy. Based on this simplicity criterion, we introduce for eachproblem instance Hτ a problem hardness measure Υ(Hτ ) ∈ [0, 1] that is based on the Hammingτ
distance d between a discretization of vminand the whole set of all k ground states {σ0,1 , ..., σ0,k }.
We define our hardness measure accordingly as
h
i
1
τ
min
d sgn(vmin), σ ,
Υ(Hτ ) =(4)N σ∈{σ0,1 ,...,σ0,k }where larger values indicate harder instances. In Fig. 3 (center), we plot the ground state probabilitypθ (σ0 |Hτ ) and ∆Erel from our model over our hardness measure Υ(Hτ ). The latter clearlycorrelates with the performance of our model, which indicates that this measure is a suitable todetermine the hardness of problem instances.

Maximum Cut: In addition, we conduct experiments on Hamiltonians representing the Max-Cut
problem of size N = 50 and with connection density 0.5 (A.1.2). We compare our method to theGoemans-Williamson (GW) algorithm [25], which is a polynomial-time algorithm that achieves anapproximation ratio of ≈ 0.878. This means that even for the worst-case Max-Cut instance GW
yields a fairly good expected solution quality [26]. In Fig. 3 (right), we compare GW with AVA bycomputing for each of 100 test Hamiltonians the Max-Cut ratio (A.1.4) of 50 sampled states. Here,
AVA often generates states with a better Max-Cut ratio than the GW algorithm. On average, ourmethod achieves a Max-Cut ratio of ≈ 0.988 which is better than the average Max-Cut ratio of GW
of ≈ 0.971 (see A.1.4 for further comparisons).

4

Conclusion and Outlook

In this paper, we use GNNs combined with LSTMs to approximate ground states for a familiy ofIsing Hamiltonians without adapting the model to individual Hamiltonians. Our experiments showthat longer annealing and more expressive GNNs improve the solution quality considerably. We
demonstrate that our method generalises to larger system sizes without any fine tuning.
For the problem sizes studied in this work, exact solutions are still routinely obtained by traditional
methods. Future work could investigate how our method can be extended to larger systems and how
it performs on real-world data sets. Additionally, it might be interesting to investigate whether a ourmodels can be used as a starting point for problem instance specific fine tuning. Furthermore, we aim
to extend the evaluation of our hardness measure on further algorithms and problem classes.

5

Broader Impact

In this work, we develop and characterize a novel approach to combinatorial optimization problems
which aims to approximate solutions for problem classes without any instance-specific retraining, i.e.
4Figure 3: Left: The estimated number of Hamiltonians namort above which it is faster to use AVA
instead of VNA in terms of annealing steps, while reaching the same relative energy difference. Error
bars indicate errors from the fitted decay rates (A.3). For VNA we use the same architecture as for
AVA (A.5). Center: ∆Erel and probability of the ground state plotted over the hardness measuredefined in Eq. 4. Here, due to spin flip symmetry of the SK model, the hardest instances have ahardness of 0.5. The model was evaluated on SK 16 and the error bars indicate the standard error
over three independent runs. Right: Comparison of the Max-Cut ratio of 50 sampled states from the
AVA and the GW algorithm on the Max-Cut-50 problem with density 0.5. Our AVA model is trainedfor 2.5 · 104 annealing steps. The evaluation is done on 100 Hamiltonians.
with a fixed set of model parameters. We expect that this method will be useful for many real-world
problems in operations research and eventually for applications in material design.

Acknowledgments and Disclosure of Funding
We thank Giuseppe Carleo for his useful comments to this work.
The ELLIS Unit Linz, the LIT AI Lab, the Institute for Machine Learning, are supported by theFederal State Upper Austria. IARAI is supported by Here Technologies. We thank the projects AIMOTION (LIT-2018-6-YOU-212), AI-SNN (LIT-2018-6-YOU-214), DeepFlood (LIT-2019-8-YOU213), Medical Cognitive Computing Center (MC3), INCONTROL-RL (FFG-881064), PRIMAL(FFG-873979), S3AI (FFG-872172), DL for GranularFlow (FFG-871302), AIRI FG 9-N (FWF36284, FWF-36235), ELISE (H2020-ICT-2019-3 ID: 951847). We thank Audi.JKU Deep Learning
Center, TGW LOGISTICS GROUP GMBH, Silicon Austria Labs (SAL), FILL Gesellschaft mbH,
Anyline GmbH, Google, ZF Friedrichshafen AG, Robert Bosch GmbH, UCB Biopharma SRL, Merck
Healthcare KGaA, Verbund AG, Software Competence Center Hagenberg GmbH, TÜV Austria,
Frauscher Sensonic and the NVIDIA Corporation.

References[1] A.J. Berlinsky and A.B. Harris. Statistical Mechanics: An Introductory Graduate Course. Graduate Texts
in Physics. Springer International Publishing, 2019. ISBN 9783030281878.
[2] Gemma De las Cuevas and Toby S Cubitt. Simple universal models capture all classical spin physics.
Science, 351(6278):1180–1183, 2016.
[3] Andrew Lucas. Ising formulations of many np problems. Frontiers in physics, page 5, 2014.
[4] Scott Kirkpatrick, C Daniel Gelatt Jr, and Mario P Vecchi. Optimization by simulated annealing. science,
220(4598):671–680, 1983.
[5] Dian Wu, Lei Wang, and Pan Zhang. Solving statistical mechanics using variational autoregressive
networks. Physical review letters, 122(8):080602, 2019.
[6] Or Sharir, Yoav Levine, Noam Wies, Giuseppe Carleo, and Amnon Shashua. Deep autoregressive models
for the efficient variational simulation of many-body quantum systems. Physical review letters, 124(2):020503, 2020.
[7] Mohamed Hibat-Allah, Martin Ganahl, Lauren E Hayward, Roger G Melko, and Juan Carrasquilla.
Recurrent neural network wave functions. Physical Review Research, 2(2):023358, 2020.
[8] Mohamed Hibat-Allah, Estelle M Inack, Roeland Wiersema, Roger G Melko, and Juan Carrasquilla.
Variational neural annealing. Nature Machine Intelligence, 3(11):952–961, 2021.

5[9] Shoummo Ahsan Khandoker, Jawaril Munshad Abedin, and Mohamed Hibat-Allah. Supplementing
recurrent neural networks with annealing to solve optimization problems. arXiv preprint arXiv:2207.08189,
2022.
[10] Tianchen Zhao, Giuseppe Carleo, James Stokes, and Shravan Veerapaneni. Natural evolution strategies
and variational monte carlo. Machine Learning: Science and Technology, 2(2):02LT01, 2020.
[11] Joseph Gomes, Keri A McKiernan, Peter Eastman, and Vijay S Pande. Classical quantum optimization
with neural network quantum states. arXiv preprint arXiv:1910.10675, 2019.
[12] Gal Yehuda, Moshe Gabel, and Assaf Schuster. It’s not what machines can learn, it’s what we cannot teach.
In International conference on machine learning, pages 10831–10841. PMLR, 2020.
[13] Yoshua Bengio, Andrea Lodi, and Antoine Prouvost. Machine learning for combinatorial optimization: amethodological tour d’horizon. European Journal of Operational Research, 290(2):405–421, 2021.
[14] Tianchen Zhao, James Stokes, Oliver Knitter, Brian Chen, and Shravan Veerapaneni. Meta variationalmonte carlo. arXiv preprint arXiv:2011.10614, 2020.
[15] Franco Scarselli, Marco Gori, Ah Chung Tsoi, Markus Hagenbuchner, and Gabriele Monfardini. The
graph neural network model. IEEE transactions on neural networks, 20(1):61–80, 2008.
[16] Quentin Cappart, Didier Chételat, Elias Khalil, Andrea Lodi, Christopher Morris, and Petar Veličković.
Combinatorial optimization and reasoning with graph neural networks. arXiv preprint arXiv:2102.09544,
2021.
[17] Nikolaos Karalias and Andreas Loukas. Erdos goes neural: an unsupervised learning framework for
combinatorial optimization on graphs. Advances in Neural Information Processing Systems, 33:6659–6672,
2020.
[18] Martin JA Schuetz, J Kyle Brubaker, and Helmut G Katzgraber. Combinatorial optimization with physicsinspired graph neural networks. Nature Machine Intelligence, 4(4):367–377, 2022.
[19] Jan Toenshoff, Martin Ritzert, Hinrikus Wolf, and Martin Grohe. Graph neural networks for maximum
constraint satisfaction. Frontiers in artificial intelligence, 3:580607, 2021.
[20] Haoran Sun, Etash K Guha, and Hanjun Dai. Annealed training for combinatorial optimization on graphs.
arXiv preprint arXiv:2207.11542, 2022.
[21] Cheng Zhang, Judith Bütepage, Hedvig Kjellström, and Stephan Mandt. Advances in variational inference.
IEEE transactions on pattern analysis and machine intelligence, 41(8):2008–2026, 2018.
[22] Peter W Battaglia, Jessica B Hamrick, Victor Bapst, Alvaro Sanchez-Gonzalez, Vinicius Zambaldi, Mateusz
Malinowski, Andrea Tacchetti, David Raposo, Adam Santoro, Ryan Faulkner, et al. Relational inductive
biases, deep learning, and graph networks. arXiv preprint arXiv:1806.01261, 2018.
[23] Sepp Hochreiter and Jürgen Schmidhuber. Long short-term memory. Neural computation, 9(8):1735–1780,
1997.
[24] Kirill P. Kalinin and Natalia G. Berloff. Computational complexity continuum within ising formulation of
np problems. Communications Physics, 5:1–10, 2022.
[25] Michel X Goemans and David P Williamson. Improved approximation algorithms for maximum cut andsatisfiability problems using semidefinite programming. Journal of the ACM (JACM), 42(6):1115–1145,
1995.
[26] Uriel Feige and Gideon Schechtman. On the optimality of the random hyperplane rounding technique formax cut. Random Structures & Algorithms, 20(3):403–440, 2002.
[27] David Sherrington and Scott Kirkpatrick. Solvable model of a spin-glass. Physical review letters, 35(26):1792, 1975.
[28] Jimmy Lei Ba, Jamie Ryan Kiros, and Geoffrey E Hinton. Layer normalization. arXiv preprint
arXiv:1607.06450, 2016.
[29] Abien Fred Agarap. Deep learning using rectified linear units (relu). arXiv preprint arXiv:1803.08375,
2018.
[30] Liyuan Liu, Haoming Jiang, Pengcheng He, Weizhu Chen, Xiaodong Liu, Jianfeng Gao, and Jiawei Han.
On the variance of the adaptive learning rate and beyond. arXiv preprint arXiv:1908.03265, 2019.

6

[31] Ilya Loshchilov and Frank Hutter. Sgdr: Stochastic gradient descent with warm restarts. arXiv preprint
arXiv:1608.03983, 2016.
[32] Filippo Vicentini, Damian Hofmann, Attila Szabó, Dian Wu, Christopher Roth, Clemens Giuliani, Gabriel
Pescia, Jannes Nys, Vladimir Vargas-Calderón, Nikita Astrakhantsev, et al. Netket 3: Machine learning
toolbox for many-body quantum systems. SciPost Physics Codebases, page 007, 2022.
[33] James Bradbury, Roy Frostig, Peter Hawkins, Matthew James Johnson, Chris Leary, Dougal Maclaurin, George Necula, Adam Paszke, Jake VanderPlas, Skye Wanderman-Milne, et al. Jax: composable
transformations of python+ numpy programs. Version 0.2, 5:14–24, 2018.

Checklist
1. For all authors...
(a) Do the main claims made in the abstract and introduction accurately reflect the paper’s contributions and scope? [Yes] We support our claims in Section 3.
(b) Did you describe the limitations of your work? [Yes] In Section 4 we discuss the limitations of
our work.
(c) Did you discuss any potential negative societal impacts of your work? [N/A](d) Have you read the ethics review guidelines and ensured that your paper conforms to them? [Yes]2. If you are including theoretical results...
(a) Did you state the full set of assumptions of all theoretical results? [N/A](b) Did you include complete proofs of all theoretical results? [N/A]
3. If you ran experiments...
(a) Did you include the code, data, and instructions needed to reproduce the main experimental
results (either in the supplemental material or as a URL)? [No] We did not yet include code
because this work is in progress. However, a public code repository is under preparation.
(b) Did you specify all the training details (e.g., data splits, hyperparameters, how they werechosen)? [Yes] See Section A.5 and also A.4 for information about hyperparameters and theircorresponding grid searches. For more information on how the data was sampled see A.1.1 andA.1.2.
(c) Did you report error bars (e.g., with respect to the random seed after running experiments
multiple times)? [Yes] We report error bars with respect to multiple random seeds in all of our
figures.
(d) Did you include the total amount of compute and the type of resources used (e.g., type of GPUs,
internal cluster, or cloud provider)? [Yes] In section A.5.5 we state on which GPUs the models
for different system sizes were be trained on. I addition we report the number of update steps in
the training for each of our results.
4. If you are using existing assets (e.g., code, data, models) or curating/releasing new assets...
(a) If your work uses existing assets, did you cite the creators? [Yes] See Section A.7.
(b) Did you mention the license of the assets? [N/A](c) Did you include any new assets either in the supplemental material or as a URL? [N/A](d) Did you discuss whether and how consent was obtained from people whose data you’re using/curating? [N/A](e) Did you discuss whether the data you are using/curating contains personally identifiable information or offensive content? [N/A]
5. If you used crowdsourcing or conducted research with human subjects...
(a) Did you include the full text of instructions given to participants and screenshots, if applicable?
[N/A](b) Did you describe any potential participant risks, with links to Institutional Review Board (IRB)
approvals, if applicable? [N/A](c) Did you include the estimated hourly wage paid to participants and the total amount spent on
participant compensation? [N/A]7

Figure 4: Comparison of the histograms of the maximum (top, left), minimum (top, right) and mean(bottom) Max-Cut ratios over 50 states on 100 test Hamiltonians between AVA and the GW algorithm.

A

Appendix

A.1

Definitions of Optimization Problems

A.1.1

Sherrington-Kirkpatrick model

The Sherrington-Kirkpatrick [27] spin-glass model (SK) is defined as
1 X
H(σ) = − √Jij σi σj ,
N i>j

(5)

where the couplings Jij are sampled independently from a standard Gaussian.

A.1.2

Max-Cut

The maximum cut (MC) of an unweighted graph is defined asn1on X (1 − σ σ ) o
i j
=
max
σT L σ ,MC = max
2
4
σ∈{−1,1}N(6)(i,j)∈E

where spins σi,j are associated with the graph nodes and E is the set of all edges. The graph Laplacian
L = D − A is defined in term of the degree matrix D and the adjacency matrix A. Since D is diagonal, σ TDσis constant and we can drop this part and just minimise σ TAσ, which is equivalent to maximizing the MC as
defined in Eq. 6. Instead of minimizing σ TAσ we minimize N1 σ TAσ, so that hyperparameters are more easilytransferable to larger system sizes. In our experiments, we sample connected graphs with a density of 0.5 andmake sure that each node can be reached by any other node.

A.1.3

Relative Energy Difference

We define the relative energy difference of a model with parameters θ with respect to a set of NH HamiltoniansasNHNσ
1 X 1 X Hτ (σi ) − Hτ (στ,0 )∆Erel (θ) =,(7)NH τ =1 Nσ i=1
Hτ (στ,0 )where Nσ is the number of sampled states for each Hamiltonian, στ,0 is the true ground state of HamiltonianHτ , and σi ∼ pθ (σ|Hτ ).

A.1.4

Max-Cut Ratio

We define the Max-Cut ratio as
RMC (τ, σ) =8MCτ (σ),MCτ (στ,0 )(8)

where the relation between the approximation ratio (AR) and the Max-Cut ratio is given by
AR = inf
RMC (τ, σ) ,Eτ σ∼pAlgo(9)

where pAlgo refers to the probability distribution stemming from an algorithm from which the configurations σ
are sampled, and the infimum is taken over the set of all possible graphs.
In Fig. 4 we provide further comparisons of Max-Cut ratios between GW and AVA on MC-50. We compare themaximum (top, left), minimum (top, right) and mean (bottom) Max-Cut ratios of both algorithms. While in mostcases our method samples better states than the GW algorithm, the maximum Max-Cut ratio averaged over 100
Hamiltonians of GW of ≈ 0.996 is better than our method which achieves a value of ≈ 0.991.

A.2

Free Energy Minimisation

The Free Energy minimisation objective for one Hamiltonian Hτ is given byhiFτ =Hτ (σ) + T log(pθ (σ|Hτ )) .
E

(10)

The gradient for the weight update is then given byihHτ (σ) + T log(pθ (σ|Hτ )) − b ∇θ log(pθ (σ|Hτ )) ,∇θ Fτ =
E(11)σ∼pθ (σ|Hτ )σ∼pθ (σ|Hτ )where we use b = Eσ∼pθ (σ|Hτ ) (Hτ + T log(pθ (σ|Hτ ))) for variance reduction (see [8; 5] for details).
In order to obtain the amortised optimisation objective, we calculate the expectation of Eq. 11 over a batch of
Hamiltonians.

A.3

Amortisation

Beyond a certain number of problem instances, denoted by namort , it pays off to use a single AVA model. We
AVA
estimate this number by dividing the number of annealing steps Nanneal, which AVA needs to obtain a certain
VNA∆Erel , and divide it by the number of annealing steps Nanneal
that VNA takes in order to obtain the same
relative energy difference.
The equation for namort is therefore given bynamort (∆Erel ) =AVA
Nanneal
.
VNA
Nanneal
∆Erel

(12)

In order to estimate this amortization number of Hamiltonians namort , we fit the relative energy differencek
as a function of annealing steps as ∆Erel = ∆E0 /Nanneal. We fit ∆E0 and k for both methods on results(VNA)
for SK-16 and obtain fit parameters of ∆E0= 10.72 ± 3.28, kVNA = 1.45 ± 0.12 for VNA, and
(AVA)
∆E0
= 1.29 ± 0.05, kAVA = 0.45 ± 0.01 for AVA. For AVA we have used the results on 100 Hamiltonians
and for VNA we have evaluated the method on 36 Hamiltonians. Since the decay rate of VNA is larger than for
AVA, VNA will achieve a lower ∆Erel for the same number of annealing steps.
The amortization number for a specific relative energy difference can then be calculated as−k 1 +k 1(AVA) k 1
(VNA) − k 1) AVA (∆E0
) VNA ∆Erel AVA VNA .

namort (∆Erel ) = (∆E0

(13)

In Fig. 3 (left), the shaded region indicates the errors from the fits, using the following upper and lower estimates
(AVA)
(VNA)
of the fit parameters. The parameters of the lower bound are estimated by ∆E0= 1.29+0.05, ∆E0=10.72 − 3.28, kAVA = 0.45 − 0.01, kVNA = 1.45 + 0.12 and the parameters of the upper bound are
(AVA)
(VNA)
∆E0
= 1.29 − 0.05, ∆E0= 10.72 + 3.28, kAVA = 0.45 + 0.01, kVNA = 1.45 − 0.12.

A.4
A.4.1

Architecture details
LSTM

In our SK experiments, we use three LSTM layers with 120 hidden neurons each. In our Max-Cut experiments
we additionally use dilated LSTMs, where the l-th (l ∈ {0, 1, 2}) LSTM layer at site i receives the hidden statesfrom the i−2l -th site. In order to obtain probabilities out of our model, we apply a Probability MLP (ProbMLP)
on the hidden state of the LSTM.

9

A.4.2Graph Embedding Network

To obtain a graph embedding of the Hamiltonian, we first pass its couplings Jij and possibly also external0fields Fi through an Encoder MLP, with which we calculate edge embeddings Eij= MLPEnc,J (Jij ) and node
embeddings Ni0 = MLPEnc,F (Fi ).
Afterwards, we use GNNs to apply several message passing updates, which can be written asl
Nil+1 = ln(Φ(Nil , □j∈N (i) Ψ(Nil , Njl , Eij)) + W Nil ).

(14)lHere, neighbouring nodes Nil , Njl and their correspoinding edges Eijat layer l are first processed by a message

MLP denoted as Φ. The resulting messages are then aggregated by a permutation invariant aggregation □j∈N (i) ,which is in our case the mean aggregator. Afterwards, the node update function Ψ is applied on the aggregatedmessages together with the current node embedding Nil . Finally, we use a skip connection, where a weightmatrix W is applied on Nil and apply layer normalization (ln) [28]. For simplicity, in this work we always
employ fully connected graphs where the edge weight denotes whether there is an edge or not.
Similarly, we update edges with
l+1
l
l
Eij
= ln(ΨE (Nil , Njl , Eij)) + WE Eij
),

(15)

where ΨE is the edge update MLP and WE is a weight matrix which is used for skip connections.
Our GNN architecture is similar to the Full GN block architecture introduced in [22] but without global graphattributes.
Afterwards, we apply our decoder MLP on the output nodes of the GNN NiL to obtain the node-wise Hamiltonian
embedding ei (Hτ ).
Our graph embedding network uses three-layer MLPs with hidden size of [128,96,64] (encoder), [64,64,64](processor) and [64, 96, 128] (decoder). The decoder MLP, encoder MLP, Φ and Ψ are always three-layer MLPs.
Within MLPs we always use Relu [29] activation functions and apply layer normalization after every layer. For
Experiments in Fig. 2 (left) and 3 (left) we determined the best learning rate with a grid search in order to obtain
good hyperparameters of each setting.

A.4.3

MLP Embedding Network

In order to obtain an MLP embedding, the couplings of the Ising Hamiltonian are concatenated to a vector with
fixed ordering. We then use an L-Layer MLP with hl × N Neurons, where N is the number of spins and hl are
the number of hidden neurons per spin in layer l. The output e(Hτ ) of the MLP is then reshaped into a matrix
eij of dimension hout × N , where the j-th column is used as the Hamiltonian embedding at node j.
For the result in Fig. 2 (left) we did a grid search over lr ∈ {10−3 , 5 · 10−4 , 10−4 }, hL ∈ {64, 80} ∀L and overL ∈ {3, 4, 5} in order to find the best hyperparameters for the MLP embedding. L = 4, h1,2,3,4 = 64 and a lrof 1 · 10−3 worked best in our experiments. The network with the MLP embedding network from Fig. 2 (left) iswith ≈ 3.8 · 106 parameters, approximately 3.4 times larger than the network with the graph embedding networkwith 10 GNN layers, which has ≈ 1.1 · 106 parameters.

A.5

Hyperparemeters

As an optimizer, we use the RAdam optimizer [30]. All other hyperparameters can be taken from Tab. 1.
In the following, we provide further details on the Learning Rate Schedule, Annealing Schedule and other
implementation choices.

A.5.1

Learning Rate Schedule

In all of our architectures we use a Cosine decay learning rate schedule with restarts [31]. During the warmup
phase we use a cosine schedule with restarts of frequency two, where lrstart is increased by a factor of two
and then decreased by the same factor. During Temperature annealing we use restarts with a frequency of 10
while decreasing the learning rate with a cosine schedule down to lrend . In our experiments we always choose
lrstart = 10−1 · lrend .

A.5.2

Gradient Clipping

In order to use LSTMs we observed that it is necessary to use gradient clipping at a gradient norm of 1.0.

10

Table 1: Table with Hyperparameters on different problem instances.
Problem Instancen GNN Layers
lrstartNwarmupNstatesNσ
FigureNH
TstartProbMLP

A.5.3

45 · 10−410450
16
2 (left)50
2.
[128,96,2]

6
5 · 10−410450
16
2 (left),50
2.
[128,96,2]

SK
8
5 · 10−410450
16
2 (left)50
2.
[128,96,2]

10
1 · 10−4
104
50
16, 212 (left), 2 (right), 3 (center)50
2.
[128,96,2]

MC
10
1 · 10−4
2 · 10350
50
3 (right), 450
0.2[128,2]Annealing Schedule

We adapt the temperature annealing schedule from [8], where first a warmup phase of Nwarmup is applied andthen the temperature is decreased linearily for Nanneal steps from T = Tstart to T = 0. Each step in Nwarmup
and Nanneal contains a cycle of Nequil = 5 gradient update steps.

A.5.4

VNA Hyperparameters

For VNA experiments for Fig. 3 (left) we have used the same architecture as for AVA and made a hyperparametersearch over the number of GNN layers L ∈ {3, 5, 10} and the lr ∈ {10−4 , 5 · 10−4 }. Here, the lr of 10−4 with
5 GNN layers performed best. In case of VNA, we run 2 · 103 warmup steps. All other hyperparameters werechosen as for SK-16 in Tab. 1.

A.5.5

Computational Resources

All Models for SK-16 can be trained on GTX 1080Ti (11GB) GPUs. For Max-Cut 50 distributed training on twoof such GPUs is necessary.

A.6

Evaluation Data

For small system sizes, we use exact diagonalization by reformulating the classical spin system into the quantum
setting. This procedure is justified as the (quantum) Ising Hamiltonian with no external fields,
X z
H = −Jijσi ⊗ σjz ,
(16)
i>j

is diagonal in the computational basis. Here, σ z =

1 0
0 −1



N

N

is the Pauli-z operator, and H ∈ C2 ×2 . The
Neigenvectors with the smallest eigenvalues are the ground states |Ψ0 ⟩ ∈ C2 and the corresponding eigenvalue
is the ground state energy. Since we work with randomly sampled real-valued couplings, there is typically only
one ground state solution |σ0 ⟩ and its spin-flipped counterpart |σ̄0 ⟩. These are products of computational basis
states, i.e., each individual spin is either up (|↑⟩ = (1, 0)T ) or down (|↓⟩ = (0, 1)T ), e.g., |σ0 ⟩ = |↑⟩|↓⟩ · · · |↓⟩.
There is a trivial equivalence between the quantum states |σ0 ⟩ and the classical vectors σ0 ∈ {−1, +1}N as thek-th spin in the product state (up or down) corresponds to the k-th entry of the classical spin vector (−1 or +1).
Any superposition c1 |σ0 ⟩ + c2 |σ̄0 ⟩ (with |c1 |2 + |c2 |2 = 1) of these ground states is a ground state as well. In
the main text and Fig. 3 (center), we use the ground state probability in brief notation p(σ0 ), which is computed
as p(|σ0 ⟩ + p(|σ̄0 ⟩).
For exact diagonalization, we make use of implementation in Netket [32]. For larger system sizes we use thespinglass server http://spinglass.uni-bonn.de/.

A.7

Packages

The code is written in Jax [33]. For the GW algorithm we use the package cvxgraphalgs https://pypi.org/
project/cvxgraphalgs/.

11