Multi-Stage Predict+Optimize for (Mixed Integer) Linear Programs

Authors: xinyi  hu
jasper c h lee
jimmy h m lee
peter j stuckey
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-13 08:37:19

Abstract: The recently-proposed framework of Predict+Optimize tackles optimization problems with parameters that are unknown at solving time, in a supervised learning setting. Prior frameworks consider only the scenario where all unknown parameters are (eventually) revealed simultaneously. In this work, we propose Multi-Stage Predict+Optimize, a novel extension catering to applications where unknown parameters are revealed in sequential stages, with optimization decisions made in between. We further develop three training algorithms for neural networks (NNs) for our framework as proof of concept, both of which handle all mixed integer linear programs. The first baseline algorithm is a natural extension of prior work, training a single NN which makes a single prediction of unknown parameters. The second and third algorithms instead leverage the possibility of updating parameter predictions between stages, and trains one NN per stage. To handle the interdependency between the neural networks, we adopt sequential and parallelized versions of coordinate descent for training. Experimentation on three benchmarks demonstrates the superior learning performance of our methods over classical approaches.

Main paper: 
 introduction: 
 Constrained optimization problems can frequently model applications in everyday life, yet, the parameters of the problem are unknown at the time of solving. Consider, for example, a real-world application where hospital administrators need to schedule shifts for nurses, so as to minimize the total costs for hiring nurses while meeting the patient load. Here, the shifts need to be decided before the (real-time) patient demand is known, which requires predicting the demand when scheduling.In the present work, we focus on the supervised learning setting, where unknown parameters can be predicted using relevant features, and historical (features, parameters) pairs serve as training data for the prediction model. The goal is to estimate the unknown parameters based on the related features, such that when we solve the optimization problem using the estimated parameters, the resulting solution is good even under the later-revealed true parameters. Classic approaches, for example by learning predictors using (regularized) ℓ 2 loss, do not necessarily work well -low prediction error in parameter space does not guarantee good performance of the estimated solution according to the optimization objective. The influential framework of Predict+Optimize proposed by Elmachtoub and Grigas [9] instead uses the more effective regret loss, which incorporates information about the optimization problem. However, the framework limits the unknown parameters to appear only in the objective and not the constraints -if uncertainty in constraints is mis-predicted, the resulting estimated solution might not even be feasible under the true parameters. Recent works by Hu et al. [15, 16] thus propose a Two-Stage adaptation of Predict+Optimize, explicitly modelling 1) a prediction stage and 2) a recourse stage which corrects infeasible solutions into feasible ones. The new two-stage framework is therefore applicable even when the optimization constraints contain uncertainty.However, the two-stage framework essentially assumes that all the unknown parameters are revealed simultaneously, excluding applications where such information is gradually released and new decisions need to be made across many stages (e.g. in a daily/weekly manner). Crucially, in these applications, predictions can also be updated between stages, in light of the new information and past committed actions. Consider again the example of scheduling shifts for nurses. A typical facility might have an appointment system, with reservations closing the day before each working day. As opposed to a two-stage modelling, a more practical approach would treat each day in a work week as its own stage where new information (the precise appointments the next day) is learned, inducing both new optimization decisions and updated predictions.The concrete contributions of this paper are three-fold.Framework We propose and formalize the new framework of Multi-Stage Predict+Optimize (Section 3), where unknown parameters are revealed across multiple stages, inducing new optimization decisions and updated parameter predictions.The flexibility to update (future) parameter predictions in each stage introduces intricate challenges to the training process, which should train a prediction model per stage. The challenges are both in predictive power and in computation time. The performance of predictors across stages are intertwined and interdependent: the "goodness" of a prediction depends on actions in other stages, which in turn depends on predictions of those other stages. Such dependency can also cause serialization issues that could drastically lengthen training time.In Section 4, we propose three neural network training algorithms for our framework, assuming that the optimization problems can be formulated as mixed integer linear programs (MILPs): 1) a baseline algorithm that directly generalizes the two-stage algorithm of Hu et al. [16] , training only a single neural network predictor, 2) a sequential coordinate descent training algorithm which trains a neural network model per stage, and each stage is considered a "coordinate", and 3) a parallel version of coordinate descent. These algorithms trade off between training time and predictive performance.We apply these methods to three benchmarks (Section 5) to empirically demonstrate their superior learning performance over classical training methods, as well as the computation/prediction tradeoff between the proposed methods.We note that there are other lines of work tackling similar settings, where unknown parameters in optimizations are also revealed in a sequential fashion. Perhaps the most well-known is multi-stage stochastic optimization [14, 30] . The main difference between our work and multi-stage stochastic optimization is supervised vs unsupervised learning. Our framework (and Predict+Optimize in general) has features that help making parameter predictions, whereas (non-contextual) stochastic optimization does not and requires different techniques to tackle. See Appendix A.3 for a detailed discussion on the connection and comparison between our framework and stochastic optimization.We include a brief literature review here in the main body. See Appendix A for a detailed exposition.Elmachtoub and Grigas proposed the influential framework of Predict+Optimize [9] , with lots of followup work in the community on improving computational efficiency [23, 24] , predictive accuracy [7, 18, 22] , on types of applicable optimization problems [12, 17, 36] , and applying to specific real-world scenarios [6, 34, 35] . More recently, Hu et al. [15, 16] proposed adaptations of the framework to handle uncertainty in optimization constraints, including the Two-Stage framework which our work is most related to.Predict+Optimize also sits in a broader line of work on decision-focused learning, including works that learn prediction models for unknown parameters but with different goals/loss measures [25, 28] , as well as other works that invent techniques for differentiating through optimization problems [1, 2, 36] .Outside of decision-focused learning, our work is also somewhat related to (multi-stage) stochastic programming [14, 30] . The main difference, again, is supervised vs unsupervised learning. See Appendix A.3 for a more detailed comparison.Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?Answer: [Yes] Justification: Our abstract and introduction state our contributions, including the new framework and training algorithms. We state clearly that the algorithms are applicable to MILPs.Guidelines:• The answer NA means that the abstract and introduction do not include the claims made in the paper. • The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or NA answer to this question will not be perceived well by the reviewers. • The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings. • It is fine to include aspirational goals as motivation as long as it is clear that these goals are not attained by the paper. 
