Phase transitions in integer linear problems

TL;DR

This paper explores phase transitions in the solvability of sparse integer linear problems, revealing a Van-Der-Waals phase diagram and increased computational difficulty near critical points.

Contribution

It introduces a phase diagram framework for understanding the computational complexity of integer linear problems based on ensemble analysis.

Findings

Feasible solution space exhibits a Van-Der-Waals phase diagram.

Computational difficulty peaks near the critical point.

Increased complexity observed in the coexistence region.

Abstract

The resolution of linear system with positive integer variables is a basic yet difficult computational problem with many applications. We consider sparse uncorrelated random systems parametrised by the density $c$ and the ratio $\alpha=N/M$ between number of variables $N$ and number of constraints $M$. By means of ensemble calculations we show that the space of feasible solutions endows a Van-Der-Waals phase diagram in the plane ($c$, $\alpha$). We give numerical evidence that the associated computational problems become more difficult across the critical point and in particular in the coexistence region.

Figures

Equations

$$\sum _{j=1}^N{\xi }_{j,\mu }{n}_j=0\mu =1,\cdots ,M.$$

$$\xi_{j,\mu}=\left\{\begin{array}[]{rlc}0&\mbox{with probability }&1-p\;,\\ 1&\mbox{with probability }&p/2\;,\\ -1&\mbox{with probability }&p/2\;.\end{array}\right.$$

$$H=\sum _{\mu }{\left(\sum _j{\xi }_{j,\mu }{n}_i\right)}^2=\sum _{j,k}{J}_{j,k}{n}_j{n}_k,$$

$$J_{j,k}=\sum _{\mu }{\xi }_{j\mu }{\xi }_{k\mu }.$$

$$Z=\sum _{\{n_j\}}{e}^{-\beta H}.$$

$$\overline{Z}=\sum _{\{n_j\}}\overline{e^{-\beta H}}$$

$$=\sum _{\{n_j\}}\prod _{\mu }\int d{m}_{\mu }{e}^{-\beta m_{\mu }^2}\overline{\delta (m_{\mu }-\sum _j{\xi }_{j,\mu }{n}_j)}$$

$$={\left(\frac{1}{2\pi }\right)}^M\sum _{\{n_j\}}\int \int \left(\prod _{\mu =1,M}d{\lambda }_{\mu }\right)\left(\prod _{\mu =1,M}d{m}_{\mu }\right)e^{-\beta |\vec{m}{|}^2+i\vec{\lambda }\cdot \vec{m}}\prod _{j,\mu }\overline{e^{-i{n}_j{\lambda }_{\mu }{\xi }_{j,\mu }}}$$

$$={\left(\frac{1}{4\pi \beta }\right)}^{\frac{M}{2}}\int \left(\prod _{\mu =1,M}d{\lambda }_{\mu }\right)e^{-\frac{1}{4\beta }|\vec{\lambda }{|}^2}\sum _{\{n_j\}}\prod _{j,\mu }(1-p+p\cos ({\lambda }_{\mu }{n}_j))$$

$$\displaystyle=\bigg{(}\frac{1}{4\pi\beta}\bigg{)}^{\frac{M}{2}}\int\left(\prod_{\mu=1,M}d\lambda_{\mu}\right)e^{-\frac{1}{4\beta}|\vec{\lambda}|^{2}}\left(1+\prod_{\mu}\left(1-p+p\cos\left(\lambda_{\mu}\right)\right)\right)^{N}\;,$$

$$\overline{Z}=\sum _{k=0}^N\left(\genfrac{}{}{0ex}{}{N}{k}\right)(1-p{)}^{kM}{\left(\frac{1}{\sqrt{4\pi \beta }}\int d\lambda e^{-\frac{1}{4\beta }{\lambda }^2}{\left(1+\frac{p}{1-p}\cos \lambda \right)}^k\right)}^M.$$

$$e^{kx\cos \lambda }=I_0(kx)+2\sum _{n=1}^{\infty }{I}_n(kx)\cos (n\lambda ),$$

$$\frac{1}{\sqrt{4\pi \beta }}\int d\lambda e^{-\frac{{\lambda }^2}{4\beta }}\cos (n\lambda )=e^{-\beta n^2},$$

$$\overline{Z}\simeq \sum _{k=0}^N\left(\genfrac{}{}{0ex}{}{N}{k}\right)(1-p{)}^{kM}(I_0(kx)+2\sum _{n=1}^{\infty }{I}_n(kx)e^{-\beta n^2}{)}^M.$$

$$\overline{Z}\equiv \sum _{k=0}^N{e}^{N{f}_{\alpha ,\beta ,c}(k/N)},$$

$$f_{\alpha ,\beta ,c}(r)=-r\log r-(1-r)\log (1-r)-\frac{c}{\alpha }r+\frac{1}{\alpha }\log (G_c(r)),$$

$$G_c(r)\equiv I_0(cr)+2\sum _n{I}_n(cr)e^{-\beta n^2}.$$

$$r=(1-r)e^{\frac{1}{\alpha }(G^{\prime }/G-c)},$$

$$f_{\alpha ,\infty ,c}(r)=-r\log r-(1-r)\log (1-r)-\frac{c}{\alpha }r+\frac{1}{\alpha }\log (I_0(cr)).$$

$$r=(1-r)e^{\frac{c}{\alpha }(I_1/I_0-1)},$$

$$\alpha (r)=c\frac{1-I_1(cr)/I_0(cr)}{\log (1/r-1)}.$$

$$\overline{Z}=\sum _{\{n_j\}}\overline{e^{-\beta H}}=$$

$$=\sum _{\{n_j\}}\prod _{\mu }\int d{m}_{\mu }{e}^{-\beta m_{\mu }^2}\overline{\delta (m_{\mu }-\sum _j{\xi }_{j\mu }{n}_j)}=$$

$$\displaystyle=\bigg{(}\frac{1}{2\pi}\bigg{)}^{M}\sum_{\{n_{j}\}}\int\int d\vec{\lambda}d\vec{m}e^{-\beta|\vec{m}|^{2}+i\vec{\lambda}\cdot\vec{m}}\prod_{j\mu}\overline{e^{-in_{j}\lambda_{\mu}\xi_{j\mu}}}=$$

$$\displaystyle=\bigg{(}\frac{1}{4\pi\beta}\bigg{)}^{\frac{M}{2}}\int d\vec{\lambda}e^{-\frac{1}{4\beta}|\vec{\lambda}|^{2}}\sum_{\{n_{j}\}}\prod_{j\mu}(1-p+p\cos(\lambda_{\mu}n_{j}))=$$

$$\displaystyle=\bigg{(}\frac{1}{4\pi\beta}\bigg{)}^{\frac{M}{2}}\int d\vec{\lambda}e^{-\frac{1}{4\beta}|\vec{\lambda}|^{2}}(1+\sum_{s=1}^{L}\prod_{\mu}(1-p+p\cos(s\lambda_{\mu})))^{N}.$$

$$(1+\sum _{s=1}^L\prod _{\mu }(1-p+p\cos (s{\lambda }_{\mu })){)}^N=$$

$$\sum _{k_1+\dots k_L\le N}\left(\genfrac{}{}{0ex}{}{N}{k_1\dots k_L}\right)\frac{1}{(N-{\sum }_i{k}_i)!}(1-p{)}^{M{\sum }_i{k}_i}\prod _{s=1}^L\prod _{\mu }(1+\frac{p}{1-p}\cos (s{\lambda }_{\mu }){)}^{k_s}$$

$$(\frac{1}{\sqrt{4\pi \beta }}\int d\lambda e^{-\frac{1}{4\beta }{\lambda }^2}(1+\frac{p}{1-p}\cos (s\lambda ){)}^{k_s}{)}^M$$

$$e^{k_{s}x\cos (s\lambda )}=I_0(k_{s}x)+2\sum _n{I}_n(k_{s}x)\cos (ns\lambda ),$$

$$\overline{Z}\simeq$$

$$\sum _{k_1+\dots k_L\le N}\left(\genfrac{}{}{0ex}{}{N}{k_1\dots k_L}\right)\frac{1}{(N-{\sum }_i{k}_i)!}(1-p{)}^{M{\sum }_i{k}_i}\prod _{s=1}^L{F}_s^M(k_{s}x)$$

$$F_s(k_{s}x)=I_0(k_{s}x)+2\sum _n{I}_n(k_{s}x)e^{-\beta s^2{n}^2}$$

$$f_L(\mathbf{r})=-\sum _s{r}_s\log (r_s)-(1-\sum _s{r}_s)\log (1-\sum _s{r}_s)+$$

$$+\frac{1}{\alpha }\sum _s\log (F_s(c{r}_s))-\frac{c}{\alpha }\sum _s{r}_s$$

$$f_L(\mathbf{r})$$

Full text

1 Introduction

The use of concepts and techniques from statistical mechanics to analyze hard combinatorial problems has been very fruitful[1, 2], with examples that range from the paradigmatic K-sat[3, 4, 5, 6] to the traveling salesman problem[7] and number partitioning[8] to cite few.

Generically, upon defining an ensemble of random instances of the problem described by some parameters we can analyze how the solutions, and their structure, depend on these parameters. A typical example is a system with constraints connecting different variables. One significant parameter will, then, be the average connectivity of a single variable and another one the ratio between the number of variables and the number of constraints. If we are able to assess the statistics of the number of solutions $N_{sol}$ and compute the average $\langle\log N_{sol}\rangle$ over different instances of the problem, we can adopt this solution entropy, possibly extensive, as a thermodynamic potential and look at singularities in the derivatives with respect to the parameters that can be interpreted as phase transitions among different phases [9, 10].

It turns out that, typically, we can distinguish a region in the parameter space in which solutions to the problem can always be found

• the SAT phase - from a region in which no solution can be found - the UNSAT phase. The associated decision problem, i.e. determining whether a solution exists or not, can be efficiently solved in regions deeply within each phase whereas it becomes effectively hard on the boundary between them [11].

A fundamental problem in combinatorics deals with the resolution of linear systems with non negative integer variables. The latter integrality constraint makes this issue computationally hard but still it inherits some useful results from linear and convex algebra[12]. Several kinds of problems can be assessed: decide if the system has non-trivial solutions (decision problem), describe the space by exhibiting a so-called Hilbert basis[13], count the lattice points inside a convex region[14] or maximize a linear function over integer vectors in a polyhedron (integer linear programming[15]). The latter problem is particularly interesting since it admits a straightforward approximation upon relaxing the integrality constraint and thus defining a linear programming problem, that is efficiently solvable in polynomial time and is connected to an important dual problem[12]. A fundamental strategy in combinatorial optimization thus consists in formulating problems as integer programs. Apart from classical applications in operations research and planning, integer-linear modeling is massively used in the field of constraint-based models of metabolic networks, in which the resolution of such systems is required for the computation of extreme pathways[16], elementary flux modules[17], conserved moieties[18, 19] and thermodynamically unfeasible cycles[20, 21]. We do point out, further, the recent surge of interest in statistical physics of the Gardner problem[22] of calculating the volume of the space of interactions of neural network models for given patterns, as it has been shown that this can be considered as a mean field model of jamming[23] of hard objects, that - in turn - can be connected to the resolution of integer linear systems. Indeed, it has been recently shown that the feasibility of the constraints defining such spaces is ruled by duality theorems defining dual spaces representing unfeasible patterns, that are represented by integer linear systems of equations[24].

Inspired by the aforementioned statistical mechanical approaches, we analyze here the solution spaces of large random sparse integer linear systems. In Sec. 2 we define the model and the ensemble, in Sec. 3 we present the thermodynamical behavior of these problems and their phase diagram obtained by means of annealed and of quenched replica calculations. In Sec. 4, we show how this framework can be used to analyze algorithms apt to solve integer programs and, finally, draw our conclusions.

2 The model

We consider a linear system of $M$ equations in $N$ unknowns $n_{j}$:

[TABLE]

We will consider both the case where the variables $n_{i}$ are non-negative integers and the case where they are Boolean variables, i.e. they can only take the values [math] and $1$. We will consider an ensemble of random systems where the coefficients $\xi_{j,\mu}$ are independent, identically distributed random variables:

[TABLE]

We study systems where $p=\frac{c}{N}$ (with $c$ a intensive constant). These systems are sparse, with a Poissonian structure, where $c$ is the average number of variables per equation. We will analyze large systems, where $N\to\infty$, $M\to\infty$, keeping $\alpha\equiv\frac{N}{M}$ fixed.

We can map the solutions of system (1) onto ground states of a model whose Hamiltonian has the form

[TABLE]

where

[TABLE]

It is interesting to consider the more general problem of calculating the partition function

[TABLE]

For example computing $Z$ in the limit $\beta\to\infty$ gives the number of solutions of (1).

The case where $n_{j}=0,1$ will turn out to be different from the general case where the $n_{j}$ are semi-positive definite integer variables (always with $\xi$ coefficients that can take the three values [math], $+1$ and $-1$).

3 The disorder average

3.1 The annealed approximation for $n_{i}=0,1$, finite T.

We will start by studying the $n_{i}=0,1$ model in the annealed approximation, i.e. by averaging the partition function (and not its logarithm, as we will do in the quenched approach) over the disorder (2). We denote by an over-line the average over the disorder. We have that

[TABLE]

and, by developing the binomial,

[TABLE]

Since $x\equiv\frac{p}{1-p}$ is small, we can set $(1+x\cos\lambda)^{k}\simeq e^{kx\cos\lambda}$ and use the relation

[TABLE]

where the $I_{n}$ are the modified Bessel functions of the first kind of order $n$ ([25]).

Now the Gaussian integral can be evaluated by noticing that

[TABLE]

and it gives

[TABLE]

By using the Stirling approximation for the binomial coefficients and the fact that $p$ is small (that allows to write $-p$ for $\log(1-p)$) we write

[TABLE]

that defines

[TABLE]

where $r\equiv\frac{k}{N}$ and

[TABLE]

The use of the Stirling approximation is justified when $M$ and $N$ are large with $c$ finite (i.e. our approach is valid for diluted models). We have also verified this fact numerically, by exact enumeration. We can now evaluate $\overline{Z}$ in our approximation by a saddle point approximation. We have to find the maxima of $f$ as a function of $r$ for $r\in[0,1]$ ($r$ emerges here as a natural order parameter).

The order parameter $r$ that emerges from the saddle point calculation is the average number of active variables per solution of system (1). Its asymptotic value for $N\gg M$, where every vector variable is a solution, is $\frac{1}{2}$.

By deriving $f_{\alpha,\beta,c}(r)$ with respect to $r$ we obtain a self consistent equation for its stationary points

[TABLE]

that can be solved numerically for each value of $c,\alpha,\beta$ in order to obtain thermodynamic quantities like the free energy ($F=-\frac{N}{\beta}f$), the energy ($E=N\frac{df}{d\beta}$) and the entropy ($S=\beta(E-F)$). In Fig. 1 we show the energy density obtained from these annealed analytic computation for different choices of $\alpha$, $\beta$ and $c$ and from a Monte Carlo simulations (of the original, quenched theory). No unexpected effects are seen, and the agreement is very reasonable.

By taking the limit $\beta\to\infty$ we have

[TABLE]

The self-consistent equation now reads (noticing that $I_{0}^{\prime}=I_{1}$)

[TABLE]

and can be simply solved in $\alpha$. In this way we obtain the $T=0$ equation of state

[TABLE]

In Fig. 2 we can see that for some values of $\alpha,c$ there are three solutions (one on an unstable branch): the curves (16) are similar to the isothermal curves of a Van Der Waals fluid and the equilibrium curve can be obtained by Maxwell construction (Fig. 2, right).

We can calculate the stationary points from $\alpha^{\prime}(r)=0$, whose solutions inserted back into (16) give a coexistence region that shrinks into a second-order critical point. The situation is summarized in the phase diagram in Fig. 2 (bottom): there are two phases with respectively low (UNSAT) and high (SAT) values of $r$ with a coexistence region that shrinks into a second order critical point akin to the Van Der Waals phase diagram of a fluid.

3.2 The annealed approximation for $n_{i}=0,\dots,L$, $T=0$.

Consider the annealed sum, where now $n_{j}=0,1,\dots,L$:

[TABLE]

We can expand

[TABLE]

and we have the integrals

[TABLE]

that, once again, can be solved upon approximating

$(1+\frac{p}{1-p}\cos(s\lambda))^{k_{s}}\simeq e^{k_{s}x\cos(s\lambda)}$, where $x=\frac{p}{1-p}$, and using the Bessel functions formula

[TABLE]

obtaining

[TABLE]

Upon defining the variables $r_{s}=k_{s}/N$ the saddle point $\overline{Z}\simeq e^{Nf_{L}^{SP}}$, where

[TABLE]

it comes from the solution of the following optimization problem:

[TABLE]

For $\beta\to\infty$, $f_{L}$ is symmetric under permutation of the variables. In this limit we search for a symmetric solution $r_{s}=r\leq 1/L\quad\forall s$, and we obtain an extended expression of $f$ for $L>1$:

[TABLE]

Upon taking the maximum, the equation of state now reads

[TABLE]

The picture is qualitatively the same as before, i.e. a Van der Waals picture with first order transition ending in a second order critical point. The values of the first order critical and spinodal lines shift to higher values of $c$ at increasing $L$. In fig 3 we plot the boundary off the coexistence region in the plane $(c,\alpha)$ for several values of $L=1,2,3$: upon increasing $L$ we witness simply to a shift of such a region towards higher values of $c$ and $\alpha$.

3.3 The quenched case for $n_{i}=0,1$, $T=0$

Let us now look to the quenched case, where we assume that the coefficients $\xi_{j,\mu}$ are fixed, and compute the expectation value of the logarithm of the partition function $Z$. We will use replicas, and by introducing $n$ copies of the system labeled by a further index $a$, $\{n_{j,a}\}$, we will compute

[TABLE]

By expanding, as before, with first order modified Bessel functions we get that

[TABLE]

The terms containing at least a cosine function upon integration produce terms proportional to $e^{-\beta}$, that go to zero when $\beta\to\infty$, that is in the limit of interest for us. So we have that

[TABLE]

where $t$ runs over the $2^{n}-1$ possible combinations of $n$ indices. This defines, with the same set of approximation we had used before in the annealed case,

[TABLE]

The saddle point approximation can be computed by solving the optimization problem (with $t=\{1\dots 2^{n}-1\}$) where one maximizes $F_{c,\alpha}({\bf r})$ under the constraints $r_{t}\geq 0$ and $\sum_{t}r_{t}\leq 1$. Both the function $F_{c,\alpha}({\bf r})$ and the domain of the search are symmetric under permutation of variables. We look for a replica symmetric solution, with $r_{t}=r$. The case $n=1$ coincides with the annealed approximation. In the limit $n\to 0$ we have:

[TABLE]

The equation of state $\frac{\partial f_{RS}}{\partial r}=0$ is

[TABLE]

Interestingly, in this case the boundary of the coexistence region can be given in parametric form from the equation $\frac{\partial\alpha(r)}{\partial r}=0$, i.e.:

[TABLE]

where $x\equiv c\,r$.

Figure 4 shows a behavior that is qualitatively very similar to the one of the annealed approximation (see Fig. 2). Upon looking at the shrinking of the coexistence region we can see that $\delta r\propto(c-c_{c})^{1/2}$, i.e. it has a mean field behavior.

We will show how the knowledge of the system thermodynamics discussed so far can be used to analyze algorithms solving integer-linear optimization problems. We will show that the difficulty of performing an optimization task is related to the position of the instance in the phase diagram. In essence, sticking to the Van-Der-Waals metaphor, we can distinguish a purely gaseous ($c<c_{c}$) from a vapor ($c>c_{c}$) phase, where the finite size scaling exponent of the execution time increases and becomes strongly parameters-dependent, in particular reaching a maximum in the coexistence region in correspondence of the transition line.

4 Resolution of single instances

The statistical mechanical analysis of the system, that we have discussed in the previous section, is of paramount importance, since it allows us to qualify the typical behavior of the system, to discuss fluctuations, to examine critical and collective behavior and phase transitions. A second side of the problem is very important, and it is based on analyzing single instances of the problem. We do not consider here statistical averages, but the detailed behavior of a system for a given realization of the coefficients. For example the work of [6] has stressed how useful this can be (in the context, in that case, of satisfiability problems), even for developing new, powerful optimization algorithms. We consider here for simplicity a linear optimization problem (integer programming) defined over the system (1), where we maximize the objective function

[TABLE]

We employ the Matlab© solver intlinprog. We consider two vertical lines in the phase diagram at fixed $c=4,12$ (on different sides with respect to the critical point $c_{c}$) and perform a sweep over $\alpha\in\left[1,5\right]$, $20$ points with fixed number of variables $N=50$ averaged over $10^{3}$ instances. In the figure below we report the relative average length of solutions (left) and the typical machine time (right, quad-core running at $2.7$ GHz) as a function of $\alpha$ for the two cases.

In both cases we see a crossover from an UNSAT to a SAT phase that seems to be steeper for $c=12$. The machine time is almost constant for $c=4$ while for $c=12$ it develops a maximum in correspondence of the coexistence region. Next we studied the machine time as a function of the system size ($N$ number of variables) upon fixing the point in the phase diagram $(c,\alpha)$ (finite size scaling): on two vertical lines ($c=4,12$) we consider $\alpha=1.5,2.5,4$.

In all cases we observe an exponential trend $t=t_{0}e^{BN}$. While for $c=4$ the exponent is small ($B\simeq 10^{-4}$, i.e. we arrive up to $N=10^{4}$) and weakly dependent on $\alpha$, for $c=12$ the exponent is big ($B\simeq 10^{-2}$ we arrive up to $N\simeq O(10^{2})$) and dependent on $\alpha$, with a peak in the coexistence region.

Conclusions

The analysis and resolution of difficult combinatorial problems has gained many insights from statistical mechanics, both in terms of a characterization of general properties of the solution space and in terms of providing new powerful resolution algorithms. A widely studied class of hard problems is the class of integer linear systems, given the wide scope of applications and the theoretical insights that are inherited from linear and convex algebra. In this note we have considered ensembles of sparse uncorrelated random systems parametrised by the average number $c$ of variables per equation and the ratio $\alpha=N/M$ between total number of variables $N$ and total number of constraints $M$. We have shown that the space of feasible solutions endows a Van-Der-Waals phase diagram in the plane ($c$, $\alpha$). We have found that the phase transition scenario depends on the integrality constraint, eventually vanishing for generic positive integers. We gave numerical evidence that associated computational problems become more difficult in regions where a SAT-UNSAT phase transition is present and, in particular, in the coexistence region.

Several further analysis could be of interest. First, we gave a picture for the whole solution space while it would be interesting to consider specific sets of optimization problems upon adding additional terms in the Hamiltonian (3). It is worthwhile noticing that the phase transition scenario that we have depicted is essentially within the universality class of the Ising model at odds with K-SAT problems that endows more complex phase transition scenario [5]. Along this line a possible development would be to test the robustness of the universality class across different ensembles in order to check whether it is related to some universal feature of the problem and at odds with the K-SAT behavior, e.g. in the convexity property. Finally, this statistical mechanical treatment could give relevant insights in algorithmic strategies like message passing and Monte Carlo, whose application in integer linear problems has been already used in a different context[19].

Appendix: Modified Bessel function of the first kind

We report here some useful formulas on modified Bessel function of the first kind from[25] that have been used in the above calculations (analytical and/or numerical).

Representations:

[TABLE]

Generating function:

[TABLE]

Approximation:

[TABLE]

Derivatives:

[TABLE]

For ratios see[26]:

[TABLE]

Acknowledgments

The research leading to these results has received funding from the People Programme (Marie Curie Actions) of the European Union’s Seventh Framework Programme ($FP7/2007-2013$) under REA grant agreement $n[291734]$ (DDM). DDM would like to thank Alexandr Kazda and Davide Raimondo for interesting discussions.