Scalable Coupling of Deep Learning with Logical Reasoning

1 Introduction

In recent years, several hybrid neural architectures have been proposed to integrate discrete reasoning or optimization within neural networks. Many of the architectures incorporate an optimization or reasoning layer in a neural network where the previous layer outputs the parameters defining the criteria of the discrete problem [Wang et al., 2019; Amos and Kolter, 2017; Mandi and Guns, 2020; Poganˇci´c et al., 2020; Mandi et al., 2022; Sahoo et al., 2023]. Learning is then achieved by back-propagating gradients from known solutions and natural inputs (such as symbols, images, text or molecule geometries). In this paper, we are more specifically interested in scalable learning when the underlying discrete reasoning problem incorporates unknown logical (deterministic) information or constraints.

Our contributions are twofold: we introduce a hybrid architecture combining arbitrary deep learning layers with a final discrete NP-hard Graphical Model (GM) reasoning layer and propose a new loss function that efficiently deals with logical information. We use discrete GMs [Cooper et al., 2020]

* An extended version of this paper with proofs and additional details can be found at http://arxiv.org/abs/2305.07617.

as the reasoning language because GMs have been used to represent both Boolean functions (in propositional logic, constraint networks) and constrained numerical functions (partial weighted MaxSAT, Cost Function Networks). Also, a probabilistic interpretation of such models is available in stochastic GMs such as Markov Random Fields (MRFs), where infinite costs represent zero probabilities (infeasibility).

Besides the discrete nature of variables that defines zero gradients almost everywhere, the Boolean nature of logic creates an additional challenge. Indeed, constraints are either satisfied or violated and the continuous relaxation of Boolean satisfiability or constraint satisfaction in weighted variants [Wang et al., 2019; Brouard et al., 2020; Palm et al., 2018] may not always account for the potentially redundant nature of constraints: in some contexts, a meaningful constraint c may be implied by a set of other constraints C, making the learning of c from examples impossible if C has already been learned.

In our architecture, during inference, an arbitrary deep learning architecture receives natural inputs and outputs a fully parameterized discrete GM. This GM can then be solved using any suitable optimization solver. To benefit from the guarantees of logical reasoning on proper input, we use the exact GM prover toulbar2 [Allouche et al., 2015]. During learning, GMs computed from natural inputs need to be gradually improved from solutions in the data set. Given the NP-hard nature of discrete GM reasoning and our target of scalable learning, using an exact optimization during learning seems inadequate: even if the finally predicted GMs may be empirically easy to solve (as is the case, e.g., for the Sudoku problem), the GMs predicted in the early epochs of learning are essentially random, defining empirically extremely hard instances that cannot be solved to optimality in a reasonable time [Zhang, 2001]. Relying instead on more scalable convex relaxations of the discrete GM optimization problem [Durante et al., 2022] would come at the cost of sacrificing the guarantees of logical reasoning on proper input.

We instead exploit the probabilistic interpretation of weighted GMs and target the optimization of the asymptotically consistent and scalable Negative Pseudo-LogLikelihood (NPLL) of the training data set [Besag, 1975]. However, this loss function is plagued with an incapability of dealing with large costs [Montanari and Pereira, 2009] and therefore with constraints. In this paper, we analyze the reasons for this incapability and show that it can be explained by the context-dependent redundancy of logical information. We propose a variant, the E-NPLL, for scalable decisionfocused learning. With this differentiable informative loss, our architecture is able to efficiently learn from natural inputs to, e.g., solve the Visual Sudoku [Wang et al., 2019; Brouard et al., 2020]. We also apply it to a many-solutions Sudoku data set [Nandwani et al., 2021] and to a data set of Protein Design instances including several hundreds of variables. With an exact prover used during inference, 100% accuracy can be reached on the symbolic Sudoku, even with small data sets. Moreover, the output of the penultimate layer being a full GM, it can be scrutinized to detect possible symmetries [Lim et al., 2022] for example. It is also possible to introduce a posteriori constraints to get solutions that satisfy desirable properties or even bias the learned criteria by adding new functions to the produced GM.

2 Preliminaries

2.1 Background