background -two-stage predict+optimize: 
 In this section, we describe Two-Stage Predict+Optimize [16] , which is prerequisite to understanding our contributions. Without loss of generality, the framework is stated for minimization problems.A parameterized optimization problem (Para-OP) P (θ) is defined as finding:x * (θ) = arg min x obj(x, θ) s.t. C(x, θ) where x is a vector of decision variables, θ is a vector of parameters, obj is a function mapping decisions x and parameters θ to a real objective value that is to be minimized, and C is a set of constraints that must be satisfied over x under parameters θ. We call x * (θ) an optimal solution under parameters θ, and obj(x * (θ), θ)) the optimal value under parameters θ. When the parameters are all known, a Para-OP is just a classical optimization problem (OP).In the Predict+Optimize setting (from the original framework [9] , to the two-stage framework [16] , and to our multi-stage framework later on), each instantiation of the true parameter vector θ has an associated feature matrix A. These features are relevant information that can help a model predict θ.In the following, we number the stages by 0 and 1 in Two-Stage Predict+Optimize [16] to make the framework look more similar to the multi-stage framework we propose later in Section 3.Stage 0 The practitioner uses a prediction model, which takes in a feature matrix A, to compute a vector of estimated parameters θ. The Stage 0 solution x0 is then computed as x(0) = arg min x obj(x, θ) s.t. C(x, θ)The Stage 0 solution x(0) should be interpreted as a form of soft commitment that can be modified in Stage 1, subject to a penalty.The true parameters θ are revealed, and the practitioner wishes to compute an updated Stage 1 solution x(1) from x(0) , subject to a (problem-specific) penalty function P en(x (0) → x(1) , θ) which depends on both the softly-committed Stage 0 x(0) , the final Stage 1 solution x(1) and the true parameters θ. More specifically, the Stage 1 solution x(1) is computed asx(1) = arg min x obj(x, θ) + P en(x (0) → x, θ) s.t. C(x, θ)The Stage 1 solution x(1) should be interpreted as a hard-committed final action, and note that it is guaranteed to be feasible under the true parameters θ.The prediction θ is evaluated using the post-hoc regret [16] , which is the sum of two terms: (a) the difference in objective between the true optimal solution x * (θ) and the final Stage 1 solution x(1) under the true parameters θ, and (b) the penalty incurred by modifying x(0) to x(1) . Formally, the post-hoc regret function P Reg( θ, θ) (for minimization problems) is:P Reg( θ, θ) = obj(x (1) , θ) -obj(x * (θ), θ) + P en(x (0) → x(1) , θ)The goal of a prediction model is to make predictions θ so as to minimize the post-hoc regret. We emphasize again that the main difference between Predict+Optimize frameworks and stochastic programming frameworks is that in Predict+Optimize, a prediction model has access to features relevant to the true parameters in order to make a prediction. Stochastic programming, on the other hand, frequently operates solely at the level of the distribution over the true parameters.In the second experiment, we showcase our framework on an MILP. The unknown parameters appear in both the objective and constraints. A person wants to make an investment plan for buying several types of financial products next year to maximize the investment profit, under limited capital. The investment profit contains 3 parts: 1) the interest gained from the products owned, 2) the market prices of the products owned at the end of the year, and 3) profits from selling products minus costs for buying products minus transaction fees. The capital for the whole year is given. However, the market price of each product in each quarter/month is revealed only at the beginning of the quarter/month, and the interest of owning each product in each quarter/month is revealed only at the end of the quarter/month. The person will estimate the market prices and interests based on past experiences, considering features such as the product type, the financial condition and operational capabilities of the company to which the product belongs, and so on.Suppose there are T quarters/months, and N investment products. The decision variables are: 1) a natural vector x ∈ N T ×N , where x i,j represents the number of product j on hand at the end of quarter/month i, 2) an integer vector y ∈ Z (T -1)×N , where y i,j represents the number of product j bought or sold in quarter/month i, and 3) a natural vector z ∈ N (T -1)×N , where z i,j represents whether the transaction fee is paid for product j in month i. Let p i,j denote the interest of product j in month i, w i,j denote the market price of product j in month i, C denote the capital for the whole year.We assume that the end of quarter/month t, i.e., the beginning of quarter/month (t + 1), is Stage t. At Stage 0, i.e., the beginning of quarter/month 1, the person can buy some products without paying a transaction fee. The market price of each product at this time is known, i.e., w 1 = (w 1,1 , . . . , w 1,N ) are given. The unknown parameters in this OP are p ∈ R T ×N and w = (w 2 , . . . , w T ) ∈ R (T -1)×N . At the beginning of each subsequent quarter/month, the person can buy more products or sell the products owned but needs to pay a transaction fee. For simplicity, we assume that the transaction fee for buying/selling product i in quarter/month j is linear in the market price of product i in quarter/month j. Here, the linearity factor is independent of the request. That is, if the person buys/sells k number of product i in quarter/month j, the person has to spend kσw ij , where σ ≥ 0 is a non-negative tunable scalar parameter, and we call it transaction factor. At Stage 0, i.e., the beginning of quarter/month 1, the person uses the predicted interests p(0) and market prices ŵ(0) to make the plan. The Stage 0 OP can be formulated as:EQUATIONEQUATIONx i,j = y i,j + x (i-1),j , ∀i ∈ {2, . . . , T }, j ∈ {1, . . . , N } (15)EQUATIONwhere obj( p0 , w 1 , ŵ0 , x, y, z)= T i=1 p(0) i,j x i,j + N j=1 ŵ(0) T,j x T,j -( N j=1 w 1,j x 1,j + T i=2 N j=1 σ ŵ(0) i,j z i,j + T i=2 N j=1 ŵ(0) i,j y i,j )represents the objective, which is to maximize the investment profit; Equations ( 13) and ( 14) ensure that the money spent on buying products and transaction fees will not exceed the capital available; Equations ( 15), ( 16), and ( 17) formulate the relationships among three decision variables x, y, and z. i represents the market price and the interest of quarter/month i estimated on quarter/month t. The investment decisions x, y, z for the prior t quarters/months have already happened and cannot be changed, which are committed variables. There is no penalty function in this scenario.Atx(t) , ŷ(t) , ẑ(t) = arg max x,y,z obj(p[1 : t] ⊕ p(t) , w 1 , w[2 : t + 1] ⊕ ŵ(t) , x, y, z) s.t. N j=1 w 1,j x 1,j ≤ C, N j=1 w 1,j x 1,j + k i=2 N j=1 σw i,j z i,j + k i=2 N j=1 w i,j y i,j ≤ C, ∀k ∈ {2, . . . , t} N j=1 w 1,j x 1,j + t+1 i=2 N j=1 αw i,j z i,j + k i=t+2 N j=1 α ŵ(t) i,j z i,j + t+1 i=2 N j=1 w i,j y i,j + k i=t+2 N j=1 ŵ(t) i,j y i,j ≤ C, ∀k ∈ {t + 1, . . . , T } x i,j = y i,j + x (i-1),j , ∀i ∈ {2, . . . , T }, j ∈ {1, . . . , N } z i,j ≥ y i,j , ∀i ∈ {2, . . . , T }, j ∈ {1, . . . , N } z i,j ≥ -y i,j , ∀i ∈ {2, . . . , T }, j ∈ {1, . . . , N } x i,j = x(t-1) i,j, ∀i ∈ {1, . . . , t}, j ∈ {1, . . . , N }y i,j = ŷ(t-1) i,j , z i,j = ẑ(t-1) i,j, ∀i ∈ {2, . . . , t}, j ∈ {1, . . . , N } whereobj(p[1 : t] ⊕ p(t) , w 1 , w[2 : t + 1] ⊕ ŵ(t) , x, y, z) = t i=1 N j=1 p i,j x i,j + T i=t+1 N j=1 p(t) i,j x i,j - N j=1 w 1j x 1j -( t+1 i=2 N j=1 αw i,j z i,j + T i=t+2 N j=1 α ŵ(t) i,j z i,j ) -( t+1 i=2 N j=1 w i,j y i,j + T i=t+2 N j=1 ŵ(t) i,j y i,j ) + N j=1 ŵ(t) T j x T jQuestion: Does the paper discuss the limitations of the work performed by the authors?Answer: [Yes]Justification: Section 3 states the "total recourse" assumption that we assume all Stage t optimization problems are always satisfiable. Section 4 states that we only work with MILP formulations. Section 5 discusses the tradeoff between training time and predictive accuracy.Guidelines:• The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper. • The authors are encouraged to create a separate "Limitations" section in their paper.• The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be. • The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on a few datasets or with a few runs. In general, empirical results often depend on implicit assumptions, which should be articulated. • The authors should reflect on the factors that influence the performance of the approach.For example, a facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or a speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon. • The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size. • If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness. • While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an important role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations. • Providing as much information as possible in supplemental material (appended to the paper) is recommended, but including URLs to data and code is permitted. 
multi-stage predict+optimize: 
 In this section, we present our new framework of Multi-Stage Predict+Optimize, which models applications where unknown parameters are revealed across T different stages.Consider again the Para-OPx * (θ) = arg min x obj(x, θ) s.t. C(x, θ)We view the true parameter vector θ as (θ 1 , . . . , θ T ), where each θ t is the sub-vector of parameters released at Stage t. Similarly, we also view the vector of decision variables x as (x 0 , . . . , x T ), where each x t is the sub-vector of decision variables that are hard-committed in Stage t (e.g. via a concrete real-world action taken at Stage t) and soft-committed in prior stages (e.g. a tentative nurse schedule).At a high level, the parameters θ t are revealed at Stage t, and a model makes a prediction θ(t) = ( θ(t) t+1 , . . . , θ(t) T ) of all the unrevealed parameters. Then, the practitioner solves the Stage t optimization problem which we define later in the section. The decision variables x t are newly hard-committed, whereas the decision variables x t+1 , . . . , x T are soft-committed with potential to be modified in future stages (at the cost of a penalty). This process is repeated until all stages are completed.In the rest of the section and paper, we will use the standard notation of θ[t 1 : t 2 ] = (θ t1 , . . . , θ t2 ) to denote sub-vectors (treated as arrays), and use the notation ⊕ for vector concatenation.Now we formally define Multi-Stage Predict+Optimize framework. In Appendix B, we also present the hospital scenario from the Introduction as a detailed example of applying this framework.Stage 0 None of the true parameters have been revealed. M odel 0 takes the feature matrix A and predicts θ(0) = ( θ(0) 1 , . . . , θ(0) T ). The practitioner then computes the Stage 0 solution x(0) asx(0) = arg min x obj(x, θ(0) ) s.t. C(x, θ(0) )The decision variables x(0) 0 are hard commitments, whereas the rest of the decision vector x(0) 1 , . . . , x(0) T are soft commitments.Stage t (for 1 ≤ t ≤ T ) The true parameters θ 1 , . . . , θ t-1 were previously revealed, and θ t is newly revealed. M odel t makes a prediction θ(t) = ( θ(t) t+1 , . . . , θ(t) T ) using 1) the feature matrix A, 2) the previous stage solution x(t-1) and 3) the revealed true parameters θ 1 , . . . , θ t . For computational efficiency reasons, M odel t may instead take any subset or derived functions of the above inputs. For example, M odel t can choose whether or not to incorporate the revealed true parameters θ 1 , . . . , θ t as input. While these revealed parameters can serve as additional features, potentially guiding and correcting current predictions more effectively, they also increase training time (and inference time to a smaller extent). The trade-off between prediction improvement and additional training time depends on the optimization problem, model structure, and training data. Therefore, whether to utilize the revealed true parameters can be considered a hyperparameter that should be tuned for each application using the available training data. See Appendix H.1 for a more detailed discussion.Afterwards, the practitioner computes the Stage t solution x(t) using the following Stage t optimization problem, which crucially modifies the original Para-OP by: 1) introducing a penalty term P en t (x (t-1) → x, θ[1 : t]) capturing the cost of changing the Stage t -1 solution x(t-1) to the new Stage t solution x, and 2) introducing the constraint that x[1 : t -1] = x(t-1) [1 : t -1], namely, hard commitments from prior stages cannot be modified in the current Stage t. This constraint is a form of a non-anticipativity constraint in the stochastic programming literature [14] .x(t) = arg min x obj(x, θ[1 : t] ⊕ θ(t) ) + P en t (x (t-1) → x, θ[1 : t]) s.t. C(x, θ[1 : t] ⊕ θ(t) ) and x[1 : t -1] = x(t-1) [1 : t -1]The Stage t solution, by construction, has x(t) [1 : t -1] being equal/compatible with the hard commitments from prior stages. The new hard commitments are x(t) t , and the rest of the decision vector x(t) [t + 1 : T ] are new soft commitments we make for the future stages.At t = T , the prediction θ(T ) is a length-0 vector since all the true parameters will have been revealed.For the rest of the paper, we make the assumption that these Stage t optimizations are always satisfiable, regardless of the prior stage solutions, prior+current predictions and revealed parameters. For practical applications, this assumption is both natural and essential. In real-world scenarios, encountering an unsatisfiable condition can lead to catastrophic outcomes. Therefore, before using the application, the domain expert should always have designed the underlying real-world system to have recourse actions to mitigate bad prior commitments (at cost/penalty) and to prevent catastrophe, and furthermore model such recourse actions in the (multi-stage) optimization problem. Any system and the corresponding formulation of multi-stage optimization problem lacking such recourse should not be used/executed. It is thus a reasonable assumption and a practical responsibility our framework asks of the practitioner, that recourse actions are always designed into the underlying system and modelled, so that our feasibility assumption is satisfied.Evaluation: The sequence of predictions θ(0) , . . . , θ(T ) , which along with the true parameters θ induces the sequence of solutions x(0) , . . . , x(T ) , is evaluated using a generalized notion of post-hoc regret, defined as follows (for a minimization problem):P Reg(( θ(0) , . . . , θ(T ) ), θ) = obj(x (T ) , θ) -obj(x * (θ), θ) + t P en t (x (t-1) → x(t) , θ[1 : t])where x * (θ) is again the optimal in-hindsight vector of decisions for the original Para-OP.We note that if a problem has only 2 stages (Stages 0 and 1), then our framework of Multi-Stage Predict+Optimize indeed captures Two-Stage Predict+Optimize described in Section 2. 
