Scaling Decision-Focused Learning to Large Problems with Lagrangian Decomposition

1 Introduction

Many real-world problems can be formulated as constrained optimization problems (COPs), where the goal is to identify an optimal decision that maximizes or minimizes an objective function subject to a set of constraints. Both the objective function and the feasible region are typically parametric, depending on decision variables as well as contextual parameters. Although decades of research have produced highly efficient solvers for scenarios in which these contextual parameters are known, such information is often unavailable in practical applications. Consequently, inmost realistic settings, the contextual parameters must be predicted from external observations, for example, using a machine learning model, before being used by an optimization solver to determine the downstream optimal decision. This framework is commonly referred to as the predict-then-optimize paradigm. In such situations, suboptimal downstream decisions generally stem not from the limitations of the COP solver, but from inaccuracies in the predicted contextual parameters. Designing a predictive model that leads to near-optimal downstream decisions is therefore a central challenge. Traditional machine learning approaches, which optimize for predictive accuracy alone, often fall short because they ignore the combinatorial structure and decision sensitivity inherent to the underlying optimization problem.

To overcome this limitation, the decision-focused learning (DFL) paradigm has emerged. Rather than training the predictive model to minimize a standard prediction loss, this approach explicitly incorporates the quality of the downstream decisions into the learning process, thereby aligning prediction and optimization objectives. Although Elmachtoub et al. [Elmachtoub et al., 2023] demonstrated that predictionfocused learning (i.e. treating estimation and decision making as independent steps) performs best when the underlying predictive model is well-specified, it is outperformed by DFL approaches in the presence of epistemic or aleatoric uncertainty (i.e. when the model is under-specified or the data are noisy). Early work on DFL proposed direct differentiation through the argmin operator, initially limited to unconstrained relaxations [Gould et al., 2016], and later extended to more general settings via differentiation of the Karush-Kuhn- Tucker conditions [Amos and Kolter, 2017]. Subsequent developments introduced several alternative formulations, including the SPO+ surrogate loss [Elmachtoub and Grigas, 2021] IMLE [Niepert et al., 2021], CaVE [Tang and Khalil, 2024a], LAVA [Berden et al., 2025], and the neuro-symbolic E-PLL [Defresne et al., 2023] for discrete graphical models, among others. For a comprehensive review and benchmark comparison, we refer the reader to the survey by Mandi et al. [Mandi et al., 2024]. While DFL effectively aligns learning with decision quality, it suffers from substantial computational overhead: at each training step and for every instance, a constrained optimization problem must be solved. Since these problems are often NP-hard, and as the number of solver calls scales with both dataset size and the number of epochs, the computational cost increases sharply with instance size. This turns the training process into a sequence of expensive combinatorial optimizations that severely limit scalability and usability to solve large-scale instances.

Separately, Lagrangian decomposition [Guignard and Kim, 1987] (LD) is a scalable optimization technique that has been successfully applied to large-scale combinatorial problems, in both mixed-integer programming and constraint programming contexts. Compared to the standard Lagrangian relaxation, LD can provide tighter dual bounds and handle arbitrary constraint structures, but at the expense of an increased computational cost: Guignard and Kim [Guignard and Kim, 1987] (see Theorem 3.1) prove that for linearly-constrained problems the LD bound is at least as tight as the Lagrangian relaxation. The key idea behind LD is to duplicate variables and partition the initial problem into independent subproblems. Ideally, these subproblems should preserve the combinatorial nature of the original formulation while being significantly faster to solve. The solutions of the subproblems can then be combined to derive a valid dual bound on the original problem. However, this bound depends on a vector of Lagrangian multipliers that coordinate the subproblems; minimizing it amounts to optimizing these multipliers, typically through subgradient methods [Shor, 2012]. Lagrangian decomposition also integrates naturally with branch-and-bound algorithms: as noted by Guignard and Kim, disagreements between variable copies directly suggest branching candidates, yielding a natural dichotomic branching rule. In the context of constraint programming, H`a et al. [H`a et al., 2015] demonstrated that LD can serve as an automatic, though computationally expensive, bounding mechanism that substantially reduces the size of the search tree. More recently, Bessa et al. [Bessa et al., 2025] showed that machine learning techniques can be leveraged to warm-start the subgradient optimization, thereby improving the efficiency of the LD process.