A discrete graphical model is a concise description of a joint function of many discrete variables as the combination of many simple functions. Depending on the nature of the output of the functions (Boolean or numerical), and how they are combined and described, GMs cover a large spectrum of AI NP-hard reasoning and optimization frameworks including Constraint Networks, Propositional Logic as well as their numerical additive variants Cost Function Networks and partial weighted MaxSAT [Cooper et al., 2020]. Following cost exponentiation and normalization, these numerical joint functions can describe joint probability distributions, as done in Markov Random Fields (MRFs). In this paper, we use Cost Function Networks for their ability to express both numerical and logical functions. We assume here that cost functions take their value in ¯R = R ∪ {∞}.

In the rest of the document, we denote sequences, vectors, tensors in bold. Variables are denoted in capitals with a given variable Yi ∈ Y being the ith variable in the sequence Y. An assignment of the variables in Y is denoted y and yi is the assignment of Yi in y. Y-i denotes the sequence of variables Y after removal of variable Yi and similarly for y-i given y. For a set of variables M ⊂ Y, we also note Y-M = Y \ M and y-M their values in y. The domain of a variable Yi is a set denoted Di with |Di | ≤ d, the maximum domain size. For a sequence of variables Y we denote as DY the Cartesian product of all Di with Yi ∈ Y. A cost function over a subset of Y described by a tensor (cost matrix) over ¯R is called an elementary cost function.

Definition 1. Given a sequence Y = {Y1, . . . , Yn} of n finite domain variables, a cost function network C is defined as a set of elementary cost functions. It defines a joint cost function, also denoted C(·) = ∑F∈C F, involving all variables in Y. The optimization problem, known as the Weighted Constraint Satisfaction Problem (WCSP), is to find an assignment y that minimizes the joint function C(y). If C(y) <∞, y is called a solution (a model in propositional logic).

In stochastic GMs such as MRFs, the joint function C(·) is used to define a joint probability distribution P C ∝ exp(-C(·)) = ∏F∈C exp(-F) requiring the expensive computation of a #-P-hard normalizing constant.

When a given function F is never larger than another function F ′ (F ≤ F ′ ), F is known as a relaxation of F ′ . A constraint is a cost function F such that F(t) ∈ {0,∞}: it exactly forbids all assignments t such that F(t) = ∞. When F ≤ F ′ are constraints, we say that F is a logical consequence of F ′ : whenever F ′ is satisfied (equal to 0), F is satisfied too. For a set of constraints C, F ∈ C is redundant w.r.t. C iff C and C \ {F} define the same function. At a finer grain, we say F is partially redundant if ∃F ′ < F such that (C \ {F}) ∪ {F ′ } and C define the same function. Consider for example Y = {Y1, Y2, Y3, Y4} with domains {0, 1} and C = {Y1 ̸= Y2, Y2 + Y3 > 1, Y3 ̸= Y4}. No constraint is redundant in C, but in the context of the assignment {Y2 = 1, Y3 = 1}, the constraint Y2 + Y3 > 1 becomes redundant w.r.t. C ′ = C ∪ {Y2 = 1, Y3 = 1} . In the context of {Y1 = 0}, Y2 + Y3 > 1 becomes partially redundant, as it could equivalently be replaced by the weaker Y2 = Y3. Because of redundancies, the observed values in a sample can create a context that makes some constraints redundant and therefore not learnable.

For n variables, a strictly pairwise graphical model C (∀F ∈ C, F involves exactly 2 variables) can be described with n(n - 1)/2 elementary cost function with tensors (matrices) of size at most d2 . We denote by C[i, j] the tensor describing the cost function Fij between variables Yi and Yj . Thus Fij(Yi = a, Yj = b) is simply C[i, j](a, b).

2.2 Problem Statement

In this work, we assume that we observe samples (ω, y) of the values y of the variables Y as low-cost solutions of an underlying constrained optimization problem with parameters influenced by natural inputs ω. From a data set S of pairs (ω, y), we want to learn a model N (in our case, a neural network) which predicts a pairwise graphical model N(ω) such that y ∈ argminy∈DY N(ω)(y). Such a graphical model C = N(ω) defines the last layer of our hybrid neural+graphical model architecture (see Figure 1).