end-to-end training algorithms on milps: 
 In this section, we give 3 training algorithms for neural network models for the Multi-Stage Pre-dict+Optimize framework, under the assumption that all Stage t optimization problems can be formulated as (mixed integer) linear programs (MILPs).The first training algorithm, our baseline (Section 4.1), is a straightforward generalization of the one proposed for the two-stage framework [16] . This algorithm only trains a single neural network and reuses the same parameter predictions across all stages. Although the approach is computationally efficient, it fails to fully exploit the power of the framework, which allows for predictions to be updated at each stage.Our second and third algorithms (Sections 4.2 and 4.3) instead train one neural network per stage, each making new parameter predictions for the corresponding stage. As mentioned in Section 1, it is delicate to train these neural networks. The quality of a prediction in one stage depends on decisions in other stages, which in turn depends on predictions made in those other stages. To handle this dependency, we employ a coordinate-descent strategy, where each stage/neural network is a coordinate. We present both a sequential version and a parallel version of this strategy as training algorithms, which trade off between training time (sequential being slower) and predictive performance.We also point out that it is technically possible to train all networks simultaneously, without using coordinate descent like in Sections 4.2 and 4.3. To do so, we would instead use ground truth parameters in place of prior and future stage predictions. However, intuitively, this simpler approach should have worse predictive ability than the proposed two methods, given the interdependency of the predictors. We show experimental comparisons in Appendix H.2 confirming this intuition.In Section 5, we show empirical results comparing these algorithms and classic non-Predict+Optimize learning algorithms (e.g. standard regression models), which demonstrate that even our baseline training approach outperforms classic non-Predict+Optimize algorithms. Also, our more sophisticated approaches yield even better learning performance, at the cost of additional training time.We first present a baseline Predict+Optimize training algorithm for our multi-stage setting, via a natural extension of the two-stage approach [16] .This baseline algorithm trains a single neural network NN , which takes a feature matrix A and returns the prediction θ = NN (A). The same predictions are then reused across all the stages. More specifically, in the language of Section 3, we choose θ(t) = θ[t + 1 : T ] for this basic approach.We use standard gradient methods to train the neural network NN , with the goal of minimizing post-hoc regret as defined in Section 3:P Reg = obj(x (T ) , θ) -obj(x * (θ), θ) + t P en t (x (t-1) → x(t) , θ[1 : t])Noting the second term is independent of θ(0) and hence NN , the gradient with respect to an edge weight w e in NN isd P Reg dw e = d obj(x (T ) , θ) dw e + T t=1 d P en t (x (t-1) → x(t) , θ[1 : t]) dw eBy the law of total derivative, we can expand this tod P Reg dw e = d obj(x (T ) , θ) dx (T ) dx (T ) d θ d θ dw e + T t=1 ∂ P en t ∂ x(t-1)x(t)dx (t-1) d θ + ∂ P en t ∂ x(t)x(t-1)dx (t) d θ d θ dw eThe term d θ dwe is calculated via standard backpropagation, while the terms d obj(x (T ) ,θ)dx (T ) , ∂ P ent ∂ x(t-1) x(t)and ∂ P ent ∂ x(t) x(t-1) are calculable when the objective and penalty functions are smooth. The only non-trivial calculation is fordx (t) d θ for all t ∈ [T ].Recall that x(t) is computed from the Stage t optimization problem, as a deterministic function of x(t-1) and θ(t) (which is a sub-vector of θ here, since we reuse predictions), while x(t-1) itself also depends on θ. We thus further decompose dx (t) d θ into the following recursive computationdx (t) d θ = ∂ x(t) ∂ x(t-1) θ dx (t-1) d θ + ∂ x(t) ∂ θ x(t-1)Calculating ∂ x(t) ∂ x(t-1) θ and ∂ x(t) ∂ θ x(t-1) requires differentiating through a MILP. So instead of directly using MILP formulations for the Stage t optimization problems, we use the convex relaxation by Hu et al. [16] , which in turn adapts the approach of Mandi and Guns [22] .We also note that it is possible to use other convex relaxations and approaches to differentiate through the MILP, for example using tools like CvxpyLayers [1] . We chose Hu et al.'s calculations because their experiments showed the computational efficiency of their approach over other tools.Training only one neural network, the baseline algorithm is efficient but fails to fully harness the power of the framework in Section 3. Each stage t makes new decisions, and the "goodness" of future decisions depends on these previously committed decisions. Thus, new predictions should be made each stage for future parameters, based on prior-stage optimization decisions, so as to yield current and future optimization decisions that work well with the already-committed ones. However, the baseline algorithm ignores this information and does not update the predictions accordingly.We thus propose our second training algorithm, which trains one neural network per stage, from Stages 0 to T -1. The neural network NN t for Stage t takes the feature matrix A as input, as well as all the prior decision vector x(t-1) , and outputs the prediction θ(t) for parameters θ[t + 1 : T ]. The astute reader might recall that the proposed framework in Section 3 allows N N t to utilize the revealed parameters from stage t -1 as additional features as input. However, preliminary experiments in Appendix H.1 indicated that including such parameters does not really enhance prediction quality, while merely increasing training time. Therefore, in our current implementation, N N t does not include these revealed parameters. In general, however, such choice should be treated as a hyperparameter and made for each application.To address the dependency between the neural networks, we employ a coordinate descent approach, where each stage/neural network is treated as a "coordinate". We train the neural networks NN 0 , . . . , NN T -1 in cyclic order, repeating until termination (e.g. convergence or timeout). See Algorithm 1 for high-level pseudocode description of this sequential coordinate descent approach, as well as the parallelized version described in the next subsection.Concretely, we train NN t by considering all the other neural networks as fixed. Here, we will focus on describing the forward pass, since the backward pass gradient computations follow essentially the same strategy described in Section 4.1 -see Appendix C for details of the gradient computations.Forward pass Consider a historical (feature, parameter) pair (A, θ). We first iteratively generate the sequence of solutions x(0) , . . . , x(t) using θ(i) = NN i (A, x(i-1) ) for i ∈ [0, t -1]. Then, we compute the Stage t prediction θ(t) = NN t (A, x(t-1) ), and generate the remaining sequence of solutions x(t) , . . . , x(T ) using θ(i) = θ(t) [i + 1 : T ] for i ∈ [t, T ].Backward pass The goal is to compute the derivative of the post-hoc regret with respect to each edge weight w e of NN t . Similar to Section 4.1, instead of directly using the MILP formulation of all the Stage t optimization problems in the forward pass, we use the convex relaxation of Hu et al. and Mandi and Guns. This allows us to differentiate the modified (due to convex relaxations) post-hoc regret with respect to each w e in NN t . As mentioned, the calculations are quite similar to those in Section 4.1, and so we defer them to Appendix C. Experimental results in Section 5 demonstrate that the sequential coordinate descent training approach outperforms both the classic non-Predict+Optimize methods and the baseline Predict+Optimize approach from Section 4.1.We note that in the above description of the training implementation of NN t , there is a lot of repeated computation that can be pre-computed and reused. Since, during coordinate descent, the prior neural networks NN 0 , . . . , NN t-1 are considered fixed, the solutions (x (0) , . . . , x(t-1) ) are also fixed for a given (features, parameters) pair (A, θ) no matter how NN t is updated during training. Thus, for each (A, θ) we always pre-compute and save the sequence of solutions (x (0) , . . . , x(t-1) ), and only recompute (x t , . . . , xT ) as we update NN t through training gradient steps.While the sequential coordinate descent training approach yields accurate predictors, the computational cost is also high. Training a single neural network already requires solving sequences of optimization problems over many iterations; the serialization from training the neural networks one at a time can make the resulting training time prohibitive for applications.We thus propose to parallelize the coordinate descent approach, slightly sacrificing prediction quality for efficiency. In each coordinate descent step, we train all the neural networks in parallel. When training a particular neural network NN t , we use copies of NN 0 , . . . , NN t-1 , NN t+1 , . . . , NN T -1 from the previous descent step, but otherwise the training implementation remain the same as in Section 4.2. See also Algorithm 1.This simple change drastically improves running time (Appendix G), while only slightly decreases predictive accuracy: the post-hoc regret of models trained using parallel coordinate descent sits between that of the baseline training algorithm and the sequential coordinate descent approach. 
