The Wishart planted ensemble: A tunably-rugged pairwise Ising model with a first-order phase transition

TL;DR

The paper introduces the Wishart planted ensemble, a tunably-hard Ising model with a first-order phase transition, useful for benchmarking optimization algorithms and understanding rugged energy landscapes.

Contribution

It presents a new class of Ising models with tunable hardness, analytical insights into their thermodynamics, and characterization of their phase transition and energy landscape.

Findings

Exhibits a first-order phase transition in temperature.

Displays a wide range of algorithmic hardness with a peak in solution time.

Analytical expressions for hardness peak and energy distribution are derived.

Abstract

We propose the Wishart planted ensemble, a class of zero-field Ising models with tunable algorithmic hardness and specifiable (or planted) ground state. The problem class arises from a simple procedure for generating a family of random integer programming problems with specific statistical symmetry properties, but turns out to have intimate connections to a sign-inverted variant of the Hopfield model. The Hamiltonian contains only 2-spin interactions, with the coupler matrix following a type of Wishart distribution. The class exhibits a classical first-order phase transition in temperature. For some parameter settings the model has a locally-stable paramagnetic state, a feature which correlates strongly with difficulty in finding the ground state and suggests an extremely rugged energy landscape. We analytically probe the ensemble thermodynamic properties by deriving the Thouless-Anderson-Palmer equations and free energy and corroborate the results with a replica and annealed approximation analysis; extensive Monte Carlo simulations confirm our predictions of the first-order transition temperature. The class exhibits a wide variation in algorithmic hardness as a generation parameter is varied, with a pronounced easy-hard-easy profile and peak in solution time towering many orders of magnitude over that of the easy regimes. By deriving the ensemble-averaged energy distribution and taking into account finite-precision representation, we propose an analytical expression for the location of the hardness peak and show that at fixed precision, the number of constraints in the integer program must increase with system size to yield truly hard problems. The Wishart planted ensemble is interesting for its peculiar physical properties and provides a useful and analytically-transparent set of problems for benchmarking optimization algorithms.

Tables

Table 1. Table 1: Parameters of the simulations. N sa subscript 𝑁 sa N_{\rm sa} is the number of samples, N sw = 2 b subscript 𝑁 sw superscript 2 𝑏 N_{\rm sw}=2^{b} is the total number of Monte Carlo sweeps for each of the N T subscript 𝑁 𝑇 N_{T} replicas for a single sample, T min subscript 𝑇 min T_{\rm min} [ T max subscript 𝑇 max T_{\rm max} ] is the lowest [highest] temperature simulated, and N T subscript 𝑁 𝑇 N_{T} is the number of temperatures used in the parallel tempering scheme for each system size N 𝑁 N and α 𝛼 \alpha .

$\alpha$	$N$	$N_{\mathrm{sa}}$	$b$	$T_{\min }$	$T_{\max }$	$N_T$
$0.50$	$48$	$1000$	$22$	$0.0650$	$1.40000$	$150$
$0.50$	$64$	$1000$	$22$	$0.0650$	$1.40000$	$150$
$0.50$	$96$	$1000$	$26$	$0.0650$	$1.40000$	$150$
$0.75$	$48$	$1000$	$21$	$0.1150$	$1.40000$	$130$
$0.75$	$64$	$1000$	$21$	$0.1150$	$1.40000$	$130$
$0.75$	$96$	$1000$	$21$	$0.1150$	$1.40000$	$130$
$0.75$	$128$	$1000$	$21$	$0.1150$	$1.40000$	$130$
$1.00$	$48$	$1000$	$21$	$0.1150$	$1.40000$	$130$
$1.00$	$64$	$1000$	$21$	$0.1150$	$1.40000$	$130$
$1.00$	$96$	$1000$	$21$	$0.1150$	$1.40000$	$130$
$1.00$	$128$	$1000$	$21$	$0.1150$	$1.40000$	$130$
$1.00$	$192$	$1000$	$21$	$0.1150$	$1.40000$	$130$

Equations

$$H(\mathbf{s})=-\frac{1}{2}\sum _{i\ne j}{J}_{ij}{s}_i{s}_j,$$

$$H(\pm \mathbf{t})=\underset{\mathbf{s}\in {\mathbb{S}}^N}{\min }H(\mathbf{s}).$$

$${\mathbf{W}}^T\mathbf{t}=0.$$

$$G(\mathbf{s})=\frac{1}{2}{\mathbf{s}}^T\mathbf{W}{\mathbf{W}}^T\mathbf{s}=\frac{1}{2}\Vert {\mathbf{W}}^T\mathbf{s}{\Vert }_2^2$$

$$\tilde{\mathbf{J}}=-\frac{1}{N}\mathbf{W}{\mathbf{W}}^T=-\frac{1}{N}\sum _{\mu =1}^M{\mathbf{w}}^{\mu }\otimes {\mathbf{w}}^{\mu }$$

$$\mathbf{J}=\tilde{\mathbf{J}}-\text{diag}(\tilde{\mathbf{J}}),$$

$$H(\mathbf{s})=-\frac{1}{2}{\mathbf{s}}^T\mathbf{J}\mathbf{s}$$

$$H(\mathbf{t})$$

$$=\frac{1}{2}\text{Tr}(\tilde{\mathbf{J}}),$$

$${\mathbf{w}}^{\mu }\sim \mathcal{N}(0,\mathbf{\Sigma }),$$

$${\boldsymbol{\Sigma}}=\frac{N}{N-1}\Big{[}{\boldsymbol{I}}_{N}-\frac{1}{N}{\boldsymbol{t}}{\boldsymbol{t}}^{T}\Big{]}$$

$$\mathbb{E}[w_i{w}_j]=-\frac{t_i{t}_j}{N-1}.$$

$${\boldsymbol{\Sigma}}^{\frac{1}{2}}=\sqrt{\frac{N}{N-1}}\Big{[}{\boldsymbol{I}}_{N}-\frac{1}{N}{\boldsymbol{t}}{\boldsymbol{t}}^{T}\Big{]}.$$

$$\text{solve }{\mathbf{W}}^T\mathbf{s}=0$$

$$\text{subject to }\mathbf{s}\in {\mathbb{S}}^N.$$

$$\tilde{\mathbf{J}}$$

$$=-\frac{\alpha }{M}\mathbf{W}{\mathbf{W}}^T$$

$$\to -\alpha \mathbf{\Sigma },$$

$$J_{ij}=\alpha \frac{t_i{t}_j}{N-1},$$

$${\mathbf{\Sigma }}^{\frac{1}{2}}=\frac{1}{\sqrt{N(N-1)}}\mathbf{A},$$

$${\boldsymbol{A}}\triangleq N\Big{[}I_{N}-\frac{1}{N}{\boldsymbol{t}}{\boldsymbol{t}}^{T}\Big{]},$$

