High dimensional optimization under non-convex excluded volume constraints

Several complex systems can be described by a large number of variables whose dynamics is trying to optimize some cost function under a set of constraints. In the most general setting, the cost function is arbitrary and the constraints may be either inequalities or equalities. In this work, we consider the generic situation in which there is a set of $N$ real variables $\underline{x}=\{x_{1},\ldots,x_{N}\}$ and a set of $M=\alpha N$ inequality constraints

each characterized by a function $g^{\mu}(\underline{x})$ and we consider a cost function $H(\underline{x})$.

An example of this situation is found in classification problems in machine learning. In the simplest setting of binary classification, one has a set of inequalities that enforce the correct separation of the training set into two classes. The parameters of the machine (that correspond to the variables $x_{i}$ here) must be adjusted to correctly satisfy these constraints. Therefore, the supervised learning task is recast into a continuous constraint satisfaction problem (CCSP). In many cases one wants to enforce that, once a solution of the problem is found, it satisfies some requirements. For example, in order to prevent numerical instabilities, one would like to have the parameters of the network be quantities that are of the good order of magnitude and do not explode to infinity. To achieve this, one typically uses regularization schemes as, for example, Ridge regularization. In other cases, one is interested in finding sparse networks in which many of the parameters are zero. This can also be achieved by imposing an appropriate cost function on the parameters Neyshabur (2020).

In all these cases, the machine performing the classification task is defined by the form of the functions $g_{\mu}(\underline{x})$ and the regularizing cost function is given by $H(\underline{x})$ which introduces some penalty to prevent that the parameters of the network behave in a way that is not the desired one. In this work, we will be rather generic and we will not focus on any practical application. We will rather focus on the characterization of the landscape of the cost function $H$ subjected to the constraints ${\cal C}$. In particular, we want to find the ground state of the cost function defined by

In order to have a well-posed problem, we will assume that there exists a region of the phase space of $\underline{x}$ such that all constraints are satisfied. We will call this region the SAT region (the terminology comes from constraint satisfaction problems in computer science Franz et al. (2017)). Examples of optimization problems as Eq. (2) can be easily constructed and classified according to the properties of the corresponding cost function. We list here just a set of general cases in which the cost function is:

• •

  separable: this means that each degree of freedom tries to optimize a cost function that is independent of the configuration of the other degrees of freedom. An example in this class is

  $$H(\underline{x})=\frac{1}{2}\sum_{i=1}^{N}(x_{i}-1)^{2}$$

  (3)

  but other non-convex potentials with local minima may be chosen.

• •

  non-separable and convex: the cost function of each degree of freedom depends on the configuration of part or the whole system and the global cost is a convex function when it is considered in the whole phase space (in absence of the excluded volume constraints). As an example one can consider:

  $$H(\underline{x})=\frac{1}{P}\sum_{\mu=1}^{P}x_{i}\eta^{\mu}_{i}\eta^{\mu}_{j}x_{j}$$

  (4)

  being $\underline{\eta}^{\mu}$ i.i.d. random vectors. Another example is the following: consider the vector $\underline{x}$ to be unconstrained. Then one wants to minimize the volume of the phase space of the vector $\underline{x}$, for example by minimizing $|\underline{x}|^{2}$, such that there is a configuration that satisfies all constraints. This problem is related to the packing problem and has been analyzed in Franz and Parisi (2016); Franz et al. (2015, 2017); Yoshino (2018); Brito et al. (2018); Franz et al. (2019a); Ikeda et al. (2019); Ikeda (2020), see also Geiger et al. (2020) for a review with the corresponding connection to deep learning.

• •

  non-separable and non-convex: we have again a non-separable cost which is also non-convex. As an example, one can think of high dimensional random Gaussian functions Castellani and Cavagna (2005)

  $$H(\underline{x})=\frac{1}{N}\sum_{i<j<k}J_{ijk}x_{i}x_{j}x_{k}$$

  (5)

  and the coupling constants are Gaussian variables of the order of $1/N$. In other words, one can superimpose a spin glass landscape on top of a continuous constraint satisfaction problem (CCSP) identified by the set of constraints.