experimental evaluation: 
 We evaluate the proposed 3 training algorithms: Baseline, Sequential Coordinate Descent (SCD), and Parallel Coordinate Descent (PCD) on 3 benchmarks described in Appendix D. We compare these algorithms to classic non-Predict+Optimize regression models [10] : ridge regression (Ridge), k-nearest neighbors (k-NN), classification and regression tree (CART), random forest (RF), and neural network (NN). The single predictions from these classical regression models are used in test time identically to our Baseline approach (Section 4.1). We tune all algorithm hyperparameters via cross-validation -Appendix E gives all the hyperparameter types and chosen values. In particular, the termination criteria for SCD and PCD in Algorithm 1 are based on a threshold that measures the difference in training set post-hoc regrets between two (outermost) iterations of the training coordinate descent. This threshold is also treated as a hyperparameter that was tuned per each application; in the experiments presented in this paper, we used a threshold of 0.1.Due to space limitations, we present only the best results obtained among all standard regression methods (BAS) as one column in the main paper. See Appendix F for full results. Furthermore, we report mainly the prediction accuracy -see Appendix G for computational setup and detailed runtime comparisons.Given the lack of datasets specific to these benchmarks, we follow a standard Predict+Optimize experimental approach [16, 24, 7] and use real data from different problems as numerical values in our experiment instances. We include details in the individual subsections. For each benchmark, we run 30 simulations, each simulation containing a 70/30 training/test data split.Our first benchmark is a linear programming (LP) problem. An oil company is developing a production and sales plan for the upcoming T quarters/months. The company aims to maximize profits -sales revenues minus production costs -under the constraint that the amount of oil sold each quarter/month cannot exceed the customer demands. The production cost and sales price for each quarter/month are known, but the demand is revealed only at the beginning of each quarter/month after the company receives the orders. See Appendix D.1 for the detailed description and LP model.We generate production costs and sales prices following the method by Ardjmand et al. [3] . Production costs are randomly generated from [50, 100]. Two groups of sales prices are considered: low-profit product prices are randomly generated from [50, 100]; high-profit product prices are from [120, 150] . Customer demands are the unknown parameters and need prediction -we use real data from a knapsack benchmark [28] as demand parameters, where each parameter is related to 4096 features.We conduct experiments on T = {4, 12}, corresponding to 4 quarters or 12 months. For NN, Baseline, SCD, and PCD, we use 5-layer fully connected networks with 512 neurons per hidden layer.Table 1 reports the mean post-hoc regrets and standard deviations across 30 simulations for the proposed 3 methods and BAS on the problem. Appendix F.1 gives a full data table (Table 5 ) with all standard regression methods. We also report the mean and standard deviations of True Optimal Values (TOV) in the last column -optimal objective under true parameters in hindsight -for readers to use as (rough) normalization for relative errors. The results demonstrate the advantage of Predict+Optimize methods. All three proposed methods, even Baseline, beat all standard regression methods. SCD consistently achieves the best performance, followed closely by PCD. Baseline falls between the two coordinate descent methods.In Appendix F.1, Table 6 shows the percentage improvements of the proposed 3 methods against BAS. From that table, we observe that the advantage of our methods increases with the number of stages.Considering the relatively large standard deviations in Table 1 , to show the substantial performance improvements clearly, we provide "win rate" results in Figure 1 and related information in Table 7 in Appendix F.1. Here, "win rate" refers to counting the number of simulations (out of 30) where a Our second benchmark shows our framework on a MILP problem. The task is to create an investment plan for trading financial products over a year, under limited capital. The net profit consists of: 1) dividends gained from holding products, 2) market price of products in the portfolio at the end of the year and 3) net profits from trading products (selling price minus buying price) minus trading transaction fees. For simplicity, we assume that transaction fees are linear in the products' market prices. We call the proportionality constant the transaction factor, which we vary across experiments.The market price of (resp. dividends from owning) each product in each quarter/month is revealed only at the beginning (resp. end) of the quarter/month, requiring predictions. As with the previous benchmark, we use real data from a knapsack benchmark [28] for parameters needing prediction.The detailed problem description and the MILP formulation are in Appendix D.2We conduct experiments on five financial products, and T = {4, 12}, corresponding to 4 quarters or 12 months. We consider two different initial capital values: 25 and 50, and three different transaction factors: 0.01, 0.05, or 0.1. For NN, Baseline, SCD and PCD, we again use a 5-layer fully connected network with 512 neurons per hidden layer.Table 2 reports the mean post-hoc regrets and standard deviations across 30 simulations for the proposed three methods and BAS on the problem when capital=50. We report full data tables (Tables 8 and 9 ) in Appendix F.2. TOVs are reported in the last column, again for reference. The results confirm the performance ordering of SCD > PCD > Baseline > classical regression methods.Table 10 in Appendix F.2 presents the percentage improvement of our proposed methods against BAS. Curiously, contrasting the previous benchmark, the advantage of our methods against BAS appears to be non-increasing with the number of stages. We hypothesize this is an LP versus MILP issuetraining for this benchmark required relaxing integrality constraints, introducing more approximations with more stages. It is interesting future work to study if the performance decay in the number of stages is a general phenomenon for MILPs and if it can be mitigated.For the investment problem benchmark, we again provide win rate results in Table 11 in Appendix F.2.The results show that the proposed 3 methods outperform BAS in more than 96.67% of the simulations.Our last benchmark is the nurse rostering MILP problem described in Section 1. The task is to minimize the costs of hiring nurses to meet the (unknown) patient demands. After each stage's patient Similarly to the first two benchmarks, we report win rate results in Table 14 and show comparisons between the proposed methods against BAS across each individual simulation in Figure 2 in Appendix F.3. Table 14 shows that SCD pretty consistently outperforms BAS, achieving win rates of 86.67% or higher in most scenarios. PCD also demonstrates competitive performance, with win rates ranging from 70% to 83.33%, making it a viable alternative as we hypothesized.Appendix G gives the training times for each method. The training times follow the order of SCD > PCD > Baseline > classic regression methods, which is the same for predictive accuracy in the benchmarks, indicating a tradeoff between training time and accuracy. Our methods take longer runtime than most other regression methods due to having to solve sequences of linear programs during training. Among our methods, the coordinate descent methods take longer time than Baseline due to having to train more neural networks, and the sequential version naturally takes longer than the parallelized version. 
conclusion and future work: 
 We propose the first Predict+Optimize framework for scenarios where unknown parameters are revealed progressively over stages. Specifically, our proposal allows better predictions and (re)optimization at each stage as more parameters are made known. Algorithmically, we focus on MILPs--a large and widely-studied class of problems--and present three training methods for our novel framework. Empirical results in three benchmarks demonstrate better predictions from our methods over classical ones. Our methods trade off between predictive accuracy and training time.Our work establishes the feasibility of Multi-Stage Predict+Optimize, and furthermore shows that even our baseline algorithm of training a single predictor outperforms classical non-Predict+Optimize approaches. Looking into the future beyond the present work, there is ample space for algorithmic improvements. As we observed in Section 5, the current experimental results suggest that our multi-stage predict+optimize methods display rather different behaviors depending on whether the optimization is a linear program or a mixed integer program. In particular, for linear programs, the advantage of our methods over classical methods increases with the number of stages, whereas the opposite happens for integer programs. We believe it is important to investigate whether this phenomenon holds more generally. If so, it is a prudent research direction to understand whether such decay is inevitable for MILPs, or if there are algorithmic techniques to mitigate this effect.Question: Does the paper specify all the training and test details (e.g., data splits, hyperparameters, how they were chosen, type of optimizer, etc.) necessary to understand the results? Answer: [Yes] Justification: We provide all data splits, how the data was generated/sourced, and hyperparameter choices. Guidelines:• The answer NA means that the paper does not include experiments.• The experimental setting should be presented in the core of the paper to a level of detail that is necessary to appreciate the results and make sense of them. • The full details can be provided either with the code, in appendix, or as supplemental material. 