In terms of supervision, we assume that the variables in Y are identified but we also want to exploit any information that would be available on elements of ω. Some of these natural inputs may be direct constraints or assignments of variables in Y that can be directly incorporated into the GMN(ω), others may be known to influence only a subset of all the variables Y. In the symbolic Sudoku problem, a partially filled grid of numbers is observed in ω and each observed value in the grid is known to be the value of its corresponding variable. Similarly, in the visual Sudoku, each image of a digit in the grid is known to represent the value of a single variable, observed in y, providing grounding information [Chang et al., 2020].

Assuming the data set S contains i.i.d. samples from an unknown probability distribution P(Y|ω), a natural loss function for the GM N(ω) is the asymptotically consistent negative logarithm of the probability NLL(S) =

-log(∏(ω,y)∈S P N(ω) (Y = y)) of the observed samples. This negative loglikelihood is however intractable because of the #P-hard normalizing constant. We instead rely on the asymptotically consistent tractable negative pseudo-loglikelihood [Besag, 1975] NPLL(S) = ∑(ω,y)∈S -log(∏i P N(ω) (yi|y-i)). The NPLL works at the level of each variable Yi, in the context of y-i, the assignment of all other variables. It requires only normalization over one variable Yi, a computationally easy task.

Property 1. For Yi, we have P(Yi|y-i) = softmax(-mi(Yi)) where mi(Yi) = ∑j̸=i Fij(Yi, yj). In a message passing interpretation, mi(·) ∈ ¯R|Di | represents the sum of all messages received from neighbor variables Yj through the incident functions Fij, given Yj = yj. Computing the NPLL is in O(n(n - 1)d) per sample and epoch. It can easily be vectorized (computed independently on each variable).

The NPLL enables scalable training from natural inputs assuming an underlying GM NP-hard optimization problem. However, the proof of asymptotic consistency of the NPLL [Besag, 1975; Geman and Graffigne, 1986] relies on identifiability assumptions that obviously do not hold in the context of constraints (zero probabilities) because of redundant constraints. Unsurprisingly, the NPLL is known to perform poorly in the presence of large costs [Montanari and Pereira, 2009]. Empirically, we observed that the resulting architecture completely fails at solving even the simplest symbolic Sudoku problem where ω contains the digits from the unsolved Sudoku grid and y is the corresponding solution.

3 The E-NPLL

To understand the incapacity of the NPLL to deal with large costs, it is interesting to look into the contribution of every pair (ω, y) to the gradient ∂NPLL ∂N(ω)[i,j](vi,vj) of the NPLL for a given pair of values (vi, vj) of a pair of variables (Yi, Yj).

Property 2 (see full paper). The contribution of a pair (ω, y) to the gradient ∂NPLL ∂N(ω)[i,j](vi,vj) is equal to

The two terms in the gradient above come from NPLL terms computed on variables Yi and Yj respectively. Consider our previous example with four variables, C = {Y1 ̸= Y2, Y2 + Y3 > 1, Y3 ̸= Y4} and y = (0, 1, 1, 0). We focus on the variables Yi=2 and Yj=3 and assume that C should hold under ω, which means that the pair of values (Y2 = 0, Y3 = 0) should be predicted as forbidden. Being forbidden under ω, 1(y2 = 0, y3 = 0) = 0. If additionally, the forbidden pairs (Y1 = 0, Y2 = 0) and (Y3 = 0, Y4 = 0) have already reached a high cost under ω, then both P N(ω) (Y2 = 0|y-2) and P N(ω) (Y3 = 0|y-3) will be close to zero, as well as the gradient itself. This will lead to a negligible (if any) change in the cost of the pair (Y2 = 0, Y3 = 0): learning will be blocked or tremendously slowed down. Said otherwise, the fact that, in the context of (ω, y), the forbidden pair (Y2 = 0, Y3 = 0) is redundant w.r.t. already identified forbidden pairs (Y1 = 0, Y2 = 0) and (Y3 = 0, Y4 = 0) effectively prevents any change in the cost N(ω)[2, 3](0, 0). The issue with the NPLL lies therefore in the dynamic of the stochastic gradient optimization: the early identification of some high costs under ω will prevent the increase of other significant costs which are redundant in the context of the observed y, but not redundant in the unconditioned problem.