Inspired by this recent line of work, and by other contributions highlighting the benefits of combining learning with Lagrangian methods [Parjadis et al., 2024; Abbas and Swoboda, 2024], we propose to integrate Lagrangian decomposition into a decision-focused learning pipeline. This integration enables DFL to scale effectively to very large combinatorial optimization problems. Our approach contrasts with most existing DFL research, which primarily focuses on developing new differentiation techniques. Instead, we introduce a general methodology that can be seamlessly combined with any existing differentiation algorithm, providing a flexible and scalable framework for decision-focused learning.

Our contributions are as follows: (1) a methodology to leverage Lagrangian decomposition into DFL, yielding a new surrogate function to optimize, (2) two loss functions to train a prediction model, and (3) the integration of our approach on two common DFL techniques (SPO+ and IMLE). The experiments are conducted on instances of the multi-dimensional knapsack problem with up to 300 variables and 10 constraints, which is roughly six times larger than what is commonly considered in prior work, and on instances of the portfolio optimization problem with up to 400 variables, representing an 8-fold increase compared to typical benchmarks. To our knowledge, this work constitutes the first attempt to incorporate Lagrangian decomposition into a DFL pipeline.

2 Problem Setting

In many industrial applications, certain parameters of the constrained optimization problem are uncertain and must be inferred from contextual data, commonly referred to as features. The settings we consider involve estimating these parameters through predictive inferences made by machine learning models, after which the final decisions are obtained by solving the corresponding optimization problems based on these predictions. In this framework, the overall decisionmaking process can be formulated as a predict-then-optimize problem, where the optimization depends on parameters inferred from observed features.

Let c ∈ Rk denote a vector of k unknown parameters, X decision variables defined over domain D, C = (C1, . . . , Cm) a set of m constraints, and f : D×Rk → R an objective function. The problem we consider is as follows:

The optimal solution is denoted by X∗ . Although the optimal solution may not be unique, we only require a consistent tiebreaking rule so that c ↦→ X∗ (c) is well-defined as a function. In practice, this is satisfied since the solver returns a single deterministic solution.

In practice, the true cost parameters c are not directly observable and must be inferred from feature vectors z ∈ Rq . This inference is performed using a predictive model Ωθ : Rq → Rk which produces an estimate ˆc = Ωθ(z), where θ denotes the parameters of the model. The decision ultimately implemented is referred to as the prescriptive solution, denoted X∗ (ˆc). This solution is generally suboptimal, and its quality can be assessed through the regret, defined as:

Our objective is to learn a predictive model that produces prescriptive solutions as close as possible to the true optimal ones over a set of historical observations, or equivalently, to minimize the expected regret. To perform this prediction task, we employ neural networks as a differentiable predictive models. Training is carried out using gradient descent, which requires a loss function whose gradient can be efficiently computed or approximated. The choice of an appropriate loss is critical, as it has a major impact on the overall performance of the model.

Prediction-focused learning is arguably the most intuitive approach to address this problem. In this paradigm, estimation and decision making are treated as two fully independent steps. The predictive model Ωθ is trained to estimate c as accurately as possible by minimizing a prediction loss, typically the mean squared error between the true parameters c and their estimates ˆc. Formally, this loss is defined as follows:

However, this approach has several limitations, particularly in under-specified settings, as demonstrated by Elmachtoub et al. [Elmachtoub et al., 2023]. This is why decision-focused learning is often more appropriate, as it directly targets the quality of the prescriptive decision. In this setting, the regret, as defined in Equation (1), is used as the loss function. However, this formulation introduces new challenges. First, computing the gradient of Equation (1) is difficult, and sometimes even meaningless, as ∂X∗ ∂ˆc can potentially be null almost everywhere. This occurs, for example, when the mapping c ↦→ X∗ (c) is piecewise constant, as is common in combinatorial optimization. Fortunately, prior work has proposed ways to obtain exact gradients in special cases, or useful gradient approximations in more general scenarios, thereby enabling DFL in principle. Second, DFL suffers from substantial computational cost. Because it relies on the prescriptive decision during training, a new COP must be solved for every instance at every epoch. When the underlying problem is NP-hard, these repeated solver calls often become the primary bottleneck in the training pipeline. As a result, classical DFL methods do not scale well to large COPs, whose complexity grows exponentially. Our contribution is an automatic mechanism based on Lagrangian decomposition to tackle this challenge.