experiment statistical significance: 
 Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments? Answer: [Yes] Justification: We give standard deviations in our results. Guidelines:• The answer NA means that the paper does not include experiments.• The authors should answer "Yes" if the results are accompanied by error bars, confidence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper. • The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions).• Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations. • The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)? Answer: [NA] Justification: Our work does not release any artifacts that can be easily misused. Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters. • Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images. • We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort. 12. Licenses for existing assets Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected? Answer: [Yes] Justification: We cite relevant data sources in Section 5. Guidelines:• The answer NA means that the paper does not use existing assets.• The authors should cite the original paper that produced the code package or dataset.• The authors should state which version of the asset is used and, if possible, include a URL. • The name of the license (e.g., CC-BY 4.0) should be included for each asset.• For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided. • If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset. 
 
 Appendix: In this section, we first summarize the related works in Predict+Optimize, and then summarize other works related to learning unknowns in optimization problems, but outside of Predict+Optimize.As mentioned in Section 1, prior works all focus on the case where all unknown parameters are revealed simultaneously. Most of them have focused on the regime where the unknown parameters only appear in the objective and use the regret proposed by Elmachtoub et al. [9] as the loss function.Since the regret loss is usually not (sub-) differentiable, and gradient-based methods do not apply, they proposed ways to overcome the non-differentiability of the regret. Elmachtoub et al. [9] propose a differentiable surrogate function for the regret function, while Wilder et al. [36] relax the integral objective in constrained optimization and solve a regularized quadratic programming problem. Mandy and Guns [22] focus on mixed integer linear programs and propose an interior point based approach.In addition to computing the (approximate) gradients of the regret function or approximations of it, another way to deal with the non-differentiability of the regret is to exploit the structure of optimization problems to train models without computing gradients. Demirović et al. [7] investigate problems amenable to tabular dynamic programming and propose a coordinate descent method to learn a linear prediction function. Hu et al. [17] extend their framework, to enable problems solvable with a recursive or iterative algorithm to be tackled in Predict+Optimize. Guler et al. [12] proposes a divide and conquer algorithm, extending the work of Demirović et al. [7] in a different manner so that the algorithm can deal with problems with the linear objective function.As for Predict+Optimize with unknown parameters also in constraints, Hu et al. [15] first propose the post-hoc regret loss function and a framework for packing and covering LPs with unknown parameters in both objectives and constraints. They [16] further advocate a conceptually simpler framework, which enable solving MILPs with unknown constraints. Besides, there are also works applying Predict+Optimize to a wide range of real-world problems, including maritime transportation [35] , last-mile delivery [6] , and trading in renewable energy [34] .Now we summarize other works related to learning unknowns in optimization problems, particularly those outside of Predict+Optimize. These works can be placed into two categories.One category considers learning unknown parameters but with very different goals and measures of loss. For example, CombOptNet [28] and Nandwani et al. [25] focus on learning parameters so as to make the predicted optimal solution (Stage 0 estimated solution in the proposed framework) as close to the true optimal solution x * as possible in the solution space/metric. However, these works also assume that all unknown parameters are revealed simultaneously, and thus cannot be applied to applications where unknown parameters are revealed progressively over several stages. Furthermore, experiments in Two-Stage Predict+Optimize [16] show that these other methods yield worse predictive performance when evaluated on the post-hoc regret. [15, 16] and Mandi and Guns [22] to perform gradient computations, we could in principle use any of the aforementioned other works. However, we point out that our main contribution is not the gradient computation method but the two training algorithms of the set of NNs. Nonetheless, experiments in Two-Stage Predict+Optimize framework [16] demonstrate that the gradient calculation method they used (which we also use) performs at least as well in post-hoc regret performance as other gradient methods, while being (significantly) faster. This is the reason we follow Hu et al.'s method and implementation.As mentioned in Section 1, while stochastic programming and Predict+Optimize are related frameworks, the technical challenges are very different. The most important difference is that Pre-dict+Optimize is a supervised learning problem, whereas stochastic programming is unsupervised learning. In Predict+Optimize frameworks, the true parameters (which need prediction) are always associated with relevant features that help prediction. On the other hand, stochastic programming frameworks have no such features, and typically assume that the entire distribution over the unknown parameters is given to the algorithm -in practice, the distribution needs to be estimated from historical data over the unknown parameters, which is an unsupervised density estimation problem.Due to the different starting assumptions, Predict+Optimize and stochastic programming formulate optimization problems rather differently. In stochastic programming, since the assumption is that the full parameter distribution is given, the optimization problem (or problems, across stages) would explicitly include the expectation operator in the objective -the goal is to solve for optimization decisions so that the expected objective, with expectation taken over the parameter distribution, is maximized/minimized. Predict+Optimize frameworks approach this rather differently: while the goal is still to optimize the expected objective, the optimization problems themselves are phrased such that they take predicted parameters, and the problem asks for the optimal decisions assuming the predicted parameters. It then becomes the goal of the learning algorithm to learn to make predictions from features, such that the expected objective is optimized overall. This is achieved via empirical risk minimization over training data, which we assume are samples from the underlying (feature,parameter) joint distribution.We also note the dimensionality of the objects being learnt in the different frameworks. In stochastic programming, the entire distribution over the unknown parameters needs to be learnt. On the other hand, in Predict+Optimize, we learn a mapping from features to predicted parameters, which, under smoothness assumptions or bounded model complexity assumptions (e.g. by restricting to using a fixed neural network architecture), can effectively be regarded as a (much) lower dimensional object than learning an entire distribution over unknown parameters.The integration of machine learning and (discrete) optimization is an increasingly popular field. In this section, we mention some works on using machine learning to help solve mixed-integer programs and especially stochastic programs.In the context of integrating machine learning with stochastic programming, several noteworthy contributions have emerged. Donti et al. [8] propose an end-to-end learning framework that directly aligns the training of probabilistic machine learning models with the ultimate task-based objective in the context of stochastic programming. Patel et al. [27] use a neural network to approximate the expected value function in two-stage stochastic programming problems, enabling more efficient solution approaches compared to traditional methods. Bae et al. [4] propose a neural network-based stagewise decomposition algorithm that can effectively approximate value functions for large-scale multistage stochastic programming problems. Rychener et al. [31] develop a principled end-to-end learning framework for training neural network decision maps that can effectively solve stochastic optimization problems, including empirical risk minimization and distributionally robust optimization formulations.The application of using ML to help mixed-integer programming (MIP) has also seen substantial progress: ML algorithms for exact solving of MIPs by branch-and-cut based algorithms [5, 11, 37] , ML algorithms for exact solving of MIPs by decomposition-based algorithms [19, 21] , ML algorithms for approximate solving MIPs by large neighborhood search based algorithm [33, 20] , and so on.In this section, we use the hospital scenario, i.e., the nurse rostering problem (NRP), mentioned in Section 1 as a running example for the Multi-Stage Predict+Optimize framework described in Section 3.1.Here we describe the NRP in detail. A hospital needs to make nurses schedule for the whole week (7 days) two weeks beforehand so that the nurses can be well prepared for the work and also plan for their leisure activities. The goal of the hospital is to minimize the total costs for hiring nurses and meet the patients' demands.There are full-time nurses in the hospital. If there are too many patients and the hospital's nurses are understaffed, the hospital can temporarily hire some extra nurses at a higher salary. Since the number of patients that will come in each shift on each day is unknown two weeks beforehand, the hospital needs to predict the number of patients to make a schedule for the full-time nurses and plan to hire extra nurses. The hospital will learn the predictor based on historical hospital records, considering features such as time of year, day of the week and temperature.To provide better service to patients, the hospital has an appointment system that requires patients to schedule an appointment in advance to receive medical care. Reservations for the next day close the night before. At this point, the hospital knows the precise number of patients for each shift of the current day. Therefore, at the night of day (t -1), i.e., Stage t (1 ≤ t ≤ 7), the true numbers of patients for each shift of the current day are revealed. Consider the time that the schedules need to be made as Stage 0. The hospital learns the predictor and uses the estimated number of patients d(0) to optimize for that week's schedule. The Stage 0 OP, the NRP using the estimations, can be formulated as:EQUATIONEQUATIONEQUATIONx i,j,k = 1, ∀i ∈ {1, . . . , n}, j ∈ {1, . . . , 7}x i,j,3 + x i,j+1,1 ≤ 1, ∀i ∈ {1, . . . , n}, j ∈ {1, . . . , 6}x i,j,4 ≤ 2, ∀i ∈ {1, . . . , n}x ∈ {0, 1}, σ ≥ 0 (6) where Equation (1) represents the objective, which is to minimize the total costs for hiring full-time nurses and extra nurses; Equation (2) ensures that the schedule can satisfy the patient demand under each shift; Equation (3) ensures that each full-time nurse will be scheduled for exactly one shift each day; Equation (4) ensures that no full-time nurse will be scheduled to work a night shift followed immediately by a morning shift; and Equation (5) ensures that each full-time nurse gets one or two day-off shifts per week.After Stage 0, the schedules for day 1 are hard commitments and cannot be changed, i.e.,x(0) 0 = {x i,1,k | ∀i ∈ {1, ..., n}, k ∈ {1, 2, 3}}, whereas the rest of the decisions are soft commitments. Example 2. (Continued) At the night of day t -1, i.e., Stage t (for 1 ≤ t ≤ 7), the reservations for the next day close, and the true numbers of patients for the three shifts of the next day θ t = (d t,1 , d t,2 , d t,3 ) ∈ N 3 are revealed. By Stage t, all the true numbers of patients for the prior t -1 days are also revealed. The number of patients for the later 7 -t days are still uncovered and are estimated as θ(t) = ( θ(t) t+1 , . . . , θ(t) T ), whereθ(t) i = ( d(t) i,1 , d(t) i,2 , d(t) i,3) ∈ N 3 represents the numbers of patients on day i estimated on day t.Hard commitments contain two parts: 1) the schedule for the day t and the prior t -1 days, and 2) the number of extra nurses hired in the prior t -1 days, i.e., here x[1 : t -1] represents {x i,j,k | ∀i ∈ {1, ..., n}, j ∈ {1, ..., t}, k ∈ {1, 2, 3}} ∪ {σ j,k | ∀j ∈ {1, ..., t -1}, k ∈ {1, 2, 3}}. The hospital may update the predictions and reschedule for the later (7 -t) days. But such rescheduled leads to extra costs for hiring full-time nurses, which are recorded by the penalty function P en(x (t-1) → x, θ[1 : t]). The more temporarily the shift is rescheduled, the larger the increase in the costs. For simplicity, we assume that the extra cost is linear in the original cost for hiring each full-time nurse. In this scenario, the penalty function can be formulated as Extra(x (t-1) → x):Extra(x (t-1) → x) = n i=1 7 j=1 3 k=1 Extra(x (t-1) → x) i,j,kwhere the (i, j, k)-th item in the penalty function is:Extra(x (t-1) → x) i,j,k = (T -j + t)ρ i c i , x i,j,k > x(t-1) i,j,k 0, x i,j,k ≤ x(t-1) i,j,kAs mentioned in Section 3.1, the Stage t optimization problem modifies the original Para-OP by: 1) introducing the penalty term capturing the cost of changing the Stage t -1 solution x(t-1) to the new Stage t solution x to the objective: Extra(x (t-1) → x) i,j,k and 2) introducing the constraint that hard commitments from prior stages cannot be modified in the current Stage t:x(t) , σ(t) =x i,j,k = x(t-1) i,j,k , ∀i ∈ {1, ..., n}, j ∈ {1, ..., t}, k ∈ {1, 2, 3}} σ j,k = σ(t-1) j,k , ∀j ∈ {1, ..., t -1}, k ∈ {1, 2, 3}Besides, for the first constraint in Equation (2), the Stage 0 predicted parameters d0 are replaced by (d The post-hoc regret used to train NN t can be written as:P Reg( θ(t) , θ[t + 1 : T ]) = obj(x (T ) , θ) -obj(x * (θ), θ) + T i=t P en i (x (i-1) → x(i) , θ[1 : i]) (7)Noting the second term is independent of θ(t) and hence NN t , the gradient with respect to an edge weight w e in NN t isd P Reg dw e = d obj(x (T ) , θ) dw e + T i=t d P en i (x (i-1) → x(i) , θ[1 : i]) dw eBy the law of total derivative, we can expand this tod P Reg dw e = d obj(x (T ) , θ) dx (T ) dx (T ) d θ(t) d θ(t) dw e + T i=t ∂ P en i ∂ x(i-1)x(i)dx (i-1) d θ(t) + ∂ P en i ∂ x(i)x(i-1)dx (i) d θ(t) d θ(t) dw eSimilar to the gradient computation in Section 4.1, the term d θ(t) dwe is calculated via standard backpropagation, while the terms d obj(x (T ) ,θ)dx (T ), ∂ P eni ∂ x(i-1) x(i) and ∂ P eni ∂ x(i) x(i-1) are calculable when the objective and penalty functions are smooth. The only non-trivial calculation is for dx (i) d θ(t) for all i ∈ [t : T ]. Recall that x(i) is computed from the Stage i optimization problem, as a deterministic function of x(i-1) and θ(t) , while x(i-1) itself also depends on θ(t) . We thus further decompose dx (i) d θ(t) into the following recursive computationdx (i) d θ(t) = ∂ x(i) ∂ x(i-1) θ(t) dx (i-1) d θ(t) + ∂ x(i) ∂ θ(t) x(i-1)Calculating ∂ x(i) ∂ x(i-1) θ(t) and ∂ x(i) ∂ θ(t) x(i-1) requires differentiating through a MILP. So instead of directly using MILP formulations for the Stage t optimization problems, we use the convex relaxation by Hu et al. [16] , which in turn adapts the approach of Mandi and Guns [22] .In this section, we give a detailed description for the other two benchmarks used in Section 5.We first demonstrate, using the example of the production and sales problem, how our framework can tackle LPs. An oil company intends to develop a production and sales plan for the upcoming T quarters/months. The goal is to maximize the profits, i.e., the sales revenues minus the production costs, under the constraints that the amount of oil product sold each quarter/month cannot exceed the customer demands. The production cost and sales price for each quarter/month are known, but the customer demand is revealed only at the beginning of each quarter/month after the company receives the orders. The company will estimate the customer demands based on historical sales records, considering features such as oil type, oil consumption of different areas, and so on.The decision variables are: 1) a real vector x ∈ R T , where x i represents the amount of product produced in month i, and 2) a real vector y ∈ R T , where y i represents the amount of product sold in month i. Let p i denote the unit profit of selling product in month i, c i denote the unit cost of producing product in month i, d i denote the customer demands in month i. The unknown parameters are d ∈ R T . At Stage 0, i.e., the time that the production and sales plan needs to be made, there is no order yet and the customer demands are unknown. The company learns the predictor and uses the predicted demands d(0) to make the plan. The Stage 0 OP can be formulated as:EQUATIONEQUATIONx ≥ 0, y ≥ 0 (11) where Equation (8) represents the objective, for maximizing the profits, i.e., the sales revenues minus the production costs; Equation (9) ensures that the amount of oil product sold each quarter/month will not exceed the customer demands; Equation (10) ensures that the amount of oil product sold each quarter/month will not exceed the inventory at that quarter/month.At the beginning of each quarter/month, the company receives orders, and the demand for that quarter/month is revealed. We assume that the beginning of each quarter/month is one stage. Then, by Stage t (1 ≤ t ≤ T ), all the true demands for the prior (t -1) quarters/months as well as the t quarter/month are revealed. The demands for the later (T -t) quarters/months are still uncovered and are estimated as θ(t) = ( θ(t) t+1 , . . . , θ(t) T ), where θ(t) i = d(t) i represents the demand of quarter/month i estimated on quarter/month t. The production and sales for the quarter/month t and the prior (t -1) quarters/months have already happened and cannot be changed, which are committed variables. There is no penalty function in this scenario. Therefore, the Stage t OP can be formulated as:x(t) , ŷ(t) = arg max x,y T i=1 p i y i - T i=1 c i x i s.t. y i ≤ d i , ∀i ∈ {1, . . . , t} y i ≤ d(t) i , ∀i ∈ {t + 1, . . . , T } y i ≤ i-1 j=1 x j - i-1 j=1 y j , ∀i ∈ {1, . . . , T } x i = x(t-1) i , y i = ŷ(t-1) i , ∀i ∈ {1, . . . , t -1} x ≥ 0, y ≥ 0 D.The methods of k-NN, RF, NN, Baseline, SCD and PCD have hyperparameters, which we tune via cross-validation: for k-NN, we try k ∈ {1, 3, 5}; for RF, we try different numbers of trees in the forest {10, 50, 100}; for NN, Baseline, and PCD, we treat the optimizer, learning rate, and epochs as hyperparameters. For Baseline, SCD and PCD, we treat the optimizer, learning rate, the early-cut-off value of log barrier regularization term (µ), and epochs as hyperparameters.Table 4 show the final hyperparameter choices for the three problems: 1) production and sales problem, 2) investment problem, and 3) nurse rostering problem. Ridge, k-NN, CART and RF are implemented using scikit-learn [29] . The neural network is implemented using PyTorch [26] . To compute the two stages of optimization at test time for our method, and to compute the optimal solution of an (MI)LP under the true parameters, we use the Gurobi MILP solver [13] .Table 5 reports the mean post-hoc regrets and standard deviations across 30 simulations for all training methods on the production and sales problem. Compared among standard regression models, NN performs well and achieves the best performance in all cases listed in Table 5 , while Ridge and RF also have decent performances.Table 6 shows improvement ratios among the proposed 3 algorithms and BAS. "A vs B" refers to the improvement in the percentage of method A over method B. Take "Baseline vs BAS" as an example, the improvement percentage of the baseline over BAS is (355.56-305.26)/355.56×100% = 14.15% when T = 4 in the low-profit price group. Comparing numbers in "SCD VS BAS", "PCD VS BAS", and "Baseline VS BAS" when stage num = 4 and 12, we can see that the advantages of the proposed methods over BAS are more distinct when the number of stages is larger. Additionally, comparing numbers in "SCD VS Baseline" and "PCD VS Baseline" when stage num = 4 and 12, we also note that the advantages of SCD and PCD over the Baseline are more distinct when the number of stages is larger. Figure 1 shows post-hoc regret comparisons between BAS and the proposed methods (Baseline, SCD, and PCD) for each run. The x-axis refers to the number of each simulation, and the y-axis refers to the ratio of the post-hoc regret achieved by BAS and the proposed methods (Baseline, SCD, and PCD) corresponding to that simulation. To easily read the comparisons, we sorted all simulations by the increasing order of the post-hoc regret ratios of BAS/SCD. The red dashed line where the post-hoc regret ratio is 1.0 represents the boundary line where (Baseline, SCD, or PCD) performs better or worse than BAS. When the point representing the post-hoc regret ratio of BAS/(Baseline, SCD, or PCD) falls above the red dashed line, it means that (Baseline, SCD, or PCD) performs better than BAS. Conversely, when the point falls below the red dashed line, it means BAS performs better than (Baseline, SCD, or PCD). Observing Figure 1 , SCD outperforms BAS across all simulations in all 4 scenarios. While not as stable as SCD, PCD and Baseline also outperform BAS in most of the simulations. Compared with Figure 1a , there are more BAS/Baseline points that fall below the red dashed line in Figure 1b , while the number of BAS/SCD points and the number of BAS/PCD points that fall below the red dashed line are similar in Figure 1a and Figure 1b . The same phenomenon can be observed when comparing Figure 1c and Figure 1d , demonstrating the advantage of SCD and PCD over Baseline.The win rate of a method against other methods can be directly computed from Figure 1 by counting the number of simulations with ratios > 1, and we provide the win rates table as a reference in Table 7 . Table 7 indicates that in more than half of the simulations, SCD achieved the best performance. PCD followed closely, achieving the best performance in the remaining 20% to 50% of the simulations. Table 8 and Table 9 report the mean post-hoc regrets and standard deviations across 30 simulations for all training methods on the investment problem. Compared among standard regression models, NN performs well and achieves the best performance in most cases, while Ridge and RF also have decent performances and obtain the smallest mean post-hoc regret in some cases.Table 10 shows improvement ratios among the proposed 3 algorithms and BAS. Comparing "SCD vs BAS", "PCD vs BAS", and "Baseline vs BAS" performance under the same capital and the same stage number, we observe that the advantages of the proposed methods (SCD, PCD, and Baseline) over the standard regression approaches become more pronounced as the transaction factor increases. Besides, comparing "SCD vs Baseline" and "PCD vs Baseline" under the same capital and the same transaction factor but different stage numbers, the advantages of SCD and PCD over Baseline are amplified as the number of stages increases. This trend is consistent with the findings from the experiments on the production and sales problem. One interesting phenomenon is that under the same capital and the same transaction factor, the advantage of the proposed methods over BAS appears to be similar or even less obvious when the number of stages is 12 compared to when it is 4. This observation differs from the pattern seen in the production and sales problem experiments. We hypothesize that this divergence may be attributable to the fundamental differences between the problem settings. While the production and sales problem is a pure LP, the investment problem is an IP with several integrality constraints. The proposed methods relax these integrality constraints and treat the problem as an LP for the purpose of forward and backward propagation. The gaps between the original IP and the relaxed LP may accumulate as the number of stages grows larger, potentially diminishing the advantages of the Predict+Optimize approaches. Table 12 reports the mean post-hoc regrets and standard deviations across 30 simulations for all training methods on the nurse rostering problem. Compared among standard regression models, NN always performs well and achieves the best performance, while Ridge and RF also have decent performances.Table 13 shows improvement ratios among the proposed 3 algorithms and BAS. Comparing "SCD vs BAS", "PCD vs BAS", and "Baseline vs BAS" performance with different extra nurse payments, we observe that the advantages of the proposed methods (SCD, PCD, and Baseline) over the standard regression approaches become more pronounced as the extra nurse payment increases.Figure 2 shows post-hoc regret comparisons between BAS and the proposed methods (Baseline, SCD, and PCD) for each run. To easily read the comparisons, we again sorted all simulations by the increasing order of the post-hoc regret ratios of BAS/SCD. Observing Figure 2 , the proposed methods outperform BAS in most of the simulations. Since the nurse rostering problem is an IP with several integrality constraints, and the proposed methods just relax these constraints and treat the problem as an LP for the purpose of forward and backward propagation. We hypothesize that the gap between the original IP and the relaxed LP may diminish the advantages of the proposed methods, and thus, BAS sometimes performs slightly better than the proposed methods.Table 14 presents the win rate results for Baseline, SCD, PCD, and BAS. Same as in the first two benchmarks, SCD demonstrates high win rates against BAS across all scenarios. SCD achieves win rates of 86.67% or higher in three out of four cases, peaking at 96.67% when the extra payment is 25. PCD also performs well, with win rates ranging from 70% to 83.33% against BAS. While not as consistently high as SCD, PCD's performance remains competitive, indicating its viability as an alternative approach. There are two stopping criteria for SCD and PCD. We set the maximum iteration number of SCD and PCD as 5. Besides, if the difference between the post-hoc regret of two consecutive iterations is less than 1, we consider that the algorithm has converged. This is the second stopping criterion.Table 15 shows the average training time across 30 simulations for the three problems. For the investment problem, since the training times under different transaction fees are similar when the capital and the number of stages are the same, we do not report them one by one but report only the average. For the nurse rostering problem, since the training times under different extra nurse payments are similar when the numbers of stages are the same, we also do not report them one by one but report only the average.Since the proposed 3 methods require solving multiple LPs when training, their training times are usually longer than standard methods. But since the production and sales problem is an LP, the solving time of which is not too long, the training time of Baseline is around double of NN. Table 15 also shows that SCD and PCD usually converge in 4 iterations in the production and sales problem.In the investment problem, the training times of Baseline are better than that of RF. The solving time of the OP, i.e., the difficulty of solving the OP, can largely affect the training times of the proposed methods. When the number of stages grows larger, the investment problem is more difficult to solve. Therefore, the training times of Baseline when there are 4 stages are quite comparable with that of NN, but the training times of Baseline when there are 12 stages are much larger than that of NN. In addition, when the OP becomes more complex, the number of iterations required for SCD and PCD convergence also increases. SCD and PCD convergence usually in 2-3 iterations when there are 4 stages, and usually in 3-4 iterations when there are 12 stages.In the NRP, since the solving time of 1 NRP is large, the training times of the proposed methods are larger than standard regression methods, which indicates that one future research direction is the speed-vs-prediction accuracy tradeoffs on Multi-Stage Predict+Optimize.As mentioned in Section 3 and Section 4, other types of training algorithms can be used within the proposed framework. Considering the performance in both post-hoc regret and training time, we proposed three training algorithms in Section 4. However, there are additional possibilities, and we will discuss two of them in this section.As mentioned in Section 3, the proposed framework allows the training algorithm to choose whether to incorporate the revealed true parameters as additional features as input. In the current implementation of the proposed SCD and PCD methods, we did not include these revealed parameters because, in our preliminary experiments, including them does not really improve prediction quality while just increasing training time.Table 16 reports the results from the experiments on the production and sales problem by incorporating the revealed parameters from stage t -1 as inputs to the stage t networks. As Table 16 shows, the performance using this expanded input did not improve over the results of the proposed three approaches. The reason that including revealed parameters does not lead to better predictive performance is due to a combination of the nature of the data and the neural network architecture. Consider one extreme: the true parameter vectors in every single stage are always equal. A neural network model (or other reasonable models) that takes in the prior-stage true parameters should be able to learn to pick up on that and use the information. Consider the other extreme: if the parameter vectors are completely independent across stages, then no model will be able to use prior-stage true parameter information, since no actual information exists. Reality (and our particular dataset) is somewhere in between: there can be some correlation between stages, but, it really depends on whether such information is extractable by the neural network architecture (or whatever other prediction model one wants to use within our multi-stage framework).For general applications, the decision whether to use the prior-stage true parameters is essentially another hyperparameter, to be tuned for that particular application using the available training data (though one can be safe and always include it, at the expense of training time).We propose two coordinate descent based methods in Section 4, but it is also possible to train all networks simultaneously, without the use of coordinate descent. In this section, we compare the proposed two methods against a simultaneous training method, to better motivate the necessity of coordinate descent.A reasonable way (and the only way we can think of) to train all networks simultaneously is replace prior and future stage predictions in SCD and PCD by the ground truth parameters. However, intuitively, this is a worse approach than the proposed methods, given the interdependency of the predictors: the performance of a predictor depends on the predictions and choices made in past and future stages. We have explored this alternative training method and here are the preliminary results on the production and sales problem in Table 17 . We observe that the solution quality achieved by the simultaneous training method was worse than the SCD and PCD methods.The coordinate descent strategy in Section 4 was needed to capture the complex interactions between the networks in different stages, and crucial for the strong performance.NeurIPS Paper Checklist 