Inspired by "dropout" in deep learning [Srivastava et al., 2014], we introduce the Emmental NPLL (E-NPLL) as an alternative to the NPLL that should still work when constraints (infeasibilities) are present in S.

Definition 2. The E-NPLL is a stochastic loss defined as E-NPLL(y) = - ∑Yi∈Y log(P N(ω) (yi|y-({i}∪Mi))) where each Mi is a random subset of {1, . . . , n} \ {i}.

The idea of the E-NPLL follows directly from the previous gradient analysis: to prevent a combination of incident functions with already-learned high cost from shrinking gradients, we mute a random fraction of the incident functions. This is used in combination with an L1 regularization on the output of the learned network N(ω) to favour the prediction of exact zero costs which makes the GM optimization problem easier to solve. Because the E-NPLL is designed to fight the side effects of redundant constraints on gradients, we expect it to learn a GM N(ω) with all redundant pairwise constraints.

3.1 Redundancy andMany Solutions

We hypothesize that existing neural architectures where an exact solver is called during training will instead be insensitive to redundant constraints and will tend to not predict them. Contrarily to the NPLL that analyzes N(ω) independently for every variable Yi in the context of y-i, a solver will deal with the global problem and the presence or absence of globally redundant constraints will not affect its output and therefore the associated loss. With L1 regularization, we hypothesize that redundant constraints will tend to disappear. We will test this using the Hinge loss [Tsochantaridis et al., 2005], a wellknown differentiable upper bound of the Hamming distance between a solver solution argminy∈DY N(ω)(y) and the observed y. Note that in our settings, the Hinge loss is equivalent (under conditions detailed in the full paper) to the recent loss of [Sahoo et al., 2023].

When a solver is called during training on (ω, y), the output of the solver is compared to the solution y. But reasoning problems usually have more than one solution and non-zero gradients get back-propagated when the predicted solution is different from the observed one, even if it is a correct one. It has been argued that this can severely hamper the training ability of solver-based layers, and Reinforcement Learning strategies have been defined to fight this [Nandwani et al., 2021]. Because it does not rely on the comparison of an arbitrary solver solution of the underlying discrete problem with the solution y in the data set, our approach should be directly able to deal with problems with many solutions.

4 RelatedWorks

As [Palm et al., 2018; Wang et al., 2019; Amos and Kolter, 2017; Brouard et al., 2020; Poganˇci´c et al., 2020; Sahoo et al., 2023], we assume we have a data set of pairs (ω, y) where y is sampled from a distribution of feasible high-quality solutions of a discrete reasoning problems whose parameters are influenced by ω. We want to be able to predict solutions for new values of ω. While [Brouard et al., 2020] proposed a non-differentiable architecture, most recent proposals, including ours, provide a differentiable DL architecture that enables learning from observables ω including natural inputs. All architectures exploit an optimization model with continuous parameters that are predicted by the previous layer. For training, a discrete differentiable reasoning layer can be used, requiring one [Tsochantaridis et al., 2005; Sahoo et al., 2023; Berthet et al., 2020; Niepert et al., 2021] or two [Poganˇci´c et al., 2020] solver calls per training sample at each epoch. On NP-hard problems, such as MaxSAT or integer linear programming, these architectures may have excruciating training costs and are therefore applied only on tiny instances (20cities traveling salesperson problems [Poganˇci´c et al., 2020] or 3 jobs/10 machines scheduling problems [Mandi et al., 2022]). For more scalability, relaxations of the underlying NP-hard problem using continuous optimization [Amos and Kolter, 2017; Wang et al., 2019; Mandi and Guns, 2020; Wilder et al., 2019] or lifted message passing [Palm et al., 2018] have been used for efficient approximate solving.

For training, the architecture we propose relies instead on a dedicated loss function (that can therefore not be easily changed). At inference, the output of the architecture is a full GM that can be optimized with any suitable optimizer. This implies that our layer is the last layer of the architecture, which is probably a reasonable assumption for most decisionfocused learning problems [Wilder et al., 2019].

