CG4AI: A Column Generation Framework for Training AI Models Under Constraints

# CG4AI: A Column Generation Framework for Training AI Models Under Constraints
Source: [https://arxiv.org/html/2608.26375](https://arxiv.org/html/2608.26375)
Youcef MagnoucheAffiliation:Huawei Technologies Ltd\., France Research Center, 18 Quai du Point du Jour, Boulogne\-Billancourt, 92100, FranceAbderrahmane DriouchEmail:[driouch\.ab\.dr@gmail\.com](mailto:[email protected])Affiliation:Huawei Technologies Ltd\., France Research Center, 18 Quai du Point du Jour, Boulogne\-Billancourt, 92100, FranceSébastien MartinAffiliation:Huawei Technologies Ltd\., France Research Center, 18 Quai du Point du Jour, Boulogne\-Billancourt, 92100, FrancePierre BauguionEmail:[firstname\.lastname@huawei\.com](mailto:[email protected])Affiliation:Huawei Technologies Ltd\., France Research Center, 18 Quai du Point du Jour, Boulogne\-Billancourt, 92100, France

###### Abstract

Standard machine\-learning training minimizes a loss function over a dataset, but does not guarantee that the resulting model will satisfy predefined rules or constraints on its outputs\. In many real\-world applications, ranging from autonomous systems to network routing, such guarantees are essential\. We propose*CG4AI*, a framework that builds a convex combination of AI models while enforcing linear constraints on the combined output\. A master linear program \(LP\) determines the optimal mixture weights, while a pricing subproblem generates new models guided by LP dual variables, focusing attention on the most violated constraints\. A cutting\-plane procedure extends feasibility guarantees beyond the training set\. We apply CG4AI to two problems: \(i\) digit classification on MNIST, where we demonstrate four distinct uses of constraints, learning from constraints alone, improving adversarial robustness, correcting misclassified examples, and enforcing output relabeling; and \(ii\) the multi\-commodity flow problem, where link capacity constraints are enforced on neural\-network routing predictors\. Experiments on MNIST and standard SNDLIB benchmark networks show that CG4AI reliably produces feasible predictors while achieving better accuracy than single\-model baselines\.

###### keywords

Column generation, Hard constraints, Neural networks, Ensemble learning, Multi\-commodity flow, Trustworthy AI, Linear programming, Safety

### 1Introduction

Artificial intelligence models are increasingly used to support or automate decisions in critical domains such as transportation, healthcare, and telecommunications\. These models are typically trained by minimizing a loss function over a labeled dataset\. This process does not provide any explicit guarantee that predictions will respect safety rules, physical laws, or regulatory requirements: a model can behave correctly on the training data and still produce infeasible or unsafe outputs on new inputs\.

This limitation becomes particularly important when predictions must satisfy*hard constraints*, strict conditions that cannot be violated under any circumstances\. Hard constraints arise naturally in many settings, as for autonomous vehicles to enforce a collision\-avoidance rule, or on radiotherapy planning where the dose delivered to healthy tissue must not exceed a given bound for any tumor configuration\. They also arise on network routing to predict as a split of traffic satisfying link capacity constraints\. Another interesting application is in regulation and compliance where a lending model must never violate legal fairness requirements, even for unusual applicant profiles\.

Enforcing such constraints in learned models is challenging\. The most common approach adds penalty terms to the training objective to discourage violations\. This transforms hard constraints into soft ones and cannot guarantee feasibility at test time\. Tuning penalty weights is also notoriously difficult\. More recent methods project the model output onto a feasible set\([Min and Azizan, 2024](https://arxiv.org/html/2608.26375#bib.bib1);[Tordesillas et al\., 2023](https://arxiv.org/html/2608.26375#bib.bib2)\)or embed convex optimization as a differentiable layer\([Amos and Kolter, 2017](https://arxiv.org/html/2608.26375#bib.bib3)\)\. Projection methods enforce constraints exactly but require solving an optimization problem at inference time, which increases latency and limits scalability\. Optimization layers work well for convex constraints but do not extend naturally to combinatorial or large\-scale settings\.

In this paper we propose CG4AI, a framework that enforces hard linear constraints*during training*using column generation\. Instead of training a single model, we build a convex ensemble: the final predictor is a weighted combination of simpler models \(the*columns*\), and the weights are determined by a master LP\. Constraints on the output of the ensemble enter the master LP directly as linear inequalities, which ensures that the combined predictor is feasible whenever the LP is feasible\. Each new column is obtained by solving a*pricing problem*that minimizes the reduced cost of the LP: a quantity combining the supervised training loss with the dual values of the most violated constraints\. The dual variables act as importance weights during training, directing the new model toward the inputs where the current ensemble fails to satisfy the constraints\. When the set of constraints is large, a cutting\-plane procedure adds only the most violated constraints at each iteration, keeping the LP tractable\.

The key advantage of CG4AI over penalty\-based and projection\-based approaches is that constraint satisfaction is enforced at the level of the LP, not through soft penalties or post\-hoc corrections\. Once the LP is feasible, the convex combination of models is guaranteed to satisfy all enforced constraints, independently of the architecture of the individual models\. The CG4AI considers any kind of AI model as sub\-problem while it accepts a custom loss function\.

#### Contributions

The proposed mechanism of adding output satisfaction constraints to a column\-generation LP allows us to make the following contributions\.

1. 1\.CG4AI: LP\-certified ensemble learning under output constraints\.We propose a general framework that builds a convex ensemble of AI models whose combined output satisfies a given set of linear constraints\. A master LP \(with the same column\-generation structure as LPBoost\([Demiriz et al\., 2002](https://arxiv.org/html/2608.26375#bib.bib4)\), but a different objective and a new family of output constraints\) determines the mixture weights; a pricing subproblem trains new models guided by LP dual variables\. This extends LPBoost in three ways: \(a\) the master LP contains*external*constraints on the combined outputs, encoding safety, physical, or user\-defined requirements; \(b\) the constraints are associated with AI models inputs that can be entirely disjoint from the training data, separating what the model*must*satisfy from what it*optimizes*; \(c\) multiple constraint types are naturally handled, including adversarial robustness, correction of misclassified examples, and output relabeling, without modifying the framework\.
2. 2\.LP\-based training vs\. inference\-time projection\.Unlike HardNet\([Min and Azizan, 2024](https://arxiv.org/html/2608.26375#bib.bib1)\)and RAYEN\([Tordesillas et al\., 2023](https://arxiv.org/html/2608.26375#bib.bib2)\), which project model outputs onto a feasible set at every forward pass, CG4AI enforces constraints during training through the LP master\. Once built, the Convex Ensemble of AI models \(CE\) inherits feasibility from the LP with zero runtime overhead at inference\.
3. 3\.Cutting\-plane extension to continuous input spaces\.A separation oracle searches for the most violated input in the constraint data set and adds it as a new LP row, extending the feasibility guarantee beyond any finite constraint set\.
4. 4\.Two applications\.On MNIST \(proof of concept\), we illustrate four constraint types: learning from constraints alone, adversarial robustness, misclassification correction, and output relabeling\. On the multi\-commodity flow problem \(main use case\), we enforce link capacity constraints on neural\-network routing predictors across 41 SNDLIB benchmark configurations\.

The paper is organized as follows\. Section[2](https://arxiv.org/html/2608.26375#S2)reviews related work\. Section[3](https://arxiv.org/html/2608.26375#S3)describes the CG4AI framework\. Sections[4](https://arxiv.org/html/2608.26375#S4)and[5](https://arxiv.org/html/2608.26375#S5)present the two applications\. Section[6](https://arxiv.org/html/2608.26375#S6)provides conclusions and new research directions\.

### 2Related Work

We review five lines of work that are directly relevant to CG4AI: penalty and Lagrangian methods, differentiable optimization layers, projection and architectural approaches, MIP\-based formulations, and the use of column generation in machine learning\.

##### Penalty and Lagrangian methods\.

The most common approach for incorporating constraints in learning adds penalty terms for violations to the training loss\. While easy to implement, these soft\-constraint methods cannot guarantee feasibility at inference time and require careful tuning of penalty weights\([Goodfellow et al\., 2016](https://arxiv.org/html/2608.26375#bib.bib5)\)\. To provide a more principled treatment,[Chamon and Ribeiro \(2020\)](https://arxiv.org/html/2608.26375#bib.bib6)developed a generalization theory for constrained learning based on the PAC framework, showing that any PAC\-learnable class is also learnable under additional constraints, using dual\-ascent algorithms\.[Elenter et al\. \(2024\)](https://arxiv.org/html/2608.26375#bib.bib7)extended this work and showed that primal iterates of dual\-ascent methods can be near\-feasible and near\-optimal in over\-parameterized regimes\. Augmented Lagrangian methods have been applied to physics\-informed neural networks to enforce hard equality constraints\([Lu et al\., 2021](https://arxiv.org/html/2608.26375#bib.bib8)\), but they do not naturally generalize to inequality or combinatorial constraints\.

##### Differentiable optimization layers\.

[Amos and Kolter \(2017\)](https://arxiv.org/html/2608.26375#bib.bib3)introduced OptNet, which embeds a convex quadratic program as a differentiable layer within a neural network, allowing constraint\-aware end\-to\-end training\. Feasibility is enforced for inputs seen during the forward pass, but not for arbitrary unseen inputs, and the approach becomes costly at scale\.[Agrawal et al\. \(2019\)](https://arxiv.org/html/2608.26375#bib.bib9)generalized this to arbitrary convex cone programs through the CVXPYlayers framework\.[Donti et al\. \(2021\)](https://arxiv.org/html/2608.26375#bib.bib10)proposed DC3, which enforces hard constraints via a differentiable completion\-and\-correction procedure: equality constraints are satisfied by implicit completion, while inequality constraints are handled by unrolled gradient steps\. DC3 does not guarantee feasibility in general and is sensitive to the number of gradient steps\.

##### Projection and architectural methods\.

[Min and Azizan \(2024\)](https://arxiv.org/html/2608.26375#bib.bib1)proposed HardNet, a projection\-based framework that guarantees constraint satisfaction by projecting network outputs onto the feasible region after each forward pass, with universal approximation guarantees\.[Tordesillas et al\. \(2023\)](https://arxiv.org/html/2608.26375#bib.bib2)introduced RAYEN, which parametrizes a feasible convex region analytically and projects efficiently onto it without iterative solvers, focusing on continuous convex constraints\.[Balestriero and LeCun \(2023\)](https://arxiv.org/html/2608.26375#bib.bib11)proposed POLICE, which enforces affine output constraints by construction through a reformulation of the network as a continuous piecewise affine mapping over a given region; this avoids any runtime overhead but is limited to a single convex region and affine constraints\.[Constante Flores et al\. \(2025\)](https://arxiv.org/html/2608.26375#bib.bib12)combine a task network \(minimizing loss\) and a safe network \(satisfying constraints\) through a fixed convex combination, where feasibility of the combination is guaranteed by robust optimization duality\. All these projection and architectural methods are limited to convex or affine constraints and do not scale naturally to combinatorial or high\-dimensional settings\.

##### MIP\-based formulations\.

[Aftabi et al\. \(2025\)](https://arxiv.org/html/2608.26375#bib.bib13)formulate trained neural networks with binary activations as mixed\-integer programs, enabling exact enforcement of combinatorial constraints\. The approach is theoretically complete but does not scale to large networks\. Related work by[Tjeng et al\. \(2019\)](https://arxiv.org/html/2608.26375#bib.bib14)uses MIP formulations to evaluate the robustness of neural networks, and[Anderson et al\. \(2020\)](https://arxiv.org/html/2608.26375#bib.bib15)develop strong relaxations for ReLU networks in the MIP framework\.

##### Column generation in machine learning\.

Column generation has been used in machine learning primarily in the context of boosting\.[Demiriz et al\. \(2002\)](https://arxiv.org/html/2608.26375#bib.bib4)proposed LPBoost, which casts the boosting problem as a linear program and solves it via column generation: the pricing problem generates new weak learners guided by dual misclassification costs\.[Aziz et al\. \(2024\)](https://arxiv.org/html/2608.26375#bib.bib16)extended this idea to deep CNN base learners on CIFAR\-10\.

##### CG4AI as an extension of LPBoost\.

CG4AI builds directly on the LPBoost framework of[Demiriz et al\. \(2002\)](https://arxiv.org/html/2608.26375#bib.bib4)\. Letn∈ℕn\\in\\mathbb\{N\}be the size of the training data set, andIIthe set of all possible AI models\. LPBoost solves the following LP over mixture weightsλi≥0\\lambda\_\{i\}\\geq 0for learnerssi:ℝd→\{−1,\+1\}s\_\{i\}:\\mathbb\{R\}^\{d\}\\to\\\{\-1,\+1\\\}:

max⁡ρ−ν​∑j=1nδjs\.t\.∑i∈Iλi=1,∑i∈Iλi​si​\(x¯j\)​y^j\+δj≥ρ,∀j=1,n¯\\max\\rho\-\\nu\\sum\_\{j=1\}^\{n\}\\delta\_\{j\}\\quad\\text\{s\.t\.\}\\quad\\sum\_\{i\\in I\}\\lambda\_\{i\}=1,\\sum\_\{i\\in I\}\\lambda\_\{i\}\\,s\_\{i\}\(\\bar\{x\}\_\{j\}\)\\,\\hat\{y\}\_\{j\}\+\\delta\_\{j\}\\;\\geq\\;\\rho,\\ \\forall j=\\overline\{1,n\}\(1\)wherex¯j∈ℝd\\bar\{x\}\_\{j\}\\in\\mathbb\{R\}^\{d\}is thejj\-th training example,y^j∈\{−1,\+1\}\\hat\{y\}\_\{j\}\\in\\\{\-1,\+1\\\}its label,ρ\\rhothe minimum margin to be maximized,δj≥0\\delta\_\{j\}\\geq 0the per\-example slack, andν∈ℝ\+\\nu\\in\\mathbb\{R\}^\{\+\}a regularization parameter\. The dual variables of the margin constraints play the role of misclassification weights \(analogous to AdaBoost\), and the column\-generation pricing finds the learnersis\_\{i\}that maximally reduces the weighted training error\.

CG4AI retains the column\-generation structure of LPBoost \(column index setII, mixture weightsλi\\lambda\_\{i\}, convexity constraint\), wheresis\_\{i\}is a general AI model rather than a binary classifier i\.e\.,si:ℝd⇒ℝrs\_\{i\}:\\mathbb\{R\}^\{d\}\\Rightarrow\\mathbb\{R\}^\{r\}\. It differs in two fundamental ways\.*\(i\) Different objective\.*CG4AI minimizes a training loss∑iℒi​λi\\sum\_\{i\}\\mathcal\{L\}\_\{i\}\\lambda\_\{i\}rather than maximizing a margin; this allows any differentiable loss \(cross\-entropy, MSE, etc\.\) and is not restricted to binary classification\.*\(ii\) New family of output constraints\.*In addition to the convexity constraint∑iλi=1\\sum\_\{i\}\\lambda\_\{i\}=1, CG4AI adds rows of the form∑i∈Iλi​∑q=1rajq​siq​\(x¯j\)≥bj\\sum\_\{i\\in I\}\\lambda\_\{i\}\\sum\_\{q=1\}^\{r\}a^\{q\}\_\{j\}s^\{q\}\_\{i\}\(\\bar\{x\}\_\{j\}\)\\geq\\;b\_\{j\}for each pointjjin a constraint set that can be entirely disjoint from the training set\. These constraints encode external requirements on the model’s output \(safety rules, physical limits, user\-defined requirements\) and are formally defined in Section[3\.2](https://arxiv.org/html/2608.26375#S3.SS2)\. The dual variable of a rowjjmeasures how difficult it is to satisfy the constraint atx¯j\\bar\{x\}\_\{j\}and guide the pricing toward the most violated inputs\.

The nature of the guarantee changes accordingly: LPBoost provides a probabilistic generalization bound via large training margins\. CG4AI provides a*deterministic feasibility guarantee*: when the LP is feasible, the ensemble satisfies all output constraints, regardless of whether those inputs appeared during training\. The empirical evidence for this gain is reported in Section[4\.4](https://arxiv.org/html/2608.26375#S4.SS4)\.

##### Distinction from HardNet and projection methods\.

[Min and Azizan \(2024\)](https://arxiv.org/html/2608.26375#bib.bib1)\(HardNet\) and[Tordesillas et al\. \(2023\)](https://arxiv.org/html/2608.26375#bib.bib2)\(RAYEN\) guarantee constraint satisfaction by*projecting*the model’s output onto the feasible set at every forward pass, during both training and inference\. This projection happens after the model computes its output, which means \(i\) a dedicated optimization problem must be solved at inference time for every new input, adding computational overhead; \(ii\) the projection is defined only over a fixed convex or affine feasible region and does not naturally extend to combinatorial or input\-dependent constraints; and \(iii\) the model is not directly trained to satisfy constraints, it is trained to minimize a loss and then corrected\. In CG4AI, constraints are enforced in the LP master, which directly determines the mixture weights that make the ensemble feasible\. No projection is needed at inference: the ensemblef=∑iλi​sif=\\sum\_\{i\}\\lambda\_\{i\}s\_\{i\}inherits feasibility from the LP, and each columnsis\_\{i\}is trained by the pricing to*actively contribute*to constraint satisfaction, guided by the dual variables\. CG4AI is therefore a training\-time rather than an inference\-time method, with no runtime overhead after the ensemble is built\.

CG4AI’s three distinguishing properties are thus: \(1\) hard constraints are enforced in the LP master and hold for any convex combination of the columns; \(2\) the pricing trains new models guided by LP dual variables, directing capacity toward the most violated constraints; and \(3\) the constraint set can differ from the training set, allowing user\-defined requirements to be imposed without retraining from scratch\.

### 3The CG4AI Framework

We consider the following general setting\. In contrast with LPBoost the size of training data and the number of hard constraints may be different\. Letn∈ℕn\\in\\mathbb\{N\}be the size of the training data set andm∈ℕm\\in\\mathbb\{N\}the size of the hard constraints set\. LetℐT​r​a​i​n=\{\(x¯k,y^k\)\}k=1n\\mathcal\{I\}^\{Train\}=\\\{\(\\bar\{x\}\_\{k\},\\hat\{y\}\_\{k\}\)\\\}\_\{k=1\}^\{n\}be a training dataset, wherex¯k∈ℝd\\bar\{x\}\_\{k\}\\in\\mathbb\{R\}^\{d\}is an input andy^k∈ℝr\\hat\{y\}\_\{k\}\\in\\mathbb\{R\}^\{r\}is the corresponding label\. LetℐC​o​n​s​t=\{x¯j\}j=1m\\mathcal\{I\}^\{Const\}=\\\{\\bar\{x\}\_\{j\}\\\}\_\{j=1\}^\{m\}be a set of inputs on which output constraints must hold;ℐC​o​n​s​t\\mathcal\{I\}^\{Const\}may overlap withℐT​r​a​i​n\\mathcal\{I\}^\{Train\}or may be disjoint\. The goal is to find a linear combination of predictorsf:ℝd→ℝrf:\\mathbb\{R\}^\{d\}\\to\\mathbb\{R\}^\{r\}that solves:

min⁡∑k=1nf⁡L⁡\(f⁡\(x¯k\),y^k\)s\.t\.∑q=1rajq​fq​\(x¯j\)≥bj,∀j=1,…,m,\\min\_\{f\}\\;\\sum\_\{k=1\}^\{n\}L\\bigl\(f\(\\bar\{x\}\_\{k\}\),\\hat\{y\}\_\{k\}\\bigr\)\\quad\\text\{s\.t\.\}\\quad\\sum\_\{q=1\}^\{r\}a^\{q\}\_\{j\}f\_\{q\}\(\\bar\{x\}\_\{j\}\)\\;\\geq\\;b\_\{j\},\\quad\\forall j=1,\\ldots,m,\(2\)
whereajq∈ℝa^\{q\}\_\{j\}\\in\\mathbb\{R\}is the coefficient associated toqt​hq^\{th\}output offf,b∈ℝmb\\in\\mathbb\{R\}^\{m\}, andLLis a loss function \(e\.g\., cross\-entropy or mean squared error\)\. This problem is generally intractable for expressive model classes such as deep neural networks because the constraints impose nonlinear requirements on the model parameters\.

#### 3\.1Convex Ensemble Structure

Letℳ\\mathcal\{M\}be a large \(potentially infinite\) set of candidate models\. Each modelsi:ℝd→ℝrs\_\{i\}:\\mathbb\{R\}^\{d\}\\to\\mathbb\{R\}^\{r\}has its own parameters \(weights, biases, etc\.\)\. The combined predictor is defined as a convex combination:

f⁡\(⋅\)=∑i∈ℳλi​si​\(⋅\),λi≥0,∑i∈ℳλi=1\.f\(\\cdot\)\\;=\\;\\sum\_\{i\\in\\mathcal\{M\}\}\\lambda\_\{i\}\\,s\_\{i\}\(\\cdot\),\\quad\\lambda\_\{i\}\\geq 0,\\quad\\sum\_\{i\\in\\mathcal\{M\}\}\\lambda\_\{i\}=1\.\(3\)This ensemble structure has two key properties, stated as a theorem\.

###### Theorem 1\(Tractability of the convex ensemble\)\.

Letf=∑i∈ℳλi​sif=\\sum\_\{i\\in\\mathcal\{M\}\}\\lambda\_\{i\}s\_\{i\}withλi≥0\\lambda\_\{i\}\\geq 0,∑iλi=1\\sum\_\{i\}\\lambda\_\{i\}=1\.

1. \[\(i\)\]
2. 1\.*\(Loss bound\.\)*IfL⁡\(⋅,y^\)L\(\\cdot,\\hat\{y\}\)is convex in its first argument, then by Jensen’s inequality[Jensen \(1906\)](https://arxiv.org/html/2608.26375#bib.bib17): L⁡\(∑i∈ℳλi​si​\(x¯\),y^\)≤∑i∈ℳλi​L​\(si​\(x¯\),y^\)\.L\\\!\\Bigl\(\\sum\_\{i\\in\\mathcal\{M\}\}\\lambda\_\{i\}\\,s\_\{i\}\(\\bar\{x\}\),\\,\\hat\{y\}\\Bigr\)\\;\\leq\\;\\sum\_\{i\\in\\mathcal\{M\}\}\\lambda\_\{i\}\\,L\\bigl\(s\_\{i\}\(\\bar\{x\}\),\\hat\{y\}\\bigr\)\.\(4\)Minimizing the right\-hand side, linear inλ\\lambdafor fixedsis\_\{i\}also reduces the actual mixture loss\.
3. 2\.*\(Constraint linearization\.\)*For anya∈ℝr×ma\\in\\mathbb\{R\}^\{r\\times m\}, and fixed modelssis\_\{i\}: ∑q=1rajq​\(∑i∈ℳλi​si​\(x¯j\)\)≥bj⇔∑i∈ℳλi​\(∑q=1rajq​si​\(x¯j\)\)≥bj\.\\sum\_\{q=1\}^\{r\}a^\{q\}\_\{j\}\\\!\\Bigl\(\\sum\_\{i\\in\\mathcal\{M\}\}\\lambda\_\{i\}\\,s\_\{i\}\(\\bar\{x\}\_\{j\}\)\\Bigr\)\\;\\geq\\;b\_\{j\}\\;\\iff\\;\\sum\_\{i\\in\\mathcal\{M\}\}\\lambda\_\{i\}\\,\\bigl\(\\sum\_\{q=1\}^\{r\}a^\{q\}\_\{j\}s\_\{i\}\(\\bar\{x\}\_\{j\}\)\\bigr\)\\;\\geq\\;b\_\{j\}\.\(5\)This equivalence follows from the linearity of matrix\-vector multiplication:∑q=1rajq​\(∑iλi​vi\)=∑iλi​\(∑q=1rajq​vi\)\\sum\_\{q=1\}^\{r\}a^\{q\}\_\{j\}\(\\sum\_\{i\}\\lambda\_\{i\}v\_\{i\}\)=\\sum\_\{i\}\\lambda\_\{i\}\(\\sum\_\{q=1\}^\{r\}a^\{q\}\_\{j\}v\_\{i\}\)for any vectorsvi∈ℝnv\_\{i\}\\in\\mathbb\{R\}^\{n\}and scalarsλi\\lambda\_\{i\}\. The output constraints thus become*linear*in the weightsλ\\lambdafor fixed models, making the master problem tractable as an LP\.

In the following table, we summarize the notations used in the remaining of the paper\.

Table 1:Summary of notationsThe key observation is that \(i\) the loss upper bound and \(ii\) the constraint linearization together reduce the original intractable problem \([2](https://arxiv.org/html/2608.26375#S3.E2)\) to a linear program in the weightsyy\.

#### 3\.2Master Linear Program

Given a set of trained models\{si:i∈ℳ\}\\\{s\_\{i\}:i\\in\\mathcal\{M\}\\\}, the mixture weights are determined by solving the following LP:

min\\displaystyle\\min\\quad∑i∈ℳℒi​λi\+Ms​∑j=1mδj\\displaystyle\\sum\_\{i\\in\\mathcal\{M\}\}\\mathcal\{L\}\_\{i\}\\,\\lambda\_\{i\}\\;\+\\;M\_\{s\}\\sum\_\{j=1\}^\{m\}\\delta\_\{j\}\(6\)s\.t\.∑i∈ℳλi=1,\[γ\]\\displaystyle\\sum\_\{i\\in\\mathcal\{M\}\}\\lambda\_\{i\}=1,\\qquad\\qquad\[\\gamma\]\(7\)∑i∈ℳλi∑q=1rajqsiq\(x¯j\)\+δj≥bj,∀j=1,…,m,\[βj\]\\displaystyle\\sum\_\{i\\in\\mathcal\{M\}\}\\lambda\_\{i\}\\sum\_\{q=1\}^\{r\}a^\{q\}\_\{j\}s^\{q\}\_\{i\}\(\\bar\{x\}\_\{j\}\)\\;\+\\;\\delta\_\{j\}\\;\\geq\\;b\_\{j\},\\qquad\\forall j=1,\\ldots,m,\\quad\[\\beta\_\{j\}\]\(8\)λi≥0,∀i∈ℳ,\\displaystyle\\lambda\_\{i\}\\geq 0,\\qquad\\qquad\\forall i\\in\\mathcal\{M\},\(9\)δj≥0,∀j∈1,…,m,\\displaystyle\\delta\_\{j\}\\geq 0,\\qquad\\qquad\\forall j\\in 1,\\dots,m,\(10\)whereℒi=∑k=1nL⁡\(si​\(x¯k\),y^k\)\\mathcal\{L\}\_\{i\}=\\sum\_\{k=1\}^\{n\}L\(s\_\{i\}\(\\bar\{x\}\_\{k\}\),\\hat\{y\}\_\{k\}\)is the training loss of modeliionℐT​r​a​i​n\\mathcal\{I\}^\{Train\},Ms≫1M\_\{s\}\\gg 1is a large feasibility penalty\. The slack variableδ∈ℝ\+m\\delta\\in\\mathbb\{R\}\_\{\+\}^\{m\}is anmm\-dimensional vector \(one component per hard constraint\); the termδ\\deltaaggregates violations into a scalar penalty\. In our applications \(Sections[4](https://arxiv.org/html/2608.26375#S4)and[5](https://arxiv.org/html/2608.26375#S5)\), each constraint is scalar \(m=1m=1\), soδj∈ℝ\+\\delta\_\{j\}\\in\\mathbb\{R\}\_\{\+\}\. Dual variablesγ∈ℝ\\gamma\\in\\mathbb\{R\}\(equality\) andβ∈ℝ\+m\\beta\\in\\mathbb\{R\}\_\{\+\}^\{m\}\(≥\\geqinequality\) are associated with constraints \([7](https://arxiv.org/html/2608.26375#S3.E7)\) and \([8](https://arxiv.org/html/2608.26375#S3.E8)\), respectively\.

Design of the slack variablesδj\\delta\_\{j\}\.When the constraint setℐC​o​n​s​t=\{x¯1,…,x¯m\}\\mathcal\{I\}^\{Const\}=\\\{\\bar\{x\}\_\{1\},\\ldots,\\bar\{x\}\_\{m\}\\\}is given explicitly, we introduce one slackδj≥0\\delta\_\{j\}\\geq 0per inputx¯j\\bar\{x\}\_\{j\}, yieldingmmslack variables\. WhenℐC​o​n​s​t\\mathcal\{I\}^\{Const\}is large or continuous and constraints are added incrementally by the cutting\-plane procedure \(Algorithm[2](https://arxiv.org/html/2608.26375#alg2)\), we use a*single*shared slack variable representing the worst\-case violation over all active constraints\. This avoids introducing a new variable for each generated constraint and keeps the master LP compact as the cutting\-plane loop progresses\. The shared slack is sufficient because the separation step \(Eq\. \([13](https://arxiv.org/html/2608.26375#S3.E13)\)\) identifies the single most violated input, so only one slack needs to be driven to zero at each outer iteration\.

Constraint \([7](https://arxiv.org/html/2608.26375#S3.E7)\) requires that the mixture weights form a valid convex combination; its dual variableγ\\gammarepresents the marginal value of adding one more model to the ensemble and appears in the reduced cost\. Constraint \([8](https://arxiv.org/html/2608.26375#S3.E8)\) requires that the ensemble’s weighted prediction satisfies the output constraints at each pointx¯j∈ℐC​o​n​s​t\\bar\{x\}\_\{j\}\\in\\mathcal\{I\}^\{Const\}; its dual variableβj∈ℝ\+m\\beta\_\{j\}\\in\\mathbb\{R\}\_\{\+\}^\{m\}quantifies how difficult it is to satisfy constraintjjgiven the current columns, a largeβj\\beta\_\{j\}signals a highly violated constraint and increases the weight placed on inputx¯j\\bar\{x\}\_\{j\}during the next pricing step\.

When all slacks are zero, the ensemble satisfies all constraints inℐC​o​n​s​t\\mathcal\{I\}^\{Const\}\. The LP is always feasible \(slacks absorb any violation\), and the penaltyMsM\_\{s\}ensures that feasibility is prioritized over loss minimization\.

#### 3\.3Column Generation

The setℳ\\mathcal\{M\}is typically exponentially large\. Column generation\([Barnhart et al\., 1998](https://arxiv.org/html/2608.26375#bib.bib18);[Demiriz et al\., 2002](https://arxiv.org/html/2608.26375#bib.bib4)\)maintains a restricted master LP \(RMP\) with a subsetℳ′⊆ℳ\\mathcal\{M\}^\{\\prime\}\\subseteq\\mathcal\{M\}of active models and iteratively adds columns \(models\) with negative reduced cost\. The*reduced cost*of modeliiis:

ri=ℒi−∑j=1mβj⊤​∑q=1rajq​siq​\(x¯j\)−γ,r\_\{i\}\\;=\\;\\mathcal\{L\}\_\{i\}\\;\-\\;\\sum\_\{j=1\}^\{m\}\\beta\_\{j\}^\{\\top\}\\,\\sum\_\{q=1\}^\{r\}a^\{q\}\_\{j\}s\_\{i\}^\{q\}\(\\bar\{x\}\_\{j\}\)\\;\-\\;\\gamma,\(11\)whereβj\\beta\_\{j\}andγ\\gammaare the dual variables of the current RMP solution\. A model withri<0r\_\{i\}<0can improve the objective and should be added to the ensemble\.

The*pricing problem*finds the model that minimizes the reduced cost:

minw∈W⁡r⁡\(w\)=ℒ⁡\(w\)−∑j=1mβj⊤​∑q=1rajq​swq​\(x¯j\)−γ,\\min\_\{w\\in W\}\\;r\(w\)\\;=\\;\\mathcal\{L\}\(w\)\\;\-\\;\\sum\_\{j=1\}^\{m\}\\beta\_\{j\}^\{\\top\}\\sum\_\{q=1\}^\{r\}a^\{q\}\_\{j\}s\_\{w\}^\{q\}\(\\bar\{x\}\_\{j\}\)\\;\-\\;\\gamma,\(12\)wherewwdenotes the model weights\. This is a training problem: we seek a model with low training loss that also helps satisfy the constraints\. The dual variableβj\\beta\_\{j\}controls how much weight is placed on constraintjj: inputs with largeβj\\beta\_\{j\}receive more attention during pricing, focusing training on the most violated constraints\.

Problem \([12](https://arxiv.org/html/2608.26375#S3.E12)\) is dependent on the considered AI model\. In practice, libraries like PyTorch and scikit\-learn provide efficient algorithms for the learning of multiple types of AI models\. Whensws\_\{w\}is differentiable, the reduced cost is differentiable with respect towwand standard backpropagation applies\. For linear \(affine\) models, problem \([12](https://arxiv.org/html/2608.26375#S3.E12)\) is convex inwwand the optimizer finds a global minimum, so the CG stopping criterion is theoretically valid\. For nonlinear models \(e\.g\., networks with ReLU activations\), the pricing is non\-convex and a local optimizer finds only a local minimum; in this case Algorithm[1](https://arxiv.org/html/2608.26375#alg1)is a*heuristic*that generates improving columns as long as the local search finds one, without the LP optimality guarantee of exact CG\. The full algorithm is given in Algorithm[1](https://arxiv.org/html/2608.26375#alg1)\.

Algorithm 1CG4AI: Column Generation for Constrained Ensembles1:Initialize

ℳ′←∅\\mathcal\{M\}^\{\\prime\}\\leftarrow\\emptyset\. Add a feasibility dummy column to

ℳ′\\mathcal\{M\}^\{\\prime\}\(large loss, zero constraint contribution\)\.

2:repeat

3:Master step\.Solve the RMP over

ℳ′\\mathcal\{M\}^\{\\prime\}to obtain dual variables

βj\\beta\_\{j\},

γ\\gammaand primal weights

λi\\lambda\_\{i\}\.

4:Pricing step\.Train a new model

snews\_\{\\mathrm\{new\}\}by minimizing \([12](https://arxiv.org/html/2608.26375#S3.E12)\), guided by

βj\\beta\_\{j\}\.

5:if

r⁡\(wnew\)<−εr\(w\_\{\\mathrm\{new\}\}\)<\-\\varepsilonthen

6:Add

snews\_\{\\mathrm\{new\}\}to

ℳ′\\mathcal\{M\}^\{\\prime\}\.

7:endif

8:until

r⁡\(wnew\)≥−εr\(w\_\{\\mathrm\{new\}\}\)\\geq\-\\varepsilonand

∑jδj=0\\sum\_\{j\}\\delta\_\{j\}=0\(all constraints satisfied\)

9:Output\.Ensemble

f=∑i∈ℳ′λi​sif=\\sum\_\{i\\in\\mathcal\{M\}^\{\\prime\}\}\\lambda\_\{i\}\\,s\_\{i\}with weights from the final RMP\.

The algorithm stops when no column with negative reduced cost exists*and*the current solution is feasible\. If the column limit is reached before feasibility, the output is still a well\-defined convex combination, but some constraints may not be satisfied; the total violation is measured by∑jδj\\sum\_\{j\}\\delta\_\{j\}\.

#### 3\.4Cutting Planes for Large Constraint Sets

WhenℐC​o​n​s​t\\mathcal\{I\}^\{Const\}is large or continuous, it is impractical to include all constraints in the RMP\. We use a cutting\-plane approach\([Kelley, 1960](https://arxiv.org/html/2608.26375#bib.bib19)\): start with a small subsetℐ0C​o​n​s​t⊆ℐC​o​n​s​t\\mathcal\{I\}^\{Const\}\_\{0\}\\subseteq\\mathcal\{I\}^\{Const\}, and iteratively add the constraint corresponding to the most violated input\.

Given the current ensemble weightsλ∗\\lambda^\{\*\}, the*separation problem*is:

x¯∗=arg⁡maxx¯∈ℐC​o​n​s​t​maxj=1,…,m⁡\(cj−∑i∈ℳ′λi∗​∑q=1rajq​siq​\(x¯\)\),\\bar\{x\}^\{\*\}\\;=\\;\\arg\\max\_\{\\bar\{x\}\\in\\mathcal\{I\}^\{Const\}\}\\;\\max\_\{j=1,\\ldots,m\}\\Bigl\(c\_\{j\}\-\\sum\_\{i\\in\\mathcal\{M\}^\{\\prime\}\}\\lambda\_\{i\}^\{\*\}\\,\\sum\_\{q=1\}^\{r\}a^\{q\}\_\{j\}s\_\{i\}^\{q\}\(\\bar\{x\}\)\\Bigr\),\(13\)This finds the input inℐC​o​n​s​t\\mathcal\{I\}^\{Const\}that is most violated under the current ensemble\. If the maximum is non\-positive, the ensemble is feasible for all ofℐC​o​n​s​t\\mathcal\{I\}^\{Const\}\. In general the separation problem is nonconvex and is solved heuristically by a black\-box optimizer\. Algorithm[2](https://arxiv.org/html/2608.26375#alg2)combines column generation with cutting planes\.

Algorithm 2CG4AI with Cutting Planes1:Initialize

ℐRMPC​o​n​s​t←ℐ0C​o​n​s​t\\mathcal\{I\}^\{Const\}\_\{\\mathrm\{RMP\}\}\\leftarrow\\mathcal\{I\}^\{Const\}\_\{0\}\(small initial constraint set\)\.

2:repeat

3:Run Algorithm[1](https://arxiv.org/html/2608.26375#alg1)with

ℐRMPC​o​n​s​t\\mathcal\{I\}^\{Const\}\_\{\\mathrm\{RMP\}\}to get mixture weights

λ∗\\lambda^\{\*\}\.

4:Separation step\.Solve \([13](https://arxiv.org/html/2608.26375#S3.E13)\) to find the most violated input

x¯∗\\bar\{x\}^\{\*\}\.

5:if

OPEN\(cj−∑i∈ℳ′λi∗​∑q=1rajq​siq​\(x¯∗\)\)\)\>ε\\Bigl\(c\_\{j\}\-\\sum\_\{i\\in\\mathcal\{M\}^\{\\prime\}\}\\lambda\_\{i\}^\{\*\}\\,\\sum\_\{q=1\}^\{r\}a^\{q\}\_\{j\}s\_\{i\}^\{q\}\(\\bar\{x\}^\{\*\}\)\\Bigr\)\)\>\\varepsilonthen

6:Add

x¯∗\\bar\{x\}^\{\*\}to

ℐRMPC​o​n​s​t\\mathcal\{I\}^\{Const\}\_\{\\mathrm\{RMP\}\}\.

7:endif

8:untilno violated input found

9:Output\.Ensemble

f=∑i∈ℳ′λi∗​sif=\\sum\_\{i\\in\\mathcal\{M\}^\{\\prime\}\}\\lambda^\{\*\}\_\{i\}\\,s\_\{i\}\.

The combined procedure alternates between pricing \(new columns\), solving the master LP \(new weights\), and separation \(new constraints\)\. When all three steps report no improvement, the solution is both LP\-optimal and feasible for the full setℐC​o​n​s​t\\mathcal\{I\}^\{Const\}\.

##### MIP\-based exact separation\.

When the columnssis\_\{i\}are single\-hidden\-layer ReLU networks, both the pricing and the separation problems can be solved to global optimality via mixed\-integer programming\([Aftabi et al\., 2025](https://arxiv.org/html/2608.26375#bib.bib13);[Anderson et al\., 2020](https://arxiv.org/html/2608.26375#bib.bib15)\)\. For separation specifically, the MIP encodes the ensemble output constraints and searches for the inputx¯∗∈ℐC​o​n​s​t\\bar\{x\}^\{\*\}\\in\\mathcal\{I\}^\{Const\}that maximally violates them\. This yields*provably optimal cuts*, tightening the master LP as much as possible at each outer iteration\. In the MCF application \(Section[5](https://arxiv.org/html/2608.26375#S5)\), the black\-box NLopt optimizer[Johnson et al\. \(2014\)](https://arxiv.org/html/2608.26375#bib.bib20)suffices empirically; the MIP approach provides a path to certified optimality for applications requiring formal guarantees on the cutting\-plane procedure\.

### 4Proof of Concept: Digit Classification

We illustrate the CG4AI mechanisms on MNIST\([LeCun et al\., 1998](https://arxiv.org/html/2608.26375#bib.bib21)\)using a deliberately*minimalist*MLP \(h=4h=4hidden neurons\)\. This choice is intentional: with only 4 neurons, a single model achieves modest accuracy \(52−77%52\{\-\}77\\%depending on training size\), making the accuracy gains from adding constraints clearly measurable\. It also represents the regime of embedded or resource\-constrained systems where model capacity is limited by design\. More expressive architectures would reduce the absolute gains but the constraint\-enforcement mechanism is architecture\-agnostic\. The section is didactic and concise; the main quantitative results are in Section[5](https://arxiv.org/html/2608.26375#S5)\.

#### 4\.1Column Model and Constraint Formulations

Let us introduceλi​sic​\(x\)\\lambda\_\{i\}s\_\{i\}^\{c\}\(x\)the probability of the classccassociated to inputxxon the modelii\. For an inputxjx\_\{j\},cjc\_\{j\}is the right class of the inputxjx\_\{j\}\.

##### Architecture

Each column mapsx∈ℝ784x\\\!\\in\\\!\\mathbb\{R\}^\{784\}to 10 logits viasi​\(x\)=Wi\(2\)​ReLU​\(Wi\(1\)​x\+bi\(1\)\)\+bi\(2\)s\_\{i\}\(x\)=W\_\{i\}^\{\(2\)\}\\mathrm\{ReLU\}\(W\_\{i\}^\{\(1\)\}x\+b\_\{i\}^\{\(1\)\}\)\+b\_\{i\}^\{\(2\)\}, withWi\(1\)∈ℝ4×784W\_\{i\}^\{\(1\)\}\\\!\\in\\\!\\mathbb\{R\}^\{4\\times 784\},Wi\(2\)∈ℝ10×4W\_\{i\}^\{\(2\)\}\\\!\\in\\\!\\mathbb\{R\}^\{10\\times 4\}\. The ensemble predictsarg⁡max⁡∑ic⁡λi​sic​\(x\)\\arg\\max\_\{c\}\\sum\_\{i\}\\lambda\_\{i\}s\_\{i\}^\{c\}\(x\)\. Columns are trained with Adam \(η=10−3\\eta=10^\{\-3\}, up to 1000 epochs, early stopping at patience 50\)\. Figure[1](https://arxiv.org/html/2608.26375#S4.F1)shows the architecture\.

⋮\\vdotsInputd=784d=784Hiddenh=4h=4, ReLU⋮\\vdotsOutputq=10q=10class 0class 1class 9⋯~~\\cdotsWi\(1\)∈ℝ4×784W\_\{i\}^\{\(1\)\}\\\!\\in\\\!\\mathbb\{R\}^\{4\\times 784\}Wi\(2\)∈ℝ10×4W\_\{i\}^\{\(2\)\}\\\!\\in\\\!\\mathbb\{R\}^\{10\\times 4\}Figure 1:Architecture of column modelsis\_\{i\}used in the MNIST experiments: one hidden layer withh=4h=4neurons and ReLU activation, mapping a784784\-dimensional input to1010class logits\. The ensemble output is∑i∈ℳλi​si​\(x\)\\sum\_\{i\\in\\mathcal\{M\}\}\\lambda\_\{i\}s\_\{i\}\(x\)\.When the constraint matrixAAencodes per\-sample classification requirements, the feasibility of the LP translates into guarantees onℐC​o​n​s​t\\mathcal\{I\}^\{Const\}\. The two formulations below each use their constraint coefficient as the ensemble prediction signal, making the guarantee self\-consistent and the proof immediate\. For a pointj∈\{1,…,m\}j\\in\\\{1,\.\.\.,m\\\}:

- •CG\-proba: require∑i∈ℳλi​\(sicj​\(x¯j\)\)≥δ\\sum\_\{i\\in\\mathcal\{M\}\}\\lambda\_\{i\}\\,\(s\_\{i\}^\{c\_\{j\}\}\(\\bar\{x\}\_\{j\}\)\)\\geq\\delta, it corresponds to an absolute probability needs\. Thus ifδ\>0\.5\\delta\>0\.5thenx¯j\\bar\{x\}\_\{j\}is classified incjc\_\{j\}\.
- •CG\-margin: require∑i∈ℳλi​\[\(sicj​\(x¯j\)\)−maxc≠cj⁡sic​\(x¯j\)\]≥δ\\sum\_\{i\\in\\mathcal\{M\}\}\\lambda\_\{i\}\[\(s\_\{i\}^\{c\_\{j\}\}\(\\bar\{x\}\_\{j\}\)\)\-\\max\_\{c\\neq c\_\{j\}\}s\_\{i\}^\{c\}\(\\bar\{x\}\_\{j\}\)\]\\geq\\delta, it corresponds to a relative probability needs\. Thus ifδ\>0\\delta\>0then the probability to classifyx¯j\\bar\{x\}\_\{j\}incjc\_\{j\}is better than all other classes\.

#### 4\.2Scenarios

##### Scenario 1 — Hard set / optimization set \(key feature\)\.

The LP enforces hard constraints on inputs in hard setℐC​o​n​s​t\\mathcal\{I\}^\{Const\}while the pricing minimizes cross\-entropy on optimization setℐT​r​a​i​n\\mathcal\{I\}^\{Train\}\. WhenℐT​r​a​i​n=∅\\mathcal\{I\}^\{Train\}=\\emptysetthe model learns entirely from the dual values\. We test four splits of the per\-class budget:\(0,10\),\(2,8\),\(5,5\),\(8,2\)\(0,10\),\(2,8\),\(5,5\),\(8,2\)where the first number is the size of the Hard Set and the second is the size of the Optimization Set\. This scenario allows considering a hierarchy between two sets of data\.

##### Scenario 2 — Adversarial robustness \(CG\-robust\)\.

After training a base column, greedy saliency\-based pixel flips \(up toPPpixels;P=100P=100in our experiments, but any black\-box attack can replace this\) generate adversarial perturbationsx~j\\tilde\{x\}\_\{j\}\. The salience score is the modification of probability when a pixel changes\. Each misclassifiedx~j\\tilde\{x\}\_\{j\}adds the constraint CG\-proba or CG\-margin toℐC​o​n​s​t\\mathcal\{I\}^\{Const\}\. A black\-box optimizer searching for the minimum perturbation that breaks the prediction can tighten these constraints and reduce cutting\-plane iterations\. Figure[2](https://arxiv.org/html/2608.26375#S4.F2)shows example perturbed digits\. Cutting plane is used to generate constraints\.

![Refer to caption](https://arxiv.org/html/2608.26375v1/1-modified.png)

![Refer to caption](https://arxiv.org/html/2608.26375v1/2-modified.png)

![Refer to caption](https://arxiv.org/html/2608.26375v1/3-modified.png)

Figure 2:Three adversarially perturbed digits \(up to 100 pixels flipped\)\. The models before adding constraints predict 2, 3, and 4 respectively, but the original labels are different\. CG4AI adds constraints to enforce correct classification of these perturbed inputs\.This scenario allows us to check if there exists an AI model able to classify digits and be robust to the greedy saliency\-based pixel flips\. Note that this approach can be considered only with optimization approachℐT​r​a​i​n\\mathcal\{I\}^\{Train\}\.

##### Scenario 3 — Correcting misclassified examples \(CG\-mis\)\.

Letℰ\\mathcal\{E\}be training images misclassified by the base column\. Adding the constraint CG\-proba or CG\-margin for each element ofℰ\\mathcal\{E\}forces the ensemble to correct them\. Figure[3](https://arxiv.org/html/2608.26375#S4.F3)shows two corrected digits\.

![Refer to caption](https://arxiv.org/html/2608.26375v1/1-digit.png)

![Refer to caption](https://arxiv.org/html/2608.26375v1/4-digit.png)

Figure 3:Two training images initially misclassified: digit 1 is predicted as 3, and digit 4 is predicted as 0\. Adding hard constraints for these images forces the ensemble to classify them correctly\.This scenario guarantees that misclassified input spotted by a user for instance must be well classified\.

##### Scenario 4 — Output relabeling \(CG\-relabel\)\.

To force classc′c^\{\\prime\}for all inputs in a target setI′⊆II^\{\\prime\}\\subseteq I, add∑i∈ℳλi​sic′​\(x¯j\)≥θ\\sum\_\{i\\in\\mathcal\{M\}\}\\lambda\_\{i\}\\,s\_\{i\}^\{c^\{\\prime\}\}\(\\bar\{x\}\_\{j\}\)\\geq\\thetafor eachj∈I′j\\in I^\{\\prime\}\. The original training labels are unchanged; the LP imposes the relabeling through new constraint rows, usually requiring only a few additional columns\. Figure[4](https://arxiv.org/html/2608.26375#S4.F4)shows digit\-3 images relabeled as digit 2\.

![Refer to caption](https://arxiv.org/html/2608.26375v1/0-digit_force.png)

![Refer to caption](https://arxiv.org/html/2608.26375v1/3-digit_force.png)

![Refer to caption](https://arxiv.org/html/2608.26375v1/3-digit_force_1.png)

![Refer to caption](https://arxiv.org/html/2608.26375v1/2-digit_force.png)

Figure 4:Relabeling example: three digits \(0, 3, 3, 2\) are shown after adding constraints forcing digit 3 to be predicted as 2\. The digit\-0 image is predicted as 0 by 50%, as 2 by 42%, and as 4 by 8%\. The first 3 is predicted as 2 by 54%, the second as 2 by 100%, and the 2 as 2 by 100%\.This last scenario guarantees some inputs are relabeled, it is interesting, for instance, when some classes must be merged\.

#### 4\.3Results

Training sizes range fromn=10n=10ton=1000n=1000\(balanced across 10 classes\)\. Test set: 2000 images \(200 per class\)\. LP solver: CPLEX 12\.6\.3[IBM Corporation \(2015\)](https://arxiv.org/html/2608.26375#bib.bib22); column limit: 400\.

Figure 5:Test accuracy on 2000 MNIST images vs\. training set size\. CG\-margin and CG\-proba \(constraints only,ko=0k\_\{o\}=0\) outperform the baseline at every scale\. The “2 hard \+ 8 train” split provides a guaranteed\-correct hard set with a small accuracy cost\.Figure 6:Effect of the hard/optimization split on test accuracy \(CG\-margin left, CG\-proba right\)\. All CG4AI variants surpass the baseline\. Increasingkhk\_\{h\}extends the guarantee coverage but reduces overall test accuracy\.Table 2:Test accuracy \(%\) on 2000 MNIST images \(MLP,h=4h=4, all\-constrained mode unless stated\)\.khk\_\{h\}: images per class in hard setℐC​o​n​s​t\\mathcal\{I\}^\{Const\}\. Single run \(seed 42\)\. Note: the comparison in Table[3](https://arxiv.org/html/2608.26375#S4.T3)uses a separate 3\-seed benchmark \(seeds 42, 43, 44\), giving a slightly different baseline \(56\.1%56\.1\\%atn=100n=100\) due to seed variability\.All CG4AI all\-constrained variants outperform the baseline at every training size\. Atn=1000n=1000, CG\-proba reaches85\.7%85\.7\\%\(\+8\.4\+8\.4pp\) of accuracy\. Withkh=2k\_\{h\}=2, CG\-margin reaches78\.7%78\.7\\%\(above baseline\) while the LP guarantees that every image in the hard setℐC​o​n​s​t\\mathcal\{I\}^\{Const\}is correctly classified by the ensemble whenever the corresponding LP constraint is satisfied a formal classification guarantee no penalty\-based or projection method can provide\. CG\-robust and CG\-mis surpass the baseline despite allocating column budget to adversarial and correction constraints\.

#### 4\.4Comparison with Constrained\-Learning Baselines

Competing methods \(3\-seed averages\):Penalty \(L2\): CE\+10​‖viol‖2\+~10\\\|\\text\{viol\}\\\|^\{2\};Aug\. Lagrangian: dual\-ascent on per\-image multipliers;DC3\([Donti et al\., 2021](https://arxiv.org/html/2608.26375#bib.bib10)\):K=5K=5unrolled gradient\-correction steps;Projection: inference\-time L2\-projection using true test labels \(oracle — shown for reference only\)\. Same MLP for all\.

Table 3:Test accuracy \(%\) atn=100n=100andn=1000n=1000, averaged over 3 seeds \(baselines\) and single C\+\+ run \(CG4AI\)\. Same MLP for all methods \(1 hidden layer, 4 neurons\)\. Best result excluding the oracle Projection isbolded\.Table 4:Training constraint satisfaction \(%\) atn=100n=100andn=1000n=1000, 3\-seed averages for baselines\. CG\-proba provides a formal LP guarantee; all others are heuristic\.Figure 7:Test accuracy \(left\) and training constraint satisfaction \(right\)\. DC3 and the baseline are nearly superimposed, confirming that constraint\-correction without an LP guarantee provides negligible accuracy gain\. Both CG4AI variants separate clearly from all competitors; the\+14\.6\+14\.6pp gap over DC3 atn=1000n=1000is the measurable effect of the LP output constraints\. Penalty and Aug\. Lagrangian satisfy fewer than17%17\\%of constraints at any scale\.DC3 gains only

\+0\.62\+0\.62pp \(

n=100n=100\) and

\+0\.75\+0\.75pp \(

n=1000n=1000\) over the baseline, our proxy for ensembling without LP output constraints \(cf\. LPBoost\)\. CG\-margin exceeds DC3 by

\+13\.0\+13\.0and

\+14\.6\+14\.6pp: this gap is the direct effect of the LP constraint rows \([8](https://arxiv.org/html/2608.26375#S3.E8)\) and their dual\-guided pricing\. Penalty methods collapse to

≤18\.8%\\leq 18\.8\\%at

n=100n=100\(below the baseline of

56\.1%56\.1\\%\), confirming that soft constraints are unreliable with limited model capacity\. Both formulations provide a formal classification guarantee when the LP constraint is satisfied, each under its consistent prediction rule\. CG\-proba reaches full LP feasibility at every training size \(100% of constraints, averaged\-probability prediction\)\. CG\-margin achieves 69\.7% constraint satisfaction at

n=1000n=1000\(400\-column limit reached\); the remaining 30\.3% are not formally guaranteed under the averaged\-logit prediction rule\. CG\-proba is recommended when complete coverage is required; CG\-margin when maximum test accuracy is the priority\.

### 5Main Use Case: Network Routing Under Capacity Constraints

Section[4](https://arxiv.org/html/2608.26375#S4)established CG4AI’s mechanisms on a simple classification task: a constraint forces correct classification of a designated input, the LP dual measures how much the constraint is violated, and the pricing trains the next column to address that violation\. The same mechanism now operates on a realistic combinatorial optimization problem where constraints arise not from user designations but from physics: link capacity limits in a telecommunication network\.

The key difference from Section[4](https://arxiv.org/html/2608.26375#S4)is scale and structure\. The constraint setℐC​o​n​s​t\\mathcal\{I\}^\{Const\}is continuous \(all possible demand vectors\), so the cutting\-plane procedure \(Algorithm[2](https://arxiv.org/html/2608.26375#alg2)\) is essential: at each outer iteration, a black\-box optimizer finds the demand vectorb∗b^\{\*\}that maximally leading to violations in capacity constraints, adds it toℐC​o​n​s​t\\mathcal\{I\}^\{Const\}, and triggers another column\-generation round\. The dual variablesβba\\beta\_\{b\}^\{a\}quantify how much arcaais overloaded under a demand scenario, directing each new network column toward routing decisions that respect capacity\.

#### 5\.1Problem Definition

We apply CG4AI to the Multi\-Commodity Flow Problem \(MCFP\), a natural testbed for hard constraints because the feasibility requirements \(capacity limits\) are well\-defined and computationally non\-trivial to satisfy with a learned model\.

We model the communication network as a directed graphG=\(V,A\)G=\(V,A\), whereVVis the set of nodes \(routers\) andAAis the set of arcs \(links\)\. Each link\(i​j\)∈A\(ij\)\\in Ahas a limited bandwidth capacityci​j∈ℝ\+c\_\{ij\}\\in\\mathbb\{R\}^\{\+\}\. We are given a set of demandsDD\. Each demandd∈Dd\\in Dis characterized by a source, a destination, and a traffic volumebdb^\{d\}\. Let𝒫d\\mathcal\{P\}\_\{d\}be a set of paths connecting the source of to the destination ofdd\. The goal is to route all demands while minimizing the*maximum link utilization*\(MLU\):

min\\displaystyle\\min\\quadM​L​U\\displaystyle MLU\(14\)s\.t\.∑p∈𝒫dfpd=1,\\displaystyle\\sum\_\{p\\in\\mathcal\{P\}\_\{d\}\}f\_\{p\}^\{d\}=1,∀d∈D,\\displaystyle\\forall d\\in D,\(15\)∑d∈D∑p∈𝒫d\(i,j\)∈pbd​fpd≤ci​j​M​L​U,\\displaystyle\\sum\_\{d\\in D\}\\;\\sum\_\{\\begin\{subarray\}\{c\}p\\in\\mathcal\{P\}\_\{d\}\\\\ \(i,j\)\\in p\\end\{subarray\}\}b^\{d\}\\,f\_\{p\}^\{d\}\\;\\leq\\;c\_\{ij\}\\,MLU,∀\(i,j\)∈A,\\displaystyle\\forall\(i,j\)\\in A,\(16\)fpd≥0,\\displaystyle f\_\{p\}^\{d\}\\geq 0,∀p∈𝒫d,∀d∈D,\\displaystyle\\forall p\\in\\mathcal\{P\}\_\{d\},\\;\\forall d\\in D,\(17\)whereM​L​UMLUis the maximum link utilization andfpdf\_\{p\}^\{d\}is the fraction of demandddrouted on pathp∈𝒫dp\\in\\mathcal\{P\}\_\{d\}\. The first constraint ensures that each demand is fully routed, the second enforces capacity relative toM​L​UMLU, and the third ensures non\-negativity\. MinimizingM​L​UMLUis equivalent to minimizing the maximum link utilization\.

#### 5\.2Learning Routing Splits for MCFP

Instead of solving \([14](https://arxiv.org/html/2608.26375#S5.E14)\)\-\([17](https://arxiv.org/html/2608.26375#S5.E17)\) for each new demand vector, one may train AI models to predict the routing splitsfpdf\_\{p\}^\{d\}directly\. For a given demand vectorb=\(bd\)d∈Db=\(b^\{d\}\)\_\{d\\in D\}, the model outputs predicted splitsf^pd≥0\\hat\{f\}\_\{p\}^\{d\}\\geq 0with∑pf^pd=1\\sum\_\{p\}\\hat\{f\}\_\{p\}^\{d\}=1for each demanddd\(the demand\-satisfaction constraint is enforced by a softmax output layer\)\. The model is trained on a dataset of\(b,f⋆\)\(b,f^\{\\star\}\)pairs, wheref⋆f^\{\\star\}is the optimal routing obtained by solving \([14](https://arxiv.org/html/2608.26375#S5.E14)\)\-\([17](https://arxiv.org/html/2608.26375#S5.E17)\)\. The loss function is:

ℒ⁡\(w\)=∑d∈D∑p∈𝒫d\(f^pd​\(b,w\)−fp⋆d​\(b\)\)2\.\\mathcal\{L\}\(w\)=\\sum\_\{d\\in D\}\\sum\_\{p\\in\\mathcal\{P\}\_\{d\}\}\\bigl\(\\hat\{f\}\_\{p\}^\{d\}\(b;w\)\-f\_\{p\}^\{\\star d\}\(b\)\\bigr\)^\{2\}\.\(18\)Once trained, the model can produce routing splits for new demand vectors in milliseconds, without solving an LP\. However, pure neural network predictions may violate capacity constraints, especially for high\-utilization scenarios \(MLU close to 1\)\.

Figure[8](https://arxiv.org/html/2608.26375#S5.F8)illustrates this issue: for instances where0\.99<MLU<10\.99<\\text\{MLU\}<1, the neural network often produces splits that exceed capacity\.

![[Uncaptioned image]](https://arxiv.org/html/2608.26375v1/small_violations.png)

Figure 8:Constraint violations near the MLU boundary: for high\-utilization scenarios \(0\.99<MLU<10\.99<\\text\{MLU\}<1\), a plain neural network predictor regularly exceeds the link capacity limit\.

#### 5\.3CG4AI Formulation for MCFP

In contrast with Section[5\.2](https://arxiv.org/html/2608.26375#S5.SS2), in the following we use the CG4AI framework, introduced in Section[3](https://arxiv.org/html/2608.26375#S3), for MCF problem\. The columnssis\_\{i\}is an AI model predicting routing splits\. The cutting\-plane procedure \(Algorithm[2](https://arxiv.org/html/2608.26375#alg2)\) extends the guarantee to continuous demand sets\. Ford∈Dd\\in D, letℳd\\mathcal\{M\}^\{d\}be the set of prediction models associated to demanddd\. Each model ofℳd\\mathcal\{M\}^\{d\}has\|D\|\|D\|inputs corresponding to the bandwidth demands and𝒫d\\mathcal\{P\}^\{d\}corresponding to the split on each path of the demanddd\.

The master problem for the MCFP is as follows:

min\\displaystyle\\min\\quad∑d∈D∑i∈ℳdℒid​λid\\displaystyle\\sum\_\{d\\in D\}\\sum\_\{i\\in\\mathcal\{M\}^\{d\}\}\\mathcal\{L\}\_\{i\}^\{d\}\\,\\lambda\_\{i\}^\{d\}\(19\)s\.t\.∑i∈ℳdλid=1,\\displaystyle\\sum\_\{i\\in\\mathcal\{M\}^\{d\}\}\\lambda\_\{i\}^\{d\}=1,∀d∈D,\\displaystyle\\forall d\\in D,\(20\)−∑d∈Dbd∑i∈ℳdλid∑p∈𝒫da∈psip\(b\)≥−ca,\\displaystyle\-\\sum\_\{d\\in D\}b^\{d\}\\sum\_\{i\\in\\mathcal\{M\}^\{d\}\}\\lambda\_\{i\}^\{d\}\\sum\_\{\\begin\{subarray\}\{c\}p\\in\\mathcal\{P\}^\{d\}\\\\ a\\in p\\end\{subarray\}\}s\_\{i\}^\{p\}\(b\)\\,\\;\\geq\\;\-c\_\{a\},∀a∈A,∀b∈ℐC​o​n​s​t,\\displaystyle\\forall a\\in A,\\;\\forall b\\in\\mathcal\{I\}^\{Const\},\(21\)λid≥0,\\displaystyle\\lambda\_\{i\}^\{d\}\\geq 0,∀d∈D,∀i∈ℳd\.\\displaystyle\\forall d\\in D,\\;\\forall i\\in\\mathcal\{M\}^\{d\}\.\(22\)Constraint \([20](https://arxiv.org/html/2608.26375#S5.E20)\) enforces a valid convex combination of models inℳd\\mathcal\{M\}^\{d\}per demand\. Constraint \([21](https://arxiv.org/html/2608.26375#S5.E21)\) enforces link capacity for all demand vectorsb∈ℐC​o​n​s​tb\\in\\mathcal\{I\}^\{Const\}\. One slack variable per demand is introduced at equation \([20](https://arxiv.org/html/2608.26375#S5.E20)\) penalized in the objective function\. It allows to keep the RMP feasible during the column generation algorithm\.

##### Pricing\.

For demandddand modeli∈ℳdi\\in\\mathcal\{M\}^\{d\}, the reduced cost is:

rid​\(wid\)=ℒid\+∑b∈ℐC​o​n​s​t∑a∈A∑p∈𝒫da∈pβba​bd​sip​\(b\)−γd,r\_\{i\}^\{d\}\(w\_\{i\}^\{d\}\)\\;=\\;\\mathcal\{L\}\_\{i\}^\{d\}\+\\sum\_\{b\\in\\mathcal\{I\}^\{Const\}\}\\sum\_\{a\\in A\}\\sum\_\{\\begin\{subarray\}\{c\}p\\in\\mathcal\{P\}^\{d\}\\\\ a\\in p\\end\{subarray\}\}\\beta\_\{b\}^\{a\}\\,b^\{d\}\\,s\_\{i\}^\{p\}\(b\)\-\\gamma^\{d\},\(23\)whereβba≥0\\beta\_\{b\}^\{a\}\\geq 0is the dual of the capacity constraint for arcaaand demand vectorbb, andγd\\gamma^\{d\}is the dual of the convex combination constraint for demanddd\. The pricing problem minimizesridr\_\{i\}^\{d\}over the network weightswidw\_\{i\}^\{d\}\.

##### Cutting\.

The setℐC​o​n​s​t\\mathcal\{I\}^\{Const\}may be infinite in size\. The separation problem searches for a demand vectorb∈ℐC​o​n​s​t∖ℐR​M​Pb\\in\\mathcal\{I\}^\{Const\}\\setminus\\mathcal\{I\}^\{RMP\}that maximizes the utilization violation on any arca∈Aa\\in A:

max∑d∈Dbd∑i∈ℳd∑p∈𝒫da∈psip\(b\)λid−ca\.\\max\\;\\sum\_\{d\\in D\}b^\{d\}\\sum\_\{i\\in\\mathcal\{M\}^\{d\}\}\\sum\_\{\\begin\{subarray\}\{c\}p\\in\\mathcal\{P\}^\{d\}\\\\ a\\in p\\end\{subarray\}\}s\_\{i\}^\{p\}\(b\)\\,\\lambda\_\{i\}^\{d\}\\;\-\\;c\_\{a\}\.\(24\)This is solved by a black\-box optimizer \(NLopt\), with one thread per arc to allow parallelization\. The most violating demand vectorb∗b^\{\*\}and arca∗a^\{\*\}are added inℐR​M​P\\mathcal\{I\}^\{RMP\}as a new capacity constraint\.

#### 5\.4Numerical Results

We evaluate CG4AI on the MCFP using benchmark instances from SNDLIB\([Orlowski et al\., 2010](https://arxiv.org/html/2608.26375#bib.bib23)\), a standard library for network design problems\. SNDLIB provides both real\-world and synthetic telecommunication topologies, making it a widely used benchmark for routing algorithms\.

##### Experimental setup\.

The number of demands per instance is chosen from\{5,10,15\}\\\{5,10,15\\\}\. Candidate paths per demand are computed using Yen’s algorithm\([Yen, 1970](https://arxiv.org/html/2608.26375#bib.bib24)\)to obtain theKK\-shortest paths; each demand has between 3 and 5 candidate paths\. The training set consists of 1000 randomly generated demand vectors with LP\-optimal routing splits \(solved by CPLEX\)\. Models are validated on 10 000 random demand vectors\.

###### Architecture\.

All column models are feedforward networks trained in PyTorch\. A hyperparameter search over hidden layer sizes\{64,128,256\}\\\{64,128,256\\\}and depths\{1,2\}\\\{1,2\\\}was conducted; the best configuration uses a single hidden layer with 256 neurons and GELU activation\. Training uses Adam\([Kingma and Ba, 2015](https://arxiv.org/html/2608.26375#bib.bib25)\)with learning rate10−310^\{\-3\}, early stopping \(patience 50\), and a learning rate decay of0\.90\.9every 10 epochs without improvement \(up to 1000 epochs\)\. The LP is solved by CPLEX on an Intel Xeon Platinum 8164; the separation problem is parallelized across arcs \(one thread per arc\)\.

##### Results\.

Table[6](https://arxiv.org/html/2608.26375#S7.T6)\(Appendix\) reports detailed results for all instances\. Here we highlight the main findings\.

##### Feasibility\.

CG4AI successfully produces feasible ensembles \(zero slack\) for most instances, as shown by the slack convergence plots in Figure[9](https://arxiv.org/html/2608.26375#S5.F9)\. The slack variable converges rapidly to zero, with occasional small increases when a new constraint is added, which is quickly resolved by the next column generation step\.

![Refer to caption](https://arxiv.org/html/2608.26375v1/slack_plot.png)\(a\)Slack evolution for nobel\-us\(\#​d,\#​p\)=\(5,3\)\(\\\#d,\\\#p\)=\(5,3\)\.
![Refer to caption](https://arxiv.org/html/2608.26375v1/slack_plot_from10.png)\(b\)Zoomed view from iteration 10\.

Figure 9:Slack variable evolution during the column\-generation and cutting\-plane iterations\. The slack converges rapidly to zero, confirming that the algorithm achieves feasibility\.
##### Prediction accuracy\.

Figure[10](https://arxiv.org/html/2608.26375#S5.F10)shows the predicted versus optimal MLU for the nobel\-germany instance with 15 demands and 3 paths\. The predicted values align well with the optimal values\. Figure[11](https://arxiv.org/html/2608.26375#S5.F11)shows the newyork instance \(5 demands, 3 paths\), where predictions are less accurate for high\-utilization scenarios; the neural ensemble tends to underestimate MLU in the most congested cases\. Importantly, all predictions remain below the critical value of 1\.0 for both instances, confirming constraint satisfaction\.

![Refer to caption](https://arxiv.org/html/2608.26375v1/mlu_vs_true_comparison-nobel-germany3.png)Figure 10:Predicted vs\. optimal MLU on the nobel\-germany instance\(\#​d,\#​p\)=\(15,3\)\(\\\#d,\\\#p\)=\(15,3\)\. The ensemble closely tracks the optimal solution\.![Refer to caption](https://arxiv.org/html/2608.26375v1/mlu_vs_true_comparison-newyork1.png)Figure 11:Predicted vs\. optimal MLU on the newyork instance\(\#​d,\#​p\)=\(5,3\)\(\\\#d,\\\#p\)=\(5,3\)\. Greater dispersion occurs for high\-utilization scenarios, but predictions remain feasible\.
##### Constraint analysis\.

Figure[12](https://arxiv.org/html/2608.26375#S5.F12)shows the violation history for Arc 13 in the Polska instance during optimization process\. Initial violations exceed the capacity by 25%\. The violations decrease over iterations, though not monotonically, reflecting the interplay between adding new columns and new constraints\.

![Refer to caption](https://arxiv.org/html/2608.26375v1/arc13_violations.png)\(a\)Utilization of Arc 13 over iterations\.
![Refer to caption](https://arxiv.org/html/2608.26375v1/arc13_demands.png)\(b\)Demand vectors generated for Arc 13\.

Figure 12:Arc violations for Arc 13 in the Polska instance\(\#​d,\#​p\)=\(5,3\)\(\\\#d,\\\#p\)=\(5,3\)\.The analysis of the demand vectors generated for Arc 13 reveals an informative pattern: three demands are consistently set to their maximum valuemm, while the remaining two are varied across\[0,m\]\[0,m\]\. Inspection of the network topology shows that those three demands have no path avoiding Arc 13, making them*captive*to this bottleneck\. The optimizer varies the non\-captive demands to find combinations where the ensemble’s routing fails to exploit alternative paths\. This insight could be used to seed the initial constraint setℐ0C​o​n​s​t\\mathcal\{I\}^\{Const\}\_\{0\}with adversarial scenarios \(captive demands atmm, others sampled strategically\), potentially reducing the number of cutting\-plane iterations\.

##### Instance\-level analysis\.

Across the SNDLIB benchmark \(Table[6](https://arxiv.org/html/2608.26375#S7.T6)\), the results confirm several consistent patterns\. The algorithm produces feasible solutions for most instances\. Harder instances \(more demands, denser topologies\) require more iterations and generate more columns; for example, di\-yuan with 10 demands requires 999 iterations and 986 columns\. Simpler instances converge in just 1–3 iterations \(ta1 with 5 demands, france with 10 demands\)\. The nobel family \(eu, germany, us\) consistently achieves very precise predictions, with validation losses on the order of10−410^\{\-4\}\. The final active model count is typically much smaller than the total number of generated columns, confirming that the LP selects a compact ensemble\. For three instances \(geant with 15 demands, di\-yuan with 15 demands and 10 demands\), the iteration limit is reached without full convergence, suggesting that these configurations require a larger column budget\.

##### Ablation: alternative objective\.

To verify that constraint satisfaction in CG4AI is driven by the dual\-guided pricing and not by properties of the MLU objective itself, we replace the MLU with an alternative cost function: the total weighted path cost∑d∈D∑p∈𝒫dwp⋅fpd⋅bd\\sum\_\{d\\in D\}\\sum\_\{p\\in\\mathcal\{P\}^\{d\}\}w\_\{p\}\\cdot f\_\{p\}^\{d\}\\cdot b^\{d\}, wherewpw\_\{p\}is the sum of arc weights along pathpp\. This function has no inherent regularizing effect on capacity constraints\.

Figure[13](https://arxiv.org/html/2608.26375#S5.F13)shows the results on the ta1 instance with 10 demands and 2 paths\. CG4AI still produces a feasible ensemble that respects all capacity constraints, although more iterations are needed compared to the MLU case\. The optimality gap on the cost function is larger \(11\.20% at the final iteration\), reflecting the harder learning problem when the objective is less aligned with the constraints\. This confirms that the dual\-guided pricing is the primary mechanism for enforcing constraint satisfaction, independent of the specific training objective\.

![Refer to caption](https://arxiv.org/html/2608.26375v1/cost.png)\(a\)Predicted vs\. optimal total path cost\.
![Refer to caption](https://arxiv.org/html/2608.26375v1/mlu_cost.png)\(b\)MLU while optimizing the total path cost\.

Figure 13:Constrained optimization with the alternative cost function on the ta1 instance\. CG4AI satisfies all capacity constraints even though the training objective is not aligned with them\.

### 6Conclusion

We have presented CG4AI, a column generation framework for training convex ensembles of AI models that satisfy hard linear constraints\. The framework combines three key ingredients: a master LP that enforces feasibility through the mixture weights; a pricing subproblem that trains new models guided by LP dual variables to focus on the most violated constraints; and a cutting\-plane procedure that extends feasibility guarantees beyond a finite training set\.

We demonstrated CG4AI on two applications\. On MNIST, we showed that hard constraints can be used as a replacement for labeled data \(constraint\-only learning\), as a mechanism for improving adversarial robustness, as a tool for correcting misclassified examples, and as a way to implement output relabeling without full retraining\. On the multi\-commodity flow problem, CG4AI produces routing predictors that satisfy all link capacity constraints on standard SNDLIB benchmark networks, with good MLU accuracy and fast inference at test time\.

The results highlight several interesting properties of the framework\. Themarginconstraint formulation is particularly efficient, requiring only 3–6 columns even for 1000 training images, with zero residual violation\. The dual\-guided pricing is effective regardless of the training objective, as shown by the ablation with an alternative cost function\. The cutting\-plane separation provides a principled mechanism for finding adversarial demand vectors that expose model weaknesses\.

Several directions remain open for future work\. First, enforcing an integer ensemble by restrictingλi∈\{0,1\}\\lambda\_\{i\}\\in\\\{0,1\\\}via a branch\-and\-bound scheme would yield a single, highly interpretable model\. Second, the model\-agnostic nature of the master LP framework permits the integration of alternative column classes, such as decision trees or gradient\-boosted models, which are well\-suited for structured constraint patterns\. Third, computational efficiency could be significantly enhanced by replacing the general black\-box separator with problem\-specific cut generation heuristics, or by employing an exact ILP\-based separation formulation for single\-layer ReLU networks to guarantee provably optimal cuts\. Finally, the framework can be extended dynamically through an active learning procedure that updates training losses using newly identified violating inputs, or scaled to distributed AI agents to enable decentralized deployment while maintaining global feasibility guarantees\.

## References

- Min and Azizan \[2024\]Min, Y\., Azizan, N\.: HardNet: Hard\-constrained neural networks with universal approximation guarantees\. arXiv preprint arXiv:2410\.10807 \(2024\)
- Tordesillas et al\. \[2023\]Tordesillas, J\., How, J\.P\., Hutter, M\.: RAYEN: Imposition of hard convex constraints on neural networks\. In: arXiv Preprint arXiv:2307\.08336 \(2023\)
- Amos and Kolter \[2017\]Amos, B\., Kolter, J\.Z\.: OptNet: Differentiable optimization as a layer in neural networks\. In: Proceedings of the 34th International Conference on Machine Learning \(ICML\)\. PMLR, vol\. 70, pp\. 136–145 \(2017\)
- Demiriz et al\. \[2002\]Demiriz, A\., Bennett, K\.P\., Shawe\-Taylor, J\.: Linear programming boosting via column generation\. Machine Learning46\(1–3\), 225–254 \(2002\)[https://doi\.org/10\.1023/A:1012470815092](https://doi.org/10.1023/A:1012470815092)
- Goodfellow et al\. \[2016\]Goodfellow, I\., Bengio, Y\., Courville, A\.: Deep Learning\. MIT Press, Cambridge, MA \(2016\)
- Chamon and Ribeiro \[2020\]Chamon, L\.F\.O\., Ribeiro, A\.: Probably approximately correct constrained learning\. In: Advances in Neural Information Processing Systems \(NeurIPS\), vol\. 33, pp\. 17–17 \(2020\)\. arXiv:2006\.05487
- Elenter et al\. \[2024\]Elenter, J\., Chamon, L\.F\.O\., Ribeiro, A\.: Near\-optimal solutions of constrained learning problems\. In: International Conference on Learning Representations \(ICLR\) \(2024\)\. arXiv:2403\.09704
- Lu et al\. \[2021\]Lu, L\., Pestourie, R\., Yao, W\., Wang, Z\., Verdugo, F\., Johnson, S\.G\.: Physics\-informed neural networks with hard constraints for inverse design\. SIAM Journal on Scientific Computing43\(6\), 1105–1132 \(2021\)[https://doi\.org/10\.1137/21M1397908](https://doi.org/10.1137/21M1397908)
- Agrawal et al\. \[2019\]Agrawal, A\., Amos, B\., Barratt, S\., Boyd, S\., Diamond, S\., Kolter, J\.Z\.: Differentiable convex optimization layers\. In: Advances in Neural Information Processing Systems \(NeurIPS\), vol\. 32 \(2019\)
- Donti et al\. \[2021\]Donti, P\.L\., Rolnick, D\., Kolter, J\.Z\.: DC3: A learning method for optimization with hard constraints\. In: International Conference on Learning Representations \(ICLR\) \(2021\)\. arXiv:2104\.12225
- Balestriero and LeCun \[2023\]Balestriero, R\., LeCun, Y\.: POLICE: Provably optimal linear constraint enforcement for deep neural networks\. In: ICASSP 2023 – IEEE International Conference on Acoustics, Speech and Signal Processing, pp\. 1–5 \(2023\)\.[https://doi\.org/10\.1109/ICASSP49357\.2023\.10096520](https://doi.org/10.1109/ICASSP49357.2023.10096520)
- Constante Flores et al\. \[2025\]Constante Flores, G\., et al\.: Enforcing hard linear constraints in deep learning models with decision rules\. arXiv preprint arXiv:2505\.13858 \(2025\)
- Aftabi et al\. \[2025\]Aftabi, N\., Moradi, N\., Mahroo, F\.: Feed\-forward neural networks as a mixed\-integer program\. Engineering with Computers41, 2411–2429 \(2025\)[https://doi\.org/10\.1007/s00366\-025\-02114\-2](https://doi.org/10.1007/s00366-025-02114-2)\. arXiv:2402\.06697
- Tjeng et al\. \[2019\]Tjeng, V\., Xiao, K\.Y\., Tedrake, R\.: Evaluating robustness of neural networks with mixed integer programming\. In: International Conference on Learning Representations \(ICLR\) \(2019\)\. arXiv:1711\.07356
- Anderson et al\. \[2020\]Anderson, R\., Huchette, J\., Ma, W\., Tjandraatmadja, C\., Vielma, J\.P\.: Strong mixed\-integer programming formulations for trained neural networks\. Mathematical Programming183, 3–39 \(2020\)[https://doi\.org/10\.1007/s10107\-020\-01474\-5](https://doi.org/10.1007/s10107-020-01474-5)\. Extended abstract at IPCO 2020\. arXiv:1811\.08359
- Aziz et al\. \[2024\]Aziz, V\., Wu, O\., Nowak, I\., Hendrix, E\.M\.T\., Kronqvist, J\.: On optimizing ensemble models using column generation\. Journal of Optimization Theory and Applications203\(2\), 1794–1819 \(2024\)[https://doi\.org/10\.1007/s10957\-024\-02391\-9](https://doi.org/10.1007/s10957-024-02391-9)
- Jensen \[1906\]Jensen, J\.L\.W\.V\.: Sur les fonctions convexes et les inégalités entre les valeurs moyennes\. Acta mathematica30\(1\), 175–193 \(1906\)
- Barnhart et al\. \[1998\]Barnhart, C\., Johnson, E\.L\., Nemhauser, G\.L\., Savelsbergh, M\.W\.P\., Vance, P\.H\.: Branch\-and\-price: Column generation for solving huge integer programs\. Operations Research46\(3\), 316–329 \(1998\)[https://doi\.org/10\.1287/opre\.46\.3\.316](https://doi.org/10.1287/opre.46.3.316)
- Kelley \[1960\]Kelley, J\.E\.: The cutting\-plane method for solving convex programs\. Journal of the Society for Industrial and Applied Mathematics8\(4\), 703–712 \(1960\)
- Johnson et al\. \[2014\]Johnson, S\.G\., et al\.: The NLopt nonlinear\-optimization package \(2014\)
- LeCun et al\. \[1998\]LeCun, Y\., Bottou, L\., Bengio, Y\., Haffner, P\.: Gradient\-based learning applied to document recognition\. Proceedings of the IEEE86\(11\), 2278–2324 \(1998\)
- IBM Corporation \[2015\]IBM Corporation: IBM ILOG CPLEX 12\.6\.3 User’s Manual\. \(2015\)\. IBM Corporation\. Armonk, NY
- Orlowski et al\. \[2010\]Orlowski, S\., Wessäly, R\., Pióro, M\., Tomaszewski, A\.: SNDlib 1\.0 – Survivable Network Design Library\.[http://sndlib\.zib\.de](http://sndlib.zib.de/)\(2010\)
- Yen \[1970\]Yen, J\.Y\.: An algorithm for finding shortest routes from all source nodes to a given destination in general networks\. Quarterly of Applied Mathematics27\(4\), 526–530 \(1970\)
- Kingma and Ba \[2015\]Kingma, D\.P\., Ba, J\.: Adam: A method for stochastic optimization\. In: International Conference on Learning Representations \(ICLR\) \(2015\)\. arXiv:1412\.6980

### 7Experimental Tables

Table 5:Column definitions for Table[6](https://arxiv.org/html/2608.26375#S7.T6)\.Table 6:Summary of CG4AI results on SNDLIB instances\.instance\#d\#p\#iter\#cols\#rowscol\_timerow\_timetot\_timeopt\_mlu\#modval\. lossmax\_diffabilene53440284\.83875\.471160\.417\.63e\-0220\.0210\.38103573692207\.117831\.8710040\.361\.57e\-01340\.0130\.46153282101336\.921215\.772553\.484\.26e\-0230\.0150\.42551184463\.051190\.411653\.699\.14e\-02100\.0270\.42atlanta534526102229\.1013427\.9015658\.016\.07e\-01430\.0100\.331036342162832\.3019045\.3021879\.009\.91e\-02530\.0170\.371531275956319\.7810198\.4016520\.902\.95e\-02260\.0270\.5755550388\.271279\.121667\.535\.16e\-0130\.0040\.34dfn\-bwin53282403582\.93774\.634358\.141\.01e\-0110\.1010\.651033531564225530\.8095600\.80121141\.192\.90e\-011040\.0370\.49153406138830845\.7022207\.0053064\.822\.73e\-01390\.0360\.45di\-yuan53155147812730\.7012305\.7025039\.903\.59e\-03330\.0030\.161039999865971904\.8052162\.70124094\.289\.15e\-031430\.0090\.23153999845068495\.300\.0068516\.256\.09e\+00580\.0580\.5055350344921822\.8014763\.5036594\.964\.45e\-03440\.0040\.17france103330273\.24241\.05514\.362\.28e\-0210\.0240\.491539701300\.781288\.172589\.152\.88e\-0210\.0290\.5255220384\.32940\.691325\.075\.97e\-021NaNNaNgeant53220393\.921058\.571452\.539\.98e\-0210\.0990\.66103161401682\.604160\.635843\.641\.68e\-0170\.0400\.481539992981355938\.1013270\.1069230\.614\.00e\+02870\.0420\.4555231741683\.582526\.604210\.791\.08e\-01160\.1020\.66newyork53121101969\.16993\.962963\.482\.39e\-0210\.0250\.491036029224803\.4649169\.5053974\.741\.12e\-01710\.0170\.34153201422502\.165126\.397629\.035\.20e\-0290\.0490\.5655110288\.411080\.681369\.115\.90e\-0210\.0560\.61nobel\-eu538868682456205\.7014245\.1070476\.87\-5\.37e\-02660\.0010\.131036345175047\.6131214\.6036263\.941\.98e\-04572\.3e\-40\.101532091802915888\.8045788\.2061682\.893\.64e\-04776\.0e\-40\.13nobel\-germany53180131947834\.6245883\.7053722\.763\.40e\-041893\.2e\-40\.09153383242753\.765227\.277982\.082\.47e\-04174\.1e\-40\.115514296795303\.8137377\.3042684\.531\.87e\-041641\.3e\-40\.07nobel\-us53463482605\.427542\.8910149\.665\.81e\-05307\.1e\-50\.05103258167758832\.3485008\.3093847\.174\.18e\-042346\.0e\-40\.14153189167478804\.4333584\.5042393\.731\.01e\-031000\.0010\.155519153974\.582233\.843208\.821\.03e\-04131\.1e\-40\.06polska532111151285319\.4755888\.3061212\.434\.91e\-032820\.0020\.22ta153110386\.521291\.101677\.642\.01e\-0210\.0200\.4910314265589850\.20151503\.00161356\.947\.31e\-021750\.0180\.42153722709054\.514745\.6213802\.147\.07e\-0260\.0210\.4255440730\.991773\.092504\.204\.75e\-0220\.0470\.66##### Note on anomalousopt\_mluvalues\.

geant \(15d\):opt\_mlu=4\.00×102=4\.00\\times 10^\{2\}— the LP becomes unbounded for this demand/capacity configuration \(demand exceeds capacity\); CG4AI still generates a feasible ensemble within the path set\.nobel\-eu \(5d\):opt\_mlu=−5\.37×10−2=\-5\.37\\times 10^\{\-2\}— floating\-point artefact of the CPLEX dual at near\-zero MLU; functionally equivalent to00\.

The results across SNDLIB instances confirm the following patterns\. The algorithm produces feasible solutions for the large majority of tested configurations\. Computational effort scales with problem difficulty: di\-yuan and polska are the hardest instances, requiring hundreds to nearly a thousand iterations, while france and ta1 converge in just a few\. The nobel family achieves the highest accuracy, with validation losses on the order of10−410^\{\-4\}and maximum differences below 0\.15 MLU units\. The final active model count \(column \#mod\) is consistently smaller than the total number of generated columns, indicating that the LP selects a compact ensemble that covers the constraint requirements without retaining all generated columns\.