$$A_{ij}=\left\{\begin{array}[]{ll}N-1&\textrm{$i=j$}\\ -1&\textrm{$i\neq j$}\end{array}\right..$$

$$\mathbf{w}=\frac{1}{\sqrt{N(N-1)}}\mathbf{A}\mathbf{z},$$

$${\mathcal{S}}_{\mathbf{w}}\triangleq \{-2(N-1),\dots ,-2,0,2,\dots ,2(N-1)\}.$$

$$J_{ij}=-\frac{1}{N}\sum _{\mu =1}^M{w}_i^{\mu }{w}_j^{\mu },$$

$$N^2(N-1)J_{ij}\in {\mathcal{S}}_{\mathbf{J}},$$

$${\mathcal{S}}_{\mathbf{J}}\triangleq \{-4M(N-1{)}^2,\dots ,-4,0,4,\dots 4M(N-1{)}^2\}.$$

$$J_{ij}=\frac{1}{N}\sum _{\mu =1}^p{\xi }_i^{\mu }{\xi }_j^{\mu },$$

$${\tilde{J}}_{ij}=-\frac{1}{N}\sum _{\mu =1}^M{w}_i^{\mu }{w}_j^{\mu }.$$

$$q=\frac{1}{N}\sum _{i=1}^N{m}_i^2$$

$$V=\frac{\alpha (1-q)}{1+\beta (1-q)}.$$

$$m_{i}=\tanh\Big{[}\beta\Big{(}\sum_{i\neq j}J_{ij}m_{j}-\beta Vm_{i}\Big{)}\Big{]}.$$

Full text

I Introduction

The interface between physics and computational complexity has yielded fruitful insights over decades of research. Hard optimization problems—which are ubiquitous throughout the natural sciences and domains such as operations research—are of significant importance to humans and are widely believed to admit no efficient algorithms for their solution over all members of their class. It was recognized that such problems show analogous features to those found in statistical mechanical systems, for example the existence of algorithmic phase transitions Hogg et al. (1996); Monasson et al. (1999), under which typical problems show a dramatic increase in the difficulty faced by known exact and heuristic algorithms.

In some cases, insights from the physics of spin glasses and disordered systems have inspired remarkable new algorithms; for example Mézard, Parisi, and Zecchina Mézard et al. (2002) studied algorithmic hardness transitions in random 3-satisfiability (3-SAT) problems using tools from statistical physics and subsequently proposed survey propagation as a promising method for solving such problems. In addition, physics-based approaches have suggested ensembles of very hard problems; for example locked constraint satisfaction problems Zdeborová and Mézard (2008) owe their difficulty to the fragmentation of the solution space into widely separated sets. The NK model Kauffman and Levin (1987); Kauffman and Weinberger (1989) is a well-studied class of tunably rugged cost functions proposed to capture the complexity of a variety of physical and biological systems.

An important special class of hard problem ensembles are those whose solutions are known to the constructor; these are often known as optimization problems with planted solutions. Aside from their theoretical interest, such problems are noteworthy for several reasons. They may, for example, serve as candidates for cryptographic one-way functions, that is functions whose outputs are cheaply computable for any input but for which determining an input yielding a given output is hard. Furthermore, they serve as useful benchmark problems for evaluating heuristic or exact algorithms. In recent years, the need for such benchmarks has increased with the advent of physical devices implementing quantum annealing Johnson et al. (2011) and related (e.g., Ref. Wang et al. (2013)) algorithms. In such situations, it is desirable to not only have access to a set of problems of tunable hardness but to also be able to compare an algorithm’s performance with the correct answer.

A physics-based approach for generating hard 3-SAT problems with planted solutions was proposed by Barthel et al. Barthel et al. (2002); more recently, Krzakala and Zdeborová presented a technique known as quiet planting Krzakala and Zdeborová (2009) for devising graph $q$-coloring problems with known solutions whose properties are indistinguishable from those of a random ensemble. This concept can be generalized to a variety of sparse problems, and has a close connection to reconstruction on trees Ricci-Tersenghi et al. (2019).

While the aforementioned techniques and analyses have yielded numerous elegant insights, they all share the undesirable property of considering problems that are structurally far-removed from contemporary physics-based optimization devices. $K$-satisfiability problems, for example, require energy functions to include terms whose value depends on groups of $K$ variables; in all the previously mentioned work over binary variables, $K\geq 3$. Realistic spin system models and optimization hardware on the other hand are typically restricted to pairwise interactions; emulating higher-order interactions on such systems can require tremendous overhead. In contrast, $q$-coloring problems directly correspond to the antiferromagetic Potts model of statistical physics and are thus expressible in terms of $2$-body interactions, but $q$ must be larger than $2$ to yield hard problems as $2$-coloring instances can be solved in linear time. On devices natively encoding problems consisting of binary variables, this can be problematic. Various techniques have recently been published Hen et al. (2015); Wang et al. (2017); Hamze et al. (2018); Hen (2019) for constructing planted Ising instances on sparse graph topologies, which in some cases appear to yield quite difficult problems Hamze et al. (2018), but in common with short-ranged disordered models in general, most of their known properties are inferred from numerical simulation and much remains unexplored about which features make them amenable as benchmarks for given algorithms.

In this paper, we propose a simple randomized procedure for generating systems of binary-constrained integer programs with a known solution Pattison et al. (2018). The system coefficients are generated according to a specific type of correlated multivariate Gaussian distribution. When translated to an effective Ising Hamiltonian, we obtain a novel type of disordered system with known ground state which we call the Wishart planted ensemble; the name is inspired by the distribution followed by the resultant random matrix of couplers.

The Hamiltonian includes interactions among all pairs of variables; this aspect makes it a less-than-perfect fit for testing on devices which implement short-range topologies. Unlike hard problem ensembles based on $K$-SAT or graph coloring however, the variable domains are binary valued and the interactions are pairwise. Vitally, the ensemble displays an rich array of thermodynamic and computational properties of considerable relevance to both classical and quantum algorithms. Computationally, the problems can be tuned to range in difficulty from very easy to extremely difficult at quite modest system sizes. We emphasize that nearly all statements about problem “hardness” in this paper refer to the empirically observed typical-case difficulty encountered by heuristics (and in all likelihood exact algorithms) and have no bearing on theoretical computer science questions concerned with worst-case difficulty.

A striking physical property of the Wishart ensemble is the existence of a first-order phase transition, a discontinuous jump in the free energy derivative at some system-dependent critical temperature. A key parameter in the generation procedure is $\alpha$, which specifies the number-of-equations–to–number-of-variables ratio; $\alpha$ is analogous to parameters such as the clause-to-variable ratio in the satisfiability problems and exerts critical influence on the physical and algorithmic complexity. By deriving the Thouless-Anderson-Palmer Thouless et al. (1977) equations for the ensemble with special care to account for the correlated couplers, we obtain the mean-field free energy, from which we find that at any finite $\alpha$, the internal energy drops discontinuously from some excited value at some $\alpha$-dependent critical temperature $T_{c}$. We verify this temperature and the nature of the transition with extensive parallel tempering Monte Carlo Geyer (1991); Hukushima and Nemoto (1996) simulations, which reveal that the system converges to its asymptotic predicted properties quite rapidly. For large $\alpha>1$, the magnitude of the discontinuity decreases monotonically and the thermodynamics smoothly change character toward a traditional second-order ferromagnetic transition. When $\alpha<1$, the paramagnetic state, i.e., the set of all configurations having no correlation with the planted solution, is stable at any nonzero temperature. This feature signals difficulty for classical heuristic algorithms as it behaves as a deceptive dynamical “trap.” More specifically, following free-energy gradients as done by methods like simulated annealing Kirkpatrick et al. (1983) will overwhelmingly lead to solutions far from the true optimum. Only by fortuitous initialization within the ground-state basin of attraction will the problem be solved with high probability. It turns out that $\alpha$ modulates the size of the ground-state basin, with larger values increasing the probability of solution by lucky initialization. When $\alpha\geq 1$, on the other hand, the paramagnetic state becomes unstable for some temperature $T_{u}<T_{c}$; at that point local algorithms can successfully find the solution by “rolling downhill” and hence such problems are typically easy. Remarkably, Barthel et al. Barthel et al. (2002) also argue that the hardness of their hard 3-SAT planted ensemble is predicated on the existence of a first-order ferromagnetic transition. We confirm the results of the Thouless-Anderson-Palmer (TAP) analysis with two alternative approaches: the replica method Edwards and Anderson (1975) and the annealed approximation.

First-order phase transitions are well known to exist in the $q$-state Potts model for $q\geq 3$ and in some Ising-like systems such as the Blume-Capel model Blume (1966); Capel (1966) in which variables nonetheless assume more than two states. They are, however, quite unusual in the pairwise Ising model in zero field (though see Ref. Yedidia and Georges (1990)); in particular neither the Sherrington-Kirkpatrick Sherrington and Kirkpatrick (1975) nor the Edwards-Anderson Edwards and Anderson (1975) spin glasses exhibit such a transition. The Hopfield model Hopfield (1982) displays a first-order transition between the spin glass and retrieval phases when a relatively small number of patterns ($\alpha<0.05$) are stored Amit et al. (1987), though as discussed in Sec. III, the Hopfield model is less appropriate as a class of problems with controllable hardness than the Wishart planted ensemble. The puzzling presence of such a transition in our system is accounted for by noting that the couplers are correlated to enforce the existence of the planted solution rather than independently disordered. Simulating systems with a first-order transition is widely known to be challenging.

The Wishart planted ensemble is of particular interest because it shares several features with models that have been shown by Nishimori and Takada Nishimori and Takada (2017) to be promising candidates for exhibiting (limited) exponential speedup when simulated using so-called nonstoquastic quantum driver Hamitonians; such systems cannot be simulated classically and hence represent a “strong” type of quantum effect. The advantage of the Wishart planted ensemble over the $p$-spin models considered in Ref. Nishimori and Takada (2017) is once again that the interactions are naturally pairwise rather than requiring terms of order $p\geq 3$.

Suitable nonstoquastic devices are not available as of the writing of this paper, but the presented model can be used to explore interesting problems on near-term stoquastic devices as well. In particular, the combination of a rough multimodal landscape (replica symmetry breaking) in the space orthogonal to the planted solution alongside tunable control of the energy of the planted solution (at leading order in $N$) opens up the possibility to explore the hard population transfer problem recently proposed by Smelyanskiy et al. Smelyanskiy et al. (2018). This model is useful for distinguishing physical quantum dynamics from classical dynamics such as quantum Monte Carlo in transverse field Ising models due to the multipath tunneling phenomena (miniband resonance). In principle one can prepare a state in a planted mode, tune the mode energy to equality (resonance) with other modes, and explore the rate of escape.

While this paper is mostly concerned with classical properties, these intriguing connections to both types of quantum devices will certainly be explored in future work.

Computationally, the Wishart planted ensemble emerges from our procedure for generating a certain type of random integer linear program Papadimitriou and Steiglitz (1998) and has numerous connections with well-studied Fu and Anderson (1989); Mertens (2001); Borgs et al. (2001); Karmarkar et al. (1986); Lagarias and Odlyzko (1985) optimization problems such as the number partitioning and subset sum Garey and Johnson (1979) problems. In common to these problems, the allowable precision over the problem parameters has important influence over combinatorial properties. There is also, however a crucial distinction from these problems; when the parameter $\alpha$ is fixed, the number of equations in the integer program scales linearly with the number of variables rather than remain constant (at unity, in the case of the previous two problems).

Random problem ensembles without a planted solution typically display a parametrized “easy-hard” difficulty transition (e.g., Refs. Gent and Walsh (1996) andMertens (2001)) in their optimization variants. The Wishart ensemble, on the other hand, shows an “easy-hard-easy” character. One of the easy regimes is due to the presence of a very large number of acceptable solutions coexisting with the planted ground state; given the task of locating any one of them, an optimization method has a relatively high likelihood of success. The other is due to the planting procedure effectively constraining the search space, providing “hints” to the algorithm toward the solution. The hard regime, however, is seen to be extremely difficult: Numerical experiments using a distributed, state-of-the art parallel tempering implementation show a dramatic hardness peak for small ($N=32$) system sizes. When $N=64$, parallel tempering Monte Carlo fails to locate even an approximate solution under lax and permissive target criteria within the allotted simulation time of around 11 h on contemporary high-speed hardware.

While we derive the Wishart planted ensemble in terms of a somewhat abstract random integer program, the resultant model turns out to have a remarkable structural similarity to the Hopfield model Hopfield (1982); Amit et al. (1987) of biological neurons and more particularly, a sign-inverted variant Nokura (1998) proposed to model neural “unlearning”. As discussed in Sec. III, there turn out to be several important differences between the models; nonetheless, it is exciting that such completely different starting points as integer programming and unlearning in neural networks result in models with close connections.

The rest of the paper is structured as follows. In Sec. II we describe our procedure to generate Wishart planted instances based on random integer programming problems. We also discuss the ensemble’s computational properties and how to represent its members, whose parameters are defined to take continuous values, with finite precision. Section III analyzes the physical properties of the class; the TAP free energy (derived in the Appendix) is analyzed and shown to have global minima along a one-dimensional subspace of the set of spin magnetizations. Furthermore, it has a locally stable paramagnetic state for all temperatures when $\alpha<1$ and for $T>\alpha-1$ when $\alpha\geq 1$. These properties give rise to the first-order transition between the paramagnetic and planted states; we determine the transition temperature in terms of $\alpha$. We show that as $\alpha$ grows, the system begins to increasingly behave like a ferromagnet, i.e., with a second-order transition. The predicted first-order transition temperature is validated with extensive Monte Carlo simulation. In Sec. IV, we turn our attention to empirical algorithmic properties; under finite-precision representation, the ensemble displays an easy-hard-easy relation with respect to parallel tempering time to solution as $\alpha$ is varied. After showing that the ensemble-averaged energies of the Wishart planted ensemble follow a gamma distribution and introducing the notion of an intrinsic search space, we analytically predict the location of the hardness peak for any target energy threshold and confirm the prediction using optimized parallel tempering simulations. The prediction is shown to be precisely accurate even for approximate solution criteria. We show that generating difficult problems under constant precision restriction requires scaling the number of constraints in the integer program approximately linearly. The Appendices VII contain most of our calculations, and confirm the TAP results through replica analysis and an annealed approximation.

II The Wishart planted ensemble

II.1 Generation procedure

Our goal is to construct an ensemble of zero-field Ising Hamiltonians over the $N$-spin complete graph with planted ground state ${\boldsymbol{t}}$, in other words, having the form

[TABLE]

where ${\boldsymbol{s}}$ and ${\boldsymbol{t}}$ refer to configurations on the $N$-spin Ising model configuration space $\mathbb{S}^{N}\triangleq\{\pm 1\}^{N}$ and such that

[TABLE]

When not explicitly stated ${\boldsymbol{t}}$ will be taken to be the ferromagnetic ground state ${\boldsymbol{t}}=(+1,+1,\ldots,+1)$ and its $\mathbb{Z}_{2}$ image, as the minimizer to such a problem can be subsequently concealed by gauge randomization.

Consider the $N\times M$ real-valued matrix ${\boldsymbol{W}}\in\mathbb{R}^{N\times M}$, whose $M\geq 1$ columns are denoted by ${\boldsymbol{w}}^{\mu}$ for $\mu=1,\ldots,M$. The value of $M$ turns out to modulate the ensemble properties such as thermodynamics and hardness; in this work we are primarily concerned with the regime in which $M$ scales linearly as a constant factor of $N$, i.e., $M=\alpha N$ for $\alpha>0$.

Given a desired ground state ${\boldsymbol{t}}$, our procedure seeks to construct a consistent homogeneous Ising-constrained linear system with ${\boldsymbol{s}}=\pm{\boldsymbol{t}}$ as a solution, in other words, to obtain ${\boldsymbol{W}}$ such that

[TABLE]

This is because the positive semidefinite quadratic form

[TABLE]

would then attain its minimum value of zero at ${\boldsymbol{s}}={\boldsymbol{t}}$, and hence if we define

[TABLE]

and zero its diagonal to form

[TABLE]

then the Hamiltonian

[TABLE]

attains its ground state at ${\boldsymbol{s}}={\boldsymbol{t}}$ with energy

[TABLE]

where the property that $s_{i}^{2}=1$ for ${\boldsymbol{s}}\in\mathbb{S}^{N}$ has been used. The scaling by $1/N$ in the definition of ${\boldsymbol{\tilde{J}}}$ is to make the energy extensive, i.e., scaling linearly with system size.

We obtain the linear system by individually generating the $M$ columns $\{{\boldsymbol{w}}^{\mu}\},i\in{1,\dots,M}$ of ${\boldsymbol{W}}$ such that $\langle{{\boldsymbol{w}}^{\mu}},{{\boldsymbol{t}}}\rangle=0$. We propose a simple projective method for efficiently generating correlated Gaussian variates satisfying the summation and other desirable properties. More precisely, the column vectors are set to be distributed as

[TABLE]

where the covariance matrix is given by

[TABLE]

with ${\boldsymbol{I}}_{N}$ the $N$-dimensional identity matrix. In other words, for each column vector ${\boldsymbol{w}}$, all elements have unit variance, and for all variable pairs $i\neq j$, the covariances are

[TABLE]

Note that $\textrm{rank}({\boldsymbol{\Sigma}})=N-1$ as expected; given any $N-1$ components of ${\boldsymbol{w}}$, the remaining one follows deterministically. To generate the column vectors, we first determine the square root of ${\boldsymbol{\Sigma}}$, i.e., ${\boldsymbol{\Sigma}}^{\frac{1}{2}}$ such that ${\boldsymbol{\Sigma}}={\boldsymbol{\Sigma}}^{\frac{1}{2}}{\boldsymbol{\Sigma}}^{\frac{1}{2}}$, to be

[TABLE]

We then iterate over the loop described in Algorithm 1:

One can readily verify that $\langle{{\boldsymbol{w}}^{\mu}},{{\boldsymbol{t}}}\rangle=0$ for all $\mu$ and that ${\boldsymbol{w}}^{\mu}$ is distributed according to Eq. (6). While the components of ${\boldsymbol{w}}$ are correlated, the Gaussian nonetheless has strong structure that simplifies the subsequent analysis. Following an appropriate gauge transformation to the ferromagnetic state, the elements of ${\boldsymbol{\Sigma}}$ imply that the distribution is exchangeable, i.e., invariant to a permutation of the components. Exchangeability is a stronger property than stationarity, which merely requires the covariance of components $i$ and $j$ to depend only on $|i-j|$. The property is used in Sec. III.1 when deriving the TAP free energy and again in Sec. IV.1 when obtaining the ensemble energy distribution.

The random matrix ${\boldsymbol{W}}{\boldsymbol{W}}^{T}$ follows a Wishart distribution, a well-studied matrix generalization of the $\chi^{2}$ distribution; in light of this we call our problem class the Wishart planted ensemble (WPE). When $M<N$, the support of the Wishart density lies on a low-dimensional subspace of $N\times N$ matrices Uhlig (1994). As we see in Sec. III.2, its spectral distribution is a key feature underlying the phase behavior. We note also that while sampling ${\boldsymbol{w}}^{\mu}$ according to the Gaussian $\mathcal{N}(\mathbf{0},{\boldsymbol{\Sigma}})$ simplifies the analysis it is by no means necessary in practice for our results to hold in the large $N$ limit. For example, if $\{{\boldsymbol{z}}^{\mu}\}$ used to obtain $\{{\boldsymbol{w}}^{\mu}\}$ were vectors of length $N$ consisting of independent and uniform $\{\pm 1\}$ variates rather than uncorrelated Gaussians, then central limit arguments show that ${\boldsymbol{W}}{\boldsymbol{W}}^{T}$ is nonetheless asymptotically Wishart. We note however, that using finite precision introduces the possibility that states other than ${\boldsymbol{t}}$ may attain the ground-state energy. The optimization problem in this paper is defined to be that of locating any ground state; the algorithmic implications of the presence of several solutions is discussed in the next section and in detail in Sec. IV.

II.2 Computational properties

Before proceeding to an examination of the WPE thermodyamics, we discuss its properties from a computational perspective in light of its interpretation as a constrained homogeneous linear system. Readers familiar with related settings such as linear error correcting codes should bear in mind that arithmetic here is over the real numbers rather than a finite field such as ${\rm GF}(2)$.

Given a matrix ${\boldsymbol{W}}$, the task of finding the ground state of Eq. (5) is equivalent to finding a solution to the following NP-hard problem called integer programming feasibility

[TABLE]

Suspending for a moment the fact that $M$ scales with $N$ in the WPE, we can obtain a sense for how it may impact problem difficulty. If ${\boldsymbol{W}}$ consists of $M<N$ independent columns, then $\textrm{dim}(\textrm{null}({\boldsymbol{W}}^{T}))=N-M$. The dimension of the nullspace of ${\boldsymbol{W}}^{T}$ implies the search space for potential Ising state solutions to Eq. (9), and so the larger it is, the more difficult the problem may be guessed to be. In particular, when ${\boldsymbol{W}}^{T}$ consists of a single row ($M=1$) and hence an $N-1$ dimensional nullspace, the problem may be surmised to be maximally hard. When ${\boldsymbol{W}}$ is specified with relatively low precision, this turns out to not be the case; exponentially many solutions other than $\pm{\boldsymbol{t}}$ overwhelmingly appear as $N$ increases and so locating any such satisfying ${\boldsymbol{s}}$ can be quite easy. This is reflected in the Ising Hamiltonian (5); when $M=1$, ${\boldsymbol{J}}$ can be verified to be fully frustrated com (a), which gives rise to tremendous low-energy degeneracy.

At the other extreme in which $M$ assumes large values, the nullspace of ${\boldsymbol{W}}^{T}$ becomes one-dimensional and hence the two Ising solutions to Eq. (9) are trivially recovered by inspection from a vector spanning the nullspace. In the Hamiltonian in Eq. (5), making $M$ large results in a ferromagnetic system. To see this, we note that

[TABLE]

where the limit follows from the law of large numbers. From the covariance matrix ${\boldsymbol{\Sigma}}$, the couplers thus uniformly approach

[TABLE]

implying that the system reduces to a gauge-transformed and rescaled Curie-Weiss ferromagnet, whose ground state is easy to find due to the lack of frustration. Hence one may reasonably guess that for a given allowable precision, the most difficult problems occur for some intermediate value of $M$. This easy-hard-easy profile is shown to indeed hold and is discussed in Sec. IV, in which we make more explicit the role of $\textrm{null}({\boldsymbol{W}}^{T})$ and develop the conjecture that the difficulty peak occurs at the value of $M$ in which the fewest number of solutions occur relative to the volume of the nullspace.

A prototypical special case of Eq. (9), corresponding to $M=1$, is the subset sum problem; its decision variant, that of establishing whether a satisfying ${\boldsymbol{s}}$ to Eq. (9) exists, was one of Karp’s Karp (1972) original 21 NP-complete problems. We note that ${\boldsymbol{W}}$ must be specifiable to arbitrary precision for the optimization variant in Eq. (9) to avoid efficient solution via the technique of dynamic programming as the complexity of this algorithm scales exponentially in the number of bits needed to specify the coefficients. This requirement holds more generally: when the value of $M$ is fixed to any integer and the matrix elements belong to a finite set, integer programming problems of the form of Eq. (9) can be solved in pseudopolynomial time using the same technique Papadimitriou (1981).

Subset sum problems with feasible solutions have been of particular interest due to their complexity underpinning the security of an early family of cryptographic systems Merkle and Hellman (1978). Remarkably, a method has been devised Lagarias and Odlyzko (1985) for solving with high probability a family of low-density random subset sum problems, in which the maximum element of the single-column ${\boldsymbol{W}}$ is large compared to $N$, based on an idea known as lattice basis reduction Lenstra et al. (1982).

A well-studied further specialization of the subset sum problem is known as the number partitioning problem, corresponding to the task of partitioning a base set of positive integers into two blocks with sums of minimum absolute difference (or “discrepancy”). When the integers are chosen independently from a uniform distribution (the problem does not respect a planted solution), number partitioning displays several interesting properties, in particular an algorithmic easy-hard phase transition Gent and Walsh (1996); Borgs et al. (2001); Mertens (2001, 2006). More specifically, if the positive integers forming the base set are bounded by $2^{\kappa N}$ for a fixed $\kappa>0$, then when $\kappa>\kappa_{c}$, partitions with a discrepancy of zero exist with vanishing probability, and typical instances become difficult for known heuristic and complete algorithms. Conversely, when $\kappa$ lies below this threshold, partitions with zero discrepancy are abundant, and the problems are easily solved.

While the subset sum and number partitioning problems bear obvious connections to the model presented in this work, it is apparent that our assumptions and focus are different. Broadly speaking, in the WPE a state must be found which now simultaneously satisfies the $M$ relations in the integer program; the fact that $M$ scales as $\alpha N$ rather than remaining fixed introduces consequences in the algorithmic hardness properties.

II.3 Representing a WPE instance

The WPE construction presented in Sec. II.1 was defined in terms of Gaussian variables; in a practical implementation, however, one must contend with finite precision and generally cannot represent the required continuous-valued parameters. In the WPE, this can be dealt with in one of two ways. The first, which will be examined in this section, is to maintain the correlation structure defined by ${\boldsymbol{\Sigma}}$ but replace the Gaussian variates ${\boldsymbol{z}}$ with a discrete zero-mean, unit variance ensemble; the simplest such choice is to have ${\boldsymbol{z}}$ uniformly and independently take the two values $\{\pm 1\}$ (sometimes called a Rademacher distribution). This turns out to allow an exact representation of the problem parameters using a logarithmic (in $N$) number of bits; further, as mentioned in Sec. II.1, the coupler matrix is nonetheless asymptotically Wishart and the same physical properties derived in this paper result.

The second and more heuristic way to represent a problem is to simply round the parameters to the closest machine-representable number. This introduces numerical errors; for example, the planted solution may no longer have its theoretically intended energy. Consequently, one must introduce a tolerance on what defines a “solution.” We use this approach in our algorithmic hardness experiments because as discussed in Sec. IV.3, it minimizes the potential discrepancy between the finite-size statistical properties and our analytical predictions, which may arise due to using non-Gaussian generator variables on the small systems we were forced to consider. Nonetheless, using the theoretically-derived energy histogram, we are able to account for the observed hardness peak under this approximate representation.

We now discuss the exact discretized WPE over a Rademacher distribution, i.e., where $\{z_{i}^{\mu}\}$ are independently and uniformly in $\{-1,1\}$ rather than drawn from $\mathcal{N}(0,1)$ for $i\in\{1,\ldots,N\}$ and $\mu\in\{1,\ldots,M\}$. We first rewrite

[TABLE]

where the integer-valued matrix

[TABLE]

i.e.,

[TABLE]

Because

[TABLE]

then up to the leading constant of $1/\sqrt{N(N-1)}$, the elements of ${\boldsymbol{w}}$ may assume values in the set of $2N-1$ equally spaced integers

[TABLE]

Thus, in the integer programming formulation, the Rademacher-discretized WPE takes approximately $\log(2N-1)$ bits to encode the parameters.

Obtaining the required precision in the Hamiltonian formulation (i.e., on ${\boldsymbol{J}}$) is more messy but analogous. Recalling that

[TABLE]

then using the previous restriction on the values of ${\boldsymbol{w}}$ we can show that

[TABLE]

where

[TABLE]

Not all elements in this set of integers are actually attainable by $J_{ij}$, but it provides a useful upper bound on the number of possible values and shows that representing $J_{ij}$ in the Rademacher-discretized WPE takes no more than on the order of $\log(M)+2\log(N)$ bits, which when $M=\alpha N$ is $O(\log N)$.

III Thermodynamic properties

Analyzing long-range disordered systems has a rich history in statistical mechanics Nishimori (2001). The Sherrington-Kirkpatrick (SK) model Sherrington and Kirkpatrick (1975), a fully connected Ising model with independently sampled Gaussian bond strengths, is a prototypical example of such systems. The replica method Edwards and Anderson (1975) is a powerful framework for performing such analyses, and has yielded great successes such as the solution to the SK model Parisi (1979) which have since been rigorously proved to be correct (see, for example, Ref. Panchenko (2013) and the references therein). This approach is pursued in Appendix VII.5 were we are able to recover the transition properties developed in this section, and through connections with the anti-Hopfield model identify some additional interesting transitions within the model.

A different and in many senses complementary approach to the analysis of weakly coupled, fully connected disordered systems is due to Thouless, Anderson, and Palmer Thouless et al. (1977). The TAP equations are a set of nonlinear relations satisfied by the local magnetizations for a given instance. They can be arrived at in one of several ways Opper and Winther (2001). For the SK model, they are often interpreted as correcting the “naïve” mean field equations with a so-called Onsager reaction term. The approach of Plefka Plefka (1982) arrives at the TAP equations by second-order expansion of the free energy at constant magnetization, which turns out to have an appealing information geometric interpretation Tanaka (2000). In this section we use another approach called the cavity method Mézard et al. (1987).

Determination of the TAP equations for the WPE is somewhat complicated by the fact that $J_{ij}$ are not independent variates as they are for the SK model. The TAP equations for systems with correlated ${\boldsymbol{J}}$ have been determined in the past: A notable example, which turns out to have a close connection to our ensemble, is the Hopfield model Hopfield (1982) of a biological neural network. The couplings of the Hopfield model are given by

[TABLE]

where the $p$ vectors $\{{\boldsymbol{\xi}}^{\mu}\}$, known as patterns to be stored for later retrieval, consist of $N$ independent zero-mean binary random variables. The physics of the Hopfield model was first studied with the replica approach by Amit et al. Amit et al. (1985, 1987); the TAP equations were derived by Mézard et al. Mézard et al. (1987) using the cavity method, but yielded results inconsistent with the replica analysis. TAP equations consistent with Ref. Amit et al. (1985) were obtained by Nakanishi and Takayama Nakanishi and Takayama (1997) using Plefka’s method, and subsequently by Shamir and Sompolinski Shamir and Sompolinsky (2000) via an elegant cavity approach.

The connection of the WPE to the Hopfield model is apparent if we express the elements of ${\boldsymbol{\tilde{J}}}$ defined in Eq. (3) as

[TABLE]

An obvious difference from the Hopfield model is the presence of the leading negation. Consequently, while the Hopfield Hamiltonian tends to favor spin configurations aligned with the patterns $\{{\boldsymbol{\xi}}^{\mu}\}$, the WPE Hamiltonian penalizes configurations overlapping with the directions $\{{\boldsymbol{w}}^{\mu}\}$. An additional distinction, however, is that while the $M$ vectors are independently drawn from the previously defined Gaussian distribution, the components of each vector are now correlated, a property which emerged due to the solution planting procedure. Remarkably, a model defined by negating the sign of the Hopfield model Hamiltonian (only the first of the two differences above) has been proposed as a model of “unlearning” paramagnetic configurations and thereby enhancing learning for biological networks. A replica analysis of this “anti-Hopfield” model was undertaken by Nokura Nokura (1998). The additional “layer” of correlations and the presence of a planted solution in the WPE turns out to lead to very different behavior from that of the anti-Hopfield model; in particular, the anti-Hopfield model has no first-order transition comparable to the transition into the planted solution, and the large-$\alpha$ regime is Sherrington-Kirkpatrick-like rather than ferromagnetic.

Another Ising ensemble related to the WPE is the random orthogonal model (ROM) proposed by Parisi and Potters Parisi and Potters (1995), where the ${\boldsymbol{J}}$ matrices are generated by uniformly sampling an orthogonal matrix ${\boldsymbol{O}}$, forming diagonal matrix ${\boldsymbol{D}}$ whose elements $D_{ii}\in\{\pm 1\}$ (often additionally assumed to have a trace of zero), and setting ${\boldsymbol{J}}={\boldsymbol{O}}{\boldsymbol{D}}{\boldsymbol{O}}^{T}$. Having been devised with different aims, this model is also quite different from the WPE. First, no planting takes place, so the ground state is uncontrolled. Further, the ground-state energy is only known if an Ising-feasible state happens to be an eigenvector of ${\boldsymbol{J}}$, which become exceedingly unlikely for even moderately sized systems. Finally, the eigenvalue distributions of the ${\boldsymbol{J}}$ matrices are not the same, nor are the TAP equations presented in Ref. Parisi and Potters (1995). The systems thus have quite differing thermodynamic properties; in particular, a first-order transition analogous to the planting transition is absent in the ROM.

It is natural to wonder about the suitability of the original Hopfield model to the task at hand. Unfortunately, it is not possible to generate sufficiently rugged energy landscapes while maintaining control over the ground states using this model. Hence, while the Hopfield model is an appealing abstraction of an associative memory, it is not as well suited to usage as a tunable planted optimization objective. The Hopfield model does not exhibit, for any setting, a persistent metastable paramagnetic state analogous to that of the $\alpha<1$ WPE. In contrast to the WPE in which the planted ground state is undisturbed while a rough energy landscape is induced in the orthogonal subspace, ruggedness in the Hopfield model arises due to interference among the stored patterns, causing undesired spurious minima Amit et al. (1987). For values of the Hopfield $\alpha$, which specifies the ratio of number of stored patterns to number of variables, of less than around 0.05, a first-order transition between a spin glass and ”retrieval” phase takes place, but finding one of the embedded patterns is quite easy in this case (as it should be for the model’s intended purpose). At somewhat larger values of $\alpha$, the model continues to function as an associative memory, though the patterns are only assured of being local optima rather than ground states. When $\alpha>0.138$, only a second-order paramagnetic to spin glass transition takes place, and control over even the local minima is completely lost.

III.1 TAP equations

We derive the TAP equations for the WPE following the two-step cavity approach of Shamir and Sompolinski Shamir and Sompolinsky (2000), making appropriate adaptations to deal with the correlations among the components of each ${\boldsymbol{w}}$. The complete calculation is shown in Sec. VII.1. Alternatives to the TAP method based on the annealed approximation and replica method are considered in Appendices VII.4 and VII.5, yielding consistent results. In these appendices the connection with the anti-Hopfield model is explored, as is the notion of a planted solution with tuned energy that might permit other uses, such as the population transfer experiment discussed in the introduction.

Let $m_{i}=\langle{s_{i}}\rangle$ and ${\boldsymbol{m}}$ be the vector of all magnetizations. Define

[TABLE]

and

[TABLE]

The TAP equations for the WPE are

[TABLE]

Solutions to Eq. (13) are stationary points of the following TAP free energy:

[TABLE]

where the local entropy terms are

[TABLE]

To determine the quenched free energy one must consider a weighted sum over minima de Dominicis and Young (1983). In this section we present a more limited analysis that is appropriate when only a single minimum dominates the free energy. The result, which may be an upper bound to the free energy, is consistent with the replica symmetric analysis of Appendix VII.5, and in good agreement with numerical results of Sec. IV. We expect the approximation to be sufficient at higher temperature, and that in other cases the consequences of multiple fixed points may be small following the replica symmetry breaking analyses in related models as discussed in Appendix VII.5 Parisi and Potters (1995); Nokura (1998).

In general, determining the global minima of $F_{\rm TAP}$ is a difficult task as ${\boldsymbol{m}}$ lies in an $N$-dimensional space subject to the bound constraints

[TABLE]

In our case, however, the existence of the planted solution considerably simplifies things by assuring that the free-energy minimum necessarily occurs along the ground-state direction:

[TABLE]

To see this, consider the restriction of ${\boldsymbol{m}}$ to the spherical shell given by

[TABLE]

for some $q_{0}\in[0,1]$ or, alternatively,

[TABLE]

for arbitrary unit vector ${\boldsymbol{e}}$. The claim is that

[TABLE]

minimizes $F_{\rm TAP}$ on the shell; from Eq. (14), the term

[TABLE]

is minimized for ${\boldsymbol{m}}=\pm\sqrt{q_{0}}{\boldsymbol{t}}$ due to the planting procedure. Further, the term

[TABLE]

is minimized on the shell at the $2^{N}$ points at which the $\{m_{i}\}$ have the same magnitude, which holds when ${\boldsymbol{m}}=\pm\sqrt{q_{0}}{\boldsymbol{t}}$, and the final term in $F_{\rm TAP}$ is independent of ${\boldsymbol{m}}$ on the sphere. The global optimum, which is the minimum over all the shells’ minimizers, thus occurs along

[TABLE]

Recall that

[TABLE]

and that

[TABLE]

Using the fact that $S(mt_{i})=S(m)$ for $t_{i}\in\{\pm 1\}$ we define the one-dimensional TAP free energy as $F(m)\triangleq F_{\rm TAP}(m{\boldsymbol{t}})$, i.e.,

[TABLE]

Consider the difference

[TABLE]

present in Eq. (15). Using the higher-order moments of correlated Gaussians, one can show that despite the dependence among the components of ${\boldsymbol{w}}^{\mu}$, the term scaling $\frac{m^{2}}{2}$ is at most $O(\sqrt{M})$. Normalizing by $N$ to obtain the free energy per spin, we obtain

[TABLE]

In the thermodynamic limit, the first term vanishes; neglecting the constant term of $-\frac{1}{2}\alpha$, the global free-energy minima of Eq. (14) can thus be found as those of the function

[TABLE]

The mimina of $\tilde{f}$ must be stationary points of Eq. (17), which are attained for solutions of the equation

[TABLE]

III.2 Thermodynamic properties of the WPE

The WPE displays phase properties quite surprising in the Ising model. Rather than undergo an SK-like spin-glass transition typical of many disordered systems, the presence of a planted solution gives rise to a thermal first-order phase transition. This transition in the WPE occurs at an $\alpha$-dependent temperature $T_{c}$ such that two stable states contribute equally to the free energy. As expected, at high temperature the disordered paramagnetic state (${\boldsymbol{m}}=\mathbf{0}$) is the unique global minimizer to $F_{\rm TAP}$. At a critical temperature $T_{c}$, this state remains globally minimizing but is no longer unique; another state ${\boldsymbol{m}}=m{\boldsymbol{t}}$ with $m\neq 0$ is also optimal. Below $T_{c}$, the paramagnet ceases to be the minimum, but whether ${\boldsymbol{m}}=0$ remains a local minimum (metastable state) to $F_{\rm TAP}$ or an unstable state (saddle point or maximizer) depends on $\alpha$ and the temperature. We find the fascinating result that when $\alpha<1$, the paramagnetic state is locally stable at all nonzero temperatures. As will be discussed, this property has considerable algorithmic implications.

The transition is first-order because of the discontinuity in the minimizing ${\boldsymbol{m}}$ and in the derivative of the free-energy minimum with respect to $T$; consequently the internal energy also undergoes a jump at $T_{c}$ as we will see. Remarkably, for any finite $\alpha$, the transition is always technically first-order, in the sense that there is a discontinuity in the log partition function derivative. The magnitude of this discontinuity decreases with $\alpha$, and the transition gradually segues from first to second order.

III.2.1 Stability of ${\boldsymbol{m}}=\mathbf{0}$

The stability of ${\boldsymbol{m}}=\mathbf{0}$ can be ascertained from the Hessian of $F_{\rm TAP}$, i.e.,

[TABLE]

The point ${\boldsymbol{m}}=\mathbf{0}$ is a local minimum of $F_{\rm TAP}$ if and only if ${\boldsymbol{H}}(\mathbf{0})$ is positive definite. It turns out that the spectral distribution of the random matrix ${\boldsymbol{J}}$ is key in determining the relation between the definiteness of ${\boldsymbol{H}}(\mathbf{0})$ and $T$ at given $\alpha$. The limiting eigenvalue distribution of ${\boldsymbol{J}}$ is calculated in Sec. VII.2; somewhat surprisingly, it is closely related to the Marchenko-Pastur Marchenko and Pastur (1967) law for the spectrum of Wishart matrices constructed with uncorrelated Gaussians. The important aspect for now is that $-\alpha$ is always (asymptotically) the smallest eigenvalue of $-{\boldsymbol{J}}$.

Computing the partial derivatives from the definition of $F_{\rm TAP}$ in Eq. (14) we obtain the Hessian matrix elements

[TABLE]

At ${\boldsymbol{m}}=\mathbf{0}$ we thus have

[TABLE]

or compactly, the Hessian at ${\boldsymbol{m}}=0$

[TABLE]

with

[TABLE]

${\boldsymbol{H}}(\mathbf{0})$ is positive definite when $c(\alpha,\beta)$ is large enough to shift the eigenvalues of $-{\boldsymbol{J}}$ to all be positive. From the spectral property mentioned above, stability hence occurs when

[TABLE]

For a given $\alpha$, this condition is equivalent to

[TABLE]

or

[TABLE]

Obviously, $T$ is non-negative; when $\alpha<1$ the relation thus holds for all $T$ while when $\alpha\geq 1$, it is satisfied when $T\geq\alpha-1$.

To summarize, when $\alpha<1$, the point ${\boldsymbol{m}}=\mathbf{0}$ is a locally stable stationary point to the TAP free energy at any temperature; on the other hand, when $\alpha\geq 1$, the Hessian becomes indefinite at temperature $T_{u}=\alpha-1$. For temperatures below $T_{u}$, note that the eigenvalue of ${\boldsymbol{H}}$ corresponding to eigenvector ${\boldsymbol{t}}$ (the planted solution) becomes negative, which implies that the paramagnetic solution becomes unstable along this direction. Interestingly, the curvature of $F_{\rm TAP}$ along ${\boldsymbol{t}}$ is of largest magnitude for $\alpha\geq 1$ at the temperature such that $c(\alpha,\beta)$ is minimal. In terms of $T$,

[TABLE]

which is smallest when $T=-1+\sqrt{\alpha}$. As local search heuristics follow free-energy gradients in their quest for the solution, it seems reasonable that this temperature is in some sense optimal to start with when searching for the ground state in the $\alpha\geq 1$ (easy) regime. A further special interpretation of this temperature relates to the anti-Hopfield model, as discussed in Appendix VII.5.

III.2.2 First-order phase transition

The first-order transition for given $\alpha$ occurs at a $T_{c}$ such that $\tilde{f}(0)=\tilde{f}(m)$ for $m\neq 0$. For any $T$, the stationary points of $\tilde{f}$ are determined numerically by iteratively solving the saddle point equation [Eq. (18)]; using binary search in $T$, we can then localize the transition temperatures at any $\alpha$. Figure 1 shows the relation of $T_{c}$ with $\alpha$, and also shows for $\alpha\geq 1$ the paramagnetic instability point $T_{u}=\alpha-1$. It is clear that $T_{c}>T_{u}$ uniformly but as $\alpha$ increases the two temperatures converge; this in turn constrains the jump magnitude between the two sides of the transition. While $T_{c}$ has no closed-form expression in $\alpha$, we can obtain a lower bound $\widehat{T_{c}}$ as

[TABLE]

which serves as a surprisingly accurate approximation for small $\alpha$ but deteriorates as $\alpha$ gets larger than $1$.

We illustrate these results by plotting $\tilde{f}(m)$ at representative temperatures for a few values of $\alpha$. Consider first the relatively small value of $\alpha=0.25$, with landscape illustrated in Fig. 2. The first-order transition at $T_{c}=0.004$ is clearly visible in the middle panel as $m=0$ is now a coexisting global minimum with the other minima close to the end points com (b). Remarkably, when $T\approx 0.003$, shown in the bottom panel, $m=0$ remains locally optimal.

In Fig. 3 we see similar behavior with an increased $T_{c}\approx 0.187$ and a persistent metastable $m=0$ state at low temperature for $\alpha=0.75$.

From an algorithmic perspective, the low-temperature stability of $m=0$ is quite intriguing as it is widely believed to correlate with genuine combinatorial hardness, sounding the death knell for heuristic approaches attempting to locally optimize trial configurations. This category of algorithms includes workhorses such as simulated annealing and parallel tempering Monte Carlo, which exploit correlations in the energy landscape to search by performing biased random walks. Paramagnetic stability suggests that a problem is hard because it implies that cooling an initially disordered configuration will overwhelmingly lead to states that are also disordered. There is no exploitable information within the landscape guiding the dynamics to the correct region of the state space, and the only hope is to begin the search from the appropriate basin. Disordered states form the vast bulk of the configuration space, however, so the probability of a correct initialization decreases exponentially with system size. While all problems with $\alpha<1$ exhibit this feature, the specific value of $\alpha$ turns out to modulate the size of the planted solution basin, i.e., all states from which the ground state can be reached at reasonable cost, which of course exerts a critical effect on the observed performance of heuristics. This aspect can be observed by comparing the intervals between the end points and maxima of $\tilde{f}$ in Figs. 2 and 3 and noting them to be wider for $\alpha=0.75$. Broadly speaking, small values of $\alpha$ lead to smaller basins and hence, one would surmise, more difficult problems. This assertion assumes, however, that the vectors $\{{\boldsymbol{w}}^{\mu}\}$ are composed of real numbers and specifiable to arbitrary numerical precision. When the precision is bounded, the interpretation is more involved and is discussed in detail in Sec. IV.

Figure 4 shows how $\tilde{f}$ is affected by $T$ for a value of $\alpha=1.5$, a regime in which we expect $m=0$ to become unstable. The first-order transition is still readily apparent at $T_{c}=0.68$. For temperatures lower than $T_{c}$ but higher than $T_{u}=0.5$, the paramagnet remains a local minimum. Finally, when cooled to $T_{u}=0.5$, $m=0$ becomes a maximum along the planted ground state direction. From a computational perspective, we anticipate that a well-designed local algorithm will be likely to succeed below this temperature as it follows the free-energy gradients leading it to the planted solution.

To further appreciate the first-order transition, we illustrate the discontinuity at $T_{c}$ of some key thermodynamic quantities for a few representative $\alpha$ values. A key observation is that the nature of the transition changes “gracefully” from first to second order as $\alpha$ increases past unity.

First, we recall (see, for example, Ref. Mezard and Montanari (2009)) that the ensemble partition function is related to the mean-field free-energy density $f(m)$ as

[TABLE]

Using Laplace’s method, we obtain the “log partition density”

[TABLE]

where an additive term arising from Laplace’s formula vanishes for large $N$. We define the minimum and minimizer of the free-energy density, respectively, as

[TABLE]

and

[TABLE]

The limiting internal energy density is

[TABLE]

so that

[TABLE]

In Fig. 5, we plot $f^{*}(\beta)$ for four values of $\alpha$, with the corresponding $\beta_{c}=1/T_{c}$ located at the vertical dashed line. For the first two values of $\alpha$ ($0.5$ and $0.75$), the log partition function is clearly a smooth function on either side of the transition but is not differentiable at $\beta_{c}$. At $\alpha=2.5$, the two sides still have different derivatives but the discrepancy is diminished, and when $\alpha=5$ it is no longer visibly discernible.

Figure 6 shows the magnetization $m^{*}(\beta)$ for the same four values of $\alpha$; the smaller $\alpha$ values show a sharp, discontinuous bifurcation at $\beta_{c}$ at which point the system abruptly shifts from being paramagnetic to being closely aligned with the planted solution. The jump is still apparent but attenuated for $\alpha=2.5$, while for $\alpha=5$ the bifurcation starts to resemble a second-order ferromagnetic transition.

Finally, the internal energy densities are plotted in Fig. 7. Once again, the energies drop discontinuously at $\beta_{c}$ from some excited values to the planted solution ground state energies (per spin) of $-\alpha/2$, which is most evident for the smaller two values of $\alpha$, diminished but still apparent when $\alpha=2.5$, and essentially invisible for $\alpha=5$. Note that the energy gap between the two sides of $\beta_{c}$ is larger when $\alpha=0.75$ than when $\alpha=0.5$, in accordance with its higher transition temperature. The gap closes again for large $\alpha$ but the energies thereafter progress continuously to their ground-state values.

The gradual reversion to a continuous transition is sensible in light of the construction procedure; when $M$ is very large, the matrix ${\boldsymbol{\tilde{J}}}$ converges to $-\alpha{\boldsymbol{\Sigma}}$, where ${\boldsymbol{\Sigma}}$ is the ${\boldsymbol{w}}$ covariance matrix. Since $\Sigma_{ij}=-\frac{1}{N-1}$, this implies that in the large $\alpha$ limit, the system becomes a (scaled) Curie-Weiss ferromagnet with a transition temperature of $T\sim\alpha$.

III.3 Monte Carlo simulation

In this section we perform finite-temperature Monte Carlo simulation of fully connected spin glasses with couplers drawn from a Wishart distribution. To detect the existence of a finite-temperature phase transition we measure the Binder cumulant Binder (1981) given by

[TABLE]

where $\langle\cdots\rangle_{T}$ represents a thermal average, $[\cdots]_{\rm av}$ represents a disorder average, and

[TABLE]

is the magnetization per spin. Although the model is disordered, it orders into a ferromagnetic phase because the ferromagnetic ground state is planted by construction. In general, the Binder ratio scales as

[TABLE]

where $\widetilde{G}(.)$ is polynomial for small values of its argument, $T$ is the temperature, and $T_{c}$ the critical temperature. In general, $g(T=T_{c})$ is independent of $N$ and, as such, on can determine $T_{c}$ by the point where data for different $N$ cross. However, because as we shall see the transition is first order, the shape of the Binder ratio as a function of temperature is somewhat different to the commonly known shape in second-order phase transitions (see, for example, Ref. Katzgraber et al. (2006)). When a first-order phase transition is present, the Binder ratio starts at $g(T\to 0)\to 1$, dips into $g(T)<0$ and then plateaus to $g(T\to\infty)\to 0$. One can determine the critical transition temperature by either extrapolating $g(T_{\rm min},L)$ for $L\to\infty$, where $T_{\rm min}$ is the minimum value of $g$. Alternatively, one can study the crossing of data for different system sizes $N$ in the range $0\leq T\leq T_{\rm min}$. The former is usually easier to use to detect a transition, because the latter is often close to $T=0$ where data for different $N$ do not splay enough to see a clean crossing. However, as shown in Refs. Binder et al. (1992); Vollmayr et al. (1993), $g(T_{c})-g(T_{\rm min},N)\sim 1/N$, whereas the corrections to scaling in Eq. (23) decrease proportional to $g(T_{c})-g(T,N)\sim 1/N^{2}$. As such, here we determine the position of the critical temperature by studying the crossing of the Binder ratio.

To further corroborate the existence of a first-order transition to a ferromagnetic phase we study the distribution $P(m)$ of the magnetization [Eq. (22)]. Close to the transition temperature where latent heat is present the order parameter should signal two competing phases, i.e., peaks at $|m|\to 1$, as well as a competing peak at $m=0$.

The simulations are done using parallel tempering Monte Carlo  Geyer (1991); Hukushima and Nemoto (1996). Thermalization is verified by ensuring that all measured observables are independent of simulation time. We do this by analyzing how the results for all observables vary when the simulation time is successively increased by factors of $2$ (logarithmic binning). We require that the last three results for all observables agree within error bars. Simulation parameters are shown in Table 1. Error bars are determined via a jackknife analysis over the disorder.

Figure 8 shows data for the Binder ratio as a function of temperature. Figures 8(a), 8(c), and 8(e) show the Binder ratio for $\alpha=0.50$, $0.75$, and $1.00$, respectively, over the whole temperature range. A negative dip signaling a first-order transition is clearly visible. Figures 8(b), 8(d), and 8(f) zoom into $T\in[0,T_{\rm min}]$. As can be seen, the data cross for all three $\alpha$ values. From the crossing points we estimate $T_{c}(\alpha=0.50)=0.065(5)$, $T_{c}(\alpha=0.75)=0.188(1)$, and $T_{c}(\alpha=1.00)=0.335(5)$, in perfect agreement with our analytical estimates. Note that $T_{c}(\alpha)\to 0$ for decreasing $\alpha$, i.e., the problems become increasingly harder numerically for smaller values of $\alpha>\alpha_{c}$.

In Fig. 9 data for the magnetization distribution are shown for $\alpha=0.50$, $0.75$, and $1.00$ and $T\approx T_{c}$. In all three cases peaks at $|m|\to 1$ are visible, signaling ferromagnetic order. However, a third peak for $m=0$ grows with increasing $N$, thus signaling a jump in the magnetization close to the transition. Note that for decreasing $\alpha$ the first-order transition is more pronounced.

IV Hardness phase transition

Having studied the phase behavior of the WPE and elicited properties that would plausibly correlate with algorithmic difficulty, we now proceed to examine some empirical results in which we probe the typical time to find the ground state using a parallel tempering (PT) Monte Carlo method. This algorithm is widely used to simulate complex physical and biological systems; contemporary implementations routinely examine spin glasses with several thousand variables and are quite competitive as optimization algorithms as well Moreno et al. (2003); Zhu et al. (2020). In essence, it is a careful stochastic local search Hoos and Stützle (2004) method such that configurations are allowed to escape metastable states by traveling to high temperatures and return afresh to probe the low-energy landscape while ensuring asymptotically correct equilibrium sampling at each temperature.

We seek to check and account for the existence of an algorithmic hardness peak, a sharp increase in the time required to find the ground state as $\alpha$ is varied. Unfortunately, our numerical studies were severely hampered by the extreme difficulty of the WPE in its hard regime, which disallowed us from empirically localizing the peak for sizes that would otherwise be considered modest by the standards of other Ising ensembles (e.g., the SK model). We systematically reduced the sizes in which PT failed, within our computational resources, to approach the planted ground-state energy for a range of $\alpha<1$ values. Reliable statistics were finally obtained for $N=32$; the plot of median time to solve the problem, or more precisely and as justified in Sec. IV.3, to find a state with energy within $\epsilon=10^{-7}$ of that of the planted ground state, is shown in Fig. 10. Other choices for $\epsilon$ are certainly sensible, and in Sec. IV.3 we consider and analyze the less stringent values of $\epsilon=10^{-5}$ and $\epsilon=10^{-3}$. Details of the simulation protocol are found in Appendix VII.6.

The hardness peak is evident at $\alpha=3/N=3/32$, where the problems typically took around $5$ minutes to solve. As expected, the instances get steadily easier for the larger $\alpha$ values, but the observation that $M=1$ is easy may initially seem surprising. As alluded to previously, this is due to the fully frustrated nature of $M=1$, giving rise to a tremendous number of low-energy states, coupled with the fact that standard double-precision floating point arithmetic was used for the computations. Consequently, the potential exists for several states that under unbounded precision would have had close but nonetheless distinct energies to be mapped to the same value. Under the low-energy degeneracy associated with $M=1$, one thus observes a huge and exponentially increasing (in $N$) number of numerically indistinguishable “ground states”; finding one such acceptable state turns out to be relatively easy com (c).

Yet, the number of solutions alone cannot be expected to predict problem difficulty; in the computationally easy large-$M$ regime, the Hamiltonian starts to look increasingly like a ferromagnet, and so is essentially assured of having a unique solution. Thus, in contrast to uncorrelated problems like number partitioning in which an easy-hard transition is observed at the parameter value such that “perfect” solutions to the problem disappear Gent and Walsh (1996); Mertens (2006), the correlations resulting from the planting procedure give rise to another factor influencing difficulty and which leads to the observed easy-hard-easy pattern. We conjecture that hardness at a given $\alpha$ is determined by two competing factors: First, the number of solutions satisfying the integer program and, second, the intrinsic search space size of the integer program. It seems plausible that the regime in which the ratio of these two quantities is minimized signals the hardness transition as this would be analogous to having the fewest needles in the biggest haystack. Our aim is to analytically predict the hardness peak for the $N$-variable WPE at given $\alpha$ and numerical tolerance $\epsilon$.

To determine the number of satisfying solutions under precision $\epsilon$, we first derive the ensemble-averaged energy histogram for the WPE, discussed in Sec. IV.1. We then formalize our notion of the intrinsic search space in Sec. IV.2; as anticipated in Sec. II.1, the size of the nullspace of ${\boldsymbol{W}}^{T}$ plays a key role. Finally, we make our quantitative conjecture regarding the transition location in Sec. IV.3 and show that it precisely predicts the location of the peak in Fig. 10. In Sec. IV.4, we present some preliminary numerically obtained features of the WPE energy landscape to complement the main analysis performed here.

IV.1 Energy histogram

We consider an ensemble of problems over $N$ spins with parameter $M$ such that

[TABLE]

with the columns $\{{\boldsymbol{w}}^{\mu}\}$ of ${\boldsymbol{W}}$ simulated as $\mathcal{N}(\mathbf{0},{\boldsymbol{\Sigma}})$ and Hamiltonian

[TABLE]

for ${\boldsymbol{s}}\in\mathbb{S}^{N}$. Note that the diagonal elements are included in this formulation as it is presently more convenient to take the ground-state energy to be zero. Without loss of generality, we assume the ferromagnetic ground state is planted so that for all $\mu$

[TABLE]

We seek the distribution of energies marginalized over the problem ensemble

[TABLE]

where $p_{E}(e|\{{\boldsymbol{w}}^{\mu}\})$ is the probability of drawing by uniform sampling a state with energy $E=e$ from the given problem.

The result is surprisingly simple; the energies follow a gamma distribution with density

[TABLE]

and cumulative distribution $P_{E}(e)\triangleq\Pr(E\leq e)$

[TABLE]

where $\Gamma(x)$ and $\gamma(x,y)$ are the gamma and incomplete gamma functions Abramowitz and Stegun (1964), respectively. The calculation of the distribution is shown in Appendix VII.3. It is quite similar to Mertens’ Mertens (2001) analysis of the number partitioning problem cost density, but here the energy is a more complicated quadratic function with a correlated coupler matrix rather than a rank-1 absolute value discrepancy over independent coefficients. Figure 11 displays $p_{E}(e)$ for the WPE when $N=64$ for a few values of $M$. When $M=1$, the density is overwhelmingly concentrated on low-energy values; combined with bounded precision, this degeneracy is responsible for the problems being easy. By the properties of the gamma distribution, the mean and standard deviation of $E$ are $M/2$ and $\sqrt{M/2}$ respectively; consequently (nonplanted) low-energy states become exponentially less likely as $M$ increases.

If the coupler matrix ${\boldsymbol{J}}$ rather than ${\boldsymbol{\tilde{J}}}$ were used, the density would simply be translated by the ground-state energy of $-\alpha N/2$, i.e.,

[TABLE]

We point out that Eq. (25) in fact describes the energy density of non-planted states, which as is shown in the Appendix, is dominated by paramagnetic states ($m=0$) far from the planted solution. In fact, the true density should be regarded as a mixture of a $\delta$ function at zero energy with weight $2^{-(N-1)}$, representing the probability of sampling the planted solution, and of the gamma distribution (25) with weight $1-2^{-(N-1)}$.

The distribution enables evaluation of the ensemble-averaged number of states whose energies are at most $\epsilon$ and which are not related by a global spin flip com (d), which turns out to be highly relevant to computational properties. This can be shown to be

[TABLE]

The leading “1” is due to the persistence of the planted solution; its presence in the expectation may initially seem strange and for small values of $\alpha$ it is indeed irrelevant, but as $\alpha$ increases, the paramagnetically contributed solutions start to vanish and its influence in fact becomes dominant.

IV.2 Intrinsic search space cardinality

We have speculated that because we aim to solve an Ising-constrained linear program or, equivalently, to solve for a state ${\boldsymbol{s}}^{*}\in\mathbb{S}^{N}$ lying in the nullspace of ${\boldsymbol{W}}^{T}$, the dimensionality of the nullspace plays some role in problem complexity. In this section, we clarify the notion and extract an $M$-dependent intrinsic search space. This is defined as a discrete set of reduced-dimensional states over which one would need to exhaustively search to solve the problem without reference to any information about the values of ${\boldsymbol{W}}$. In some sense, this is analogous to the reduced complexity enjoyed by various optimization problems on graphs of low treewidth. The result is that when ${\boldsymbol{W}}^{T}$ has $M$ independent rows, one only needs to consider $O(2^{N-M-1})$ states in brute-force search; the subtraction by one is simply to remind that the problem is invariant to global spin flips. We note that the following procedure is not meant for literal implementation because we are disregarding, for example, numerical issues that may be important in practice, but simply to show that up to a polynomial prefactor the exhaustive search complexity is exponential in $K$, where

[TABLE]

We remind the reader that the integer programs are always feasible due to the planting construction. Let the vectors

[TABLE]

span $\textrm{null}({\boldsymbol{W}}^{T})$. Such a basis can be obtained in polynomial time using, for example, the singular value decomposition of ${\boldsymbol{W}}$. The integer program can be equivalently phrased as

[TABLE]

where ${\boldsymbol{V}}$ is the $N\times K$ matrix whose columns are $\{{\boldsymbol{v}}_{i}\}$ and ${\boldsymbol{x}}\in\mathbb{R}^{K}$. For any state ${\boldsymbol{s}}$, one can in polynomial time either solve for satisfying coefficients ${\boldsymbol{x}}$ or show that no such coefficients exist using standard linear algebra techniques. While a solution will eventually be found when the method is repeated for all possible states, we now show that not all ${\boldsymbol{s}}\in\mathbb{S}^{N}$ need be checked.

We first obtain a $K\times K$ matrix ${\boldsymbol{\tilde{V}}}$ composed of a subset of $K$ independent rows from ${\boldsymbol{V}}$. Such a matrix can be determined using for example the QR decomposition with pivoting Golub and Van Loan (2012). We apply this decomposition to the columns of ${\boldsymbol{V}}^{T}$:

[TABLE]

where ${\boldsymbol{P}}$ is an $N\times N$ column-permuting matrix, $Q$ is $K\times K$ and orthogonal, and $R$ is $K\times N$ and upper-triangular. This ensures that the first $K$ columns of ${\boldsymbol{V}}^{T}{\boldsymbol{P}}$ are independent, so we take

[TABLE]

where ${\boldsymbol{P}}_{K}$ is an $N\times K$ submatrix composed of the first $K$ columns of ${\boldsymbol{P}}$ and encodes which rows of ${\boldsymbol{V}}$ were selected for ${\boldsymbol{\tilde{V}}}$. These steps are computed only once for a given problem.

Now let $\tilde{{\boldsymbol{s}}}$ be an Ising state of length $K$. The system

[TABLE]

with respect to partial state $\tilde{{\boldsymbol{s}}}$ always has a solution because ${\boldsymbol{\tilde{V}}}$ is full rank. After solving for $\tilde{{\boldsymbol{x}}}$, we verify whether this generalizes to a feasible solution over all Ising variables by computing

[TABLE]

The first $K$ components of ${\boldsymbol{y}}$ will be $\tilde{{\boldsymbol{s}}}$ by construction; if the rest are also in $\{\pm 1\}$, then ${\boldsymbol{y}}$ is a feasible solution. If not, then we choose another $\tilde{{\boldsymbol{s}}}$ and repeat the process starting with solving for $\tilde{{\boldsymbol{x}}}$. This shows that up to the polynomial overhead involved in the various steps, one only needs to consider the number of states in $\mathbb{S}^{K}$; in fact due to the spin-flip symmetry, the intrinsic search space contains $2^{K-1}$ elements. In the case where $M\geq N$ and so the nullspace is one dimensional, only a single state need be checked.

IV.3 Predicting the hardness transition

We now synthesize the two preceding antithetical factors into an expression that predicts when problems are likely to become difficult. For illustration, we first consider the case of degenerate ground states, in which the reasoning does not in fact make any assumptions about the specific generative process or family from which problems are drawn. The number of solutions not related by global flip symmetry to a specific problem is denoted by $N_{G}$, a random variable whose distribution may correspond to some parametric setting.

In the exhaustive algorithm described in the previous section, a set of no more than $2^{N-M-1}$ partial states $\{\tilde{{\boldsymbol{s}}}\}$ need be checked for feasibility. Given a traversal order over the partial states, the expected time to locate a solution to a problem is simply the average number of iterations until a specific $\tilde{{\boldsymbol{s}}}$ mapping to a feasible completion is found. Noting that each partial state corresponds to at most one solution and assuming that an adversary uniformly allocates the solutions among the partial states, by searching among $\{\tilde{{\boldsymbol{s}}}\}$ in arbitrary order, the solution time can be identified with the number of draws of a partial state (without replacement) from a population of $2^{N-M-1}$ such states, $N_{G}$ of which are solutions, until a single solution is found. This quantity follows a specific type of negative hypergeometric distribution, where for convenience we obtain the inverse of its mean $\tau$ for our particular values as

[TABLE]

When averaged over all problems in the ensemble, the expected inverse solution time is thus approximately

[TABLE]

Equation (29) makes explicit the competing effects of the number of solutions (numerator) and search space size (denominator) on problem difficulty.

Returning our attention to the WPE, we presently determine how to compute the expected number of solutions, or equivalently the mean number of ground states in the Ising formulation, defined as

[TABLE]

where the degeneracy results from the discreteness of the generator variables $\{{\boldsymbol{z}}^{\mu}\}$. We recall from Sec. II.3 that in the WPE, all independent and standardized $\{{\boldsymbol{z}}^{\mu}\}$ yield the same limiting properties as those resulting from using a Gaussian and presented in this work. Consequently, a good estimate for the ground-state degeneracy can be obtained for sufficiently large systems using the expected count in Eq. (27). For example the number of ground states in the Rademacher-discretized WPE presented in Sec. II.3 can be approximated by

[TABLE]

where $\epsilon_{1}=O(1/N^{3})$ represents the smallest attainable excited-state energy for the problem size.

As mentioned, the computational expense associated with simulating the WPE in its hard regime forced us to simulate systems of relatively small size. Such a restriction can, however, potentially introduce finite-size effects; in other words, the properties of the discretized systems may be far from those predicted by their Gaussian asymptotics. To minimize such artifacts, we opted to sample the parameters using Gaussian generators represented with the maximum available precision and consequently, to contend with numerical errors. The task was thus translated to that of approximate solution.

Bounded precision means we cannot represent all the “true” values of the Gaussian variables $\{{\boldsymbol{z}}^{\mu}\}$ used to generate ${\boldsymbol{W}}$ nor can we in general exactly represent the values of ${\boldsymbol{W}}$ resulting from the planting procedure using the rounded $\{{\boldsymbol{z}}\}$ variables or even the exact energy of a state relative to the rounded ${\boldsymbol{W}}$. Rounding and discretization cause errors in the representation of some energies; in particular, states that had distinct energies under arbitrary precision may be mapped to numerically indistinguishable values, and the relative ordering among closely spaced true energies may change. Let $\widehat{E}$ be the random variable resulting from finite-precision representation of some true continuous-valued energy $E$; its corresponding CDF $P_{\widehat{E}}(\epsilon)\triangleq\Pr(\widehat{E}\leq\epsilon)$ now has step discontinuities. If we interpret the mapping of a state’s exact energy $E$ to its representable approximation $\widehat{E}$ as a pseudo-random walk process leaving invariant the probability of having energy of at most $\epsilon$, then provided $\epsilon$ is large enough to be representable we can approximate the finite-precision energy CDF with the continuous distribution introduced in Sec. IV.1:

[TABLE]

where the latter equality follows from the nonnegativity of $E$. We can then approximate the expected number of states with observed energy of at most $\epsilon$ using Eq. (27):

[TABLE]

The value of $\epsilon$ is in principle arbitrary but must be large enough to reflect the approximation error of the planted ground-state energy. We observe that under the aggregated discretization effects of the construction procedure, the true planted ground-state energy of zero is typically distorted to a value on the order of the so-called machine epsilon $\approx 10^{-16}$ for double-precision arithmetic. A sensible strategy is thus to take $\epsilon$ to be some value not much larger than this, which would conservatively ensure that no planted ground states are missed while accepting a certain number of “false positives,” i.e., that some nonplanted states may be considered successful solutions. Excessively small values of $\epsilon$, however, resulted in prohibitively long computational times for the problems in the hard regimes. Our methodology to analyze the empirical solution time settled on using three target energy values: $\epsilon\in\{10^{-7},10^{-5},10^{-3}\}$. We hence strive to yield accurate predictions of the hardness transition for the important task of approximate (or relaxed) solution to the problem.

Finally, we define the function trading off number of solutions to effective search space size

[TABLE]

and anticipate that problems become most difficult at

[TABLE]

as this coincides with the fewest number of solutions relative to the effectively constrained search space. In Fig. 12, we plot $\log Q_{N}(M)$ for system sizes $N=32$ and $N=128$, when $\epsilon\approx 10^{-7}$ in accordance with our most stringent target. When $N=32$, $M^{*}=3$ in perfect accordance with the location of the hardness peak shown in Fig. 10. The latter size was far too large to solve in the hard regime within our constraints, but the plot suggests that to obtain the most difficult problems for this larger size, $M$ needs to increase. It is instructive to predict the required scaling of $M$ with $N$ for maximally hard problems at this level of precision.

Figure 13 shows $M^{*}$ as a function of $N$, where a clear linear scaling is apparent. By linear regression, we find this relation to be

[TABLE]

suggesting that the hardest problems at $\epsilon=10^{-7}$ occur when $\alpha\approx 0.073$. The general message is clear: To have truly difficult problems under precision constraints, the number of equations $M$ in the integer program cannot be constant.

The reader may notice that the definition of $Q_{N}$ resembles that of the mean inverse solution time (29) with respect to the exact partial state traversal algorithm. Directly relating $Q_{N}$, which considers all energies in a specified range, to the solution time of an exact solver is not straightforward however. The issue is that a partial state $\tilde{{\boldsymbol{s}}}$ may now have several completions to full Ising states whose energies match the target but which are not ground states. While we have shown that finding the ground-state completion of a given $\tilde{{\boldsymbol{s}}}$ or verifying that no such completion exists is straightforward, locating a completion guaranteed to match a generic target $\epsilon>0$ is nontrivial. We may nonetheless heuristically justify $Q_{N}$ as follows. Suppose that fixing the $N-M$ partial state variables $\tilde{{\boldsymbol{s}}}$ sufficiently constrains the remaining free variables that locating the extension of minimum energy can typically be done rapidly and with high probability with a heuristic algorithm. An “exhaustive-approximate” algorithm for finding a state whose energy is bounded by $\epsilon$ can thus proceed as follows: Traverse the $\{\tilde{{\boldsymbol{s}}}\}$ in some order, where for each $\tilde{{\boldsymbol{s}}}$, the minimizing extension is heuristically obtained, and stop if a resultant state’s energy is less than $\epsilon$. Under the assumption that the target states are uniformly distributed among the partial state “bins,” $Q_{N}$ can once again have an interpretation as an expected inverse solution time. To justify our conjecture, we now show that $Q_{N}(M)$ can serve to localize the hardness peak for generic values of $\epsilon$.

Figure 14 shows the PT median times to solution for the same $N=32$ WPE ensemble considered so far, but for three values of $\epsilon$ defining an acceptable solution: $\epsilon=10^{-7}$, $10^{-5}$, and $10^{-3}$. Unsurprisingly, the typical solution times decrease as the energy criterion becomes more permissive. Additionally, we note that the location of the hardness peak shifts to larger values of $M$ as $\epsilon$ is relaxed. Most importantly, in Fig. 15 we observe the predictive power of $Q_{N}(M)$ at these values of $\epsilon$: In all three cases, the minimizer $M^{*}$ of $Q_{N}(M)$ at the respective values of $\epsilon$ precisely corresponds with empirically observed PT solution time. We have hence proposed a robust, theoretically motivated framework for generating tunably difficult problems over a wide range of approximate solution targets.

IV.4 Properties of locally optimal states

As we have seen, the energy histogram derived in Sec. IV.1 has been useful in predicting the algorithmic properties of WPE instances. Nonetheless, this distribution does not provide information about topological aspects of the local minima, i.e., states that are energetically stable with respect to a single spin flip.

In this section, we briefly probe the properties of local minima using exhaustive search on small instances; we save analytic examination of these properties, along the lines of Bray and Moore’s Bray and Moore (1980) analysis for the SK model, for later work.

A natural statistic to analyze is the expected number of local minima as $\alpha$ is varied. Furthermore, it is instructive to define a residual energy histogram restricted to stable states. These are shown for a system with $N=24$ variables in Figs. 16 and 17, respectively. As expected, the number of minima decreases monotonically in $\alpha$ as the ensemble tends toward a ferromagnet. At small $\alpha$, we observe a large number of minima, but also that the residual energy itself is likely to be small. This is consistent with the observation that at small $\alpha$, the problems with restricted precision are easy.

To further illustrate properties of the WPE, Fig. 18 shows examples of disconnectivity graphs for four specific instances at their respective values of $M$. Disconnectivity graphs are two-dimensional representations of high-dimensional energy landscapes. In their simplest form they depict the minima of the system and the lowest energy barrier connecting any two minima, where an energy barrier is defined as the highest energy value encountered along a specific pathway. The barriers represent the minimum increase in energy necessary to transition from one minimum to another. In this work, the minima of the system were obtained via complete enumeration and the barriers were calculated using a search over all possible pathways identified with flipping spins that are misaligned between pairs of minima. This method was outlined by Garstecki, Hoang and Cieplak Garstecki et al. (1999). To deal with the large number of minima in a computationally efficient manner, we further made use of an approximation based on the relative proximity of minima; specifically, for each minimum only the barriers to the $50$ minima closest to it in Hamming distance were obtained. This approximation is based on the fact that transitions between two minima can also happen through basins of intermediate minima, and hence if two minima are separated by a large Hamming distance it is likely that the lowest energy barrier between these two minima will be already represented via transitions between intermediate minima and their corresponding barriers.

In Fig. 18 the minima are represented by vertical bars whose lowest points denote their energies. On top, they are connected by lines converging to a common point representing the height of the barrier that needs to be crossed in order to transition between the connected minima. Due to the continuous nature of the energy values, we sort the minima into a hierarchical cluster structure whose end points comprise the intervals $[E_{b}-\Delta\ell,E_{b}+\Delta\ell)$, where $E_{b}$ is the energy of the barrier and $2\Delta\ell$ denotes the length of the interval. Minima whose connecting barriers fall within the same energy interval are sorted into a common cluster. Within an individual cluster the minima are arranged based on the number of spins in the up or corresponding down states. Minima which have a high number of up-state spins are sorted toward the left and, correspondingly, minima with a high number of down-state spins are sorted toward the right. Note, that this order strictly only applies within the individual cluster; the order of the individual clusters relative to each other is determined by the hierarchical structure. In this work we set $\Delta\ell=0.075$. Figure 18 shows a clear progression in the energy landscape for which small values of $M$, e.g., $M=4$ in Fig. 18 (a), are characterized by a very large number of almost degenerate metastable states, and larger values of $M$ [Figs. 18(b) –18(d)] tend to break the degeneracy, emphasize the planted ground state, and make the landscape more funnel-like.

Figure 19 shows the distribution of the average barrier height within the first Hamming distance of its individual minima states. It represents the distribution of the average increase in energy necessary for the system to escape its minima via the shortest route, i.e., to transition to adjacent minima. For a given instance, the average $\langle\cdots\rangle$ was taken over the barriers to the minima which lie within the shortest Hamming distance to a given minimum, say minimum $i$. $E_{b}$ are the energies of the individual barriers and $E_{\mathrm{min}}$ is the energy of minimum $i$. The overall distribution is then obtained by sampling over all minima $i$ of $N=100$ sample systems for each of the values of $\alpha$.

As can be seen in Fig. 19, at small values of $\alpha$ the distribution is more spread than at large values, where additionally it is dominated by relatively small values of the average energy barrier. From an energy landscape perspective this indicates that transitions from the minima at large $\alpha$ are more likely to be achievable with less energy cost and therefore might be more probable than at small values of $\alpha$ com (e).

V Discussion

We have proposed a planted Ising ensemble with several noteworthy physical and algorithmic properties. The model exhibits a first-order temperature transition, a persistent locally stable paramagnetic state when $\alpha<1$, and when represented with finite precision, an easy-hard-easy algorithmic difficulty profile. Its physical properties are consistent with the observed hardness of finding its ground state; moderately sized problems are extremely difficult in the hard regime. This meshes well with the intuition that the transition and paramagnetic stability give rise to a golf-course-like energy landscape. After deriving the instance-averaged energy distribution—which turns out to follow a gamma law—we compare the expected number of states matching a solution criterion with a quantity we introduced quantifying the intrinsic search space size at given $\alpha$ to analytically predict the location of the hardness peak. The prediction is validated using solution times obtained with a highly optimized implementation of parallel tempering Monte Carlo.

The first-order transition between the planted and paramagnetic phases is furthermore established by the TAP analysis in Sec. III.1, with alternative derivations supported by the replica method (Appendix VII.5) and annealed approximations (Appendix VII.4). Careful Monte Carlo simulations demonstrate the correctness of our analytically predicted transition temperatures.

A connection is made with the anti-Hopfield model and this is developed analytically in Appendix VII.5 Nokura (1998). This analysis indicates that we might expect to see the impact of replica symmetry breaking if we focus on the region orthogonal to the planting space. In this paper we have not focused on this for two reasons. First, the emphasis has been on intermediate scale problems, where it is difficult to numerically establish such phenomena, finite size effects may dominate, and practical issues such as finite-precision representation may have greater impact. Second, at large $N$ replica symmetry breaking has little impact on the free-energy barrier separating the bulk of the space from the planted solution, which is the primary driver of hardness for heuristic optimization in the planted case.

Nevertheless, we believe the qualitative description of the space orthogonal to the planted solution, that of a rough energy landscape with deep solutions almost orthogonal to the planted solution, is in effect. Roughness is apparent at the small scales we have worked with empirically. With this in mind, a modification of the planting procedure is pursued in Appendices VII.4 and VII.5 whereby the planted solution is partially penalized. Its leading-order energy is then tunable at fixed $\alpha$, at the cost of losing strict guarantees that it is a ground state. This modification, it is believed, might be the basis for interesting tests of dynamics (escape into or out of the planted solution). In particular, an interesting suggestion has been made that for this type of problem in the presence of a transverse field, quantum dynamics may be differentiated from classical counterparts Smelyanskiy et al. (2018). The features of our model with this modified planting procedure represent a practical realization of many of the abstract model features underlying the population transfer hypothesis.

Future work will explore this direction more deeply, and consider as well whether the features of the WPE energy landscape lend themselves to the demonstration of fundamental speedup using emerging quantum annealing devices, in particular those whose classical simulation is known to be intractable. It is our hope that the insights into the ensemble’s physical and algorithmic classical properties presented in this work will solidly underpin these future directions.

In addition, we are interested in pursuing the relation between the WPE and the important class of low-rank estimation problems Lesieur et al. (2017) used in unsupervised feature extraction and dimensionality reduction, in further examining whether any deeper connections can be made between the ensemble and biological unlearning Nokura (1998), and in studying potential connections with models for analyzing large-system code-division multiple-access (CDMA) multiuser detectors Tanaka (2002).

VI Acknowledgements

We acknowledge useful discussions with Bill Macready, Federico Ricci-Tersenghi, Hidetoshi Nishimori, Jon Machta, Cris Moore, Koji Hukushima, Brad Lackey, Lenka Zdeborová, Salvatore Mandrà, Mohammad Amin, Andrew King, Cathy McGeogh, and Aidan Roy. H.G.K. thanks the Department of Poultry Science at Texas A&M University for providing support. H.G.K. and C.P. are supported in part by the Office of the Director of National Intelligence (ODNI), Intelligence Advanced Research Projects Activity (IARPA), via MIT Lincoln Laboratory Air Force Contract No. FA8721-05-C-0002. The views and conclusions contained herein are those of the authors and should not be interpreted as necessarily representing the official policies or endorsements, either expressed or implied, of ODNI, IARPA, or the U.S. Government. The U.S. Government is authorized to reproduce and distribute reprints for Governmental purpose notwithstanding any copyright annotation thereon.

VII Appendix

VII.1 Cavity calculation of the TAP equations

For notational convenience, we redefine the Hamiltonian of the WPE to include the constant terms corresponding to the diagonal elements of ${\boldsymbol{\tilde{J}}}$; this will not affect the final result. Further, without loss of generality we assume that the ferromagnetic solution ${\boldsymbol{t}}=(1,1,\ldots,1)$ and its reversal are to be planted as the ground state. The Hamiltonian of the WPE is

[TABLE]

where

[TABLE]

and ${\boldsymbol{w}}^{\mu}\sim\mathcal{N}(\mathbf{0},{\boldsymbol{\Sigma}})$. As before, $M=\alpha N$ for $\alpha>0$ and $N\to\infty$. The elements of the covariance matrix ${\boldsymbol{\Sigma}}$ are as follows:

[TABLE]

This covariance structure implies that the Gaussian process generating ${\boldsymbol{w}}$ is not only stationary but exchangeable. As discussed in Sec. III, this model bears many resemblances to the Hopfield network, but in addition to the negation of the interactions, we have the additional complication of the “antipatterns” $\{{\boldsymbol{w}}^{\mu}\}$ consisting of correlated elements. Nonetheless, by closely following the two-step cavity-based method presented by Shamir and Sompolinksi Shamir and Sompolinsky (2000) for the Hopfield model, taking special care to account for the correlations, we may derive the TAP equations for the WPE.

The cavity approach Mézard et al. (1987) derives the self-consistent relation for each local spin magnetization by first considering removal of the spin from the $N+1$-spin system and defining a state distribution on the $N$-spin subsystem. Remarkably, the joint distribution of the field and spin in the original system can be expressed in terms of the field distribution resulting from the subsystem, allowing the spin and field statistics for the full system to be related to those of the subsystem. On their own, these exact relations do not give much insight because they are intractable to compute, but when the subsystem field distribution can be justified to be Gaussian, the TAP equations for the magnetizations may be obtained. A substantial amount of the work is in deriving the correct parameters for the field distribution. The details for the WPE follow.

VII.1.1 Cavitating a spin

Consider an $(N+1)\times(N+1)$ WPE matrix ${\boldsymbol{\tilde{J}}}$ coupling spins $\{0,\ldots,N\}$ through $M$ vectors $\{{\boldsymbol{w}}^{\mu}\}$, i.e., $\widetilde{J}_{ij}=-\frac{1}{N+1}\sum_{\mu=1}^{M}w_{i}^{\mu}w_{j}^{\mu}$. We can decompose this Hamiltonian into a sub-Hamiltonian $H^{(N)}$ consisting only of interactions among spins $\{1,\ldots,N\}$, denoted here by $s_{1:N}$, and a term accounting for the interaction between spin zero and the others:

[TABLE]

with

[TABLE]

and

[TABLE]

The final constant in $H^{(N+1)}$ is irrelevant and will be dropped. The exact joint distribution over $(s_{0},h_{0})$ can be shown to be

[TABLE]

where

[TABLE]

and $\langle{\cdots}\rangle_{N}$ refers to thermal averaging with respect to $P^{(N)}(s_{1:N})=\frac{1}{Z_{N}}\exp(-\beta H^{(N)}(s_{1:N}))$. From this, we obtain the (intractable) relations

[TABLE]

and

[TABLE]

The next step is to compute the field statistics to be used following a Gaussian assumption for $P^{(N)}(h_{0})$.

As for the Hopfield model, define the field mean and variance as

[TABLE]

Note that

[TABLE]

where

[TABLE]

Using the definition of $\widetilde{J}_{0i}$, we obtain

[TABLE]

If we define the overlap of spins $\{1,\ldots,N\}$ with the last $N$ components of ${\boldsymbol{w}}^{\mu}$ as

[TABLE]

and its covariance under the cavitated spin distribution

[TABLE]

then we obtain

[TABLE]

and hence

[TABLE]

Noting $\chi_{ii}^{(N)}=1-\langle{s_{i}^{2}}\rangle_{N}=O(1)$ while $\chi_{ij}^{(N)}=O(\frac{1}{\sqrt{N}})$ for $i\neq j$, we proceed to determine the magnitude of $\langle{\delta\eta_{\mu}\delta\eta_{\nu}}\rangle_{N}$ in order to simplify the field variance, bearing in mind that while ${\boldsymbol{w}}^{\mu}$ is independent of ${\boldsymbol{w}}^{\nu}$, there are componentwise correlations within each vector not present in the Hopfield model. The conclusion will be that just as for the Hopfield model, $\langle{\delta\eta_{\mu}\delta\eta_{\nu}}\rangle_{N}$ is $O(\frac{1}{N^{3/2}})$ when $\mu\neq\nu$ and $O(\frac{1}{N})$ and when $\mu=\nu$.

Define

[TABLE]

where the superscript on $\chi$ has been dropped. We seek $\mathbb{E}[S_{\mu\nu}]$ and the fluctuations $\sqrt{\mathbb{E}[S_{\mu\nu}^{2}]-\mathbb{E}[S_{\mu\nu}]^{2}}$, where the expectations are taken over $\{{\boldsymbol{w}}\}$.

The linear expectations are straightforward to compute; when $\mu\neq\nu$, we have

[TABLE]

while, recalling the covariance structure of ${\boldsymbol{w}}$, when $\mu=\nu$,

[TABLE]

which is different from the Hopfield model, in which this quantity is zero.

Direct computation of the quadratic expectation $\mathbb{E}[S_{\mu\nu}^{2}]$ is tedious but straightforward. In the expansion, there will be a total of $N^{4}$ terms of the form $\chi_{ij}\chi_{kl}w_{i}^{\mu}w_{j}^{\nu}w_{k}^{\mu}w_{l}^{\nu}$. We first make a relevant partitioning of these terms for generic $\{\mu,\nu\}$ and then compute the expectations for the cases where they are equal and different.

• •

Case $i=k$, $j=l$

If $i=j$, then there are $N$ terms of the form $\chi_{ii}^{2}(w_{i}^{\mu})^{2}(w_{i}^{\nu})^{2}$

If $i\neq j$, then there are $N(N-1)$ terms like $\chi_{ij}^{2}(w_{i}^{\mu})^{2}(w_{j}^{\nu})^{2}$

• •

Case $i=k$, $j\neq l$

If $i=j$, then there are $N(N-1)$ terms as $\chi_{ii}\chi_{il}(w_{i}^{\mu})^{2}w_{i}^{\nu}w_{l}^{\nu}$

If $i\neq j$, then we have $N$ terms like $\chi_{ij}\chi_{ii}(w_{i}^{\mu})^{2}w_{i}^{\nu}w_{j}^{\nu}$ and $N(N-2)$ terms as $\chi_{ij}\chi_{il}(w_{i}^{\mu})^{2}w_{j}^{\nu}w_{l}^{\nu}$

• •

Case $i\neq k$, $j=l$

This is a rotated version of the previous case. When $i=j$, there are $N(N-1)$ terms of the form $\chi_{jj}\chi_{kj}w_{j}^{\mu}w_{k}^{\mu}(w_{j}^{\nu})^{2}$

When $i\neq j$, there are $N$ terms like $\chi_{ij}\chi_{jj}w_{i}^{\mu}w_{j}^{\mu}(w_{j}^{\nu})^{2}$ and $N(N-2)$ terms as $\chi_{ij}\chi_{kj}w_{i}^{\mu}w_{k}^{\mu}(w_{j}^{\nu})^{2}$

• •

Case $i\neq k$, $j\neq l$

When $i=j$, there are $N(N-1)$ terms like $\chi_{ii}\chi_{kk}w_{i}^{\mu}w_{k}^{\mu}w_{i}^{\nu}w_{k}^{\nu}$ and $N(N-1)(N-2)$ terms as $\chi_{ii}\chi_{kl}w_{i}^{\mu}w_{k}^{\mu}w_{i}^{\nu}w_{l}^{\nu}$

When $i\neq j$, there are $N(N-1)(N-2)$ terms as $\chi_{ij}\chi_{kk}w_{i}^{\mu}w_{k}^{\mu}w_{j}^{\nu}w_{l}^{\nu}$ and $N(N-1)(N^{2}-3N+3)$ terms like $\chi_{ij}\chi_{kl}w_{i}^{\mu}w_{k}^{\mu}w_{j}^{\nu}w_{l}^{\nu}$

Consider now the case of $\mu\neq\nu$. Adding the terms in the expansion, recalling that ${\boldsymbol{w}}^{\mu}$ is independent of ${\boldsymbol{w}}^{\nu}$ and that $\mathbb{E}[w_{i}w_{j}]=-\frac{1}{N}$ when $i\neq j$, we find that

[TABLE]

implying that $S_{\mu\nu}$ is O\Big{(}\frac{1}{N^{3/2}}\Big{)}. This quantity is of identical order in the Hopfield model despite the correlations in $\{{\boldsymbol{w}}\}$.

When $\mu=\nu$ on the other hand, we find that

[TABLE]

where the $O(\cdots)$ terms in Eq. (36) refer to the effects of the $\{\chi_{ij}\}$ products. Using the properties of the higher-order moments of a correlated Gaussian distribution (e.g., Ref. Tong (2012)), we have by the exchangeability of our distribution for any indices such that different letters refer to different values,

[TABLE]

Substituting these into Eq. (36), we obtain the result that $\mathbb{E}[S_{\mu\mu}^{2}]=O(1/N^{2})$ and hence that $S_{\mu\mu}$ is typically

[TABLE]

namely O\Big{(}1/N\Big{)}. These results imply that the field variance in Eq. (35)

[TABLE]

is

[TABLE]

where the order of the second summation follows from the independence of $\mu$ and $\nu$; by self-averaging it can hence be approximated by

[TABLE]

Finally, recall that $h_{0}=-\frac{1}{N+1}\sum_{\mu=1}^{M}w_{0}^{\mu}\sum_{i=1}^{N}w_{i}^{\mu}s_{i}\approx-\sum_{\mu=1}^{M}w_{0}^{\mu}\eta_{\mu}$. The $\{\eta_{\mu}\}$ decorrelate at the same rate as they do in the Hopfield model; further, $w_{0}^{\mu}$ follows deterministically from $\{w_{1:N}^{\mu}\}$ and adds no information about the state distribution on sites $\{1,\ldots,N\}$ and hence on the $\{\eta_{\mu}\}$. This suggests that the field can be approximated by a limiting Gaussian distribution:

[TABLE]

where $V\triangleq\lim\limits_{N\to\infty}V_{N}$. This Gaussian approximation dramatically simplifies relations in Eqs. (33) and (34), which reduce to

[TABLE]

Considering deletion of any spin $i$ rather than zero, the TAP relation

[TABLE]

follows. Note that the distinction between ${\boldsymbol{J}}$ and ${\boldsymbol{\tilde{J}}}$ disappears at this point; a consequence of the cavity method is that $\{\widetilde{J}_{ii}\}$ terms are disregarded as one expects. To fully specify the TAP relation, $V$ must be determined for the WPE; we turn to this task next.

VII.1.2 Cavitating a ${\boldsymbol{w}}$

Consider now a system of $N$ spins but with $M+1$ $\{{\boldsymbol{w}}^{\mu}\}$. The corresponding Hamiltonian can be related to that of a system with a single ${\boldsymbol{w}}^{\mu}$ removed. Specifically, if

[TABLE]

then the full Hamiltonian is

[TABLE]

where again an irrelevant constant has been dropped. Now following the cavity procedure, we obtain the distribution of $\eta_{0}$ relative to the Boltzmann distribution corresponding to $H^{(M+1)}$ in terms of that corresponding to $H^{(M)}$:

[TABLE]

The expectation $\langle{\eta_{0}}\rangle_{M}=\frac{1}{N}\sum_{i=1}^{N}w_{i}^{0}\langle{s_{i}}\rangle_{M}=0$ by self-averaging due to the stationarity of ${\boldsymbol{w}}^{0}$. The variance $\langle{(\delta\eta_{0})^{2}}\rangle_{M}$ was shown to be O\Big{(}\frac{1}{N}\Big{)} previously; for large $N$, its typical value is

[TABLE]

where

[TABLE]

Finally, assuming that $P^{(M)}(\eta_{0})$ is Gaussian and using the relation in Eq. (41), we obtain

[TABLE]

which is a Gaussian with mean zero and variance

[TABLE]

Note that this value would result for any ${\boldsymbol{w}}^{\mu}$ removed. To obtain the field variance, we use Eq. (38) to obtain

[TABLE]

This resembles that of the Hopfield model but with a changed sign in the denominator.

VII.2 Limiting spectral distribution of ${\boldsymbol{J}}$

We determine the limiting eigenvalue distribution of the ${\boldsymbol{J}}$ matrix. To do so, it suffices to determine the spectral distribution of ${\boldsymbol{\tilde{J}}}$. To see why, first note that $\widetilde{J}_{ii}\to-\alpha$ for all $i$ by the law of large numbers. Recalling that ${\boldsymbol{J}}={\boldsymbol{\tilde{J}}}-\textrm{diag}({\boldsymbol{\tilde{J}}})$, this implies that ${\boldsymbol{J}}$ is asymptotically related to ${\boldsymbol{\tilde{J}}}$ by addition of a uniform quantity to the diagonal, namely

[TABLE]

which in turn means that the eigenvalues of ${\boldsymbol{J}}$ are simply those of ${\boldsymbol{\tilde{J}}}$ translated by $\alpha$.

For Wishart matrices of the form $\frac{1}{M}{\boldsymbol{X}}{\boldsymbol{X}}^{T}$, where ${\boldsymbol{X}}$ are $N\times M$ matrices whose elements are independent zero-mean unit-variance Gaussian variates and $M=\alpha N$, the limiting spectral distribution is known as the Marchenko-Pastur Marchenko and Pastur (1967) law. The fact that the columns of ${\boldsymbol{W}}$ are correlated Gaussian variates seems at first to complicate the determination for $\frac{1}{N}{\boldsymbol{W}}{\boldsymbol{W}}^{T}$. The structure of the specific covariance matrix ${\boldsymbol{\Sigma}}$ considerably simplifies matters, however. Recalling that

[TABLE]

it is apparent that first, ${\boldsymbol{t}}$ is an eigenvector with null eigenvalue, and second, that any vector in the subspace orthogonal to ${\boldsymbol{t}}$ is an eigenvector with eigenvalue $N/(N-1)\to 1$. Hence, any orthonormal set of $N-1$ vectors orthogonal to ${\boldsymbol{t}}$ can be used to represent ${\boldsymbol{\Sigma}}$, and the variation of ${\boldsymbol{w}}$ along each of these eigenvectors is asymptotically of unit magnitude. The procedure for generating ${\boldsymbol{w}}$ is for large $N$ thus equivalent to first generating vector ${\boldsymbol{x}}$ whose first $N-1$ elements are independently $\sim\mathcal{N}(0,1)$ and whose last element is zero and next, transforming ${\boldsymbol{x}}$ by some unitary matrix ${\boldsymbol{U}}$ rotating the $N^{\textrm{th}}$ coordinate vector ${\boldsymbol{e}}_{N}=(0,\ldots,0,1)$ to ${\boldsymbol{t}}$. This implies that the spectral distribution of ${\boldsymbol{W}}{\boldsymbol{W}}^{T}$ approaches that of ${\boldsymbol{X}}{\boldsymbol{X}}^{T}$, where

[TABLE]

consists of an $(N-1)\times M$ matrix $\widetilde{{\boldsymbol{X}}}$ composed of iid normal variates and a final row of zeros; the eigenvalues of ${\boldsymbol{X}}{\boldsymbol{X}}^{T}$ are thus those of $\widetilde{{\boldsymbol{X}}}\widetilde{{\boldsymbol{X}}}^{T}$ with an extra zero added. The limiting spectral distribution of $\frac{1}{N-1}\widetilde{{\boldsymbol{X}}}\widetilde{{\boldsymbol{X}}}^{T}$ can be straightforwardly obtained from the Marchenko-Pastur law by appropriate change of variables. We then obtain the limiting eigenvalue distribution of the matrix $-{\boldsymbol{\tilde{J}}}=\frac{1}{N}{\boldsymbol{W}}{\boldsymbol{W}}^{T}$ to be

[TABLE]

where

[TABLE]

and

[TABLE]

Note that the $\delta$ spike at zero never disappears, a feature that turns out to crucially influence the phase behavior for large $\alpha$. Finally, the spectral distribution of ${\boldsymbol{J}}={\boldsymbol{\tilde{J}}}+\alpha{\boldsymbol{I}}$ follows by reflection and translation:

[TABLE]

where

[TABLE]

and

[TABLE]

The distribution consists of a continuous component over the support $[\lambda_{-},\lambda_{+}]$ and a persistent $\delta$ function at $\alpha$, where of course $\alpha\geq\lambda_{+}$.

VII.3 Energy histogram of the Wishart ensemble

Define the length $M$ vector of normalized state overlaps with the $\{{\boldsymbol{w}}^{\mu}\}$

[TABLE]

so that $H({\boldsymbol{s}})={\boldsymbol{\eta}}^{T}{\boldsymbol{\eta}}$ and

[TABLE]

The derivation of the marginal energy

[TABLE]

is simplified by decomposing the integration into sums of expectations over subsets of states with a constant number of positive elements and exploiting the exchangeability of $f({\boldsymbol{w}})$. Let $N_{+}({\boldsymbol{s}})$ be the number of elements in ${\boldsymbol{s}}$ with value $+1$; for a magnetization $m\in[-1,1]$, $N_{+}=\frac{(1+m)N}{2}$. On the $N_{+}$ constrained subset, the marginal density of $E$ is

[TABLE]

where the sums are over states with $N_{+}$ positive entries. The joint distribution

[TABLE]

by the independence of the columns, but the components of each ${\boldsymbol{w}}^{\mu}$ are correlated. They are however exchangeable, and since the same ${\boldsymbol{s}}$ appears in each term in the sum this implies that

[TABLE]

only depends on $N_{+}$. Consider the specific case of ${\boldsymbol{s}}$ whose first $N_{+}$ elements are 1. Define the sums (one for each column of ${\boldsymbol{W}}$)

[TABLE]

Because $\sum_{i=1}^{N}w_{i}^{\mu}=0$, the sum

[TABLE]

deterministically. By the independence of the ${\boldsymbol{w}}^{\mu}$ the sums $\{A_{+}^{\mu}\}$ are independent random variables and from the properties of linear transformations of Gaussian variables each is distributed according to a zero-mean Gaussian with variance

[TABLE]

When $N_{+}=0$ or $N_{+}=N$, the zero-variance Gaussian is defined to be a $\delta$ function as we expect. We express the expectation (44) at fixed $N_{+}$ as

[TABLE]

Now let

[TABLE]

with density $f_{S}$. We then have

[TABLE]

Because $S$ is the sum of squares of $M$ zero-mean iid Gaussians with variance $\sigma_{N_{+}}^{2}$, it is gamma-distributed, i.e., with density

[TABLE]

with parameters $(k,\theta)=(M/2,2\sigma_{N_{+}}^{2})$. For large $N$ and $m$ relating $N_{+}$ to $N$, $\sigma_{N_{+}}^{2}\approx\frac{1-m^{2}}{8}$. We then obtain the constrained-$m$ energy histogram

[TABLE]

for $e\geq 0$ and zero otherwise and the overall energy distribution

[TABLE]

where the sum runs over the $N+1$ values of $m$ mapping to $N_{+}\in\{0,\ldots,N\}$. As described, for example, in Refs. Mertens (2006); Mezard and Montanari (2009), applying Stirling’s large-$N$ approximation to the binomial coefficients, replacing the sum with an integral and evaluating it with Laplace’s method, we obtain

[TABLE]

or the marginal gamma energy density for the WPE

[TABLE]

VII.4 Annealed approximation

The purpose of the annealed approximation is to provide a simple bound on the typical case behavior of the free energy at leading order in $N$. The free-energy density is defined

[TABLE]

where $\langle\cdot\rangle$ denotes an average over instances. Discontinuities in the free energy describe the phase transitions, and derivatives describe the order parameter(s) and other statistically significant quantities at thermal equilibrium.

It is convenient for this section to consider a model defined by the Hamiltonian

[TABLE]

where $Z_{i,\mu}$ are independent and normally distributed random variables. For the case $\kappa=1$ this Hamiltonian is identical to the main text Hamiltonian up to the choice of the embedding solution ($t_{i}=1$, $\forall i$), inclusion of a nonzero diagonal term in the coupling matrix (which adds a constant offset to the energy), and corrections of $O(1/N)$. These restrictions are for the convenience of analysis, and have no significant impact on the analysis method or conclusions of the section.

The parameter $\kappa$ is useful in making the connection to the anti-Hopfield model (obtained for the case $\kappa=0$) Nokura (1998), and in identifying a variation on the principle of embedding discussed in the main text. Tuning of this parameter allows one to control the energy level of the planted solution at leading order in $N$, allowing an embedding of a ground state with control over the gap at fixed $\alpha$, or implanting an excited state well separated from other stable and metastable states.

The annealed approximation may be used to obtain a lower bound, $f\geq f_{A}$, on the free-energy density

[TABLE]

The physical interpretation for this approximation is that the quenched degrees of freedom ($Z$) are treated on an equal footing with the dynamical degrees of freedom (${\boldsymbol{s}}$). This means that models of lower energy can be selected disproportionately, since spin and model variables can become correlated to lower the free energy. By this process it is possible that atypical models can dominate the free energy so that typical case is not reflected. However, in this ensemble we show that the free energy is correct through much of the phase space in agreement with the TAP analysis of Sec. III.1.

Owing to the quadratic form of the Hamiltonian in the order parameter, it is rather straightforward to take the disorder average explicitly, which yields:

[TABLE]

where $X_{ij}(s)=\frac{1}{N}(s_{i}-\kappa\sum_{i^{\prime}}s_{i^{\prime}})(s_{j}-\kappa\sum_{j^{\prime}}s_{j^{\prime}})$. Following this, we notice that the eigenvalues of ${\boldsymbol{X}}$ are a function only of the sum of spin variables, thus the trace log can be simplified, defining $m=\frac{1}{N}\sum s_{i}$ we can write

[TABLE]

Using Stirlings’s approximation, and using a continuum approximation for $m$ we can write

[TABLE]

Finally, applying Laplace’s approximation yields the result

[TABLE]

The quantity $m$ can be readily associated with the planted state overlap (magnetization, for the case of $t_{i}=1\;\forall i$). The first term is the standard mean-field entropy term, the latter term being an energy term.

Maximizing this equation involves solving for $df/dm=0$, also called the saddle-point equation. At high temperature there is only the solution $m=0$, whereas at low temperature (for large enough $\kappa$) there are two additional solutions. For the case $\kappa=1$ of the main text [Eq. (55)] is identical to the dominant minima of the TAP equation free energy [Eq. (14)].

At low temperature it is interesting to observe that the crystallization transition occurs due to a competition between the energy term and the entropy term, one dominating in each regime. For the paramagnetic solution is defined by $m=0$, whereas for the planted solution $m\approx 1$. Equating these two terms for the case $\kappa=1$ we can find a simple approximation to the first-order transition, the approximation in Eq. (20) proves a very accurate lower bound for the critical temperature shown in Fig. 1 for small values of $\alpha$.

In this section we have found that the annealed approximation recovers the more general TAP result derived in Sec. III.1. For the special case of $O(N)$ planting $\kappa<1$ we find that the embedded state remains thermodynamically dominant at low-enough temperature for large-enough $\kappa$, at a diminished critical temperature (and increased ground-state energy). For the case $\kappa\approx 0$, and at small $\alpha$ the annealed approximation is insufficient to predict all features of the phase diagram, for this we can use the replica method of Sec. VII.5.

VII.5 The replica method

The annealed approximation of Sec. VII.4, and TAP analysis presented (Sec. III.1), are insufficient to describe all the features of our model phase space. To go beyond these we here introduce a replica method. The replica method when solved is able to capture non trivial properties of the exact free energy, at leading order in $N$, in a variety of disordered models closely related to our proposed model Amit et al. (1987); Parisi and Potters (1995); Nokura (1998). As presented here, the replica method is a nonrigorous method for purposes of obtaining insight through related models.

The replica method for the anti-Hopfield model ($\kappa=0$) is demonstrated by Nokura et al. Nokura (1998), and may be derived along similar lines as the Hopfield model Amit et al. (1987). Only minor variations to the form of the Hopfield replica method (free energy) are necessary to analyze the WPE, but these changes have significant consequences.

The free-energy density [Eq. (49)] can be rewritten for purposes of the replica method in terms of a replicated partition function

[TABLE]

where $\mathcal{Z}$ is the partition function, and the trick is to solve for the case of positive integer $n$ and analytically continue to real $n$. For integer $n$ we can write

[TABLE]

using Eq. (50), ${\boldsymbol{s}}^{\rho}$ being a vector of dimension $N$. From here we can follow closely the method in Section 2 of Ref. Amit et al. (1987), which is considered a standard approach. This being understood, we are sparse in our derivation. The main difference in the method is the introduction of an additional order parameter $m$ through the identity

[TABLE]

The integral identity, and scaling with $N$, being appropriate for the large system limit saddle-point, to be later identified. A similar identity is introduced for the overlap $q^{\rho,\rho^{\prime}}=\frac{1}{N}\sum_{i}s_{i}^{\rho}s_{i}^{\rho^{\prime}}$.

The replicated partition function following these manipulations be written

[TABLE]

where ${\boldsymbol{I}}$ is the $n\times n$ identity matrix, $\mathbf{q}$ and $\mathbf{{\hat{q}}}$ are matrices of the same dimension with zero on diagonal (overlap order parameters), and ${\boldsymbol{m}}$, ${\hat{\boldsymbol{m}}}$ and ${\boldsymbol{s}}$ are $n\times 1$ vectors (alignment with planted solution order parameters).

For the special case $m^{\rho}=0$ (or $\kappa=0$) this free energy is identical to that of the anti-Hopfield model Nokura (1998). The interpretation of this result is as follows: since $m$ describes the degree of alignment with the planted solution, we can argue that when there is no extensive alignment with the planted solution, such as at high temperature (above the first-order planting transition) or in general for the space orthogonal to the planted solution, the model phenomena will be identical to the anti-Hopfield model. This is interesting, since in the anti-Hopfield model, and in related classes of problems, there is understood to be a replica symmetry breaking phenomena Parisi and Potters (1995). The space becomes divided into modes separated by extensive barriers at sufficiently low temperature. We thus expect this same phenomena to carry over to our model either as a stable or metastable solution. To determine which we must analyze the free energy within an approximation. In this Appendix we will develop only the replica symmetric theory, this being understood as sufficient for the paramagnetic phase, and the crystal phase discussed in the main text.

In the replica symmetric solution we take $m^{\rho}=m,\forall\rho$ and that $q^{\rho,\rho^{\prime}}=q$ for $\rho\neq\rho^{\prime}$ (the diagonal term is 0). The free energy at leading order can be determined from the saddle-point free energy

[TABLE]

We can attempt to solve this equation as in the annealed case by setting derivatives with respect to $m$, ${\hat{m}}$, $q$ and ${\hat{q}}$ to zero. These equations can be written

[TABLE]

First we should note the solution $m={\hat{m}}=q={\hat{q}}=0$, which is the same (paramagnetic) free energy as in the annealed and TAP analysis, and describes the high-temperature solution. The local stability of this solution at $\kappa=1$, which determines the critical temperature, is identical to that determined by alternative analysis methods (Secs. VII.4 and III.1).

Next consider the subspace with ${\hat{q}}=0$ and $m\neq 0$. In this case we find the relationship $q=(2\kappa-\kappa^{2})m^{2}$, and $q=0$, and recover the planted phase with identical properties to the annealed and TAP analysis. The energy of this solution being controlled by $\kappa$. This solution, where it exists is again locally stable.

Finally, for $\alpha>1$ a solution with nonzero ${\hat{q}}$ is possible in the subspace with $m=0$, the critical temperature given by $T=-1+\sqrt{\alpha}$. Nevertheless, for $\kappa=1$ the embedding-aligned phase ($|m|>0$) is dominant. Thus the replica symmetric solution is in agreement with the understanding presented in the main text. For $\rho<1$ this solution can become relevant for a limited region $\alpha\gtrsim 1$, and describes a transition similar to that in the SK model.

However, it is known that the replica symmetric solution can be unstable, and it is necessary to go beyond these approximations to first and higher orders of replica symmetry breaking to describe the phase space correctly, this being most important at small $\alpha$. We have not probed the interesting consequences in this part of the phase diagram, but the results of the anti-Hopfield model are understood to apply and there is expected to be dynamical and static transitions in the space orthogonal to the planted solution Nokura (1998); Parisi (1995). It is the existence of these solutions, whether they are stable or metastable, alongside the planted one, that offers the interesting possibility to undertake a population transfer measurement in the context of this model Smelyanskiy et al. (2018). Being able to plant a deep stable or metastable solution may also have other interesting applications, particularly in sampling and inference applications where one must provide more than a single ground-state certificate.

The aim of this paper is to provide practical intermediate scale benchmarks, the replica method describes only the typical case properties at leading order in $N$. A tension therefore exists between the results found in this section, and practical limitations of the ensemble such as precision, and instance to instance fluctuations at finite size. These facts must always be borne in mind when considering such analysis results.

VII.6 Time to solution measurements

In this work the problem hardness is quantified via time to solution (TTS) of parallel tempering Monte Carlo. We do emphasize that we expect similar results using other heuristics. For a single disorder realization, the time to solution is defined as the run time such that there is a $99\%$ success probability to have found the solution at the end of the run.

Many solvers have parameters that affect the success probability which we denote by a set $\{\phi\}$. For parallel tempering Monte Carlo, the parameters considered are the lowest temperature and number of replicas. Often, it is faster to carry out many short attempts taking run time $R$ with a lower success probability $p(R,\{\phi\})$ which motivates the definition Boixo et al. (2014); Rønnow et al. (2014); Albash and Lidar (2018)

[TABLE]

The number of attempts is taken to be a real number for convenience.

In the case of an ensemble, averages may be poorly defined and so the median is used Boixo et al. (2014). For the TTS to be well defined, the minimum value with respect to run time and parameters must be computed Rønnow et al. (2014). The ensemble TTS is defined as

[TABLE]

where the index $i$ refers to the TTS of individual disorder samples as a function of run time and parameters.

Due to the nonsequential nature of parallel tempering Monte Carlo, the success probability as a function of run time can be efficiently measured with only a small number of runs Mandrá and Katzgraber (2018): The algorithm is run $W$ times until the solution is found. For a given run time $R$, the success probability is estimated as the percentage of runs where the solution is found in time less than $R$.

Due to finite numerical precision, solutions are considered to be any state with energy $E<E_{\rm GS}+\epsilon$ where $E_{\rm GS}$ is the planted solution energy and $\epsilon$ is a numerical constant. $\epsilon=10^{-7}$ was used for the TTS study in Fig. 10. The effects of different values of $\epsilon$ on the location of the TTS maximum are shown in Fig. 14.