In the Predict-and-optimize framework [Elmachtoub and Grigas, 2022; Mandi et al., 2020], a known optimization problem needs to be solved but some parameters v in the criterion must be predicted using historical records of pairs (ω, v). v being available at training, the optimization problem can be solved and the regret (the loss in criteria generated by using predicted instead of true values of v) can be computed and used as a loss. In our case, we only have an unlabeled solution y instead of the optimization parameters v. The reader is referred to the review of [Kotary et al., 2021] for more details on end-to-end constrained optimization learning.

Learning constraint is also central in constraint acquisition [Bessiere et al., 2017], where guaranteed learning of constraints is achieved using positive/negative examples, in interaction with an agent. We instead assume a fixed data set with just positive examples annotated by features ω, a setting where guarantees are essentially unreachable. Our approach is closer to [Beldiceanu and Simonis, 2016], where global (instead of elementary) constraints are learned from solutions of highly-structured problems, without differentiability.

5 Experiments

We test our architecture on logical (feasibility) problems with one or many solutions [Nandwani et al., 2021], ω being purely symbolic or containing images. We also apply it to a real, purely data-defined, discrete optimization problem to check the ability of the E-NPLL to estimate a criteria.

Unless specified otherwise, all experiments use a Nvidia RTX-6000 with 24GB of VRAM and a 2.2 GHz CPU with 128 GB of RAM. Our code is written in Python using PyTorch version 11.10.2 and PyToulbar2 version 0.0.0.2. We use the Adam optimizer with a weight decay of 10-4 and a learning rate of 10-3 (other parameters take default values). An L1 regularization with multiplier 2.10-4 is applied on the cost matrices N(ω)[i, j]. Code and data are available at https://forgemia.inra.fr/marianne.defresne/emmental-pll.

5.1 Learning To Play the Sudoku

The NP-complete Sudoku problem is a classical logical reasoning problem that has been repeatedly used as a benchmark in a "learning to reason" context [Palm et al., 2018; Amos and Kolter, 2017; Wang et al., 2019; Brouard et al., 2020]. The task is to learn how to solve new Sudoku grids from a set of solved grids, without knowing the game rules.

Task

Given samples (ωl , yl ) of initial and solved Sudoku grids, we want to learn how to solve new grids. Sudoku players know that Sudoku grids can be more or less challenging. As one could expect, it is also harder to train how to solve hard grids than easy grids [Brouard et al., 2020]. We use the number of initially filled cells (hints) as a proxy to the problem hardness, a grid with few hints being hard. The minimal number of hints required to define a single solution is 17, defining the hardest single-solution Sudoku grids. We use an existing data set [Palm et al., 2018], composed of single-solution grids with 17 to 34 hints. We use 1, 000 grids for training, and 256 for validation (all hardness). As in [Palm et al., 2018], we test on the hardest 17-hints instances, 1, 000 in total.

A 9 × 9 Sudoku grid is represented as 81 cell coordinates with a possible hint when available. Each cell is represented by a GM variable with domain {1, . . . , 9}. For N, we use a Multi-Layer Perceptron (MLP) with 10 hidden layers of 128 neurons and residual connections [He et al., 2016] every 2 layers. It receives the pairs of coordinates of pairs of cells (Yi, Yj) and predicts all pairwise cost matrices N(ω)[i, j]. Hints are used to set the values of their corresponding variable in N(ω). Performances are measured by the percentage of correctly filled grids; no partial credit is given for individual digits.

The resulting architecture learns how to solve the Sudoku and provides rules in N(ω)[i, j]. For pairs of cells on the same row, column or 3 sub-square, we expect soft differencelike cost function to be predicted (a matrix with a strictly positive diagonal and zeros) that prevents the use of the same value for the two variables. With the L1 regularization on the output of the neural net, other cost matrices should be 0, indicating the absence of a pairwise constraint.

Grids with a Unique Solution

We first train our network with the regular NPLL loss. As expected, it learns only a subset of the rules that suffices to make all other rules redundant: for a cell Yi in the context of y-i, a single clique of difference constraints for every row (or column, or square) is sufficient to determine the value of the cell Yi from y-i, creating vanishing gradients for all other constraints that are instead estimated as constant 0 matrices. On the test set, inference completely fails.