3 Lagrangian Decomposition for DFL

The key idea of Lagrangian decomposition is to provide a dual bound to a given multi-constrained optimization problem by solving mono-constrained ones. Let us start from (Pc) defined earlier. Let X1 refer to the initial variables. We introduce m - 1 new variables by duplicating X1, obtaining a vector ˜X = (X1, . . . ,Xm) where each Xi keeps the same domain as X1. We first observe that, by adding equality constraints, the original COP (Pc) is equivalent to:

The notation Ci(Xi) means that constraint Ci is applied to Xi. This does not require any structural assumption on the constraints since each Xi is a full copy of the original variable vector X (same dimension and domain). We then relax the new equality constraints by adding penalty terms into the objective function with multipliers µ = (µ2, . . . , µm) ∈

R|X1|×(m-1) . By relaxing the equality constraints, we obtain a dual bound:

We note that this bound is parametrized by the multipliers µ. Those are commonly referred to as Lagrangian multipliers. Since each constraint Ci acts on a distinct variable and the objective function is separable on the variables, the computation of B(µ, c) can be decomposed into m subproblems featuring only a single constraint:

In practice, we are interested in finding the tightest possible dual bound, which requires computing the best multipliers i.e. µ∗ (c) = argminµ (B(µ, c)). This can be done with a subgradient descent method to optimize the argmin, such as Shor's algorithm [Shor, 2012]. Let ˜X∗ (µ, c) = (X∗1 (µ, c), . . . , X∗m(µ, c)) be the solution of Equation (4) for a given µ. A standard subgradient expression of B(µ, c) with regard to µi is (X∗1 (µ, c) -X∗i (µ, c)) for all i ∈ [[2,m]] [Bessa et al., 2025].

3.1 Defining Appropriate Loss Functions

Lagrangian decomposition provides a fast way to compute an upper bound, although the bound's tightness depends on the chosen multipliers. Ideally, we would replace the COPsolving phase in the traditional DFL pipeline with Lagrangian decomposition using optimal multipliers, which yield the tightest dual bound. Taking a step back, we note that Equation (4), when parameterized and evaluated with optimal multipliers, becomes:

This expression is a new COP parameterized by (µ∗ (c), c). Therefore, we can reuse the tools from DFL, except that we need to handle the multipliers µ∗ (c) as well, which act as parameters in Equation (8) along c, but are not given by the predictive model upstream. In contrast, they are obtained by a third-party optimization algorithm. The next step is to design an appropriate DFL-based loss function. Let ˜X∗ (µ∗ , c) = (X∗1 (µ∗ , c), . . . ,X∗m(µ∗ , c)) denote the solution to Equation (8) and let b : (˜X, µ, c) ↦→ f(X1, c) +

∑mi=2 µi · (X1 -Xi) be a function mapping for Lagrangian decomposition. We propose two loss functions that integrate Lagrangian decomposition into a DFL pipeline:

The first loss (L1) corresponds exactly to the regret associated with Equation (8). It involves the objective function b(˜X, µ, c) where the Lagrangian multipliers appear explicitly. In contrast, in the second loss (L2), explicit penalty terms with the multipliers disappear since f(X, c) is used instead of b(˜X, µ, c). However, the multipliers are still involved in the resolution of X∗1 , so they are still ultimately considered. Intuitively, L2 mirrors the regret defined from (Pc), but it is evaluated at X∗1 (µ∗ ,ˆc) rather than X∗ (ˆc). We use X∗1 (instead of the duplicated variables ˜X∗ ) because these variables are obtained from Equation (5), which directly incorporates both the objective f and the multipliers. Since the bottleneck in DFL is the solving time of the COP at each iteration, we expect the two losses based on Lagrangian decomposition to be significantly faster. As a counterpart, as they are based on a relaxed problem (Eq. (8)), the accuracy of the estimator ˆc on the primal (Pc) is not guaranteed. The choice to use f(X, c) in L2 is motivated by the idea that it could strengthen the performance during evaluation on (Pc). Figure 1 illustrates how these two losses differ from the standard regret loss.