Decision: This paper considers the problem of learning unknown parameters for multi-stage optimization problems which are gradually known at each stage. The latent parameters are predicted to compute soft decisions; the stage's unknown parameters become known; practitioner can modify the soft decision to a hard decision for some penalty; and once the hard decision is fixed, the next stage begin. The post hoc regret is computed once all the stages are done. The overall goal is to train models to predict the latent parameters by minimizing the overall regret. The paper generalizes previous end-to-end training methods for two-stage problems. 

This is a borderline paper based on the reviews and discussion. One of the reviewer who gave a initial negative review could not participate in the discussion due to emergency situation (the author response seems to mostly address their critical comments). Here are strengths and weaknesses based on my reading and assessment. 

(++) Well-motivated new problem setting

(+) The proposed training approaches are reasonable but not particularly innovative

(-) The evaluation is on synthetic problems. It would have been nice to have at least one real-world problem to demonstrate the practical utility

I'm recommending accepting the paper with the hope that this new problem setting may inspire new research by both AI researchers and practitioners to expand the knowledge of decision-focused learning paradigm.

Reviews: review: paper_summary: paper_summary: The paper proposes a new framework called Multi-Stage Predict+Optimize for tackling optimization problems with parameters revealed in multiple stages, rather than simultaneously. The authors develop three training algorithms for neural networks within this framework, particularly for mixed integer linear programs (MILPs). These algorithms include a baseline extension of prior work and two novel algorithms leveraging sequential and parallel coordinate descent methods. The paper demonstrates the efficacy of these methods through experiments on three benchmarks, showing superior learning performance over classical approaches.
strength_weakness: strength_weakness: 1. The proposed Multi-Stage Predict+Optimize framework, while an extension of existing two-stage methods, may not significantly differentiate itself from prior frameworks in practical applications. The core idea of updating parameter predictions and decisions in multiple stages is not sufficiently innovative compared to existing work in multi-stage stochastic optimization. To improve, the authors could more clearly highlight unique aspects and potential new applications of their framework that are not covered by existing methods.