We replaced the NPLL by the E-NPLL, ignoring messages from k randomly chosen other variables. In terms of accuracy, the training is largely insensitive to the value of the hyper-parameter k (see Table 1) as long as it is neither 0 (regular NPLL) nor close to n - 1 (no information). However, larger values of k tend to lead to longer training. We set k = 10 for all Sudoku experiments. In this case, training takes less than 15 minutes. At inference, the predicted N(ω) leads to 100% accuracy on the 1, 000 hard test-set grids.

In Table 2, we compare our results with previous approaches that learn how to solve Sudoku. The Graph Neural Net approach of Recurrent Relational Network (RRN) [Palm et al., 2018], the convex optimization layer SATNet [Wang et al., 2019] and a hybrid ML/CP method [Brouard et al., 2020]. It should be noted that SATNet's accuracy is measured on a test set of easy Sudoku grids (avg. 36.2 hints). All the methods that rely on an exact solver are able to reach 100% accuracy but the modeling flexibility of Deep Learning offers the ability to describe the grid geometry in ω, leading to better inductive bias and far better data efficiency. In our case, with k = 10, we still obtain 100% accuracy on 17-hints grids using a training set of just 200 grids. More epochs are necessary (less than 200) but the training time does not increase (587 ± 32s over 10 training using 10 different seeds).

Interpreting the neural output. In our settings, besides a solution, we obtain a full GM that can be scrutinized and interpreted in terms of which rules have been learned. In the benchmarking situation of the Sudoku problem, we can be confident that the accuracy of 100% observed on the test set actually extends to any Sudoku instance. More interestingly, as shown in [Lim et al., 2022], this output could also be monitored during training to detect symmetries that could speed up training and provide more human-readable information.

Redundant constraints. The rules of Sudoku are naturally described as 810 pairwise difference constraints. It is known that 162 of these (at least) are redundant [Demoen and de la Banda, 2014]. As we expected, the E-NPLL produces the full set of 810 pairwise rules, including redundant constraints. We tried to learn how to solve the Sudoku problem using an embedded exact prover and the Hinge loss [Tsochantaridis et al., 2005]. Embedding an exact prover for an NP-hard problem as a neural layer proved to be challenging even on this small problem. In particular, in the first epochs, the predicted GMs N(ω) are dense random pairwise GMs that are impossible to solve in a reasonable time [Zhang, 2001]. Therefore, we adopted a progressive learning scheme where additional hints are extracted from y, leaving an increasing number of unassigned variables: we used 20 variables initially, increasing this at each epoch until, eventually, only the initial hints were left. This schedule was difficult to adjust. Each training took 2 to 3 days on 1, 000 training grids. Still, it reaches 100% accuracy on the Sudoku test set. Depending on initialization, we observe that several redundant constraints are not learned (see Figure 2), always preserving more than 648 soft difference cost matrices, the conjectured number of non-redundant constraints for Sudoku [Demoen and de la Banda, 2014].

Multiple-solution Grids

Published Sudoku grids have only one solution. We now consider a Sudoku benchmark used in [Nandwani et al., 2021], where each grid has more than one solution. For each grid, the set of solutions is only partially accessible during training: at most 5 of them are present in the training set. The aim is to be able to predict any solution. We use 1, 000, 64 and 256 grids of the data set from [Nandwani et al., 2021] respectively for training, validating, and testing. All hyper-parameters are set as previously, the validation set is only used to stop training. The testing criteria is to be able to retrieve one of the feasible solutions (all of them are known for testing).

Since the E-NPLL never compares a solver-produced solution to the provided solution y, it is not sensitive to the existence of many solutions. With the same training procedure as previously, selecting one of the 5 provided solutions randomly at each epoch, training takes 723.4 ± 64.9 seconds (21.4±1.9 epochs) on average, leading to one of the expected solutions for 100% of the test grids. In fact, since the correct rules are identified, we could verify that thresholding the learned costs into Boolean enables a complete enumeration of all feasible solutions for all instances in the test set.

5.2 Visual Sudoku

Our previous examples show the benefits of exploiting inputs ω. To explore this capacity more deeply, we tackle the visual Sudoku problem in which hints are MNIST images. The goal is to simultaneously learn how to recognize digits and how to play Sudoku. We add a LeNet network [Lecun et al., 1998] to our previous architecture. The GM N(ω) produced is composed of the pairwise cost functions produced by the same Residual MLP with the addition of the negated logit output of LeNet as an elementary cost function on each variable Yi with a hint. This GM is fed to our regularized E-NPLL loss for back-propagating solutions. No new tuning is necessary (same learning rate for both networks).