This class of optimization problems is quite vast. In the following, we will consider the simplest case of having a simple separable convex optimization cost and the simplest set of non-convex excluded volume constraints. We will show that the non-convexity in the constraints brings about glassiness and marginal stability even when the cost function is convex by itself.

To be concrete we consider a simple toy model. We assume that each variable $x_{i}$ wants to minimize a simple convex cost function but it is constrained to be in a possibly non-convex region of phase space. We consider the simplest cost function that is

and we want to minimize it under the constraints that

with

where $\underline{\xi}^{\mu}=\{\xi^{\mu}_{1},\ldots,\xi^{\mu}_{N}\}$ is a random vector whose components are Gaussian random variables with zero mean and unit variance. Furthermore, we constrain $\underline{x}$ to be on the $N$-dimensional sphere $|\underline{x}|^{2}=N$. The constraints in Eqs. (7) and (8) define the perceptron continuous constraint satisfaction problem (CCSP). Eqs. (6), (7) and (8) define a quadratic programming problem with inequality constraints Boyd and Vandenberghe (2004). However, as we will see, the topology of the phase space and the nature of the constraints are such that the global optimization problem can become non-convex.

In Gardner (1988); Gardner and Derrida (1988); Franz and Parisi (2016); Franz et al. (2017) the properties of the space of solutions of this CCSP have been extensively studied. As a function of $\sigma$, one has a critical value $\alpha_{J}(\sigma)$ of the density of constraints below which the CCSP is SAT, meaning that with probability one (for $N\rightarrow\infty$) the phase space defined by the constraints has a finite positive volume, and above which the problem is UNSAT, meaning that there is no configuration of $\underline{x}$ that satisfies all constraints.

In Franz et al. (2017) the UNSAT phase of the perceptron has been studied by associating a harmonic energy cost to the unsatisfied constraints. In Franz and Parisi (2016) the SAT phase of the perceptron CCSP has been characterized. The present setting, instead, is quite different: we still consider the SAT phase of the perceptron CCSP, corresponding to $\alpha<\alpha_{J}(\sigma)$, but on top of it we want to optimize the cost function $H(\underline{x})$ of Eq. (6). It is therefore a problem of constrained optimization.

The nature of the constraints changes completely when $\sigma$ passes from being positive, where the CCSP is convex Franz and Parisi (2016), to a negative value, where the solution of the CCSP lies at the intersection of non-convex domains which may be in general non-convex.

In order to study the properties of the landscape of this optimization problem, we consider the partition function defined as

where $\beta=1/T$ is the inverse of the temperature $T$ and $\theta(x)$ is the Heaviside step function. Given the partition function, one constructs the average free energy as

where the overline denotes the average over the random vectors $\xi$s. To compute the properties of the landscape of the cost function, we need to study the behavior of the free energy in the zero-temperature limit $\beta\rightarrow\infty$.

The average of the logarithm of the partition function can be computed using the replica method Mézard et al. (1987)

The average over the $n$-th power of the partition function can be obtained by considering $n$ to be integer and then taking the analytic continuation down to $n\rightarrow 0$. Performing the integral over the random vectors $\xi$s (the details are essentially very close to what is reported in Franz et al. (2017) and we will omit them) we obtain

where $A(Q)$ is an action and $Q$ is an $(n+1)\times(n+1)$ integration matrix. For $N\rightarrow\infty$ we can evaluate the integral by saddle point. The expression of the function $A(\hat{Q})$ is given by

where we have assumed that $\hat{Q}_{1a}=\hat{Q}_{a1}=m$ independently on $a$¹¹1The fact that $Q_{1a}$ is independent on the replica index is a property that holds at saddle point. Indeed one can show that at the saddle point $Q_{1a}=\frac{1}{N}\sum_{i=1}^{N}\overline{\langle x_{i}\rangle}$ which is therefore independent of $a$.. At the saddle point we have the following interpretation for the overlap matrix $Q$:

