GLinSAT: The General Linear Satisfiability Neural Network Layer By Accelerated Gradient Descent
Authors: Hongtai Zeng
Chao Yang
Yanzhen Zhou
Cheng Yang
Qinglai Guo
Venue: NeurIPS 2024 poster (2024)
Keywords: Differentiable general linear satisfiability neural network layer
Constraint Satisfaction
Accelerated gradient descent
Primary area: optimization
TL;DR: We design GLinSAT based on accelerated gradient descent method, which is the first general linear satisfiability neural network layer where all the operations are differentiable and matrix-factorization-free.

Abstract: Ensuring that the outputs of neural networks satisfy specific constraints is crucial for applying neural networks to real-life decision-making problems. In this paper, we consider making a batch of neural network outputs satisfy bounded and general linear constraints. We first reformulate the neural network output projection problem as an entropy-regularized linear programming problem. We show that such a problem can be equivalently transformed into an unconstrained convex optimization problem with Lipschitz continuous gradient according to the duality theorem. Then, based on an accelerated gradient descent algorithm with numerical performance enhancement, we present our architecture, GLinSAT, to solve the problem. To the best of our knowledge, this is the first general linear satisfiability layer in which all the operations are differentiable and matrix-factorization-free. Despite the fact that we can explicitly perform backpropagation based on automatic differentiation mechanism, we also provide an alternative approach in GLinSAT to calculate the derivatives based on implicit differentiation of the optimality condition. Experimental results on constrained traveling salesman problems, partial graph matching with outliers, predictive portfolio allocation and power system unit commitment demonstrate the advantages of GLinSAT over existing satisfiability layers. Our implementation is available at https://github.com/HunterTracer/GLinSAT.

GLinSAT: The General Linear Satisfiability Neural Network Layer By Accelerated Gradient Descent
Hongtai Zeng'
Chao Yang² Yanzhen Zhou' Cheng Yang²
Qinglai Guo
State Key Laboratory of Power Systems, Department of
Electrical Engineering, Tsinghua University Decision Intelligence Lab, Alibaba DAMO Academy
zenght200mails. tsinghua, edu, cn, xiuxin. www.ii.cooaliii. rocaalibaba-. inc. com, zhouyzh0126.com com, charis yanaeaalannoaaabaaa.e com, Stsigghuu edu. cn
Abstract
Ensuring that the outputs of neural networks satisfy specific constraints is crucial for applying neural networks to real-life decision-making problems. In this paper, we consider making a batch of neural network outputs satisfy bounded and general linear constraints. We first reformulate the neural network output projection prob- lem as an entropy-regulatized linear programming problem. We show that such a problem can be equivalently transformed into an unconstrained convex optimization problem with Lipschitz continuous gradient according to the duality theorem. Then, based on an accelerated gradient descent algorithm with numerical performance enhancement, we present our architecture, GLinSAT, to solve the problem. To the best of our knowledge, this is the first general linear satisfiability layer in which the operations are differentiable and cattii-faattrizztion-rrre that can explicitly perform backpropagation based on automatic differentiation mechanism, we also provide an alternative approach in GLinSAT to calculate the derivatives based on implicit differentiation of the optimality condition. Experi- mental results on constrained traveling salesman problems, partial graph matching with outliers, predictive portfolio allocation and power system unit commitment demonstrate the advantages of GLinSAT over existing satisfiability layers. Our implementation is available at https:: /tte::///tttbt www/HunterTracee......o wra.errrreerrrrrATTT GLinSAT
Despite the fact
Introduction
Constrained decision-maken problems are pervasive across various disciplines. For example, logis- tics companies need to arrange delivery routes minimize transportation costs while ensuring that all orders are delivered on time. Power system operators need to decide how to allocate electricity pro- duction between different power plants to meet ever-changing electricity demand while maintaining system stability. Unfortunately, directly solving these complex constrained decision-making problems via commercial optimization solvers requires a large amount of time. As result, in scenarios that require rapid response, traditional solvers may not be suitable due to their long computation time. With the development of deep learning, it is hopeful that neural networks can capture the domain characteristics and complex relationships involved constrained decision-maie makinn problems through their powerful expressive capability and the solution time can be thus reduced. In recent years, research on how to use neural networks to solve constrained Decision-making problems has become a topic of general interest. Despite the great success of neural networks classification and regression tasks, making the outputs of neural networks satisfy specific constraints is not straightforward, which still needs be further investigated.
Conresponding author.
38th Conference on Neural Information Processing Systems (NcurIPS 2024).

---

A natural idea to impose constraints on the neural network outputs is to penalize the constraint violation in the training stage of supervised learning or reinforcement learning [ 3. However, such an approach requires a careful selection of each penalty coefficient to achieve a balance between decision objectives and constraint violations. As the complexity of constraints increases, choosing appropriate penalty factors may require a large number of attempts, which is time-consuming. Moreover, it is often difficult to theoretically guarantee the boundedness of constraint violations [4], which makes benalty-baset methods less attractive. Ref. [5] managed to determine the width of neural networks required for ensuring feasibility by modeling these networks using binary variables and solving complex mixed-integer bilevel programming. However, this approach necessitates shrinking the original feasible region and can only handle inequality constraints, which limits its broader application. There are also methods in the literature that are better suited for inequality constraint satisfaction. Ref. [6] uses gauge function to map the neural network outputs from a e unit ball into a given polyhedral. Despite its success in the field of control, this method may encounter difficulties in handling equality constraints since the polyhedral need to contain the origin an interior point. Ref. [7] first calculates a reference point within the interior of convex region using convex programming in the offline stage, and then computes a feasible point based on this reference point through simple arithmetic operations in the online stage. Despite its efficiency potentially being affected when constraints are not fixed, the method may encounter difficulties in satisfying equality constraints, as it requires computing the null space of the equality constraints. Since the basis for the null space of a matrix is typically dense, calculating the null space for large matrices may present both efficiency and memory challenges. Another way to encode constraints in neural networks is to reformulate the original problem as a Markov decision problem [8] 2). During the solution process, the decision variables are generated one by one and the value range of the next variable is determined by the value of the current variable so that constraints can be satisfied compulsorily, However, not all decision- making problems can be equivalently converted to Markov decision problems which limits the application of such an approach.
Due to the limitations of the above approaches, many researchers want to use more reliable way to ensure that the outputs of neural networks satisfy specific constraints. A promising way is to integrate optimization solvers as neural network layers. When we embed solver for end-to-end learning, we need to pay special attention to the following two issues: the first one is the supported constraint types, and the second one is the efficiency.
As for the issue of supported constraint types, some frameworks can directly impose combinato- rial constraints on neural network outputs through integrating black-box commercial mixed-integer programming solvers at the cost of inexact gradient estimates and poor utilization of GPUs (since modern commercial solvers are CPU-based) [TO (2). These approaches need to solve com- binatorial optimization problems in both training and inference stages, which is time-consuming. Instead of directly handling the combinatorial constraints, some researchers manage to make neural network outputs satisfy constraints obtained from the continuous relaxation of the original problem, e.g. the widely used double stochastic matrix constraint in Ref. [13 14 [5] solved by Sinkhorn algorithm [6 17 18). Another example is the positive semi-definite (PSD) matrix constraint with unit diagonals as a continuous relaxation of the original MAXSAT problem 19. However, both of the above methods can only express specific constraints, which limits their application. Recently, LinSAT, which is based on a generalized Sinkhorn algorithm, is proposed to impose positive linear constraints on neural network outputs 220 However, the requirement for all elements in constraints to be nen-negative limits the application of LinSAT. For example, even for a simple constraint x S y. namely y < O, LinSAT cannot be used due to the negative coefficient in front of y, which shows the limited expressiveness of LinSAT. Decision variables with negative constraint coefficients occur a lot in real-life peeeisionnmakingg problems, such as the bin packing problem (21), chemical process scheduling [22], power system unit commitment (23). etc.
To deal with general linear constraints in a differentiable way, currently, there are two main approaches, CaxpyLayers [24 and OptNet (25). However, when solving batch of problems, both of them may encounter efficiency issues. Although CvxpyLayers can achieve parallelism through multiprocessing on the CPU, there are only a dozen cores in one CPU, leading to limited parallelism performance. On the other hand, OptNet presents GPU-based batch quadratic programming interior point vorver where batch matrix factorization are performed to accelerate the solution process. Unfortunately, batch matrix factorization may be still a computational bottleneck even when GPU is used. Although some scholars have also studied how to parallelize parts of operations in matrix factorization on the
2

---

GPU [26] 27], the degree of parallelism still highly depend on the structure of the matrix and its elimination tree. Two nodes in the elimination tree can be computed parallel only when there is no direct branch connecting them. As a result, matrix factorization cannot fully utilize the parallel computing ability of the GPU due to the sequential characteristics in the elimination tree (28].
In this paper, we investigate how to apply bounded and general linear constraints to neural network outputs in differentiable way while ensuring batch processing can be performed efficiently on the GPU.
The contributions of this paper include:
To impose general linear and bounded constraints on neural network outputs in a differentiable way, we formulate the corresponding projection problem as entropy-regularizee linear programming where negative logistic entropic regularization terms are added into the objective function. We show that such an entropy-regularee
linear programming problem can be transformed into an unconstrained convex optimization problem with Lipschitz continuous gradient, and thus can be solved by gradient descent based algorithms where matrix factorization operation is required.
2) We design GLinSAT, general linear satisfiability layer to impose linear constraints neural network outputs based on state-of-the-ee accelerated gradient descent algorithm with numerical performance enhancement. Since the main operations GLinSAT is matrix veetor product and no matrix factorization is involved, it is convenient to execute these operations in parallel on the GPU. Although all the operations involved in AT are differentiable which means that we can directly use the automatic differentiation mechanism to perform back propagation, we also provide an alternative way for derivative calculation based on the optimality condition to reduce the memory consumption.
3) We then provide experimental results to demonstrate the capabilities of our proposed method. Experiments on constrained traveling salesman problems, partial graph matching with outliers, predictive portfolio allocation and power system unit commitment show the efficacy of our proposed method and advantages over existing methods. A comparison of methodologies for imposing constraints on neural networks outputs presented in Table A pipeline that shows how our approach works is provided in Fig.
Table Comparison with existing optimizer layers for imposing constraints on the outputs of neural networks
GPU parallel
Explicit back-
Implicit
Matrix factoriza- tion free
Exact gradi- ent
Abbr. of approach
Constraint type
back-
Ref.
comput- ing
propaga- tion
propaga- tion
TO
Perturbed optimizer
combina- torial
X
double
Sinkhorn
stochastic
14
15
matrix
PSD matrix
19
SATNet
with unit diagonals
positive
20
LinSAT
linear
linear and cone
24
CvxpyLayers
X
X
25
OptNet
linear
X
This article
GLinSAT
linear
Note: Explici//Implici backpropagation means that this algorithm performs backward propagation based on automatic differentiation mechanism/inplici differentiation. means that this algorithm feature is dependent on the implementation of the backend solver.

---