2. Some sections of the paper, particularly those explaining the algorithms and theoretical foundations, may be dense and challenging for readers not deeply familiar with the subject. This can hinder the accessibility and broader understanding of the proposed methods. Simplifying the explanations and including more intuitive examples or visual aids can make the paper more accessible to a wider audience.

3. The benchmarks used for experimentation, while useful, may not fully represent the diversity of real-world applications. Additionally, the performance metrics and evaluation scenarios could be more comprehensive to cover various practical constraints and conditions. Expanding the experimental evaluation to include a wider variety of benchmarks and more complex real-world scenarios would strengthen the paper's claims.
questions: questions: Please see the weakness. As I am unfamiliar with this topic, I am unsure of my judgment and may need to discuss it with other reviewers.
limitations: limitations: NA
paper_summary: paper_summary: The authors present an approach for learning hidden parameters for multi-stage optimization problems where parameters are gradually revealed at each stage. In this setting, latent parameters are predicted, then soft committed decisions are made based on those predictions, that stage’s parameters are revealed, and the practitioner can modify the soft committed solution to a hard committed solution for a penalty. Once that hard committed solution is fixed, the next stage begins with new predictions. After all stages are done, the post hoc regret is determined by the difference between the quality of decisions made and the optimal hindsight sequence of decisions, minus incurred solution modification penalties.

The goal is to train models that predict the next time step’s parameters, such that the overall regret is minimized. To solve this problem, the authors propose generalizing previous work on end-to-end training for the same post hoc regret developed for two stage problems. Here, the authors consider unrolling the multiple stages and propose three methods of training: baseline, coordinate descent, and parallel coordinate descent. The baseline approach considers the same model to be predicting at each stage. The other two approaches consider training independent networks for the different stages. Coordinate descent trains each network individually by fixing the other network weights and caching intermediate solutions if they don’t need to be recomputed. Parallel coordinate descent trains the different networks in parallel disregarding that the gradients for a given network may not be representative of the gradients obtained given the state of the other models.

The authors run experiments on three synthetic domains motivated by real world problems and show that their approaches outperform optimization-agnostic learning methods.
strength_weakness: strength_weakness: 1 the problem settings while well motivated are synthetic and seem to be pulling data from sources that are quite unrelated to the problem domain. For instance, the ICON challenge used for nurse scheduling represents energy data. Additionally, it is unclear what the knapsack data is supposed to represent and why it would be relevant for the oil shipment or investment settings. For the investment settings in particular, previous work [1,2,3,4] has used real world public stock data to evaluate performance. Given that the contribution is mostly empirical, it would be helpful to demonstrate performance on real-world settings with relevant data.