Therefore $\hat{Q}_{1a}=\hat{Q}_{a1}=m$ is an order parameter that quantifies the overlap between the configuration $\underline{x}_{a}$ and the minimum $[1,1,...,1]$ of $H(\underline{x})$, while $\hat{Q}_{ab}=\frac{1}{N}\underline{x}_{a}\cdot\underline{x}_{b}$, for $a,b>1$, is the overlap between the configurations of different replicas of the system. Note that the spherical constraint on $\underline{x}$ imposes that $\hat{Q}_{aa}=1$. The local partition function ${\cal Z}$ is given by

where we have defined $Q$ as the matrix obtained from $\hat{Q}$ by removing the first row and column. The saddle point equations for the matrix $\hat{Q}$ can be obtained in full generality following the same steps as in Franz et al. (2017). Here we will not repeat the computation and we will consider the replica symmetric solution and its stability.

In order to solve constrained optimization problems, one can rely on the so called Karush-Kuhn-Tucker (KKT) conditions. We consider a dynamical version of such conditions given by the following dynamical system

where $\zeta(t)$ is a Lagrange multiplier needed to enforce the spherical constraint on the variables $x_{i}$. It is simple to show that a fixed point of these equations is a solution of the optimization problem. Indeed if at long times $g_{\mu}>0$, we have $f_{\mu}=0$. Instead if we have $g_{\mu}=0$ we have $f_{\mu}>0$. The variables $f_{\mu}$, which are the Lagrange multipliers of the KKT conditions, are nothing but the forces that the constraints are exerting in order to fix the corresponding variables $g$s to zero. Therefore this system of equations, when it reaches a fixed point, directly gives the statistics of forces corresponding to the $g_{\mu}$ that are identically equal to zero.

This set of dynamical equations can be analyzed exactly in the $N\rightarrow\infty$ limit through dynamical mean-field theory. We first start by considering the case in which the dynamical system is initialized at random for the $x_{i}$ while we will fix $f_{\mu}(0)=1$ which is a convenient choice. Then, to analyze the dynamics in the $N\rightarrow\infty$ limit, we consider the Martin-Siggia-Rose-Jensenn-De-Dominicis dynamical path integral Castellani and Cavagna (2005) defined as

where the dynamical action is given by

and the overline stands for the average over the random patterns. Taking the average over these vectors gives

where we have introduced the dynamical order parameters

As usual, at the saddle point level $D(t,t^{\prime})=0$ by causality. The single gap effective process then gives the equation

where the statistics of the noise $\eta(t)$ is Gaussian and given by

Finally we have

Note that Eq.(33) must be understood as a distributional equation. Furthermore, we have used that at the initial time $t=0$ forces are initialized to $f_{\mu}(0)=1$. This initialization does not affect the endpoint of the dynamics and provides a way to enforce that all forces are either zero or positive. The single-site effective process is instead

where the statistics of the Gaussian noise $\Xi(t)$ is given by

and the kernel $K(t,s)$ is defined as

where the perturbation $P(t)$ is an additive perturbation on the right hand side of Eq. (33). Therefore we obtain the following equations for correlation, response and magnetization $m(t)=\sum_{i}x_{i}(t)/N$

Note that the Lagrange multiplier $\zeta(t)$ can be obtained directly from the equation for $C(t,t^{\prime})$.

From these dynamical equations it follows that

which is a relation connecting the response function $R(t,s)$ to the magnetization $m(t)$.

We have considered the generic setting of optimization problems under non-convex excluded volume constraints. We have analyzed a simple example of this kind in which the cost function is separable and convex and in which the non-convex excluded volume constraints are modeled by a negative perceptron.

We have found that, when the number of constraints is low enough, the cost function has a simple ground state which is captured by the replica symmetric solution of the model. At high density of constraints, instead, the optimization landscape undergoes an RSB transition where minima become marginally stable. Remarkably, we find that the RSB transition point happens before the point in which the accessible phase space defined by the constraints splits into glassy regions. We have also shown how to recover these results from the dynamical mean-field theory of the KKT algorithm. We leave the analysis of more complex cost functions for future work.