# (
arg min f fc,
1xx,xx)
d
al &
Ōx
Dataset
Neural Networks
GLinSAT layer
Loss function
Figure 1: A pipeline that shows how GLinSAT layer works.
Methodology
Sec. 2.1 formulates the neural network output projection problem as an entropy-regularized linear programming by introducing logistic entropy regularization terms the objective function. Based on duality theorem, the original problem can be transformed into unconstrained convex optimization problem with Lipschitz continuous gradient. Then, in Sec. [.2 based on variant accelerated gradient descent method, we design GLinSAT, which solves the projection problem using GPU- friendly algorithm with several numerical enhancements. The corresponding time complexity is also provided in Sec. 2.2 Morcover, although all the operations in the forward pass of GLinSAT are differentiable, in Sec. 2.3 we provide an alternative way based on the optimality condition to calculate the derivatives for memory saving.
2.1 Reformulation of the neural network output projection problem
Here, we want to use a differentiable way to project the output of the neural network c' variables x' E R= that are as similar as possible but satisfy the following constraints (1)
R into
Afa¹ s bj
(la)
A²ax 2 b2
(1b)
AG.² b/2
(1c)
U' x u'
(1d)
where A/ E xn', Ag E Rm''x"', A4 E xn' b E g
b/2 E Rm''s b/3 E 1m1 ul e ga
Moreover, we also suppose that the feasible region in (I) is non-empty.
Apparently, any general linear constraints with bounded variables like I can be converted to the standard form like (2 by shifting bounds and introducing slack variables (see Appendix A.4)
Ax b 0 < x z u
(2a)
(2b)
where m² + mlg 4mg. n m² m A RR b E RT R Here, we denote
the vector obtained from padding (m² + mm/) zeros after the original vector c' as c. Now, the original problem is transformed into problem of projecting E Ra onto E RD that satisfy constraints 2 In the following sections, we mainly focus on such a transformed problem in standard form.
In this paper, we aim for the vector after projection to be as close as possible to the vector prior projection, while adhering to specified constraints. Here, we use the dot product vector similarity. Consequently, our objective function becomes that of maximizing the dot product of vectors and æ. which is equivalently described as minimizing the dot product of c and æ. Besides, as pointed by [29), the optimal solution to an linear programming may not be differentiable (or even continuous) with respect its parameters. Therefore, additional regularization terms need to be included in the objective to make the optimization problem differentiable. Inspired by entropy- regularized optimal transport, here we formulate the projection problem as an entropy-regulariee linear programming to make the entire problem differentiable. Logistic entropy regularization terms
measure of

---

3b as y. the Lagrange dual are added into the objective as follows:
min 05 OSæ5u
(x²
min
c
1
Σ
Ij
log
Zj új
Ij Uj
log 1 Ti Wi
(3a)
0
Uj
s.t. Ax
(3b)
where 0 0 is the inverse temperature parameter that controls the approximation degree between the entropy-ree cetropy-regularized problem and the original linear programming. The regularization coefficient 1/00 controls the smoothness of the outputs. The smaller 1/0 is, the more the outputs tend to be at the extreme point of the feasible region. As 0 + +oo, the optimal solution of the antropy-regularin problem should approach that of the original linear programming.
Remark 1. It is noteworthy that unlike entropy- optimal transport problems where only regularization terms in the form of x log x are involved in the objective, here logistic entropy regularization terms with respect to both z and its complement 1 - x u are added into the objective. Actually, additionally incorporating the complementary entropy regularization terms is the most important part for the derivation of the Lagrange dual problem. Otherwise, we cannot obtain simple mosed-form expression of the dual objective in the following derivation. Similarly, if we use the common e fo-ormm as the regularization term, we cannot obtain an analytical expression either.
we denote the dual variables with respect to the equality constraints function for 3 can be expressed as follows:
9 ( = inf osa5u
(-c" x 1T x log + +1-1-0 w o log 11- #
u =)) y"
Ax
+ b² y (4)
where a b, b represents the element-wise multiplication and division of vector and b respectively.
Since the derivative magnitude log + (1 - x) log (1 - x) tends to infinity when r + 0+ or - 17, the infimum in (4) can be attained only stationary point instead of boundary point. When the infimum in is attained, by making the derivative of the inner function equal to zero, we have:
1 Ou
x
c - A
log
0
(5)
After simplifying the above formula, we can get that when the infimum in 44) attained, the optimal value of (y) can be expressed as:
(v) u (-Ou (-c A A²y))
(6)
where a (.) is the sigmoid function.
Substituting equation 6 into equation
we have:
g (y) = 1 log a (Ou (-c- +by
(7)
Since log a (-) is a strictly concave function, by minimizing the opposite of y (z), we can obtain the following Lagrange dual problem (8) which is exactly an unconstrained convex optimization problem.
min -g (3) = min "ER"
1 lo o (fu (-c - A=) b"
(8)
MER
We can easily show that f (x) is strongly convex function and -9 (3) has Lipschitz continuous gradient (see explanations Appendix(A.5) Therefore, gradient descent based algorithms can be directly applied solve such problem.
2.2 Forward pass in GLinSAT
In the previous section, we have shown that the original entropy-regularizzed linear programming problem (3) can be equivalently converted into an unconstrained convex optimization problem ( with Lipschitz continuous gradient. Theoretically, can be solved readily through gradient descent based method. However, in the actual calculation process, it will be hard to choose suitable step
5

---

size if we just use vanilla gradient descent method. If the step size is much greater than the local Lipschitz constant, the algorithm may diverge. Otherwise, the convergence may too slow.
Considering the strong convexity property of the entropy regularization terms, here we use a variant of accelerated gradient descent method, adaptive primal-dual accelerated gradient descent (APDAGD), which can adaptively approximate the local Lipschitz constant [30]. The detailed procedure of solving the entropy-regularried linear programming problem (3 in GLinSAT is provided in Algorithm |
Compared with the original version of APDAGD, here we improve the numerical performance of smoother way update the approximation
Algorithm T from the following two aspects. First, of the local Lipschitz constant M in GT. In Algorithm I M decreased only when the decrease of the dual objective satisfies the corresponding condition for at least two consecutive times. As result, when M is already a good estimate of the local Lipschitz constant, the frequency of needless updates be reduced, which will lead less computation time. Second, handle the round-off error, we also use a small number á to relax the criterion for the decrease of the objective function. Otherwise, due the existence of numerical error, the criterion of sufficient decrease in objective may be never satisfied. If we do not relax the criterion, M may become a large number and the algorithm will get stuck.
In addition, it is noteworthy that most of the operations involved in Algorithm 1 are calculation of matrix-vector products, vector-vector element wise products and unary functions. Therefore, it is convenient to execute these operations in parallel on the GPU for solving a batch of entropy- regularized linear programming problems.
As for the time complexity of Algorithm[] material of 30.. it can be easily proved that the number of iterations required by Algorithm is roughly proportional to vō and inversely proportional to VE. The corresponding result is given in Corollary | and the detailed discussions can found in Appendix A.6
| based on Theorem and Theorem in the supplementary
Algorithm 1: Solving the entropy hattopy-regularized linear programming problem in GLinSAT
Input: A E Rmxn b Im c R" RT. inverse temerature 0 0, tolerance >>
initial estimate Lipschitz constant (,()), initial estimate of dual variables y (0),
numerical precision 5 > 0
Set k = 0, ()) = (0) 17(0) = (0)) ( x = uoo (-Ou (omoo-ee (-c- y (0)).
g(0) (0) 0, f False;
while ||)x((k) -bl2 > do
Set a (1++) [1+ /1 +1++ 4M(()) (2M(());;
Set (+++1 = 3(k) + (k+1);
Set (A-1) = (k++))/30 /3(+t)))(++1);
Set A 17 (k) + ((+1]) (c(k) 17 (k))
Set (A (1+ = u 1oo (-Ou (-c A (
Set ç (k++1) = (()) - -a (b++) (Ax (x+++)) -b):
Set n (k+1) (A) ++++1] (c(++)) - (1(1));
if (-g g (n) (++1)) (-g (A (1+1)))) - 5 Ax (A(++1)) - -b1l2/ ((2M(x) then
if f = True then
Set (k+11 = (kk))( /2:
else
Set ((++1) = (A);
end
Set x (l++1 (bk) + (k++1) (x (A(++1)) - (A)) f =True;
Set k = k ++11 1;
else
Set (A) (k), = False;
end
end
Output: Optimal primal variables (+), Optimal dual variables (A))

---

Corollary 1. Assume that the optimal dual solution y of problem (3) satisfies 11y*|12 S R. Then, for given tolerance >> 0, the number of required iterations is O
0 max (u) VoR/e).
2.3 Backward pass in GLinSAT
Since all the operations involved in Algorithm are differentiable with respect c, natural idea is to directly use the auto differential mechanism to calculate the derivatives in the backward pass. However, directly backward propagation may require ever growing memory to store computational graphs and may cost much time when the forward pass requires lot of iteration steps. To save the memory usage and accelerate the derivative calculation, we also provide an alternative way based on the optimality condition to calculate the derivatives in GLinSAT.
First, by calculating the derivative of -9 (y), we can obtain the optimality condition as follows:
h (ÿ) = (u o a (-Ou o (-c - y)) b
According to implicit differentiation and chain rule, differentiating equation ( we can get:
dy
0h
äh
(10)
dc
dy
dc
Furthermore, according to equation (6). the derivative of loss function I with respect to a can be calculated as:
0l Øl Oæ Øl dæ dy Ōc dx dc Ex dy dc dy dc dx dc
Ol dy Øl dæ 0l dx
+ ol dy
Oh
-1
Oh
(11)
Zx dy
ay
Ōc
In the actual implementation of GLinSAT, do not explicitly form these Jacobian matrices dr On [hh de .y de
In ad- (see derivations Appendix A.7 we can use conjugate gradient method to calculate the inverse-matrix-eest -vecor products in GLinSAT. There- fore, only matrix-vector product operations are involved in the calculation of derivatives. Moreover, for the sake completeness, GLinSAT, we also implement derivatives with respect to A, b, u for future potential usage. The detailed derivation process of all derivatives is provided in Appendix A.7
0h dc
Instead, we directly form the matrix-vector products du de Sw oy
äl dæ at Oæ
dition, since the jacobian matrix
@h is positive By
3 Experimental Results
In this section, experiments on constrained traveling salesman problems, partial graph matching with outliers, predictive portfolio allocation and power system unit commitment used to demon- strate the advantages of GLinSAT through comparison with the state-of-the-art linear satisfiability layers LinSAT [20]. CvxpyLayers
[4 and OptNet [25]. For OptNet, LinSAT and GLinSAT, the regularization coefficients of nonlinear terms are all set to 1/0. For CvxpyLayers, the projection problem to be solved is set the same 33. The first three experiments originate from Ref. 200 The last experiment is the unit commitment problem in actual power systems. In the fol- lowing sections,
means that GLinSAT is used with dense/sparse matrix and backpropagation is performed using automatic differeettillimpllic differ- ential.
means that LinSAT is used with dense/sparse matrix and maximum iteration number is set to 100/500. The reason we cannot set the maximum iteration number in LinSAT to +00, as we do in GLinSAT, is that LinSAT may iterate endlessly and get stuck in such case. All the experiments are conducted on computer with 24-core Intel(R) Xeon(R) Platinum 8360H CPU and a NVIDIA Tesla A100 GPU through Pytorch 2.2. Our code is provided in https://gitee thu..b........ wwwm//huntee/rracee/l.....
(Denss//parsee
3.1 Constrained traveling salesman problem
Using the traveling salesman problem (TSP) dataset in 20), here we test the performance of each satisfiability layer through experiments on TSP with starting and ending cities constraint and pri- ority constraint respectively. The mathematical formulation TSP with starting and ending cities constraint Eta-SttrtEnd) and TSP with priority constraint (TSP- (TSP-Priority) is provided in Appendix
7

---

Table 2: Average allocated GPU memory and solution time of different satisfiability layers during batch processing of projection and backpropagation when set to 0.1 in TSP training phase
TSP-StartEnd
TSP-Priority
GPU Mem./MB
Time/s
GPU Mem./BB
Time/s
Layer
Proj.
Backprop.
Proj.
Backprop.
Proj.
Backprop.
Proj.
Backprop-
CvxpyLayers
112.1
18.39
116.5
19.94
OptNet
14305
5005
18.92
0.929
14333
5005
20.26
1.136
LinSAT-Dense-100
14977
181.2
0.278
0.417
15009
180.9
0.276
0.418
LinS
74108
181.2
1.323
1.927
74272
180.9
1.317
1.929
GLinnAAT-Dense-FFyyy
4380
4.868
0.382
0.240
4898
4.880
0.440
0.270
GLinSAAT-DDnnee-IIllict
13.35
53.23
0.306
0.143
13.36
53.22
0.349
0.146
LinSAT- AT--ATT--asse0000
7971
132.0
0.281
0.358
8026
132.6
0.286
0.368
LinSATT-Sparse-500
39020
133.5
1.356
1.614
39309
133
1.397
1.645
GLinSAT--Spars--Exppa
2787
5.426
0.603
0.326
3127
5.983
0.652
0.355
62.95
24.71
0.454
0.158
63.27
24.51
0.495
0.165
Note: The GPU memory used by CvxpyLayers is not counted since CvxpyLayers the CPU parallel
mechanism. Statisties of CvxpyLayers and OptNet are based on the first epoch since we cannot obtain well-trained model reasonable time.
A.8
The detailed experimental settings are provided in Appendix A.8
We report the average batch
processing performance in Table 2 where z is set to 0.1. The results when is set to 10-2 are similar therefore we display the results in Table A. and Table A.2 in Appendix AA.8]
method which outperforms all the other methods with minimum total storage and shortest total computation time. We can also find that our proposed Even though we choose the L.iiSSAT-Dens--Expaa
From Table 2 can be seen that GLinSAT-Dense--mme
GLinSAT is memory-efficient
will cost the most memory among all versions of GLinSAT, the total memory usage is still less than 100. We also attempted to set the maximum number of iterations for LinSAT
that of LinSAT-Spanta to +oc, but at this point Li will get stuck and we cannot obtain a reasonable result. Contrarily, for our proposed GT, the convergence is guaranteed, setting the maximum iteration number to +0c will not affect the result.
To obtain feasible tours, we exploit two kinds of post-processing methods in the validation stage 20) The first one is rounding and the second one is beam search where the width of the beam is set to 2048. Table shows the average tour length and feasibility ratio of each method. Since CvxpyLayers and OptNet are hundreds of times slower than LinSAT and GL,nSAT, we are unable to obtain trained models within reasonable time and thus the results are not included.
Table 3: Mean tour length and feasibility ratio obtained from using different and post-processing methods in TSP validation stage
TSP-StarEEnd
TSP-Priority
Rounding
Beamsearch with à = 10"
Rounding
Beamsearch with = 10
with = 1
with = 1
Mean
Feas. Ratio
Mean Length
Feas. Ratio
Mean Length
Feas. Ratio
Mean Length
Feas. Ratio
Layer
Length
LinSAT-Dense-100
4.007
15.3%
3.843
100%
4.114
41.6%
3.952
100%
LinSATT-ensee-500
3.926
93.6%
3.823
100%
4.098
91.5%
3.947
100 %
GLinSAAT-Demss--XXpicc
3,939
94.2%
3.817
100%
4.079
93.5%
3.934
100%
GLinSAT-
3.922
94.2%
3,811
100%
4.068
93.4%
3.927
100%
LinSAT- ATASA--Spars--000
3,843
100%
4.567
100 %
SinSAT-Sparree-555
3,818
100%
4,400
100 %
GLinSAT-Sparsee S-anSATTT--rrr--xpppii
3.939
94.2%
3.817
100% %
4.078 4.073
92.6%
3.935
100 %
BLLinSATT-Sarre-IIplle
3.929
94.6%
3.818
100%
93.4%
3.933
100%
Note: The output of LinSAT-Sparss when 1/0 10 is not real number so that the results are not shown.
8

---

From T we can see that GLinSAT-Dense -enp----eee---iiiii results in the shortest mean tour length when beamsearch is used. It is also noteworthy that LinSAT will produce poor solution when we apply rounding to the results and the max iteration number is set to 100. Although setting the maximum number of iterations to 500 can improve LinSAT's performance, LinSAT"s performance still not as good as GL.nSAT. Considering that the total computation time of LinSAI AT is more than five times that of GLinSAT this time, it makes LinSAT less competitive.
3.2 Partial graph matching with outliers
The detailed mathematical formulation of partial graph matching with outliers is provided in A.9
We earry out experiments on Pascal VOC Keypoint dataset [3]] with Berkeley annotations 32 under the unfiltered setting [20] 33).
Considering there graphs with different sizes in one batch, stack constraints block diagonal matrices and forward them to LinSAT and GLinSAT. However, Cvxpy and OptNet currently cannot handle large block diagonal matrices. Disciplined parameterized programming compilation in CvxpyLayers and matrix factorization of large matrices in OptNet will cost a significant amount of time. Therefore, we can only use for-loop handle batch with different sizes separately. The average GPU memory usage and solution time across different satisfiability layers is provided in Table A.3 of Appendix A.9 In the validation stage, for obtaining feasible integer solutions [20]. We regard the cost of matching pair of nodes as the outputs of satisfiability layers, then use Hungarian algorithm to obtain maximum matching. Finally, we use greedy strategy to preserve pairs with top-p matching scores for constraint satisfaction. The matching Fl scores between graph pairs across various satisfiability layers shown in Table 4 The result of LinSAT-Dense--05 5e0 is not given due to out-of-memory (OOM) issues. According to Table we can find GLinSAT yields the highest F1 scores across all satisfiability layers.
Hungarian algorithm and greedy strategy
Table 4: Mean Fl scores across different satisfiability layers in partial graph matching problem
10
= 10
= 10
10
#
Layer
Mcan F1
Mcan F1
Layer
Mcan F1
Mcan F1
CvxpyLayers
0.616
0.605
OptNet
0.619
0.613
LinSAT-Dense-
0.619
0.614
LinSAT-Dense-500
X
X
0.620
0.620
GLinSAT- Implicit
0.619
0.620
LinSAT-Sparse 100
0.620
0.618
LinSAT-Sparse-500
0.619
0.611
GLinSAT-SSarre-BB Explicit
0.619
0.620
GLinSAT GLinSAT-Sparse
0.620
0.620
3.3 Predictive portfolio allocation
In this section, we use the predictive portfolio allocation dataset in 200 Denote E [0, 1] as the predicted portfolio decision variable of asset i, S as the preferred portfolio asset. Our portfolio allocation needs to maximize the Sharpe ratio 34 while ensuring decision variables satisfy constraints Σi= = 1, Eies Z; where q a pre-defined positive constant. The details of experiments are provided in Appendix A. The average memory usage and solution time is shown in [..leeA..4oo Appendix 10 According Table A..4 can find that our proposed GLinSAT the fastest layer among all layers. In Table we show the mean Sharpe ratio obtained from different satisfiability layers. Our proposed method always yields high Sharpe ratio whether 0 takes 10 10-2.
Table 5: Mean Sharpe ratio obtained from different satisfiability layers in portfolio allocation problem
# 1
1
= 10
$ 10
Layer
S. Ratio
S. Ratio
Layer
S. Ratio
S. Ratio
CvxpyLayers
2.535
2.600
OptNet
2,553
2,381
LinSAT-Dense-
2.245
2.409
LinSAT-Demse-500
2.245
2.409
2,535
2.608
GLinSAT-Dense-Impa
2.535
2.608
LinSAT-
SiSAA--Spare-100
2.245
2.409
LinSAT-Sparse-500
2.245
2.409
GLinSAT-S A-nSAAASSpasse-Eppllle
2.535
2.608
GLinSAT
2.535
2.608

---

3,4 Power system unit commitment
In this section, we carry out experiments about the unit commitment problem on a real provincial power system. In the unit commitment problem, there are hard constraints and soft constraints. Generally, constraints directly related generators are regarded as hard constraints, e.g. the generator logical constraints, generator minimum up-time and down-time constraints. Constraints related to the section power and load balance are usually regarded as soft constraints where violation with large penalty coefficient is introduced in the objective [55 The unit commitment problem can be formulated into mixed-integs integer linear programming (MILP) problem, which is detailed in Appendix
A.II Based one-year power system load data, we first use Gurobi [36] to solve the MILP within 0.1% optimality gap.
After obtaining the integer solution of the unit commitment problem, we then use supervised learning to train neural networks with satisfiability layers so that they can predict the optimal state of a unit while satisfying logical constraints, minimum up-time and down- time constraints. As pointed by [37], when we consider logical constraints and minimum up-time and down-time constraints, these constraints formulate convex hull so that the extreme points of the corresponding feasible region are binary. As a result, we could expect that the outputs of satisfiability layers tend to be binary when 1/6 + 0, thereby making all constraints, including integer constraints, more likely to be satisfied after rounding operations. Once we obtain the predicted integer commitment status of the generators, we can fix the integer variables in the unit commitment problem and solve the corresponding linear programming problem, thereby providing a good initial point for the original mixed-integer programming problem.
Since negative coefficients occur in constraints, LinSAT cannot be used. In Table A.5 of Appendix A.III we compare the performance of batch processing with different layers. When we stack constraints into block diagonal form to exploit parallelism, there are about 1000000 rows and 2000000 columns in the matrix. GLinSAT-Spane Sparse-Implieii is the only way that will not report out-of-memory issues when we use GLinSAT, Both CvxpyLayer and OptNet cannot directly handle such a giant matrix within reasonable time thus we can only use sequential way instead.
We train neural networks with 1 0.1. Table 6 shows the feasibility ratio and average gap on feasible solutions obtained from fixing unit state variables to rounded outputs of neural networks and then solving the continuous unit commitment problem validation stage. Since CvxpyLayers and OptNet are significantly slower than GLinSAT, we are unable to obtain trained models within reasonable time and the results are not included in Table Table shows that if we use sigmoid function to replace the satisfiability layer in training and validation, we cannot obtain any feasible solution. As - O, the feasibility ratio increases. When using GLinSAT with } < 0.0005, the feasibility ratio reaches 100%. In addition, when set ă to exactly zero and solve the resulted projection problem in the form of linear programming (LP) via Gurobi, 100% feasible solutions are also found.
Table 6: Feasibility ratio and average gap obtained from using different 1/0 in validation
Training final layer
Validation final layer
1/0
Feasibility Ratio
Average Gap
Sigmoid
Sigmoid
0%
SLinSAT--Spars--Imme
0.01
86.23%
0.1119%
MLLnnSATTSppars--IImlii
GL.inS Imp-slia
0.005
95.41 %
0.13819
SLinSAAT-Spaase-Implle
GI AT- STpa-s---asisiiiiicci
0.001
98.17%.
0.1109 %
GLiinSAT-Sparse-
0.0005
100 %
0.1114%
GLiuSAT-Sparss -aa-Spaee---piiiiii
GLinSAT-Sparse-Innn
0.0001
100 %
0.1114%
SiiiSAAT-Sparse-Imiii
Gurobi-LP
0
100%
0.1114%
4 Conclusion
In this paper, we reformulate the neural network output projection problem into convex optimization problem with Lipschitz continuous gradient. We then propose GLinSAT, general linear satisfiability layer to impose linear constraints on neural network outputs where all the operations are differentiable and matrix-factorization- GPU. We showcase four applications of GLinSAT and the advantages of our proposed framework over existing satisfiability layers are illustrated.
GLinSAT can fully leverage the parallel computing capabilities of the
10

---

Acknowledgments and Disclosure of Funding
This work was supported in part by the National Natural Science Foundation of China under Grant U22B2097, 52321004 and in part by Alibaba Innovative Research Program. We would like to express our sincerest gratitude to the anonymous reviewers for their insightful feedback on our work. We are also immensely thankful to Wotao Yin for his invaluable support throughout the research process.
References
[1] Nikolaos Karalias and Andreas Loukas. Erdos goes neural: an unsupervised learning framework for combinatorial optimization on graphs. Advances in Neural Information Processing Systems, 33:6659-6672, 2020.
[2] Steven Bohez, Abbas Abdolmaleki, Michael Neunert, Jonas Buchli, Nicolas Heess, and Raia Hadsell. Value constrained model-free continuous control. arXiv preprint arXiv:1902.04623, 2019.
[3] Abhinav Bhatia, Pradeep Varakantham, and Akshat Kumar. Resource constrained deep rein- forcement learning. In Proceedings of the International Conference on Automated Planning and Scheduling, volume 29, pages 610-620, 2019.
[4] Runzhong Wang, Li Shen, Yiting Chen, Xiaokang Yang, Dacheng Tao, and Junchi Yan. Towards one-shot neural combinatorial solvers: Theoretical and empirical notes on the cardinality- constrained case. In The Eleventh International Conference on Learning Representations, 2022.
[5] Tianyu Zhao, Xiang Pan, Minghua Chen, and Steven H Low, Ensuring dnn solution feasibility for optimization problems with convex constraints and its application de optimal power flow problems. The Eleventh International Conference on Learning Representations, 2023.
[6] Daniel Tabas and Baosen Zhang. Computationally efficient safe reinforcement learning for power systems. In 2022 American Control Conference (ACC), pages 3303-3310. IEEE, 2022.
[7] Jesus Tordesillas, Jonathan P How, and Marco Hutter. Rayen: Imposition of hard convex constraints on neural networks. arXiv preprint arXiv:2307.08336, 2023.
[8] Irwan Bello, Hieu Pham, Quoc Le, Mohammad Norouzi, and Samy Bengio. Neural combina- torial optimization with reinforcement learning. arXiv preprint arXiv: 1e..v::611.009000 (1111.099900...11 2016.
[9] Elias Khalil, Hanjun Dai, Yuyu Zhang, Bistra Dilkina, and Song. Learning combinatorial optimization algorithms over graphs. Advances in neural information processing systems, 30, 2017.
[10] Marin Vlastelica Pogančić, Anselm Paulus, Vit Musil, Georg Martius, and Michal Rolinek. Differentiation of blackbox combinatorial solvers. In International Conference on Learning Representations, 2019.
[11] Quentin Berthet, Mathieu Blondel, Olivier Teboul, Marco Cuturi, Jean PhinlPhilippe Vert, and Francis Bach. Learning with differentiable pertubed optimizers. Advances in neural information processing systems, 33:9508-955 2020.
[12] Anselm Paulus, Michal Rolínek, Vit Musil, Brandon Amos, and Georg Martius. Comboptnet: Fit the right np-hard problem by learning integer programming constraints. In International Conference on Machine Learning, pages 1443-8453. PMLR, 2021.
[13] Rodrigo Santa Cruz, Basura Fernando, Anoop Cherian, and Stephen Gould. Deeppermnet: Visual permutation learning. In Proceedings of the IEEE Conference Computer Vision and Pattern Recognition, pages 3949-3957, 2017.
[14] Runzhong Wang. Junchi Yan, and Xiaokang Yang. Learning combinatorial embedding networks for deep graph matching. In Proceedings of the IEEF/CVF international conference on computer vision, pages 3056-3065, 2019.
11

---

[15] Marco Cuturi, Olivier Teboul, and Jean-Philipee Vert. Differentiable ranking and sorting using
optimal transport. Advances in neural information processing systems, 32, 2019.
[16] Richard Sinkhorn and Paul Knopp. Concerning nonnegative matrices and doubly stochastic matrices. Pacific Journal of Mathematicx, 21 2132)33333348 1967.
[17] Marco Cuturi. Sinkhorn distances: Lightspeed computation optimal transport. Advances in neural information processing systems, 26, 2013.
[18] Deeparnab Chakrabarty and Sanjeev Khanna. Better and simpler error analysis of the sinkhorn- knopp algorithm for matrix scaling. Mathematical Programming. 188(1):395 407, 2021.
[19] fo-Wei Wang, Priya Donti, Bryan Wilder, and Zico Kolter. Satnet: Bridging deep learning and logical reasoning using a differentiable satisfiability solver. In International Conference on Machine Learning, pages 6545-6554. PMLR, 2019.
[20] Runzhong Wang, Yunhao Zhang, Ziao Guo, Tianyi Chen, Xiaokang Yang, and Junchi Yan. Linsatnet: the positive linear satisfiability neural networks. In International Conference Machine Learming, pages 06605-36625. PMLR, 2023.
[21] Salma Mezghani, Boukthir Haddar, and Habib Chabchoub. The evolution of the rectangular bin packing problem-a review of research topics, applications, and cited papers. 2023.
[22] Christodoulos A Floudas and Xiaoxia Lin. Mixed integer linear programming in process scheduling: Modeling, algorithms, and applications. Annals of Operations Research, 139: 131-162, 2005.
[23] Bernard Knueven, James Ostrowski, and Jean-Paul Watson. On mixed-integer programming formulations for the unit commitment problem. INFORMS Journal on Computing, 32(4): 857-876, 2020.
[24] Akshay Agrawal, Brandon Amos, Shane Barratt, Stephen Boyd, Steven Diamond, and J Zico Kolter. Differentiable convex optimization layers. Advances in neural information processing systems, 32, 2019.
[25] Brandon Amos and J Zico Kolter. Optnet:
Differentiable optimization as layer in neural networks. International Conference on Machine Learning, pages 136-145. PMLR, 2017.
[26] Steven C Rennich, Darko Stosic, and Timothy A Davis. Accelerating sparse cholesky factoriza- tion on gpus. Parallel Computing, 59: 140-100, 2016.
[27] Meng Tang, Mohamed Gadou, and Sanjay Ranka. A multithreaded algorithm for sparse cholesky factorization hybrid multicore architectures. Procedia Computer Science, 108:616-625, 2017.
[28] Timothy A Davis, Sivasankaran Rajamanickam, and Wissam M Lah-dar A survey of direct methods for sparse linear systems. Acta Numerica, 25:383-566. 2016.
[29] Bryan Wilder, Bistra Dilkina, and Milind Tambe. Melding the data-decisions pipeline: Decision- focused learning for combinatorial optimization. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 33, pages 1658-1665, 2019.
[30] Pavel Dvurechensky, Alexander Gasnikov, and Alexey Kroshnin. Computational optimal transport: Complexity by accelerated gradient descent hetter than by sinkhorn's 's algorithm. In International conference on machine learning, pages 1367-1376. PMLR, 2018.
[31] Mark Everingham, Luc Van Gool, Christopher KI Williams. John Winn, and Andrew Zisserman. The pascal visual object classes (voc) challenge. International journal of computer vision, 88: 303-338, 2010.
[32] Lubomir Bourdev and Jitendra Malik. Poselets: Body part detectors trained using 3d human pose annotations. In 2009 IEFE 12th international conference on computer vision, pages 1365-1372. IEEE, 2009.
12

---

[33] Michal Rolinek, Paul Swoboda, Dominik Zietlow, Anselm Paulus, Vít Musil, and Georg Martius. Deep graph matching via blackbox differentiation of combinatorial solvers. In Computer Vision- ECCV 2020: 16th European Conference, Glasgow, UK, August 23-28, 2020, Proceedinggs, Part XXVIII 16, pages 407-424. Springer, 2020.
[34] William F Sharpe. The sharpe ratio. Streetwise-the the Best of the Journal of Portfolio Management, 3:169-185, 1998.
[35] Jianghua Wu, Peter B Luh, Yonghong Chen, Mikhail A Bragin, and Bing Yan. A novel optimization approach for sub-hourly unit commitment with large numbers of units and virtual transactions. IEEE Transactions Power Systems, 37(5):3716-3722 2021.
[36] Gurobi optimization. [Online]. Available:
com
[37] Deepak Rajan, Samer Takriti, et al. Minimum up/down polytopes of the unit commitment problem with start-up costs. IBM Res. Rep, 23628:1-14, 2005.
[38] Martin WP Savelsbergh. Preprocessing and probing techniques for mixed integer programming problems. ORSA Journal on Computing, 6(4):445-44 44444444,,4, 1994.
[39] Armin Fügenschuh and Alexander Martin. Computational integer programming and cutting planes. Handbooks in Operations Research and Management Science, 12:69-121, 2005.
[40] Tobias Achterberg. Constraint integer programming. 2007.
[41] Tobias Achterberg, Robert E Bixby, Zonghao Gu, Edward Rothberg, and Dieter Weninger. Presolve reductions in mixed integer programming. INFORMS Journal Computing, 32(2):
473-506, 2020.
[42] Ambros Gleixner, Leona Gouwald, and Alexander Hoen. Papilo: A parallel presolving brary for integer and linear optimization with multiprecision support. INFORMS Journal Computing, 35(6):1329 (166:1329-1133 1341,111,,1, 2023.
[43] Runzhong Wang, Junchi Yan, and Xiaokang Yang. Neural graph matching network: Learn- ing lawler's quadratic assignment problem with extension hypergraph and multiple-graph matching. IEEE Transactions on Pultern Analysis and Machine Intelligence, 44(9):526155299, 2021.
[44] Defu Cao, Yujing Wang, Juanyong Duan, Ce Zhang. Xia Zhu, Congrui Huang, Yunhai Tong, Bixiong Xu, Jing Bai, Jie Tong, et al. Spectral temporal graph neural network for multivariate time-series forecasting. Advances in neural information processing systems, 33:17766-17778, 2020.
13

---

A Appendix
A.1 Related Work
Constraints handling paradigm in neural networks for decision making. Some simple constraints can be directly encoded by neural network activation functions, e.g., ReLU function for non negative constraints, sigmoid function for bounded constraints, softmax function for sum-to-one constraints. However, it is difficult for neural network outputs satisfy complicated constraints by only using these activation functions. For a few problems with special structures, the constraints can be directly handled by the well-designed action space of sequential decisions in reinforcement learning [8] ( However, as the complexity of constraints increases, it will be hard to design a suitable action space. Considering the diffieulty for agents find feasible solutions through random exploration, the resulting sparse rewards may also lead to slow convergence. Another common way to deal with constraint violation is to penalize such violation in the training stage. Ref. incorporates the augment the reward function in reinforcement learning with the sum of the constraint violation penalty weighted by the Lagrange multipliers. However, it cannot be guaranteed that the constraints can be always satisfied by directly penalizing the constraint violation, see Ref. 4 Ref. [ 7 are better suited for inequality constraint satisfaction. These methods may encounter difficulties in satisfying equality constraints, either in terms of efficiency or expressiveness. Another way to handle constraints is incorporate optimization solvers into neural network layers, which is detailed in the next paragraph.
constraint violation in the loss function of supervised learning. Ref. (2) 3
Optimizers neural network layers for constraint satisfaction. To integrate optimizers into neural networks, it is necessary to calculate derivatives with respect to the parameters and perform batch processing efficiently. Ref. [10 [ [2] exploit black-box solvers impose combinatorial constraints on decision variables. However, only approximated gradients can be obtained through perturbations on the problem. Ref [19) relax the original combinatorial constraint into the positive semi-definite matrix constraint with unit diagonals and integrates GPU-based MAXSAT solver into neural network layers. Ref. 13 [4 15), Sinkhorn algorithm [6 17 18 is used to make neural network outputs satisfy the double stochastic matrix constraint, which is linear relaxation of permutation, matching and sorting constraints. Since Sinkhorn algorithm only involves iterative normalization of rows and columns in matrices, it is straightforward to parallel the operations on GPUs and calculate the derivatives based on automatic differential mechanism. Ref. [25] presents OptNet, neural network layer that integrates a GPU-based batched quadratic programming solver, qpth. The forward pass exploits a primal-dual interior point method to find the solution and the derivatives are calculated based on the matrix factorization obtained from the forward propagation. Ref. (24) further presents Cvxpyl Layers to incorporate convex programming into neural network layers. Although CvxpyLayees can represent more optimization problems, it relies CPUs for parallelism which may lead to efficiency issues if a large batch of optimization problems needs solving. Ref. [20) designs LinSAT, differentiable layer encode the positive linear constraints based on a multi-set Sinkhorn algorithm. Although such an algorithm is easy to be parallel on GPUs, it can only imposing positive linear constraints on neural network outputs.
A.2 Broader Impacts
This paper is aimed at making the outputs of neural networks satisfy bounded and general linear constraints. The proposed framework GLinSAT can be used for end-to-end neural network training while ensuring the feasibility of neural network outputs, providing promising approach for applying neural networks to decision-makinn makinng problems. There no foresceable negative societal consequence that is a direct result of the proposed method.
A.3 Limitations
In this paper, we manage to impose bounded and linear constraints on neural network outputs. Considering that real-life decision variables often have finite upper and lower bounds simultancously. our proposed method can actually be applied to a lot of decision-making problems. However, for variables with one-sided boundary or explicit boundary, our method cannot be directly used. A possible workaround is to manually calculate the implicit bounds of these variables through domain 42). In addition, it should be noted that currently our algorithm can propagation, see [38 39 40
14

---

only deal with linear constraints. In the future, we need to conduct further research on neural network layers that can efficiently handle cone constraints on GPUs.
A.4 Reformulation of general linear constraints with bounded variables into standard form
Denote c' RR" as the output of an neural network. In this section, consider projecting the output of the neural network x' E TR' into variables x' E R that satisfy the following linear constraints and bounded constraints:
A/ - bj
(A.la)
A/a' 2 b/2
(A.1b)
A/j²' b/2
(A.lc)
1' S x S u'
(A.1d)
where A E R xn by RT,, 4 E Rm Xn' b/2 RT A/ E Rm/x"" b/S E RT, 1 R" u _2"".
Obviously, we can convert all inequality constraints into equality constraints by introducing bounded slack variables as follows:
A/ s b1
(A.2a)
AGao' 8/2 b2
(A.2b)
A/z²' b/3
(A.2c)
I' S x < u
(A.2d)
0 < 6 8
(A.2e)
0 0 8 $2
(A.2f)
where s bj -
A u' + A 1' positive and negative parts of matrix , Al respectively.
A
-
A
u.
8
by, A A,2 and Ai A²_ are the
We use the following notation:
1'
A
I
b'
A ALi
u'
A
A
I
b
by
0
A
b/g
Then, the original problem is transformed into problem of projecting Ra'+ +m²s onto variables nR''-m; +mg that satisfy the following linear constraints and bounded constraints:
Ax b
(A.3a)
0 < x: S u
(A.3b)
A.5
Property of the primal and dual objective function
We first show that the primal objective function f (z) is strongly convex, where f (x) is:
- = - cTx
n
Inj
log
Jj + (1 Inj uj
log (1 Fj
U, ))
(A.4)
Uj
The second order derivative of (x) can be expressed as follows:
v2 (x: diag
(A.5)
x (u -x)
where diag (-) maps vector to its corresponding diagonal matrix, represents the element-wise product, = represents the element-wise reciprocal of vector z.
Since we have 0 < x < u. we have x o (u - x) 2 2 (xz) > @ max I, which means f is a m strongly convex function.
fu u. As a result, we have
15

---

We next show that the gradient of the opposite dual objective function V (-9 (3)) is Lipschitz continuous, where -9 (y) is:
9 (3) = 1 -01'"oog ( (-c - A²))) b"
(A.6)
The first order derivative of -g (y) can be expressed follows:
V ( (v)) Ax (z) b
(A.7)
where x (y) a (~Ou (-c y))
Using the first order derivative of -9 (y). we have:
V (-9 v (-9 (12)) = || (ar (y1) - x (92)) S || | (31) - x (32)l12 (A.8)
According to equation the strong convexity of f (a:), we have:
> (y1) (32))12 - (Vf (x (y1)) Vf (122))) (o (3/1 x (y22)
i max (u)
(-c+
1 x (v1)
u x (y1) (-c+ 1
x (y2)
(or (3/1) x (322) (A.9)
log
ou u x (22)
log
Ou
= (1111 (x: (81) (y2)) - || ll 1
Vall2 x (311 x (22)ll2
which implies:
1 max (u)² ||All2
| (yi) x (32))l2
| y2ll2
(A.10)
Combining equation
A.8 and A. 10 we have:
max (u)² |
(- (vi) - V (J2))ll2 S
11 32212
(A.11)
4
which means that V (-g (v)) is
||) Lipschitz continuous.
A.6 Time Complexity Algorithm
According Theorem 2 in the supplementary material of [30], the no-n of (k b at the k-th iteration is bounded by 16LR/k², where L is the Lipschitz constant of the gradient of the dual objective function, R is the upper bound of l2-norm of the optimal dual variables. Since have
shown that V (-)( is m)x(()" | Lipschitz in Appendix A.5
we can see that:
40 max (u)²
2 A | R
Ax (N)
b11
(A.12)
k2
Let Axx(k) -bl12 =E, We can obtain the upper bound of outer cycle iteration number as follows:
OR
t: S 2 max (w) | V
(A.13)
According to Theorem in the supplementary material of (30 the number of all the iterations after an iteration & is O (k). Therefore, the time complexity of Algorithm is O maxx (u) ||A]|, v AR/E As result, for optimization problems that with similar upper bounds of dual variables, the required iteration number is approximately proportional v0/e. Such time complexity is better than the complexity of sinkhorn-based algorithms, which the iteration number is approximately proportional to 0/22, see [18 20
16

---

A.7 Derivative Calculation in the backward pass of GLinSAT
The optimality condition can be expressed as follows:
h (y) = A (u (-Ou o (-c - A y))) b 0 Let A or b or u. according chain rule, have:
(A.14)
01 Øl dæ
Øl dx dy
01 dy
Øl dæ
01 Ōæ
01
[hh
0h
(A.15)
du dæ dv dæ dy du
dy dr
dæ de
Üx dy
dy
Ky
where æ (y) uoo Ou (-c- 4)y)).
We can first calculate
8z, as follows:
By,
Fzg Oz (Ug Ig) Apg
(A.16)
Oyp
By writing the above equation into matrix form, we have:
dx
= diag o (u x))A
(A.17)
dy
In the actual computation process, the matrix Instead, directly formulate the derivatives follows:
dn explicitly formulated for saving GPU memory. ay 81 Ow Th TApa"
81 the Dw Bbp"
81
On
al
8z. in A.15 as
Dcg"
Oug
Tw Dyp
01 dx
Øl 0zg
Ø1
Oun (wn - x.g) Up
(A.18a)
dz 0A pq
0xg AAp
0xg
Øl Ox
(A.18b)
Ox Obp
Øl dæ
Ol dzg
Ø1 0xg
0rg (ug xzg)
(A.18c)
Ox dcg
Deg dcg
(-ca- Σ
n
5q Org (ug - Izg)
Yi Aig
Øl Ōx
Ø1 Oxg
Øl
(A.18d)
Dz 0ug
0zg 0rg
D..
Uy
al Ox
Øl 0xg
"
Ol Grg (ug
Σ Dxg
Apo
(A.18e)
dx Our
Ozg Oyp
By writing the above equations into matrix form, we have:
Øl dæ
y ( Øl
fx (u x;)
(A.19a)
dx 0A
dæ
Øl dx
=0 Dæ Ub
(A.19b)
0l dx
0x o (u - x))
(A.19c)
Oa dc
dr
0l dx Øl
@ Ox (u 1 x) - A²
(A.19d)
Oz du
dr
0l dx
Øl
A
Ox o (u - x))
(A.19e)
dz dy
da.
We then calculate
0hg Zy.
as follows:
0h p Oyg
A Ozz / (uj xz) Agi
(A.20)
17

---

By writing the above equations into matrix form, we have:
Dh Oy
Adiag (Ox (u (u)+))= A
(A.21)
where diag (.) maps a vector its corresponding diagonal matrix. When is the optimal solution, we have < x < u. Therefore, Adiag (fx o (u - æ)) A is positive semi-definite matrix.
Here, we denote
At dr ch By By 8t
Hh
as
at oh
After we finish the calculation of 5h by conjugate
di
gradient method, we can calculate
8h
øt Sbp
ah
dh
ah Brg.
follows:
Sh
JApq
Deg"
Ol öh
0l Oh,
Ø1
m
Øl
Iq + Up
Aig
Bxg (uy - xzg)
(A.22a)
Ōh 0.pp
th, 0Ap
ahy
==1
ol 0h
01 Sh,
0l
(A.22b)
Öh dbv
Thp 0bp
@hp
Ol 0h
m
0l ah,
m
0l
Aig
Oxg (wg - xzg
(A.22c)
@h. dcg
Th, dcg
Ohg
Iq
Oxg (ug - Ia0)
(-co
m
YiA
Øl 0h
Øl Oh
Ol
Aig
(A.22d)
Ōh dug
3h dug
Oh,
2g
By writing the above equations into matrix form, we have:
Øl Öh
Øl xx + y ((Ar Øl
(fx (u (u=x))
(A.23a)
Ōh 0A
0h
0h
Ol 0h
Øl
(A.23b)
0h 0b
Oh
0l 0h
AT Ø1
(fx (u x²)
(A.23c)
0h dc
0h
ØL 0h
T Ø1 0h
x Ax o (u - x (-c ATy)
(A.23d)
0h. du
Finally, by substituting equations (A.19) and (A.33) gradient as follows:
into (.155)
we can obtain the corresponding
0A = y
0l
((
01
A 0l
o (fx o (u - =))
T 0l
T
(A.24a)
8x
0h
Oh
01
Øl
(A.24b)
0b
Ōh
01
0l
AT 0l
(u
o (Ox o (u - x))
(A.24c)
Dc
Ox
Th
01
Øl
Øl
Ox (u x) (-c A Ayy)
AT
(A.24d)
Ou
Fx
0h
A.8 Experimental details about constrained traveling salesman problem
Here, we first provide the mathematical formulation TSP with starting and ending cities constraint (TSP-StartEnd) and TSP with priority constraint (TSP- (TSP-Priority). We first focus on Startt,r,End, which can be modeled as follows:
n
n
min (0.1)))
Dis
X
(A.25a)
1=1
A=1
s.t.
Xx11 = 1, Xe,n = 1,
(A.25b)
x 1, = 11,
X1, = 1n
(A.25c)
18

---

where x, = 1 means city i is the j-th visited city in a tour, D..j refers to the distance between city 2 and city J.
On the basis of problem (.255). the TSP-Priority problem can be obtained by introducing the following priority constraints:
m++1
Σ
x
(A.26)
j=1
where p is the city that needs to be visited in the first m steps.
Due to the focus of this experiment on comparing the performance of each satisfiability layer, designing better network structure is beyond the scope of this paper. Therefore, we directly used the SOTA network structure in solving TSP, which consists of a three-layers Transformer, followed by three-layer MLP with ReLU activation
20) The hidden sizes all set to 256 and the number of multi-head attention is set to 8. All the training and test data is consisted of 20 nodes that are uniformly sampled from unit square. For all satisfiability layers, the constraints are the continuous relaxation of the original TSP problem and all the common settings are the same as follows. The learning rate is set 10a4. The batch size is set to 1024. When we stack constraints into block diagonal form to exploit parallelism, there are about 40,000 rows and 400,000 columns in the whole matrix. The training epoch number is set 50. The constraint tolerance is set 10 For AT, the initial estimate of Lipschitz constant is set to 0, the initial estimate of the dual variable is set to zero vector, the numerical precision is set to 10emachine. where Emachine is the machine epsilon. In each epoch, we generate 256000 random cases as the training set. In the training stage, the objective function in (A.25a) is used as the loss function. In the validation stage, we use two kinds of post-pnn procpescsssing methods. The first one is rounding. The second one is beam search where the width of the beam is set to 2048. All the experiments are conducted on a computer with a Intel(R) Xeon(R) Platinum 8360H CPU and NVIDIA Tesla A100 GPU with 80GB memory through Pytorch 2.2.
Here, we additionally show the training performance when 1/0 is set to 10 follows:
Table A. Average allocated GPU memory and solution time of different satisfiability layers during batch processing of projection and backpropagation when 1/0 is set to 10 in TSP training phase
TSP-StartEnd
TSP-Priority
GPU Mcm./BB
Time/s
GPU Mem./MB
Time/s
Satisfiability layer
Proj.
Backprop.
Proj.
Backprop.
Proj.
Backprop.
Proj.
Backprop.
CaxpyLayers
107.3
31.12
107.4
30.22
OptNet
14305
5005
18.89
0.930
14333
5005
20.40
1.125
L.nSAT--ense-100
14977
181.2
0.278
0.417
15009
180.9
0.277
0.419
LinSAT-Dense-500
74108
181.2
1.339
1,930
74272
180.9
1.318
1.932
GLinSATT-Dese-EXplicc
11095
4.852
0.976
0.552
11155
4.857
0.965
0.555
Gl GlinSSAT--Dense-Immliict
13.34
53.23
0.888
0.072
13.36
53.22
1.096
0.078
LinSAT--Spars--100
LipS-AT-Sparse-500
GLinSATT-Spans--EXppii
7085
5,051
1.552
0.774
7007
5.519
1.508
0.812
SLinSSATTSpaars--Impiic
62.95
24.71
1.355
0.078
63.27
24.51
1.721
0.083
Note: LinSATT-Deene/SSprree000000533
means that LinSAT is used with dense/sparse matrix and max (DensneseDDeeepppssse(eeEepp
means that GLinnST is used with dense/sparse matrix and backpropagation performed using automatic differential/mpplici differential. The GPU memory used by CvxpyLayers not counted since CvxpyLayers CPU-based.
iteration number is set to 100/500.
Note: The output of LinSAT-Sparse when 1/0 = 0.01 is not a real number so that the results are not shown.
From Table A.I we can see GLinS -----ense-Imlieee is the most memory efficient satisfiability layer among all these layers. As to the solution time, we can find that GLinSAT is slightly slower than LinS - terminates due to reaching the maximum number iterations rather than because of con- vergence. Actually find that LinSAT often reports warnings like "non-zero constraint violation
The reason LinSAT has a fast calculation speed is that the algorithm
19

---

within max iterations". which indicates the algorithm has not converged. After we increase the max iteration number 500, the number of warnings has decreased, but there are still some warnings that indicate the algorithm has not converged. When we set maximum iteration number to +00, the algorithm progress will stuck. Compared with LinSAT, our proposed GIT is more reliable and is guaranteed to converge. Table shows the corresponding validation results about the mean tour length and feasibility ratio. From Table A.2 to that of LinSAT.
we can see that the performance of GLinSAT is superior
Table A.2: Mean tour length and feasibility ratio obtained from using different 1/0 and post- processing methods in TSP validation stage
TSP-StartEnd
TSP-Priority
Rounding
Beamsearch with 1/0
Rounding with 1/0 10
Beamsearch with 1/0
with 1/0 10x
0.1
0.1
Mean Length
Feas. Ratio
Mean Length
Feas. Ratio
Mean Length
Feas. Ratio
Mean Length
Feas.
Layer
Ratio
LinSAT-Denns--10 100
8.090
100% %
7.031
100%
LinSAT-Densee5500
3.873
100%
4,011
100%
GlinSAAT-Denss-Exppiccc
4.030
92.7%
3.869
100 %
4.171
91.5%
3,983
100%
GLinSAT--Dnns-Immplii
3,930
94.1%
3.790
100%
4.053
94.6% %
3,924
100 %
100
LinSATTSpparse-500
LLLinSAT-SSaars--Explee
4.022
92.6%
3.859
100%
4.078
92.6%
3,935
100%
GLinSAT-SSarse-IIppal 3.926
94.8%
3,803
100%
4.073
93.4%
3.933
100%
Note: The output LinSAT-Sparse when 10 not a real number so that we cannot obtain any trained
model. The output of LinSAT- Aenee when } 10 not real number so that results are not shown.
A.9 Experimental details about partial graph matching with outliers
Here, we first provide the mathematical formulation of partial graph matching with outliers. Denote m, n as the number nodes of two graphs respectively. The partial graph matching problem with p inliers can be expressed as follows:
X S 1,
(A.27a)
X1, S 1m
(A.27b)
1s X1, =
(A.27c)
E (0, 13 muxn
(A.27d)
where X1 means the i-th node in the left graph matches the j-th node in the right graph.
When we use AT as the satisfiability layer, it is necessary to canonize the original inequality constraints. By introducing bounded slack variables into equations (A.27a) and (A.27b). we can
reformulate constraints (A.27 into the standard form as follows:
X'1n + sn 1,
(A.28a)
X1, + tm 1m 1 X1, =p
(A.28b)
(A.28c)
E (0, 13 mxn
(A.28d)
where 0,, < s, 5 On < tm S 1m.
In the training stage, we follow the experimental codes in Ref. 20 where the neural net- works trained on the basis of pretrained SOTA graph matching NGMv2 model 43 named "pretrained params_vgga 2arp11____m______ afat---vocc. Different satisfiability layers are used to make the outputs satisfy the continuous relaxation of constraints (A.27).
is noteworthy that since the sizes of graphs differ a lot one batch, we stack constraints into block diagonal forms in LinSAT and GL to exploit parallelism of the GPU. However, it is difficult
20

---

for CvxpyLayers and OptNet to directly handle large block diagonal matrices since disciplined parameterized programming compilation and matrix factorization of large matrices will cost large amount of time. Therefore, we can only use sequential way to handle batched graphs with different sizes for CvxpyL.aerrs and OptNet. The batch size is set to 128 across all the experiments. When we stack constraints into block diagonal form to exploit parallelism, there are about 2,500 rows and 13,000 columns in the whole matrix. The constraint tolerance is set to 10 In the training stage, binary cross entropy loss is used as the Ioss function. For experiments about OptNet and GLinSAT-Explicit, we use double-precision precision floating-poin point numbers pass while reporting warnings like "Returning an inaccurate and potentially incorrect solution". GLinSAT-Expln Explicit will not encounter numerical issues in its forward pass. However, sometimes the gradient calculated by auto differential may be not a real number. We believe the problem is that single-precision floating-point numbers amplify the cumulative error of backpropagation. Noating-poin numbers, everything works fine. As a result, we use
caating-point numbers during projection. If single- used, OptNet will encounter numerical issues in its forward
When we use kingleprreission floating-poin numbers on all the other layers except OptNet and Gl.inAAT-Explictt.
In the validation stage, we use Hungarian algorithm and greedy strategy for post-processing [20]. We can regard the cost of matching pair of nodes as the outputs of satisfiability layers, then use Hungarian algorithm to obtain a maximum matching. Finally, we can use greedy strategy to preserve pairs with p-highest matching scores and obtain the solution.
All the experiments are conducted on computer with a Intel(R) Xeon(R) Platinum 8360H CPU and NVIDIA Tesla A A100 GPU with 80GB memory through Pytorch 2.2. For GLinSAT, the initial estimate of Lipschitz constant set to 0, the initial estimate of the dual variable set to zero vector, the numerical precision is set to 10c where macchine is the machine epsilon. Here, we show the average memory usage and the solution time of different satisfiability layers in Table A.3
Table A.3: Average allocated GPU memory and solution time of different satisfiability layers during batch processing of projection and backpropagation in the training phase of partial graph matching
10
10
GPU Mcm./MB
Time/s
GPU Mcm
Timo/s
Layer
Proj.
Backprop.
Proj.
Backprop.
Proj.
Backprop.
Proj.
Backprop-
CvxpyLayers
64.26
16.09
12.76
15.07
OptNet
167.3
826.2
4.290
3.712
167.2
826.3
4.437
3.726
LinSAT-Dense-1000
23944
322.1
0.611
3.887
23944 X
322.1
0.611
3.887
LinSAT-Dense-500
X
X
X
X
X
GLinSATT-Dense-EEplicc
1249
0.997
1.891
4.418
1311
1.052
1.991
4.147
GLinSATT-Dense-Impliici
129.0
736.5
1.333
3.995
129.0
736.4
1.302
3.515
LinSATT-Sparse-100
803.0
397.7
1.593
3,906
803.6
397.7
1.615
3,877
ASnSSAT-Sparss-50000
2772
304.3
2.922
5.534
2772
304.3
2.923
5.292
648.7
172.9
1.794
4.435
687.6
139.8
1.986
4.620
Gl GLLinSAT-SSprrse-Impiiii
0.001
862.1
1.544
3.884
0.001
862.1
1.449
3.865
Note: The GPU memory used by CvxpyLayers is not counted since CvxpyLayers is CPU-based.
Note: The symbol "*" means the outputs of this satisfiability layer cannot meet constraints. The symbol "x" means the layer leads to out-of-memory (OOM) issucs.
In Table A.3 although LinSAT-100 seems to be the fastest method, the price is that the required constraints are not satisfied at all. When we set the maximum iteration number of LinSAT to 100, almost every batch LinSAT will report warnings like "non-zero constraint violation within max iterations". When we set the maximum iteration number to 500, LinSAT-Dense will soon run out of memory while LinSAT-Sparss will still display the warning message in almost every epoch. Notably, the computational speed of LinSAT this point has already fallen behind that of GLinSAT. If we set the maximum iteration number to +oo in LinSAT-Sparse the algorithm will get stuck. Compared with our proposed GL,nSAT, when the maximum iteration number is set to +20, the convergence of GI GL is guaranteed. In summary, we can conclude that our proposed GI.nn is the fastest satisfiability layer while ensuring the outputs satisfy the linear and bounded constraints from the results shown in Table
21

---

A.10
Experimental details about predictive portfolio allocation
Here, we first restate the mathematical formulation of predictive portfolio allocation. We denote 2s as the predicted portfolio decision variable of asset i and 8 as the preferred portfolio asset. The portfolio allocation needs to maximize the Sharpe ratio [34] while satisfying the following constraints:
Σ I4 1, Σ 24 29
(A.29a)
iES
c T E 1 zi<n
(A.29b)
where is a pre-defined positive constant. Following the codes in 20 here we set as 0.5 and set C [AAPL,, MSFT, AMZN, TSLA, GOOGL). We also the first 120-day historical data train neural network that could maximize the Sharpe ratio for the future 120 days where StemGNN 44 are used as the network backbone to extract the features.
When we use GLinSAT as the satisfiability layer, it is necessary to canonize the original inequality constraints, By introducing bounded slack variables into constraints (A.29a), can reformulate constraints A.29 into the standard form as follows:
Σ T3 = 1, Σ
Tq - w--
(A.30a)
SES
0 w S || - g. 0.05 S x S 1i
(A.30b)
where (S) refers to the number of elements the preference set S.
To conduct fair comparison, we train 50 epochs with batch size of 128, learning rate of 10 and a constraint tolerance of 1 across all satisfiability layers. When we stack constraints into block diagonal form to exploit parallelism, there are about 250 rows and 60,000 columns in the whole matrix. the training stage, weighted sum of prediction MSE error on future asset prices and the opposite of Sharpe ratio is used the loss function. All the experiments computer with a Intel(R) Xcon(R) Platinum 8360H CPU and a NVIDIA Tesla A100 GPU with 80GB memory through Pytorch 2.2. For AT, the initial estimate of Lipschitz constant is set to 0, the initial estimate of the dual variable is set to zero vector, the numerical precision set to 10mmachine- where Cuachine is the machine epsilon. In the main text, have reported the mean Sharpe ratio of each method in Table 5 Here, we provide the results related to the training performance in Table A.4 From Table A.4 can that the total memory usage is similar between each variant of AT and GLinSAT since there are only two constraints in the original problem. The projection time of LinSAT and GL is significantly less than that of CvxpyLayers and OptNet while AT-eSAATDDesse-mmlllie use the shortest total calculation time among all satisfiability layers.
conducted
Table A.4: Average allocated GPU memory and solution time of different satisfiability layers during batch processing of projection and backpropagation in the training phase of predictive portfolio allocation
1
1
GPU Mem./MB
Time/s
GPU Mem./BB
Time/s
Layer
Proj.
Backprop.
Proj.
Backprop.
Proj.
Backprop.
Proj.
Backprop.
CvxpyLayers
12.85
1,476
12.27
1.423
OptNet
2012
765.0
5.599
1.035
2012
765.0
4.606
1.279
LinSAT-Dense-100
81.34
434.9
0.091
0.393
105.8
410.5
0.127
0.413
LinSAT--Dense-500
81.34
434.9
0.266
0.386
136.6
385.6
0.315
0.413
GLinSSAT--Demse-EXpplee
112.9
402.6
0.143
0.376
107.9
407.2
0.148
0.377
GliinSAT--Dense-Impliict
1.565
514.1
0.090
0.320
1.565
513.6
0.090
0.322
LinSAT-Sparse
SEAATT-Sparse-100
223.6
295.0
0.113
0.400
290.7
227.6
0.161
0.421
LinS STASAA--SaaseSS0
223.6
295.0
0.114
0.398
371.1
224.1
0.225
0.442
GLLinSAATSparre--EXpiict
78.24
437.7
0.211
0.391
74.80
440.8
0.213
0.392
GLinSAT-SSarse-Immgp
4.501
512.0
0.143
0.316
4.501
511.7
0.139
0.315
22

---

A.11 Experimental details about power system unit commitment
In this section, we carry out experiments on power system unit commitment where the data comes from a real provincial power system.
We first briefly introduce the unit commitment problem. Unit commitment problem is optimization problem in power system operation and planning. It mainly involves deciding how to most economically and safely arrange the startup-shutdown status and output power of generators while ensuring constraints related to equipment and grid operation can be satisfied. Generally, constraints can be divided into soft constraints and hard constraints. In general, constraints directly related to generators are often regarded as hard constraints, e.g. the generator minimum up-time and down-time constraints. Constraints related the section power and load balance are usually regarded as soft constraints, For these constraints, we often penalize the corresponding violation in the objective [35]
Next, we will use the common three- brree-binary formulation 23 to model the unit commitment problem. Let T denote the total number of time steps that considered in unit commitment problem. G is the total number of generators. We denote T = (1, (1.... T and G (1, ... GG.. Let ug (t) denote whether the generator g is time Vgg denote whether the generator is turned at time t, and
we (t) denote if the generator is turned off at time t. Then, Ug (t). 0gg (t), wy (t) satisfy the following logical constraints:
Ug (1) - "g (t - 1) 0g (1) wi (1) ET
(A.31)
where T is the total number of time steps that considered in the unit commitment problem.
Units need to satisfy minimum up-time constraints and down-time constraints as follows:
"g (1) S Ug (1), E 9,1 ........
(A.32a)
i=t-U +1
Wy (i) S 1 Ug (t) g,t D
.......
(A.32b)
DTT-DT,+1
where g.g. DTg are the minimum up-time and minimum down-time for generator g respectively. Let Pg (f) denote the power produced by generator g at time t. Let Pa. Py denote the lower bound and upper bound of generator g. Then, we have constraints related to the bound of generator output as follows:
ug (t) S Pg (t) S uy (t) g GTEET
(A.33)
We also need to consider the ramping capability of each generator. The ramping constraints can be formulated as follows:
Py (t) - Py (t - 1) S -RU (t) + Pg + RUg) Ug (t) - u (t g)teE (A.34a)
Pg (t - 1) Pg (t) < RDDD (t) + Pg + RD,) Ug (t 1) Py (t) ,9 E G.E (A.34b)
where RUg, R.D,, denote the ramp-up rate and ramp-down rate of generator g.
The load generation balance constraint is modeled as the following soft constraint form:
Σ Σp, Po (t) + (1) (t) Σ Ela())tT deD
(A.35)
gEg
where D is the set of all loads, la (1) represents the d-th load at time l, s+ (1) 2 0,8 (1) 2 0 are the non-negative slack variables at time t which will be later penalized in the objective.
We denote the set of sections as K. The soft constraints on section power are provided follows:
Fx
HEgPy (t)
Σ
Heala (1) + s (t) - sp (t) < Fe, EKEET
(A.36)
9EG
dED
where FR, Fx are the lower bound and upper bound of the k-thh section power, Hsg is the generation shift factor that indicates the change of the k-th section power with respect to change in injection
23

---

at generator g. Hkg is the load shift factor that indicates the change of the k-th section power with respect to change in injection at load d. $x (t) 2 0, & (t) 2 0 are the non-negative slack variables related to the k-th section at time which will be later penalized in the objective.
Finally, we want to minimize the system operation cost and we can obtain the optimization problem as follows:
min
ΣΣ LET 9Eg
SU (coPa (1) + "g u1 (1) +
Σ IET
M (s+ (1)+ (0)+0000) (1) + ΣΣ
MR (5z (1) + 5x (1)
LET KER
(A.37a)
s.t.
A.31
A.36
" (t) E [0. 1y (t) {0, 13. my (t) (0.1).
(A.37b)
$ (1) 2 0, 0,1ET,
(A.37c)
8" (t) 0,s 0, (t) 2 20
(A.37d)
where M, MM are pre-defined penalty coefficients, co is the generator cost coefficient, g is the generator sp cost.
The power system we use in this article contains about 360 units. We set the total number of time steps T as 96 where the interval between two adjacent time steps is set to 15 minutes. More than 1,400 sections need to be considered the unit commitment problem. The penalty coefficient for load imbalance is set to 10 and the penalty coefficient for section power violation is set to 10". Based on the one-year load data, we solve the unit commitment problem via Gurobi within a 0. 1. 1% optimality gap. We can then obtain the optimal unit states for further supervised learning.
In supervised learning, we want to predict the optimal value of ug as accurately as possible so that we can fix these binary variables and quickly obtain a high-quality solution from solving a linear programming problem. Therefore, we should require the predicted variables Ug to satisfy the hard constraints (A.3I) and (A.32). namely the logical constraint and the minimum up-time and down-time constraints. In addition, Ug (t), vg (t) Wg (t) should be in the range between 0 and We use satisfiability layers to ensure the above constraints can be satisfied.
When we use GLinSAT the satisfiability layer, it is necessary canonize the original inequality constraints, By introducing bounded slack variables into constraints (A.29a), we can reformulate constraints (A.29) into the standard form as follows:
Un (i) + su, (t) = ug (n(() (t),9 E +.= -UG.
(A.38a)
i= UT
Wy (1) + swy (t) = 1 (g), E g,t DTgt T
(A.38b)
im DT. +1
Since we have Mg (t) E [o, 1), swy (t) can be quickly inferred, which is exactly sug (t) 10 swy (t) E [0. 1].
E ( wy (f) E [0, 11. the value range of variables svy (t) and
We use MLP-based neural network to learn the optimal unit states. The loads at each time-step forwarded to a 2-layer MLP where the hidden sizes are
32. Then, we concatenate the embeddings at all time steps. and forward the concatenated embedding to a 1-layer MLP where the hidden size is set to 3072. Finally, we use 2-layer MLP to read out the optimal unit states where the hidden size is set to 512. We use ReL. as the activation function in the hidden layers. We set the learning rate 3 and train neural networks for 100 epochs. We use 70% of one-year data as the training set and the other as the validation set. In the training stage, binary cross entropy loss is used as the loss function. For all satisfiability layers, the batch size is set 16 and the regularization parameter a is set to 0.1. It is worth noting that although the batch size is not that large, the scale of the optimization problem in each batch is still quite large. Each instance involves nearly 360 projection problems, and each projection problem involves about 160 constraints and 350 variables. When we stack constraints into block diagonal form to exploit parallelism, there are about 1,000,000 rows and 2,000,000 columns in the whole matrix. GLinSAT- will not report out-of-memor issues when we use AT. For GL,nSAT, the initial estimate of Lipschitz constant is set to 0, the initial estimate of the dual variable is set to zero vector, the
Implicit is the only way that
24

---

numerical precision is set to 10cmanenne. where Emachine is the machine epsilon. As to CvxpyLayers and OptNet, neither of them can directly handle such giant matrix within reasonable time thus we can only use sequential way instead. The performance of each satisfibility layer is provided in Table A..5 The results of LinSAT layers are not reported in Table A.55 since LinSAT layers only support positive linear constraints. All the experiments are conducted on computer with a Intel(R) Xeon(R) Platinum 8360H CPU and a NVIDIA Tesla A100 GPU with 80GB memory through Pytorch 2.2.
Table A.5: Average allocated GPU memory and solution time of different satisfiability layers during batch processing of projection and backpropagation in the training stage of predicting unit states
GPU Mcm./MB
Time/s
Layer
Proj.
Backprop.
Proj.
Backprop.
CvxpyLayers
2771
684.0
OptNet
33012
96.64
257.4
23.60
Gl GL.inSAT-Sparee-Immma
930.7
76.86
26.78
1.636
Note: Statisties of CvxpyLayers and OptNet are based on the first epoch since we cannot obtain welll-raiied model in reasonable time.
different values of test their performance. When σ is set to 10 we use the couble-precision floating-poin numbers during projection to avoid potential numerical issues. When solve it via Gurobi. After we have obtained the outputs of satisfiability layers, we round the outputs of final layers to 0 and then fix the unit state variables Ug using these rounded outputs. Once we fix all unit state variables to 0 or 1, the original mixed-intese turns to linear programming (LP) problem. We use Gurobi to solve such an LP problem and record whether the problem with fixed unit states is feasible and record the optimal solution if exists. We also compare the optimal objective obtained from fixing variables with the original optimal objective. The corresponding average gaps between them are shown in Table 6
In the validation stage,
set to O, the projection problem turns into linear programming problem and we
linear programming (MILP) problem
To illustrate the importance of the satisfiable layer, we additionally substitute the original satisfiability layer with a simple sigmoid activation function in both training and validation stages and report the corresponding result in Table We can easily see that satisfiability layers are essential to produce feasible unit states.
25

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: In this paper, we consider making a batch of neural network outputs satisfy bounded and general linear constraints. We present Gl GLinSAT, which is the first general linear satisfiability layer which all the operations are differentiable and
Experimental results demonstrate the advantages of GL over existing methods.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made the paper and important assumptions and limitations. A No or NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize other settings. is fine to include aspirational goals as motivation as long is clear that these goals are not attained by the paper.
2. Limitations
Question:
Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification: We have discussed our limitations in Appendix A.3 For variables with one- sided boundary or explicit boundary, our method cannot be directly used. A possible workaround is manually calculate the implicit bounds of these variables through domain propagation but we have not implemented such an algorithm in GLinSAT. Also, currently our proposed framework cannot deal with conic constraints.
Guidelines:
The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors encouraged to create separate "I "Limitations" section in their paper.
The paper should point out any strong assumptions and how robust the results to violations of these assumptions (e.g.. independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be.
The authors should reflect the scope of the claims made, e.g., if the approach was only tested on a few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution low or images are taken in low lighting. Or speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness.
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, worse outcome might be that reviewers discover limitations that aren't acknowledged the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
26

---

3. Theory Assumptions and Proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof?
Answer: [Yes]
Justification: In the derivation process, we always assume the feasible region is non-empty as
shown in Sec. 2.1
The derivation process of our results can be found in Sec. 2.1 Appendix
A.5
A.6 and A.7
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced in the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but they appear in the supplemental material, the authors are encouraged to provide short proof sketch to provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental Result Reproducibility
Does the paper fully disclose all the information needed to reproduce the main ex- perimental results the paper to the extent that affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Question:
Answer: [Yes]
Justification: We have detailed our experimental settings in Appendix A.8 A. and the code released in github can be also be used as reference.
A.9
A.
and
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, a No answer to this question will not perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is novel architecture, describing the architecture fully might suffice, or if the contribution a specific model and empirical evaluation, it may be necessary either make it possible for others replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g., in the case of large language model), releasing of a model checkpoint, or other means that appropriate to the research performed.
While NeurIPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) the contribution is primarily a new algorithm, the paper should make clear how reproduce that algorithm. (b) If the contribution is primarily the architecture clearly and fully.
model architecture, the paper should describe
(c) the contribution a new model (e.g., a large language model), then there should either be way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset instructions for how to construct the dataset).
27

---

(d)
We recognize that reproducibility may be tricky in some cases, in which case authors are welcome describe the particular way they provide for reproducibility. In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should be possible for other researchers to have some path to reproducing or verifying the results.
5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instruc- tions to faithfully reproduce the main experimental results, as described in supplemental material?
Answer: [Yes]
Justification: We provide our codes and instructions the first three experiments, but currently we are unable to open source the data and code about the unit commitment problem in real power system due to data security issues and confidentiality agreements.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurIPS code and data submission guidelines (https://nips. cc/ Sublic/guides CoodeSubmissionool.
for more details.
While we encourage the release of code and data, we understand that this might be possible, so "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurlPS code and data submission guidelines (https
7 //nips /nips.cccpublle... http/uublccggguueee
for more details.
The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible in supplemental material (appended to the paper) recommended, but including URLs to data and code permirted.
6. Experimental Setting/Details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary understand the results?
Answer: [Yes]
Justification: We have detailed our experimental settings in Appendix A. and the code released in github can be also be used as reference.
A.9 A.10 and
Guidelines:
The answer NA means that the paper does include experiments.
The experimental setting should be presented the core the paper to level of detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, supplemental material.
7. Experiment Statistical Significance
Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Question:
Answer: [No]
Justification: Error bars are not reported because it would be computationally expensive.
Guidelines:
28

---

The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by enror bars, confi- dence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper.
The factors of variability that the error bars capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g., Normally distributed errors).
should be clear whether the error bar is the standard deviation or the standard error of the mean.
It is OK to report I-sigma error bars, but should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not show in tables figures symmetric error bars that would yield results that are out of range (c.g. negative error rates).
error bars are reported in tables or plots, The authors should explain the text how they were calculated and reference the corresponding figures or tables in the text.
8. Experiments Compute Resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?
Question:
Answer: [Yes]
Justification: All the experiments are conducted on computer with a Intel(R) Xeon(R) Platinum 8360H CPU and NVIDIA Tesla A100 GPU with 80GB memory through Pytorch 2.2. We also provide the batch processing performance in Table 2 Table A.1 Table A.3
Table A.4 Table A.5
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU or GPU, internal cluster, cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental rums as well as estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported the paper (e.g., preliminary or failed experiments that didn't 't make it into the paper).
9. Code Of Ethics
Question: Does the research conducted the paper conform, in every respect, with the NeurlPS Code of Ethics https: www.erripp.ccppulliiccccpc http:///aeaaeeeccccppp////aaaaaa
Answer: [Yes]
Justification:
The research conducted in the paper conform with the NeurIPS Code of Ethics.
Guidelines:
The answer NA means that the authors have not reviewed the NeurlPS Code of Ethics.
If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (c.g.. if there special consid- eration due to laws or regulations in their jurisdiction).
10. Broader Impacts
Question: Does the paper discuss both potential positive societal impacts and negative socictal impacts of the work performed?
29

---

Answer:
[Yes]
Justification:
We have discussed potential impacts in Appendix A.2
Guidelines:
The answer NA means that there no societal impact of the work performed.
If the authors answer NA or No. they should explain why their work has no societal impact why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious or unintended uses (e.g.. disinformation, generating fake profiles, surveillance), fairesss considerations (e.g., dep.lymmlo of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and tied to particular applications, let alone deployments. However, if there is direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.
there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML).
11. Safeguards
Question:
Does the paper describe safeguards that have been put in place for responsible release of data models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Answer: [NA]
Justification:
This paper poses no risk for misusing our proposed method.
Guidelines:
The answer NA means that the paper poses no such risks.
Released models that have high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Are the creators or original owners of assets (e.g., code, data, models), used in the paper. properly credited and are the license and terms of use explicitly mentioned and properly respected?
Question:
Answer: [Yes]
Justification: We have properly credited the used Pascal VOC Keypoint dataset [1] with
Berkeley annotations 32 under the unfiltered setting 20 33
We also cite the TSP dataset
and portfolio allocation dataset provided in Ref. 2
Guidelines:
The answer NA means that the paper does not existing assets.
30

---

The authors should cite the original paper that produced the code package or dataset.
The authors should state which version of the asset is used and, if possible, include URL.
The name of the license (e.g., CC-BY 4.0) should be included for each asset. For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided.
If assets are released, the license, copyright information, and terms of use in the datasets..paerssttt paperswi thcode. com/datasets package should be provided. For popular datasets,
has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided.
If this information is not available online, the authors are encouraged to reach out to the asset's s creators.
13. New Assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [NA]
Justification: This paper does not release new assets.
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details of the dataset/code/model part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset is used.
At submission time, remember anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsourcing and Research with Human Subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well as details about compensation (if any)?
Answer: [NA]
Justification: This paper does involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing research with human subjects.
Including this information in the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included in the main paper.
According the NeurIPS Code of Bthies, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional Review Board (IRB) Approvals or Equivalent for Research with Human Subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?
Answer: [NA]
Justification: This paper does not involve crowdsourcing research with human subjects.
31

---

Guidelines:
The answer NA means that the paper does not involve crowdsoureing nor research with
human subjects.
Depending on the country in which research is conducted, IRB approval (or equivalent)
may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such the institution conducting the review.
32


Reviews:
[{"review_id": "UgwWuZbXSl", "reviewer": "Reviewer_6WNJ", "summary": "The main contribution of this paper is a method for enforcing arbitrary linear constraints on the outputs of a neural network. This is achieved in a differentiable way so that the entire pipeline can be trained end to end to solve constrained optimization problems. This is achieved by viewing the problem of projecting the outputs of the neural network to the feasible region as an entropy-regularized LP.  This in turn can be reformulated as a convex unconstrained optimization problem with Lipschitz gradients. The problem is solved via adaptive primal-dual gradient descent and the authors provide a more efficient way to perform the backward pass than by just calling autodiff. The method is significantly more efficient than previous work and it works well on several problems.", "questions": "see above", "limitations": "n/a", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "- I think this is an interesting approach to satisfying linear constraints on the outputs of a neural network and it builds nicely on ideas from the OT literature. \n- The proposed method is considerably more efficient than previous works and it improves nicely over previous work on positive linear satisfiability via OT.\n- Empirically, the proposed method seems to improve over previous work.", "weaknesses": "- In terms of contribution, I would argue that this is somewhat incremental, given that previous works can handle some of the same problems. The method here is not limited by positivity like Linsatnet,which is nice, but I am not convinced it's a sufficiently strong contribution. However, there is another stated benefit of the method which is its computational efficiency. This leads me to my next point.\n- My main concern regarding the experimental evaluation is the scale of those experiments. For example, the TSP results are on quite small instances (20 nodes) which is quite far from anything that could be considered practical. How far can the method scale and how well does it do compared to a solver in those cases? The introduction of the paper argues that solvers can be quite expensive in order to motivate a neural approach but how does Gurobi do on those same problems. How much time does it take to solve them? Is there any setting in which this approach is competitive with Gurobi?  \n\nI am asking this because for some combinatorial optimization problems (including classical TSP)  there are neural network-based approaches that can have competitive results vs Gurobi given a certain time budget. How does the method here do against Gurobi for different budgets?\n\nOverall, I like the direction of the paper and the approach it takes but I am not convinced by the empirical results. I start with a tentative score which I will update after the rebuttal.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thanks for the thorough review and valuable comments. We are encouraged with your acknowledgement to our idea, experiments and contribution. Below we respond to your specific comments.\n\n**W1: Is the contribution compared with previous works incremental?**\n\nAlthough existing satisfiability layers, such as LinSAT, CvxpyLayers and OptNet, can handle a few problems, they have significant drawbacks in either applicable scenarios or computational efficiency. **Our proposed GLinSAT is the first satisfiability layer that can fully use GPU's parallel computing capability through matrix-factorization-free operations while supporting general linear constraints at the same time.** We formulate the linear satisfaction problem as entropy-regularized linear programming and prove that it can be transformed into unconstrained convex optimization with Lipschitz continuous gradient. Thus it can be solved by gradient descent based algorithms.\n\nCompared with LinSAT, our proposed GLinSAT supports general linear constraints and thus enhances the expressiveness of satisfiability layers. **Even for a simple constraint $x\\\\le y$, namely $x-y\\\\le0$, LinSAT cannot be used due to the negative coefficient in front of $y$, which shows the limited expressiveness of LinSAT.** Moreover, we want to emphasize that **making neural network outputs satisfy general linear constraints is of great significance for practical application, since negative coefficients occur frequently in real-life constrained decision-making problems**, such as startup shutdown constraints in unit commitment [RW1.1], packing constraints in bin packing [RW1.2], flow conservation constraints in network flow [RW1.3], etc.\n\nAlthough other satisfiability layers, such as CvxpyLayers and OptNet, can be used for general linear constraint satisfaction, they may encounter efficiency issues. Take the unit commitment case as an example. In Table A.5, GLinSAT is ~10 times faster than OptNet and ~100 times faster than CvxpyLayers, which achieves significant acceleration. Considering we need ~1 day to train with GLinSAT, then we will need ~10 days with OptNet and ~100 days with CvxpyLayers, which is totally unacceptable. **Our proposed GLinSAT have dramatically shortened the training time for neural networks with satisfiability layers, making what was previously impractical now viable and efficient.** In addition, compared with LinSAT, our proposed GLinSAT can achieve a lower time complexity. **In Section A.6, we further show the time complexity of GLinSAT is $O(\\\\sqrt{\\\\theta/\\\\epsilon})$, which is lower than LinSAT's $O(\\\\theta/\\\\epsilon^2)$. Moreover, our proposed GLinSAT is guaranteed to converge while LinSAT sometimes is not.** From empirical results in Section A.8, we found LinSAT often reports warnings like \"non-zero constraint violation within max iterations\" when maximum iteration number is 100. If we make a fair comparison, namely setting the maximum iteration number in LinSAT to $\\\\infty$ just like GLinSAT, then LinSAT sometimes cannot converge. Therefore, our proposed GLinSAT is more reliable.\n\n[RW1.1] On mixed-integer programming formulations for the unit commitment problem. INFORMS Journal on Computing, 2020.\n\n[RW1.2] Machine Learning for the Multi-Dimensional Bin Packing Problem: Literature Review and Empirical Evaluation. arXiv 2023.\n\n[RW1.3] Machine learning based approaches to solve the maximum flow network interdiction problem. Computers & Industrial Engineering, 2022.\n\n**W2: Scalability issues**\n\nWe conduct experiments on 500-node TSP to show the ability of GLinSAT on large-scale problems. Despite the supplementary TSP experiments, in the original paper, the unit commitment experiments can also demonstrate such an ability.\n\n**As the number of nodes increases, the number of edges and decision variables can grow quadratically. In the 500-node TSP, there are about 250,000 decision variables, leading to large-scale problems.** We train neural networks using GLinSAT-Sparse with $1/\\\\theta=0.1$. In the validation stage, beamsearch post-processing is used for obtaining feasible tours. We also use Gurobi with 32 threads to solve the MTZ formulation of TSP. The results of Gurobi under different time budgets and our approach are shown as follows:\n\n|$\\\\hspace{11em}$TSP-StartEnd$\\\\hspace{8em}$TSP-Priority$\\\\hspace{4.4em}$|\n|:----------------------------------------------------------:|\n\n|$\\\\hspace{2em}$Method$\\\\hspace{2em}$|Mean Time|Mean Tour Length|Mean Time|Mean Tour Length|\n|:--------------------:|:------------:|:--------------:|:------------:|:--------------:|\n|GLinSAT|**3.225s**|35.323|**3.275s**|**36.005**|\n|Gurobi-1200s|1200s|52.820|1200s|82.242|\n|Gurobi-2400s|2400s|**35.266**|2400s|45.700|\n\nFrom the table above, we can find that, solving 500-node TSP via Gurobi is extremely time-consuming while our method can generate feasible solution in few seconds. Since solving via Gurobi requires a large amount of time and computer threads, we can only compare the performance on 100 randomly generated cases. In TSP-StartEnd, the performance of our method is similar to Gurobi using 32 threads within 2400 seconds. As for the more complicated TSP-Priority, the solution from Gurobi within 2400s cannot outperform GLinSAT's solution in average.\n\nMoreover, the empirical results from the unit commitment experiment can also demonstrate the scalability of the proposed method. In this experiment, we use a real-life bulk power system which consists of about 360 units while the total number of considered time steps $T$ is set to 96. **Consequently, when we stack constraints into block diagonal form to exploit parallelism, there are about 1,000,000 rows and 2,000,000 columns in the coefficient matrix.** Experimental results in Table A.5 of the Appendix also demonstrate the applicability and efficiency of our proposed method in such a large case.\n\n------\n\nWe hope this response could help address your concerns, and wish to receive your further feedback soon.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "sRj2MJ86R5", "reviewer": "Reviewer_MtC5", "summary": "This paper proposes a new architecture named GLinSAT that projects the output of a neural network to the feasible region of bounded and general linear constraints. The approach is based on a gradient-based approach for solving an entropy-regularized linear program and can be implemented with backpropagation/differentiation. Numerical results are reported for several tasks.", "questions": "How do you choose $\\theta$ in practice?", "limitations": "The authors clearly discuss limitations in Appendix A.3.", "rating": 5, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. Outputing solutions satisfying linear constraints is an important topic and the idea in this paper is reasonable and smooth.\n2. The method works for general linear constraints and can be implemented efficiently (free of matrix factorization and can be differentiated with backpropagation).\n3. The numerical results look nice compared with previous works.", "weaknesses": "The authors do not discuss how integer constraints are satisfied when some entries of $x$ are constrained to be integers. Note that several tasks in the numerical part are optimization over integer variables.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thanks for the valuable comments, constructive advice, and recognition of our idea and empirical results. Below we respond to your specific comments.\n\n**W1: How to make it possible for integer constraint satisfaction?**\n\nThanks for the comment. As shown in Section A.8, A.9, A.11 in Appendix, when it comes to integer constraints, we need to do some post-processing since GLinSAT is originally designed for linear constraint satisfaction. **We restate the details about integer constraint satisfaction in the general response.**\n\n**Q1: How to choose $θ$ in practice?**\n\nActually, the regularization coefficient $1/\\\\theta$  controls the smoothness of the outputs. The smaller $1/\\\\theta$  is, the more the outputs tend to be at the extreme point of the feasible region. A recommended way is to start with $0.1$, but users can also use methods like grid search to tune the parameter for their specific tasks. From the empirical results shown in Table 3-6, using a small value of $1/\\\\theta$ usually yields better results in these tasks. In addition, during the inference phase, using different values for $1/\\\\theta$  than what was used in the training phase is also worth a try. Users may also choose the best parameter value through grid search.\n\n------\n\nWe hope this response could help address your concerns, and wish to receive your further feedback soon.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "TGU9vtTA1T", "reviewer": "Reviewer_TMC5", "summary": "The paper introduces \"GLinSAT,\" a novel neural network layer designed to enforce general linear and bounded constraints on neural network outputs using a differentiable approach. The method leverages entropy-regularized linear programming, transforming it into an unconstrained convex optimization problem that avoids matrix factorization and is efficiently solvable on GPUs using accelerated gradient descent. The authors present experimental results that demonstrate the effectiveness of GLinSAT compared to existing solutions across several constrained decision-making applications.", "questions": "1. For L2 projection layer, the projection distance will be no larger than the approximation error; thus, the universal approximation ability of NN+L2-projection will be maintained. Meanwhile, with the dot-product-based projection layer, will the approximation ability of NN+Dot-Product-projection be affected? \n2. A deeper discussion on the motivation of this entropy regularization design for linear programming and its related works would be appreciated. For example, is it the first time to formulate such a dot-product projection problem with entropy regularization in literature?\n3. Providing more detailed information about the problem sizes used in experiments (e.g., dimensions of m and n) would help in assessing the scalability and practical applicability of GLinSAT.\n4. Since GLinSAT is designed for continuous problems, the experiments contain 3 discrete problems. It would be beneficial for the authors to explain in detail how GLinSAT is applied to these discrete problems and how feasibility over discrete decision variables is recovered.\n5. Why not compare GLinSAT with CvxpyLayers/OptNet with the L2 projection?\n6. Why can the feasibility of GLinSAT not be guaranteed in Table 6?", "limitations": "The authors have discussed limitations in Appendix A3.", "rating": 7, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1.  The paper addresses the challenge of imposing general linear constraints on neural network outputs in a differentiable manner. The approach of using entropy-regularized linear programming and transforming it into an unconstrained optimization problem is novel and technically sound.\n2. The reformulation from a projection problem to an unconstrained optimization problem not only ensures differentiability but also boosts computational efficiency, as demonstrated in the forward and backward processes through dual formulation-based optimization and implicit differentiation.\n3.  The experiments cover a diverse set of applications, including the traveling salesman problem, graph matching, portfolio allocation, and power system unit commitment. The proposed GLinSAT outperforms several baselines over those problems in terms of constraint satisfaction and computational efficiency.", "weaknesses": "1. Besides differentiable layer-related works, the paper could benefit from a more extensive review of other works about NN feasibility over linear constraints, such as [1], [2], and [3].\n2. The motivation for using a dot product in the objective function (Eq. 3) as opposed to a standard L2/L1 norm minimization for projection problems is unclear. Clarification on the rationale behind this choice, its advantages, or specific scenarios where it is particularly beneficial would be valuable. Additionally, is dot-product-based projection the only one that admits the entropy-regularized formulation and the unconstrained dual problem?\n3. The comparison between the entropy-regularized and the unregularized projection problems could be more explicitly detailed, especially how the optimality gap between the two problems is affected by $\\theta$.\n\n[1] \"Tordesillas, J., How, J. P., & Hutter, M. (2023). Rayen: Imposition of hard convex constraints on neural networks. arXiv preprint arXiv:2307.08336.\", \n[2] \"Zhao, T., Pan, X., Chen, M., & Low, S. (2023). Ensuring DNN solution feasibility for optimization problems with linear constraints. In The Eleventh International Conference on Learning Representations.\"\n[3] \"Tabas, D., & Zhang, B. (2022, June). Computationally efficient safe reinforcement learning for power systems. In 2022 American Control Conference (ACC) (pp. 3303-3310). IEEE.\"", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thanks for the thorough review and valuable comments. We are encouraged with your acknowledgment of our methodology, empirical results, and contribution. Below we respond to your specific comments.\n\n**W1: Some related works need to be reviewed**\n\nThanks for sharing these works other than differentiable-optimizer-based layers for linear constraint satisfaction.\n\nThe first work, Rayen, calculates null space for equality constraints in the offline stage and then use the basis for convex constraint satisfaction in the online stage. However, its efficiency may be affected when constraints are not fixed. As the basis for the null space of a matrix is usually dense, calculating null space for large matrices may face efficiency and memory challenges. \n\nThe second work shrinks the feasible region of inequality constraints and decides DNN size for guaranteeing feasibility through solving a complicated bilevel programming. However, as it can only handle inequality constraints, equality constraints must be removed through Gaussian elimination, which may lead to denser matrices.\n\nThe third work uses gauge function to map the neural network outputs from a $\\\\infty$-norm unit ball into a given polyhedral. Despite its success in the field of control, this method may encounter difficulties in handling equality constraints since the polyhedral need to contain the origin as an interior point.\n\nAll these works are very insightful and provide excellent ideas for constraint satisfaction. We will include these works in our new version of literature review.\n\n**W2 & Q2: Motivation for entropic regularization and dot-product based projection instead of L2/L1 normalization**\n\nThanks for the comment. **We provide the details in the general response.**\n\n**W3: Comparison between entropy-regularized and unregularized problems**\n\nHere, based on existing results in Table 6, we conducted more detailed tests and comparisons on unit commitment problems under different $1/\\\\theta$ values in the validation stage. The results after rounding operations are shown as follows:\n\n|$\\\\;1/\\\\theta$|Feasible Ratio|Average Gap for Feasible Cases|\n|:----------:|:------------:|:----------------------------:|\n|0.01|86.23%|0.1119%|\n|0.005|95.41%|0.1381%|\n|0.001|98.17%|0.1109%|\n|0.0005|100%|0.1114%|\n|0.0001|100%|0.1114%|\n|0|100%|0.1114%|\n\nFrom the table, we can find that no matter what the values of $1/\\\\theta$ are, the average optimality gaps are all around 0.1%. We can also find that in this case, as $1/\\\\theta$ decreases, the feasible ratio will increase. When $1/\\\\theta\\\\le0.0005$, our method can reach a 100% feasible ratio and the average gaps all tend to 0.1114%.\n\n**Q1: Universal approximation ability of NN+L2-projection and NN+Dot-Production-projection**\n\nAs pointed in the answer to W2 & Q2 in the general response, L2-projection and softmax function can be also regarded as dot-product projection methods. We believe the key to the universal approximation ability is not the projection type, but the nonlinearity of the input-output mapping. Since the optimal solution should satisfy the KKT condition with respect to the input parameters, according to implicit function theorem, the optimal solution can be viewed as a nonlinear function of inputs. In this way, universal approximation theorem still holds.\n\n**Q3: Detailed Problem Size**\n\nThe table below shows the size of constraint matrices when we stack constraints into block diagonal form to exploit parallelism. We will also include this information in our new version of paper.\n\n||$\\\\hspace{1.5em}$TSP$\\\\hspace{1.5em}$|Partial Graph Matching|Portfolio Allocation|Unit Commitment|\n|:--:|:--------:|:--------------------:|:------------------:|:-------------:|\n|$m$|~40,000|~2,500|~250|~1,000,000|\n|$n$|~400,000|~13,000|~60,000|~2,000,000|\n\n**Q4 & Q6: Integer constraint satisfaction & Feasibility issues in Table 6**\n\nThanks for the comment. As shown in Section A.8, A.9, A.11 in Appendix, when it comes to integer constraints, we need to do some post-processing since GLinSAT is originally designed for linear constraint satisfaction. **We restate the details in the general response.** Moreover, from our response to W3, we can find that in this case, when $1/\\\\theta\\\\le0.0005$, our method can reach a 100% feasible ratio.\n\n**Q5: Comparison with L2-projection**\n\nAs shown in the answer to W2 & Q2 in the general response, L2-projection is just one special case of dot-product-based projection where the quadratic regularization coefficient need to be $0.5$. In Section 3, we have already used OptNet where the quadratic regularization coefficient is set to $0.1$. Here, we also provide the result when using L2-projection. Since CvxpyLayers and OptNet is time-consuming, we only provide the results of partial graph matching and portfolio allocation.\n\n|$\\\\hspace{2em}$Problem Type$\\\\hspace{2em}$|Mean F1 using L2-projection|Mean F1 when Quadratic-$1/\\\\theta=0.1$|\n|:--------------------:|:-------------------------:|:-----------------------------------:|\n|Partial graph matching|0.614|0.619|\n\n|$\\\\hspace{2em}$Problem Type$\\\\hspace{2em}$|S. Ratio using L2-projection|S. Ratio when Quadratic-$1/\\\\theta=0.1$|\n|:------------------:|:--------------------------:|:------------------------------------:|\n|Portfolio allocation|2.504|2.553|\n\nFrom the tables above, we can find that using a larger regularization coefficient in these two cases didn't produce better results than before.\n\n------\n\nWe hope this response could help address your concerns, and wish to receive your further feedback soon.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "EzyZF3qPSI", "reviewer": "Reviewer_U49f", "summary": "This paper presents a new approach to enforce constraints on the output of neural networks. Termed **G**LinSAT, this paper extends the application scope of LinSAT from positive linear constraints to general linear constraints. The authors also adopt ideas from OptNet/Cvxpylayers to derive implicit gradients in order to save GPU memory. Experiment study on all 3 experiments involved in LinSAT paper and one more case study with negative constraint terms show the proposed method works as expected to enforce constraints, while at the same time, offering better efficiency in GPU memory and (in many cases) running time.", "questions": "* In the discussion of the limitation of LinSAT, the authors mentioned the bin packing problem as an example. Can you elaborate on why the bin-packing constraint has negative components?\n* Can we derive a similar form by maximizing $\\mathbf{c}^\\top\\mathbf{x}$? I think it is more intuitive that larger values are preserved to be larger after an activation layer (like GLinSAT).\n* How to compute $\\partial x/\\partial y$ in Eq (11)?\n* Please make sure the code is released and easily accessible to the research community if this paper is accepted.", "limitations": "N/A", "rating": 8, "confidence": 4, "soundness": 4, "presentation": 3, "contribution": 4, "strengths": "* Enforcing constraints to the output of neural networks is an important research direction given the fact that deep neural networks are being extended from \"traditional\" classification and regression tasks to more complicated real-world decision-making with complicated constraints.\n* The proposed approach seems technically sound. It was well-grounded with existing work (including LinSAT and OptNet) and the derivation of Lagrange duals and adoption of ADPAGD sounds novel and interesting.\n* The experiment study shows that the proposed GLinSAT is able to encode the constraint as expected, reaching similar (and marginally better) accuracy than LinSAT. At the same time, it is more GPU memory efficient, which is important for scaling up the applications.", "weaknesses": "* The major drawback of GLinSAT over LinSAT is that its inference time is longer if LinSAT is set with its default 100 iterations. It will be helpful for the readers if there is a \"discussion\" section about what are the practical tips and suggestions for choosing the method for encoding constraints.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thanks for the thorough review and valuable comments. We are encouraged with your acknowledgment of our motivation, methodology, empirical results, and contribution. Below we respond to your specific comments.\n\n**W1: Inference time of GLinSAT is sometimes longer than LinSAT with its default 100 iterations and required for discussion on practical tips**\n\nAs shown in Section A.8 in the Appendix, we have illustrated the reason for the speed difference. **When the maximum iteration number of LinSAT is set to 100, LinSAT often reports warnings like \"non-zero constraint violation within max iterations\".** After we increase the max iteration number to 500, the number of warnings has decreased, but there are still some warnings that indicate LinSAT has not converged. At this time, we already find that GLinSAT-Dense has already become faster than LinSAT-Dense from Table A.1. If we make a fair comparison, namely setting the maximum number of iterations of LinSAT to $\\\\infty$ just like GLinSAT, then LinSAT sometimes will not converge and will be always iterating.\n\nIn addition, thanks for your suggestion on adding a discussion section about practical tips. In our new manuscript, we will supplement the following tips for using GLinSAT. First, the regularization coefficient $1/\\\\theta$  controls the smoothness of the outputs. The smaller $1/\\\\theta$  is, the more the outputs tend to be at the extreme point of the feasible region. A recommended way is to start with $0.1$. Using different values for $1/\\\\theta$ in the inference stage is also worth a try. Users can also use grid search to tune the parameter for their specific tasks. Second, for large-scale problems, it will be better to use the sparse version and implicit differential for saving GPU memory.\n\n**Q1: Is there any negative component in bin-packing problem?**\n\nIn the bin-packing problem, there are indeed negative components in the constraint. Let's take the following classic bin-packing problem as an example [RQ1.1]:\n$$\n\\\\begin{array}{c}\n\\\\min\\\\limits_{y_j,x_{ij}} \\\\sum\\\\limits_{j = 1}^n{y_j}\\\\\\\\\ns.t.\\\\sum\\\\limits_{i\\\\in I}{s(i)x_{ij}}\\\\le B{y_j}\\\\\\\\\n\\\\sum\\\\limits_{j=1}^n{x_{ij}}=1\\\\\\\\\ny_j\\\\in\\\\{0,1\\\\},x_{ij}\\\\in\\\\{0,1\\\\}\n\\\\end{array}\n$$\nwhere $I$ is the set of items,  $s(i)\\\\in\\\\mathbb{Z}^+$ is item size, $B\\\\in\\\\mathbb{Z}^+$ is bin capacity, $y_j=1$ if bin $j$ is used and $x_{ij}=1$ if item $i$ is put into bin $j$. To see the negative coefficients clearly, we can reformulate the second linear constraint into the standard form $Ax\\\\leq b$ as follows:\n$$\n\\\\sum\\\\limits_{i\\\\in I}{s(i)x_{ij}}-B{y_j}\\\\le0\n$$\nDue to the negative coefficient in front of $y_j$, LinSAT cannot be directly used for such a constraint. Similarly, even for a simple constraint $x\\\\leq y$, there will be negative coefficients in front of $y$ in the canonical form. Moreover, for the more complicated 2D [RQ1.2] or 3D [RQ1.3] bin-packing problems, there will be more negative coefficients in the constraints.\n\n[RQ1.1] Bin-packing problem, Chapeter 8 in Knapsack Problems: Algorithms and Computer Implementations, John Wiley & Sons, Inc., 1990.\n\n[RQ1.2] The two-dimensional bin packing problem with variable bin sizes and costs. Discrete Optimization, 2005.\n\n[RQ1.3] Machine Learning for the Multi-Dimensional Bin Packing Problem: Literature Review and Empirical Evaluation. arXiv 2023.\n\n**Q2: Can we derive a similar form by maximizing $\\\\boldsymbol{c}^⊤\\\\boldsymbol{x}$?**\n\nYes, of course! Making larger values to be larger after an activation layer sounds more reasonable. Actually, the reason why we initially use $\\\\min \\\\boldsymbol{c}^⊤\\\\boldsymbol{x}$ is affected by optimal transport where the objective is to minimize the total transportation cost. But in the field of neural networks, maximizing $\\\\boldsymbol{c}^⊤\\\\boldsymbol{x}$ is indeed a more suitable choice.\n\nTo maximize $\\\\boldsymbol{c}^⊤\\\\boldsymbol{x}$, all we need to do is to flip the signs with respect to $c$ in the derivative and function expressions. In our new version of this paper, we will use the maximization syntax and revise the corresponding formulas to provide a more intuitive understanding. Thanks for your suggestion!\n\n**Q3: How to compute $\\\\partial\\\\boldsymbol{x}/\\\\partial\\\\boldsymbol{y}$ in Eq (11)?**\n\nIn the actual computation process, the matrix $\\\\frac{\\\\partial\\\\boldsymbol{x}}{\\\\partial\\\\boldsymbol{y}}$ is not explicitly formulated for saving GPU memory. Therefore, in the Appendix of the original manuscript, we only provide the formula of $\\\\frac{\\\\partial l}{\\\\partial\\\\boldsymbol{x}}\\\\frac{\\\\partial \\\\boldsymbol{x}}{\\\\partial\\\\boldsymbol{y}}$ in Eq (A.17e) instead of $\\\\frac{{\\\\partial \\\\boldsymbol{x}}}{{\\\\partial \\\\boldsymbol{y}}}$. In our revised manuscript, we will also supplement how $\\\\frac{\\\\partial\\\\boldsymbol{x}}{\\\\partial\\\\boldsymbol{y}}$ is calculated to provide a better understanding on the derivative calculation.\n\nWhen  $\\\\boldsymbol{x(y)}=\\\\boldsymbol{u}\\\\circ\\\\boldsymbol{\\\\sigma(-\\\\theta\\\\boldsymbol{u}\\\\circ(\\\\boldsymbol{c-A}^T\\\\boldsymbol{y}))}$, we can calculate the derivative with respect to $\\\\boldsymbol{y}$ as follows:\n$$\n\\\\frac{\\\\partial x_q}{\\\\partial y_p}=\\\\theta x_q(u_q-x_q)A_{pq}\n$$\nBy writing the above equation into matrix form, we have:\n$$\n\\\\frac{\\\\partial\\\\boldsymbol{x}}{\\\\partial\\\\boldsymbol{y}}={\\\\rm{\\\\boldsymbol{diag}}}(\\\\theta\\\\boldsymbol{x}\\\\circ (\\\\boldsymbol{u-x}))\\\\boldsymbol{A}^T\n$$\nwhere $\\\\circ$ represents the element-wise multiplication, ${\\\\rm{\\\\boldsymbol{diag}}}(·)$ maps a vector to its corresponding diagonal matrix.\n\n**Q4: Please make sure the code is released and easily accessible to the research community if this paper is accepted.**\n\nThanks for your recognition for our contribution. We are currently reorganizing our codes so that the research community can easily use them. Once this paper is accepted, our source code will be released in Github.\n\n------\n\nWe hope this response could help address your concerns, and wish to receive your further feedback soon.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:
Dear Chairs and Reviewers,

We greatly appreciate the reviewers' time, valuable comments, and constructive suggestions. We are delighted that all the reviews have expressed a positive inclination towards accepting our submission. Overall, the reviewers acknowledge our methodology as "interesting" (6WNJ, U49f), "nice" (6WNJ, MtC5), "novel" and "sound" (TMC5, U49f), with empirical results showing better efficiency than existing methods (6WNJ, MtC5, TMC5, U49f).

In the author response period, we make every effort to address reviewers’ concerns. Below are answers to some common questions.

**(1) Motivation for entropic regularization and dot-product based projection in Eq. (3) instead of L2/L1 normalization**

First, if L1-norm is used as the objective, then the optimization problem will become a linear programming (LP). As pointed by [R1.1],  the optimal solution to an LP may not be differentiable (or even continuous) with respect to its parameters. **As a result, the non-differentiability of L1-norm may make neural networks unable to train.**

**As for L2-norm, we can show that it is also a dot-product based projection, but with an additional quadratic regularization term** since $\\min\\frac{1}{2}\\left\\|\\boldsymbol{y-x}\\right\\|_2^2$ $\\Leftrightarrow$ $\\min-{\\boldsymbol{y}^T}\\boldsymbol{x}+\\frac{1}{2}\\sum{x_i^2}$, where $\\boldsymbol{y}$ is the output of the previous layer, $\\boldsymbol{x}$ is the decision variable and the dot product can be regarded as a measure of vector similarity. From above we can see that the main difference between Eq. (3) and L2-norm lies in the regularization term.

In Section 2.1 of our paper, we have already mentioned the motivation for using dot-product based projection and entropic regularization, which is inspired by optimal transport. The formulation of optimal transport is given as follows [R1.2]:
$$
\\begin{array}{c}
\\min\\limits_{X_{ij}\\ge0}\\sum\\limits_{i,j}{(c_{ij}X_{ij}+\\eta X_{ij}\\log X_{ij})}\\\\
{\\rm{s.t.}}\\boldsymbol{X1=r},{\\boldsymbol{X}^T}\\boldsymbol{1=c}
\\end{array}
$$
A similar example is the softmax function since it is the solution to the following optimization:
$$
\\begin{array}{c}
\\min\\limits_{\\boldsymbol{x\\ge0}}-\\boldsymbol{y}^T\\boldsymbol{x}+\\boldsymbol{1}^T(\\boldsymbol{x}\\circ\\boldsymbol{\\log x})\\\\
{\\rm{s.t.}}\\;\\boldsymbol{1}^T\\boldsymbol{x}=1
\\end{array}
$$
In Remark 1 of Section 2.1, we also show that **the main benefit of using entropic regularization terms instead of quadratic regularization terms mainly comes from $\\boldsymbol{x\\log x+(1-x)\\log(1-x)}$ being an essentially smooth function [R1.3].** Therefore, the infimum in Eq. (4) can be attained only on a stationary point and we can derive an unconstrained dual optimization problem.

Moreover, as pointed by reviewer U49f, making larger values to be larger after an activation layer sounds more reasonable. Therefore, in our new version of paper, we will use the $\\min-\\boldsymbol{c}^T\\boldsymbol{x}$ syntax to provide a more intuitive understanding.

As for whether the proposed approach is the only way to formulate unconstrained dual optimization problems, so far we are not sure. Perhaps this could become a direction for future research.

[R1.1] Melding the data-decisions pipeline: Decision-focused learning for combinatorial optimization, AAAI 2019.

[R1.2] Sinkhorn distances: Lightspeed computation of optimal transport. NeurIPS 2013.

[R1.3] Convex analysis, Chapter 26. Princeton university press, 1997.

**(2) How to make it possible to satisfy integer constraints?**

**As shown in Section A.8, A.9, A.11 in Appendix, when it comes to integer constraint, we need to do some post-processing since GLinSAT is originally designed for linear constraint satisfaction.** Here, we restate these methods more clearly.

**For the TSP problem, we have used two post-processing techniques, one is rounding and the other one is beam search.** In the validation stage, direct rounding on the outputs of GLinSAT can satisfy all constraints in about 94% of cases, while beam search can satisfy all constraints in 100% cases. The results can be found in Table 3 and Table A.2.

**For the partial graph matching problem, we use Hungarian algorithm and greedy strategy for post-processing.** We regard the cost of matching a pair of nodes as the outputs of GLinSAT, then use Hungarian algorithm to obtain a maximum matching. Finally, we use greedy strategy to preserve pairs with $p$-highest matching scores for constraint satisfaction. The mean F1 scores can be found in Table 4.

For the unit commitment problem, we want to predict the optimal values of unit status while satisfying the logical constraint and the minimum up-time and down-time constraints. **After we obtain the outputs of GLinSAT, we round the outputs to 0 or 1.** Although GLinSAT can guarantee the feasibility of linear constraints prior to rounding operations, the inevitable rounding off errors may lead to infeasibility in few cases after rounding operations. However, as shown in Table 6, the feasible ratio will increase as $1/\\theta$ decreases and reach 100% when $1/\\theta=0$. The reason is that when $1/\\theta$ gets smaller, the output of GLinSAT gets closer to the extreme point of the feasible region. **As pointed by [R2.1], when we consider the logical constraint and the minimum up-time and down-time constraints, these constraints formulate a convex hull so that the extreme points of the corresponding feasible region are binary.** In this case, the output of GLinSAT will tend to be binary when $1/\\theta\\to0$.

[R2.1] Minimum up/down polytopes of the unit commitment problem with start-up costs. IBM Research Report RC23628, 2005.

------

**In our individual responses, we provide detailed answers to all the specific questions raised by the reviewers.** We hope these responses could help address the reviewers' concerns, and further discussions are welcomed towards a comprehensive evaluation of our work.

Meta-review: Thank you for your valuable contribution to Neurips and the ML community. Your submitted paper has undergone a rigorous review process, and I have carefully read and considered the feedback provided by the reviewers.

This paper proposes a method to enforce arbitrary linear constraints on the outputs of a neural network. This is achieved by viewing the problem as an entropy-regularized linear program and formulating it as an unconstrained convex problem. Overall, the paper received mostly positive response from the reviewers. Some issues raised by the reviewers were successfully addressed in the rebuttal. These issues include (i) comparison with the previous work (ii) inference time (iii) certain mathematical details of the approach.

Given this positive assessment, I am willing to recommend the acceptance of your paper for publication.

I would like to remind you to carefully review the reviewer feedback and the resulting discussion. While most reviews were positive, the reviewers have offered valuable suggestions that can further strengthen the quality of the paper. Please take another careful look a the 'weaknesses' section of each reviewer comment and the issues (i)-(iii) listed above. I encourage you to use this feedback to make any necessary improvements and refinements before submitting the final version of your paper.

Once again, thank you for submitting your work to Neurips.

Best,
Area Chair
Decision: Accept (poster)