Our data set is obtained from the symbolic Sudoku data set by replacing hints with corresponding MNIST images, as in [Brouard et al., 2020]. We use the same 1, 000 grids for training, the validation set contains 64 grids (only used to stop training). To check for sensitivity to initialization (an 80% training failure rate was observed for SATNet in [Chang et al., 2020]), 10 runs with different seeds were performed.

After training, we extract the LeNet network alone: it reaches a 97.6% accuracy on MNIST, being indirectly supervised by the provided hints, through the E-NPLL (directly in SATNet [Chang et al., 2020]). We test again on 1, 000 hard grids (see Table 3). When all the hints are correctly predicted, grids are correctly filled, as proper rules have been learned. Moreover, in 8.7% of cases, the solver is able to correct LeNet's errors, leading to an overall accuracy of 76% on hard grids. Training takes an average time of 1150 ± 13s.

We also compare our architecture with SATNet, using their data set of 9K training and 1K easy test grids (average of 36.2 hints [Wang et al., 2019]). We use the exact same parameters as in the previous experiment and reuse SATNet's ConvNet to process MNIST digits. We train for at most 20 epochs (100 for SATNet) using 64 of the training grids for validation (to decide when to stop training). On this data set, SATNet's accuracy is 63.2%. SATNet's authors compared this to a theoretical maximum accuracy of 74.7%, using a 99.2% accuracy MNIST classifier and a perfect Sudoku solver. Integrating LeNet's uncertainty on classification as negated logits in the GM pushes accuracy well beyond this theoretical limit (see Table 4). Since the ConvNet is trained through the ENPLL loss, its weights are automatically adjusted to optimize the joint Pseudo-Loglikelihood that includes also the Sudoku rules being learned. This automatically calibrates its output for the task (similarly to what is done, a posteriori and with known hard Sudoku rules, in [Mulamba et al., 2020]).

5.3 Learning To Design Proteins

The problem of designing proteins has similarities with solving Sudoku [Strokach et al., 2020]. Proteins are linear macromolecules defined by a sequence of natural amino acids. Proteins usually fold into a specific 3D structure which enables specific biological or biochemical functions. Designed proteins have applications in health and green chemistry, among others [Kuhlman and Bradley, 2019]. To design new proteins, with new functions or enhanced properties, one usually starts from an input backbone structure, matching the target functions, and predicts a sequence of amino acids that will fold onto the target structure.

Considering an input protein structure as a Sudoku grid, each amino acid corresponds to a cell and must be chosen among the 20 natural amino acids instead of 9 digits. The structure is predominantly determined by inter-atomic forces, which drive folding into a minimum energy geometry: the most usual approach of the protein design problem is as an NP-hard unconstrained energy minimization problem [Pierce and Winfree, 2002; Allouche et al., 2014]. Inter-atomic forces are influenced by relative distances and atomic natures which implies that the final interactions inside a protein depend on the geometry of the input structure: we will therefore use this geometry in the natural inputs ω, as we did with the (fixed) Sudoku grid geometry. Compared to Sudoku, protein design instances have variable geometry, they may contain several hundreds of amino acids or more and they are subjected to many-bodies interactions that cannot be directly represented in a pairwise GM.

When designing proteins, the Hamming distance between the predicted and observed (native) sequences, called the Native Sequence Recovery rate (NSR), is often used for evaluation. Protein design is a multiple-solution problem: above 30% of similarity, two protein sequences are considered as having the same fold (geometry). So, one given structure can adopt many sequences and a 100% NSR cannot be reached.