2 The approaches themselves are somewhat incremental considering previous work. Mainly, the approaches consider unrolling the two-stage approach considered in previous work.

4 It is unclear how substantial the performance improvements are given that the variability is quite large and induce overlapping confidence intervals. It might be helpful to give win rate results to determine how often a given approach is the top-performing approach.

5 It is unclear why a coordinate descent approach is needed. Would it be possible to train all networks for all timesteps simultaneously? It would help to better motivate the necessity of coordinate descent by comparing against a simultaneous training approach.

Minor comments

Formal framework definition - it would be good to keep the revealed (and predicted) indices consistent. The off-by-one indices appear to be frequently mixed up. For instance, line 142, true parameters theta1,…, theta t-1 are revealed but then line 144 - theta t is revealed. Additionally, soft solutions \hat{x} seem to start at 0 but are indexed from 1 to t-1 in line 152. It is fairly clear from the writing what the terms are supposed to represent, but making things more cohesive would improve the paper.

Additionally, in regards to notation it might be helpful to disambiguate the hard vs soft commitments with different variables, since it seems that \hat{x} is considered as both a soft and hard committed solution which gets “overwritten”. Mainly, this just changes that the problem in 152 refers to \hat{x} as a soft solution in the objective but hard solutions in the constraints.

[1] Wang, Kai, et al. "Automatically learning compact quality-aware surrogates for optimization problems." Advances in Neural Information Processing Systems 33 (2020): 9586-9596.

[2] Ferber, Aaron, et al. "Mipaal: Mixed integer program as a layer." Proceedings of the AAAI Conference on Artificial Intelligence. Vol. 34. No. 02. 2020.

[3] Shah, Sanket, et al. "Decision-focused learning without decision-making: Learning locally optimized decision losses." Advances in Neural Information Processing Systems 35 (2022): 1320-1332.

[4] Zharmagambetov, Arman, et al. "Landscape surrogate: Learning decision losses for mathematical optimization under partial information." Advances in Neural Information Processing Systems 36 (2024).
questions: questions: 1. Do the networks at stage t take in as input the revealed parameters from stage t-1 in addition to the features? Are those already considered as features? Are the standard baselines in BAS also given the previous step information as input?
limitations: limitations: The limitations are adequately addressed
paper_summary: paper_summary: The paper proposes a Multi-Stage Predict+Optimize framework that addresses optimization problems where parameters are revealed in sequential stages. It introduces three neural network training algorithms tailored for mixed integer linear programs (MILPs) under this framework. The methodologies are empirically evaluated using three benchmarks, demonstrating enhanced predictive performance over classic methods.
strength_weakness: strength_weakness: 1. **Assumptions on Constraints**: The paper assumes that stage-wise optimizations are always feasible regardless of previous stages' outcomes and current parameter estimates. This assumption may not hold in practical scenarios, potentially limiting the framework's applicability.

2. **Limited Benchmarks**: The choice of benchmarks does not include complex problems like the Alloy Production Problem or the 0-1 knapsack[1], which could benefit from a multi-stage approach. This omission raises questions about the framework's applicability to a broader range of problems.

[1] X. Hu, J. C. H. Lee, and J. H. M. Lee. Two-Stage Predict+Optimize for mixed integer linear programs with unknown parameters in constraints. In Advances in Neural Information Processing Systems, 2023.
questions: questions: 1. **Parameter Inclusion in Predictions**: Why aren't revealed parameters included in subsequent predictions? Would incorporating these parameters enhance prediction accuracy? I hope the author will share more of your thoughts on this with us.

2. **Clarification on Terminology**: What does TOV stand for in the last columns of Tables 1 and 2? I could not find the full form of TOV in the main text.
limitations: limitations: 1. **Applicability Concerns**: Despite listing applications like nurse scheduling, the practical applicability across diverse real-world scenarios remains uncertain. More concrete examples and a broader range of applications would strengthen the claims.

2. **Innovativeness**: While the paper extends the two-stage predict+optimize framework to multi-stage, the core ideas closely resemble the existing models, which might diminish the perceived novelty.

**Additional Note on Writing Corrections**:
- Correct the notation for $\hat{x}_0$ to $\hat{x}^{(0)}$ where applicable in line 101.
- Ensure consistency in table font sizes, especially noting that Table 3’s font size is noticeably smaller than that in Tables 1 and 2.

This is an interesting field, and if the authors can provide sufficiently detailed responses to the questions raised, I would be willing to consider revising my score.
paper_summary: paper_summary: This paper addresses optimization problems with parameter values unknown at solving time. Specifically, the work takes a supervised learning approach to predicting these parameters, which are revealed over multiple “stages.” The authors extend the Predict+Optimize framework, which uses an optimization-inspired loss function rather than typical mean squared error. This paper addresses multi-stage MILP problems, where unknown parameters of the optimization problem are revealed over several sequential stages. Three algorithms are proposed for training the prediction model, including the joint training of individual neural network predictors for each stage of the problem. Finally, experiments are conducted on three benchmark problems showing the advantage of the new framework over classical training methods.
strength_weakness: strength_weakness: - While this work is very interesting, the fit to this venue is questionable. Firstly, the value of this work to the ML community is unclear; it seems the main impact is in operations research. Secondly, many crucial aspects are moved to the appendices, e.g., all of the optimization problems and most of the computational results. For these reasons, perhaps a leading operations research journal would provide more space and be more appropriate.
- While the integration of machine learning and optimization is an increasingly popular field, this paper does not reference many works. Adding references around using ML to help mixed-integer programming and especially stochastic programming, e.g., in Appendix A.3 or elsewhere, would better help readers place this work in the context. 
- The termination criteria for the algorithms in Algorithm 1 are omitted. They are also not discussed in Section 4.2-4.3, i.e., definition of convergence.
- Training the neural networks for these problems requires differentiating through a MILP, which is approximated by differentiating its convex relaxation. The authors do not discuss the quality of this approximation, such as if the strength of the convex relaxation impacts the performance. Often, multiple MILP formulations are possible, with different convex relaxations.
questions: questions: - Line 109: how common is the assumption that all the variables can be changed in Stage 1 (i.e., when does this occur)? This is what enables the guarantee of feasibility, unlike the hard commitments introduced later in the paper.
- Line 157: likewise, what are the implications of the assumption that future decisions are always feasible given any choice of hard commitments? This should be discussed as a limitation.
- Line 163: please check this. The two-stage framework in Section 2 includes no hard commitments.
limitations: limitations: Yes

Comments: title: Response to Rebuttal
Response to rebuttal
Response to Rebuttal
comment: We thank the reviewer again for the useful discussion, and we're glad that you appreciate our paper and our careful and thorough rebuttals/responses.
Thank you for your excellent response. Your explanation has addressed my concerns and clarified the rationale behind the decision regarding parameter inclusion. Although I am not deeply familiar with this specific area, I appreciate the thoroughness of your reply and the contributions you have made. Based on your clarifications, I will be raising my score.

Thank you again for your hard work and dedication to this research.
Thank you for taking the time to read and respond to our rebuttal. We are glad that our clarifications do address your concerns.
Thanks for clarifying; I have raised my score to 5 but with less confidence since I am unfamiliar with this field.
Thank you for appreciating that our work is very interesting and well-motivated. We definitely agree with you as well that the effect of the tightness of convex relaxation is important future work.
I have read the authors responses and thank the authors for their answers. I maintain my score and overall positive impression of this work. Note that it would indeed make this work more comprehensive to investigate the effect of the tightness of the convex relaxation, as this can directly affect the accuracy of the trained surrogate models.
We thank the reviewer for appreciating our rebuttal explanations, and for following up.

Re: parameter inclusion in next stage prediction. Thank you for pushing us further on this, this is worth adding discussion in the revised paper. For your question, the outcome necessarily has to be due to a combination of the nature of the data and the neural network architecture.

Consider one extreme: the true parameter vectors in every single stage are always equal. A neural network model (or other reasonable models) that takes in the prior-stage true parameters should be able to learn to pick up on that and use the information.

Consider the other extreme: if the parameter vectors are completely independent across stages, then no model will be able to use prior-stage true parameter information, since no actual information exists.

Reality is probably somewhere in between: there can be some correlation between stages, but, it really depends on whether such information is extractable by the neural network architecture (or whatever other prediction model one wants to use within our multi-stage framework). So the decision whether to use the prior-stage true parameters is essentially another hyperparameter, to be tuned per application using the available training data (though one can be safe and always include it, at the expense of training time).

In order to address this concern, in the revised manuscript we will:
- Explicitly include in our framework the possibility of using the prior-stage parameters as input to the models.
- Add the above discussion, which can help guide practitioner intuition.

Please let us know if you have any further questions! Thank you again for engaging with us.
Firstly, I would like to thank the authors for their detailed responses to my review comments. I appreciate the effort made to address the concerns raised.

Regarding the decision not to include the revealed parameters from stage $t−1$ as inputs for predictions at stage $t$, the authors mentioned that preliminary experiments showed that including these parameters did not significantly improve prediction quality but did increase training time. This result seems counterintuitive, as incorporating more accurate revealed information would typically be expected to guide and correct current predictions more effectively. Could the authors provide further analysis on why including these revealed parameters did not lead to better predictive performance? Is it possible that this outcome is due to specific experimental settings, the complexity of the model architecture, or the nature of the data itself? A deeper exploration of this phenomenon could provide valuable insights and help in further understand the model's performance.

Thank you again for your thoughtful responses. I look forward to your further clarifications.