3.2 Dealing withMultipliers under Uncertainty

Unfortunately, minimizing these two loss functions is challenging, mainly because the multipliers depend on the uncertain parameters, i.e., µ∗ (ˆc). A first issue is that computing µ∗ (ˆc) for each instance would significantly increase computational cost. Moreover, differentiating the losses (i.e., ∂L ∂ˆc ) requires computing both ∂X∗1 ∂ˆc and ∂µ∗ ∂ˆc . As a result, optimizing these losses becomes a bilevel optimization problem. We address this challenge by fixing the Lagrangian multipliers rather than optimizing over µ∗ (ˆc). With fixed multipliers, the optimal solutions to the subproblems (X∗2 , . . . ,X∗m) no longer depend on ˆc (Eq. (7)). The next step is to decide how to set these multipliers. A natural choice is to fix them to their optimal values under the true cost parameters, i.e., we set µ∗ (ˆc) = µ∗ (c). Intuitively, this choice reflects the multipliers we would like to obtain at the end of training. The two corresponding losses can then be reformulated as follows:

where σ : (X1, µ, c) ↦→ f(X1, c) + ∑mi=2 µi · X1. The resulting training algorithm, referred to as trainingSingle since it relies on a single fixed Lagrangian decomposition, is depicted in Algorithm 1. It illustrates how the predictive model Ωθ is trained using a dataset S =

Algorithm 1 trainingSingle(S, Ωθ , P, L, α, E)

Require: Dataset S = {(zi, ci, X∗1 (ci), µ∗ (ci))}ni=1 .

Require: Predictive model Ωθ.

Require: Combinatorial problem P = ⟨X, D, C, f⟩.

Require: Loss function L ∈ {Lnew

2 }.

1 ,Lnew

Require: learning rate α, number of training epochs E.

• 1: θ := randomInitialization()
• 2: for each epoch from 1 to E do

3:

for each sample i ∈ S do

4:

µ := µ∗ (ci)

• 5: ˆ ci := Ωθ(zi)
• 6: Φ := LD(P) {Eq. (6)}
• 7: X∗ 1 := solve(Φ, µ, ˆ ci) {Subproblem solving}

8:

θ := θ -α( ∂L

1

• 9: end for
• 10: end for
• 11: return Ωθ

ci ,µ)

1 (ˆ

∂ˆ

ci ci

∂θ

· ∂ ˆ

) {Chain rule}

· ∂X∗

∂X∗