For training, we use the data set of [Ingraham et al., 2019], already split into train/validation/test sets of respectively 17, 000, 600 and 1, 200 proteins, in such a way that proteins with similar structures or sequences are in the same set. Similarly to Sudoku, a protein is described by features ω computed on each pair of positions i, j in the protein. They include inter-atomic distance features encoded with Gaussian radial basis function [Dauparas et al., 2022], and a positional encoding of the sequence distance (|i - j|) [Vaswani et al., 2017]. Each backbone geometry ω is associated with y, the sequence of the corresponding known protein. All pair features in ω are processed by a neural network composed of a gated MLP [Liu et al., 2021] that learns an environment embedding from a central amino acid and its neighbors within 10 ˚A, fed to a ResMLP (as for Sudoku) that takes pairs of features and environment embeddings to predict 20 × 20 cost matrices. Training and architecture details are in the arxiv version.

We train the same model with the same initialization using either the NPLL or the E-NPLL loss. To adapt to variable protein sizes, the E-NPLL eliminates k% of incoming neighbor messages. With up to 500 amino acids, the optimization task is challenging at inference and we used a recent GM convex solver [Durante et al., 2022]. As shown in Table 5, we observe that the E-NPLL not only preserves the good properties of the NPLL but actually improves the NSR. While protein design is often stated as an unconstrained optimization problem, we hypothesize that this improvement results from the existence of infeasibilities: when local environments are very tight, they absolutely forbid large amino acids. Such infeasibilities could not be properly estimated by the NPLL alone.

Our architecture provides a decomposable scoring function, such as those used for protein design in Rosetta [Park et al., 2016]. We compared both approaches on the data set of [Park et al., 2016]. Rosetta's results, extracted from [Ingraham et al., 2019], give a 17.9% NSR on this dataset of small proteins. Instead, our E-NPLL learned decomposable scoring function shows an NSR of 27.8%, outperforming Rosetta's energy function. This is all the more satisfying as Rosetta's full-atom function considers all atoms of the protein while we just use the backbone geometry and amino acids identities (as in coarse-grained scoring functions [Kmiecik et al., 2016]).

6 Conclusion

In this paper, we introduce a hybrid neural+graphical model architecture and a dedicated loss function for learning how to solve discrete reasoning problems. It is differentiable and as such, allows natural inputs to participate in the definition of discrete reasoning/optimization problems, providing the ability to inject suitable inductive biases that can also enhance data efficiency, as shown in the case of Sudoku. While most discrete/relaxed optimization layers [Poganˇci´c et al., 2020; Sahoo et al., 2023; Wang et al., 2019] can be inserted in an arbitrary position in a neural net, our final GM layer with the E-NPLL loss offers scalable training, avoiding calls to exact solvers that quickly struggle with the noisy instances that are predicted in early training epochs. It is able to benefit from exact or relaxed solvers during inference. Thanks to the ENPLL, it can simultaneously identify a criterion to optimize as well as constraints. Finally, its output can be scrutinized to check properties and can be a posteriori completed with side constraints or additional criteria, in order to inject further instance-dependent information, that may have been learned, be available as knowledge or as user requirements.

On various NP-hard Sudoku benchmarks, it is able to produce correct solutions from natural input (including images), while being data-efficient and capable of generalization in incomplete multiple-solution settings. The use of an exact prover during inference results in robust prediction, allowing for the correction of otherwise noisy neural predictions. On the real-world protein design problem, our approach is also quickly able to learn a geometry-dependent pairwise decomposed function that outperforms the most recent full-atom pairwise decomposable energy functions for predicting the sequence of natural proteins.

Much remains to be done around this architecture. As for SATNet [Lim et al., 2022], the ultimate N(ω) GM layer of our architecture could be analyzed during training to identify emerging hypothetical global properties such as symmetries or global decomposable constraints, allowing for more efficient learning and improved human understanding. For memory and computational efficiency, we limited ourselves to pairwise models but the use of other languages (e.g. weighted clauses) in replacement of, or addition to, pairwise functions would enhance the capacity of the architecture to capture many-bodies interactions. Another possibility is the use of latent/hidden variables [Stergiou and Walsh, 1999].

Acknowledgements

This work was performed using HPC resources from CALMIP (Grant 2022-P21025), and Jean-Zay GENCIIDRIS (Grant 2022-AD011013779) and has been supported by the French ANR through grant ANR-19-PIA3-0004 (ANITI). M. Defresne PhD studies are funded by the EUR BioEco (grant ANR-18-EURE-0021). We thank B. Savchynskyy for referring us to the Hinge loss, and the proof-reading team from MIAT.