{(zi, ci,X∗1 (ci), µ∗ (ci)}ni=1 . Recall that zi denotes the feature vector of data sample i. We note that constructing this training set requires computing the Lagrangian multipliers, via a subgradient method, for all instances in the training set and for multiple cost parameter values, which can be computationally expensive. Instead, we propose to terminate the subgradient procedure after a predefined number of iterations and to use the resulting potentially suboptimal multipliers. As we will show in Section 5, this early termination still leads to competitive results in practice. The model weights are initialized randomly (Line 1), and training proceeds for a fixed number of epochs. Each data sample is processed once per epoch. For every sample, we first predict the cost parameters from the feature vector (Line 5), then execute the Lagrangian decomposition of the problem (Line 6), and finally solve the resulting subproblems (Line 7). As explained previously, only the variables X∗ 1 are required to compute the loss. The weights are then updated by backpropagating the loss using the chain rule (Line 8). The term ∂L ∂X∗ 1 is obtained directly from Eq. (11) or Eq. (12) and ∂ ˆ ci ∂θ is computed by differentiating through Ωθ, e.g., a neural network. Concerning the term ∂X∗ 1 (ˆ ci ,µ) ∂ˆ ci , we emphasize that the main advantage of our approach is that we do not need anymore to consider all the variables as in standard DFL. The LD subproblem is easier to solve and differentiate because it is smaller and has a simpler structure. This gradient can be obtained using classical DFL differentiation techniques or, when the subproblem has exploitable structure, via dedicated solver-free differentiation methods. These two situations are later considered in our experiments. Since our approach can be combined with any differentiation technique, it provides flexibility in selecting the most suitable one for a given use case. Note that some DFL differentiation techniques, e.g. SPO+, provide the straightforward gradient ∂L ∂ˆ ci instead of ∂X∗ 1 (ˆ ci ,µ) ∂ˆ ci , which is still compatible with Algorithm 1.

Algorithm 2 trainingMultiple()

Require: Number of decompositions m.

Require: S (d) = {(zi, ci,X∗1 (ci)(d) , µ∗ (ci)(d) )}ni=1 .

Require: Predictive model Ωθ.

Require: Combinatorial problem P = ⟨X, D, C, f⟩.

Require: Loss function L(d) ∈ {L(d)

1 ,L(d) 2 } ∀d ∈ [[1,m]].

Require: learning rate α, number of training epochs E.

• 1: θ := randomInitialization()
• 2: for each epoch from 1 to E do
• 3: d := randomSelection(1,m)
• 4: Ωθ := trainingSingle(S(d) ,Ωθ,P,L(d) , α, 1)
• 5: end for
• 6: return Ωθ

are chosen to enable direct comparison with the work of Elmachtoub and Grigas [Elmachtoub and Grigas, 2021] and Mandi et al. [Mandi et al., 2020].

Multi-dimensional knapsack. It is an extension of the classic knapsack problem. In the multidimensional knapsack problem, we are given m constraints and n items, each associated with a profit c ∈ Rn and multiple weights (w1 , . . . ,wm) ∈ Rn×m. Each constraint k ∈ {1, . . . ,m} represents a capacity limit bk for a different dimension. The objective is to select a subset of items that maximizes total profit while ensuring that the total weight in each dimension does not exceed the corresponding capacity constraint. All benchmark instances follow the protocol explained by [Tang and Khalil, 2024b] and summarized in the survey of Mandi et al. [Mandi et al., 2024]. For each data point i, shared features zi ∈ Rp are generated and then transformed into item profits (ci,j)j=1,...,n via the following projection:

Here, ϵi,j is sampled uniformly from the interval [1-¯ ϵ, 1+¯ ϵ], and B denotes a random projection matrix. Value t is used to parametrize the degree of the projection. Item weights are sampled uniformly from 3 to 8 and remain fixed across the dataset, ensuring that uncertainty affects only the objective coefficients. The capacity limits bk equal the half of the sum of each constraint weights (wk 1 , . . . ,wkn ) and also remain fixed. We use m = 10, n ∈ {100, 200, 300}, t = 8, p = 12 and ϵ = 0.5, corresponding to substantially larger

3.3 HandlingMultiple Decompositions

Our current approach considers only the subset X1 ∈ ˜X when computing the loss. Although this choice improves overall computational efficiency, it raises an important question: the selection of variables to serve as the basis for duplication (see Eq. (3)) is inherently arbitrary. This nontrivial design choice can have a significant impact on the performance of our method. We propose to mitigate this issue by dynamically changing, during training, which variables Xi ∈ ˜X are used for duplication and by leveraging multiple decompositions. A problem with m constraints then admits m possible decompositions. Let us denote by X(d) 1 and X∗(d) 1 the variables and their optimal values in the corresponding subproblem for decomposition d ∈ [[1,m]], respectively. Their associated losses L(d) 1 , and L(d) 2 are defined as follows.

Intuitively, these updated losses convey more information about the underlying problem. We propose to select which decomposition to use at random at each training iteration. As a drawback, since each decomposition d selects a different variable and constraint for the main subproblem, it yields different Lagrangian multipliers µ∗(d) (ci) and subproblem solutions X∗1 (d) (ci). We therefore need to construct a separate training set for each decomposition, S (d) = {(zi, ci,X∗1 (d) (ci), µ∗(d) (ci))}ni=1 , which is more time-consuming, but can still be performed offline. Moreover, these m constructions are fully independent and thus easily parallelizable. The resulting training algorithm is depicted in Algorithm 2. It simply consists in executing Algorithm 1 (without the random initialization of the model parameters θ) with a randomly selected decomposition at each epoch (Line 4).

4 Case Studies

This paper addresses two challenging problems in the DFL literature: the multi-dimensional knapsack problem and the quadratic portfolio optimization problem. These benchmarks problem instances than those typically considered in the literature. For example, the recent survey by Mandi et al. reports knapsack instances with only a single constraint and at most 45 items. The predictive model Ωθ is a simple linear regression. The sub-problems are one-dimensional knapsack problems and are solved using Gurobi [Gurobi Optimization, LLC, 2024] (version 12.0.0). Finally, for each size n, 1300 instances are generated and split as follows: 200 for training, 100 for validation and 1000 for the final evaluation.

Quadratic portfolio optimization. Given predicted asset returns c ∈ Rn and a positive semidefinite covariance matrix Σ ∈ Rn×n with average entry ¯Σ, the decision maker chooses portfolio weights x ∈ Rn to maximize expected return, c⊤x, while respecting a risk budget, x⊤Σx ≤ γ¯Σ, and fully investing the available capital, 1⊤x = 1. Again, the instances are generated following the procedure detailed by Mandi et al. using the hyperparameters they recommend. In particular, five features and a simple linear regression for the predictive model are used. We consider instances with n = {50, 200, 400} assets. For comparison, [Elmachtoub and Grigas, 2021] consider instances with at most 50 assets. The sub-problems are solved using Gurobi [Gurobi Optimization, LLC, 2024] (version 12.0.0). Finally, 10125 problems are generated and split as follows: 100 for training, 25 for validation and 10000 for the final evaluation. Interestingly, when choosing the risk constraint as the main constraint, the resulting LD main subproblem:

admits a closed-form optimal solution once we relax the nonnegativity constraint on the decision variable, given by

We refer the readers to Appendix B for the derivation of this exact closed-form optimal solution. This allows us to design an exact differentiation scheme to compute this solution. During training, we replace the solver-provided solution X∗1 (c, µ) with the analytical expression ¯X∗1 (c, µ) (Eq. (17)). This eliminates the need for a solver and a generic DFL technique because its gradient can be obtained directly via automatic differentiation tools.

5 Computational Experiments

This section presents the results of our computational experiments. All tests were conducted on a cluster of CPU nodes equipped with AMD EPYC 9654 and Intel Xeon Gold 6148 (Skylake) processors, with each run restricted to a single core and 4 GB of RAM. The optimization problems are solved using Gurobi [Gurobi Optimization, LLC, 2024] (version 12.0.0). We compare our approach against a standard mean squared error loss (MSE) and the DFL algorithms SPO+ [Elmachtoub and Grigas, 2021; Mandi et al., 2020] and IMLE [Niepert et al., 2021]. For these baselines, we rely on the implementations provided in PyEPO [Tang and Khalil, 2024b]. Predictive models are implemented and trained using PyTorch [Paszke et al., 2019]. We also compare our approach with other methods focussed on scaling: CaVE and CaVE+ [Tang and Khalil, 2024a] as well as LAVA [Berden et al., 2025]. These methods are less generic than ours, and cannot be applied to quadratic portfolio optimization. We do not compare with E-PLL [Defresne et al., 2023], a solver-free neuro-symbolic approach for discrete graphical models (cost function networks), as it does not directly apply to our continuous quadratic portfolio nor to the global linear capacity constraints of the multi-dimensional knapsack; a detailed comparison on suitable discrete benchmarks is left for future work. Our Lagrangian decomposition framework supports two modes, single and multiple decomposition (see the corresponding algorithms), with two loss functions available for each mode. In addition, it can be integrated into any DFL method. Finally, because our approach introduces additional overhead for optimizing the multipliers when constructing the training set, we consider two trainingbudget settings. In the first, the total budget (including dataset construction) is matched to that of the corresponding DFL baseline (=). In the second, we allow the multipliers to be trained until convergence (∞). We denote our approaches as LD(η, γ, L,B), where η ∈ {SPO+, IMLE}, γ ∈ {static, multiple, exact}, L ∈ {L1,L2}, and B ∈ {=,∞}. Method exact refers to the exact differentiation we designed for the portfolio problem. As evaluation metrics, we report the regret achieved by each approach over the test dataset, as commonly done in related works. Each model is trained with ten random seeds to assess the stability of the training procedure. The implementation is available at https://github.com/corail-research/DFL-LD.

5.1 Results: Performances with Unlimited Budget

This first experiment compares our approaches with prior baselines under an unrestricted training budget (denoted by the ∞ variant), meaning that all methods are trained to convergence, which is a commonly adopted assumption in the field. The results, reported as averages with 95% confidence intervals over 10 random seeds, are summarized in Table 1. In all settings, we use early stopping and retain the checkpoint achieving the lowest validation relative regret, the training time reported in Table 1 corresponds to the wall-clock time at which this best validation model is reached. We note that for the quadratic portfolio problem, the multipledecomposition strategy is less relevant because one of the possible choices for the main subproblem (keeping the linear budget constraint) leads to a weakly informative optimization layer. We therefore focus on the decomposition that retains the quadratic risk structure (the subproblem forms are provided in Appendix A). Similarly, our L2 variant is less relevant for SPO+, as this differentiation technique is designed as a surrogate for regret and is therefore naturally aligned with the regret-like loss L1, whereas L2 does not correspond to a regret. Overall, we observe that the best performance is consistently achieved by one of our proposed approaches. Interestingly, neither IMLE nor SPO+ dominates across all settings, which is consistent with the findings of Mandi et al. [Mandi et al., 2024]. We view this as an important strength of our framework, as it can be combined with the most suitable decision-focused learning technique depending

Quadratic portfolio optimization (2 constraints)

on the specific problem. Moreover, we observe that baselines designed to achieve low training times (CaVE, CaVE+ and LAVA) perform worse than other DFL algorithms IMLE or SPO+. This observation may be explained by the noisy nature of the datasets. Although CaVE and CaVE+ converge significantly faster than our approaches, our methods converge to models with superior final performance. Finally, the exact LD-based differentiation achieves the best results for the quadratic portfolio, showing the versatility of LD in DFL.

5.2 Analysis: Training with a Restricted Budget

This second set of experiments analyzes the behavior of the methods under a fixed training-time budget, where all approaches are allocated the same execution time. For our methods, this budget includes the construction of the training set (denoted by the = variant). We restrict the analysis to the largest instances and to our best-performing approaches. The results are shown in Figures 2a and 2b. We report the relative regret on the test set for different training-time budgets. Hyperparameters are tuned independently for each budget (via validation), so each point corresponds to a different trained model rather than successive checkpoints of a single run. As expected, thanks to the Lagrangian decomposition, our approach solves the underlying optimization problems more efficiently at each training iteration, as the resulting subproblems are easier to solve. This allows training epochs to be executed faster, which ultimately leads to improved performance compared to the baselines under the same trainingtime budget. This results suggest that the Lagrangian multipliers optimization overhead can be controlled while remaining competitive. To confirm the scalability of our approach, we report the number of training epochs completed within each training-time budget in Figures 3a and 3b. The results directly corroborate this hypothesis, as at least one LD-based approach consistently achieves the highest number of epochs for a given time budget, among DFL baselines.

5.3 Analysis: Benefits of Parallelization

Another benefit of the Lagrangian decomposition is that, in the multiple decomposition setting, each decomposition yields an independent sub-problem with its own set of Lagrangian multipliers. As a result, the offline datasetconstruction step (i.e. computing these multipliers via subgradient iterations) is naturally amenable to parallelization. We evaluate the benefits of parallelization on the multidimensional knapsack problem with 200 items. In all runs, we allocate 10 CPU cores, we then vary how many parallel worker processes are used to distribute the independent multiplier optimizations across these cores (from 1 worker, i.e., fully sequential, up to 10 workers, i.e., one per decomposition). Starting from an execution time of 4,218 seconds to perform 1,000 subgradient iterations per multiplier, we reduce the runtime to 882 seconds with 5 workers and to 462 seconds with 10 workers. This corresponds to a speed-up of approximately 9.12×, indicating that the overhead induced by parallelization is negligible. We note that this parallelization opportunity is specific to the multiple-decomposition variant, where several independent multiplier optimizations must be carried out. In particular, it provides marginal benefit for the quadratic portfolio problem in which the multiple decomposition setting is not used. Finally, these results show that our approach is well suited to scaling decisionfocused learning to large-scale problems, provided they admit a Lagrangian decomposition, as is the case for the multidimensional knapsack and the quadratic portfolio problem.

6 Conclusion

Decision-focused learning is a promising paradigm for addressing predict-then-optimize problems. However, its practical deployment is often hindered by high computational costs, as it requires solving a constrained optimization problem for each training instance at every iteration. This significantly limits its ability to scale to large combinatorial problems. In this paper, we introduced a new methodology based on Lagrangian decomposition to substantially improve the scalability of DFL algorithms. We proposed a new surrogate objective and two associated loss functions for training the prediction model. Experiments on large-scale instances, up to eight times larger than those commonly considered in the literature, demonstrate that our approach can be seamlessly integrated into standard DFL methods, leading to improved performance. Scalability can be further enhanced through parallelization. As future work, we plan to extend our approach to other combinatorial optimization problems in more realistic settings, such as shortest path and shift scheduling. Another promising direction is to investigate alternative strategies for updating the Lagrangian multipliers during training. Finally, establishing a theoretical relationship between the regret on the decomposed problem and the accuracy on the original primal problem is also an interesting direction to pursue.

A Portfolio: Forms of the Main Subproblem

The quadratic portfolio problem involves two constraints: the linear budget constraint 1⊤x = 1 and the quadratic risk constraint x⊤Σx ≤ γ¯Σ. In our Lagrangian-decomposition framework, we may choose either constraint to define the main subproblem. We detail the two resulting cases below.

Linear budget constraint as main constraint. If the budget constraint is selected, the main subproblem reads

This subproblem admits a closed-form solution: it places all the mass on an asset attaining argmaxi(ˆ ci + µ2,i) (and zero elsewhere). As a result, it is computationally trivial and does not capture the quadratic risk structure of the original formulation, which makes it less informative for training.

Quadratic risk constraint as main constraint. If instead the quadratic constraint is selected, the main subproblem becomesmax

Unlike the previous case, this formulation retains the covariance matrix Σ in the optimization layer, and is therefore more faithful to the structure of the original quadratic portfolio problem.

B The Exact Method

In the Lagrangian Decomposition, we choose the quadratic constraint in the main sub-problem. Thus, this sub-problem 1 is

Analysis of another problem. We consider the problem

We define J(w) = c⊤w, F(w) = w⊤Σw - γ ¯Σ and K the space of feasible solutions.

We note that the J and F are continuous and convex with respect to w. Moreover, F ′ (w) is nonzero for all w ∈ K with w ̸= 0, which ensures that the constraint qualification condition is satisfied.

We can thus apply the Kuhn and Tucker Theorem. For any nonzero w ∈ K, w is a global maximizer of J over K if and only if

We have

1 Since there are only two constraints in the Portfolio problem, we have µ = (µ2).

Since c ̸= 0, we must have p ̸= 0. It follows that

We conclude that the optimal solution of (21) is w∗ = √ γ¯Σ c⊤Σ-1 c Σ-1c.

Back to our sub-problem. Once the constraint X1 ≥ 0 is removed, the solution to

is given by

We observe that the closed-form solution is differentiable with respect to c, and its gradient can be obtained directly via PyTorch's automatic differentiation. These findings motivate a problem-specific differentiation scheme for the portfolio model. During training we replace the solver-provided solution X∗1 (c, µ) with the analytical expression ¯X∗1 (c, µ). This eliminates the need for a solver and a generic DFL technique.

Acknowledgments

This research has been supported by the CRM and the Simons Foundation, thematic program Combinatorial Optimization and Data Science. It has been partly funded by the European Research Council (ERC) under the EU Horizon 2020 research and innovation program (Grant No 101002802, CHAT-Opt) and by the Canada Research Chairs Program in Healthcare Analytics. We also acknowledge the support of the Natural Sciences and Engineering Research Council of Canada (NSERC RGPIN-2025-04995 and RGPIN-202203964).