End-to-End Probabilistic Framework for Learning with Hard Constraints
Authors: Utkarsh Utkarsh
Danielle C. Maddix
Ruijun Ma
Michael W. Mahoney
Bernie Wang
Venue: ICLR 2026 Poster (2026)
Keywords: scientific machine learning
conservation laws
physically constrained machine learning
partial differential equations
time series forecasting
uncertainty quantification
Primary area: applications to physical sciences (physics, chemistry, biology, etc.)

Abstract: We present ProbHardE2E, a probabilistic forecasting framework that incorporates hard operational/physical constraints and provides uncertainty quantification. Our methodology uses a novel differentiable probabilistic projection layer (DPPL) that can be combined with a wide range of neural network architectures. DPPL allows the model to learn the system in an end-to-end manner, compared to other approaches where constraints are satisfied either through a post-processing step or at inference. ProbHardE2E optimizes a strictly proper scoring rule, without making any distributional assumptions on the target, which enables it to obtain robust distributional estimates (in contrast to existing approaches that generally optimize likelihood-based objectives, which can be biased by their distributional assumptions and model choices); and it can incorporate a range of non-linear constraints (increasing the power of modeling and flexibility). We apply ProbHardE2E in learning partial differential equations with uncertainty estimates and to probabilistic time-series forecasting, showcasing it as a broadly applicable general framework that connects these seemingly disparate domains.  Our code is available at https://github.com/amazon-science/probharde2e.

Published as a conference paper at ICLR 2026
END-TO-END PROBABILISTIC FRAMEWORK FOR LEARNING WITH HARD CONSTRAINTS
Utkarsh" MIT CSAIL utkarsh5emit.ed edu
Danielle C. Maddix
AWS
dmmaddix wew.iiiiaaaaaaza.ac.ooooo com
Ruijun Ma
Michael W. Mahoney
Yuyang Wang
Amazon SCOT
Amazon SCOT zmahmi chwamazo. com
AWS yuyawangeamazon..oo
rui junma@ amazon. com
ABSTRACT
We present PdE2H, a probabilistic forecasting framework that incorpo- rates hard operational/physical constraints, and provides uncertainty quantification. Our methodology uses novel differentiable probabilistic projection layer (DPPL) that can be combined with wide range of neural network architectures. DPPL allows the model to learn the system in an end-to-end manner, compared to other approaches where constraints are satisfied either through a post-procsssing step or at inference. ProbHardE2E optimizes strictly proper scoring rule, without mak- ing any distributional assumptions on the target, which enables it to obtain robust distributional estimates (in contrast existing approaches that generally optimize likelihood-bate -ased objectives, which can be biased by their distributional assump- tions and model choices); and can incorporate range of non-linear constraints (increasing the power of modeling and flexibility). We apply ProbHa r dÆ2E in learning partial differential equations with uncertainty estimates and to prob- abilistic time-series forecasting, showcasing it as broadly applicable general framework that connects these seemingly disparate domains. Our code available at https: /ii thub. com/ ama zon- sci wciencceproobbie prw.hardebb.......
INTRODUCTION
Recently, machine learning (ML) models have been applied to a variety of engineering and scientific tasks, including probabilistic time series forecasting (Rangapuram al., 2021; Hyndman al., 2011; Taieb et al., 2017; Olivares et al., 2024b) and scientific applications (Krishnapriyan et al., 2021; Hansen et al., 2023: Négiar al., 2023; Mouli et al., 2024). Exact enforcement of hard constraints can be essential in domains where any violation of operational or physical requirements (e.g., coherency in hierarchical forecasting, conservation laws in physics, and non-negativity in economics and robotics) is unacceptable (Gould et al., 2022; Hansen et al., 2023; Donti et al., 2021). Limitations of data-driven MI. approaches arise in various disciplines where constraints need to be satisfied exactly (Rangapuram et al., 2021; Hansen et al., 2023). Within ML, constraints are typically incorporated as soft penalties, e.g., with a regularization term added to the loss function (Raissi et al., 2019; Li et al., 2024); but they sometimes incorporated via post-training correction mechanisms, e.g., to enforce a hard constraint (Hansen et al., 2023; Mouli et al., 2024; Cheng et al., 2025). Some methods have managed to enforce hard constraints "end-to-end" in general framework as a differentiable solver (Négiar et al., 2023; Chalapathi et al., 2024; Rackauckas et al., 2020), as a differentiable optimization layer (Amos & Kolter, 2017; Agrawal et al., 2019; Min et al., 2024). as an auxiliary procedure (Donti et al., 2021).
The aforementioned hard-constrained models typically provide point estimates without uncertainty quantification (UQ), limiting their use cases in operational and physical domains requiring probabilis- tic forecasts. Generating output distribution statistics under hard constraints often computationally
Work done during internship at AWS.
1 "Corrsponndence Danielle C. Maddix <dmmadd xx@ma com?.

---

Published as a conference paper at ICLR 2026
expensive or yields only approximate solutions (Robert al., 1999; Szechtman & Glynn, 2001; Girolami & Calderhead, 2011). There have been nomain-speciffe works in hierarchical probabilistic time series forecasting, which enforce coherency constraints using end-to-end deep learning models (Olivares et al., 2024b; Rangapuram et al., 2021). However, these works either apply only to linear constraints, or they require computationally expensive sampling procedure in training. Similar approaches have been proposed for computing probabilistic solutions to partial differential equations (PDEs) that satisfy constraints (Hansen et al., 2023; Mouli et al., 2024; Cheng et al., 2025; Gao et al., 2023; Utkarsh et al., 2025). In these works, however, the constraints are only applied as proc-psog step, and they do not lead to an end-to-end solution (of interest in the common situation that one wants to incorporate a hard-constrained model within larger model, and then optimize the larger model) that optimizes the evaluation accuracy. In both the forecasting and PDE application domains, none of this prior UQ work can handle complex nonlinear (hard) constraints.
In this work, we propose a novel probabilistic framework, FobHerdEE2 that integrates broad class of hard constraints (including non-linear constraints) in end-to-end fashion, while incorporat- ing UQ. By leveraging key results from statistics and optimization in novel way, we predict both the mean and covariance of the output data, moving beyond point estimate predictions. ProbHar Hataarddi rdE22 enforces nonlinear constraints with an efficient sampling- free method to generate distribution statis- ties. Our probabilistic approach enables the effective handling of exogenous spikes and jumps (or other discontinuities) by leveraging data heteroscedasticity. enhancing the model's robustness and flexibility under varying data conditions.
We summarize our key contributions as follows.
We introduce ProbHardE2E, as a general framework to learn a function in an end-to- end manner by optimizing an objective under hard constraints. The framework enables UQ by learning parameters of a multivariate probabilistic distribution. We show that P ProbHardEEE can incorporate broad class of deep learning backbone models.
The key technical novelty of P ProbHarde rdEEE is a differentiable probabilistic projection layer (DPPL) that extends standard projection methods to accommodate UQ while enforcing hard constraints. ProbHardE22 E can handle constraints ranging from linear equality to general nonlinear equality convex inequality constraints.
We use the DPPL to impose constraints directly on the marginals of the multivariate distri- bution for an efficient sampling- -free approach for posterior distribution estimation, which reduces the computational overhead by up to 3-5x during training.
We show that ProbHa ProbHardE2E related) tasks, where hard constraints are important: probabilistic time series forecasting; and solving challenging PDEs in scientific machine learning (SciML). We provide an extensive empirical analysis demonstrating that ProbHardE2E results in up to 15x lower mean- squared error (MSE) in mean forecast and 2.5x improved uncertainty estimates, measured by the Continuous Ranked Probability Score (CRPS), compared to the baseline methods. We show that training with the continuous-ranked probability score (CRPS), rather than negative log- likelihooo (NLL. leads better predictive performance. While the need for this well-known in, e.g., time series forecasting, previous PDE learning works commonly use NLL-based metrics for UQ
effective in two sveemingly-unreltted but technically-
2 RELATED WORK
There is a large body of related work from various communities, ranging from imposing constraints on neural networks for point estimates (Min et al., 2024; Donti et al., 2021), to probabilistic time series forecasting with constraints (Rangapuram et al., 2021; 2023; Olivares et al., 2024b), to imposing constraints on deep learning solutions to PDEs (Négiar et al., 2023; Hansen et al., 2023). Table 5 in Appendix A summarizes some advantages and disadvantages of these methods that are motivated by enforcing hard constraints in these domains. (See Appendix A for additional details.)
3 PROBHARDE2E: A UNIFIED PROBABILISTIC OPTIMIZATION FRAMEWORK
In this section, we introduce dEZ. See Algorithm for a summary. (See also Appendix for universal approximation guarantee.) In Section 3.1, we discuss the proper evaluation metric
2

---

Published as a conference paper at ICLR 2026
Algorithm 1 ProbHardE2E: Training and Inference
Require: Training data {(6(1) =)+)) ~ D, test data & and constraints g(-) - 0, h(-) = 0.
Ensure: Learnable function fs & + )' that outputs constrained distribution parameters.
Pick a model class θ, initialize weights A O for probabilistic unconstrained model fo : - Z.
2: while 0 not converged do
3: Predict unconstrained distribution parameters (n (c(()) (d()))
4: Training Mode: Project parameters ()o()())) (x(())) (PPL
5: Update 0 E θ by minimizing the CRPS loss { (()((())) u (+))
6: end while
7: Inference Mode: Project random variable Y (Φ) = ( (0), g()). h(.)), where Zœ( and
Yo(ő) denote the unconstrained and unconstrained random variables, respectively.
8: Return Feasible predicted sample u (z)())) ~ Yo(ő), where ze z0(o) Zo(
for a constrained probabilistic learner, and we define our objective function that corresponds that evaluation metric. In Section 3.2, we propose our differentiable probabilistic projection layer (DPPL) that enforces the hard constraints. In Section 3.3, we describe how compute the parameters of the resulting constrained posterior distribution. In Section 3.4, we discuss update rules for various types of constraints (linear equality, nonlinear equality, and convex inequality constraints). In Section 3.5, we propose sample-free formulation for satisfying the constraints while optimizing for the objective.
3.1 PROBABILISTIC EVALUATION METRICS AND OBJECTIVE FUNCTION
We formulate the problem of probabilistic learning under constraints. The goal of this problem is to learn a function fa a V, where & C II'" denotes the input space, 0 E θ C θ denotes the feasible parameter space, and y C R* denotes the space of predicted distribution parameters that meet the constraints. Given multivariate distribution class, these learned parameters induce a predictive multivariate random variable Y (o())), where (g5(9) ()) D, where g(+) T u(0) R and D denotes training data from a distribution D. Each realization of a(p(1) ~ Yo(d() E R" is and h(ā()(0)) 0. We
required satisfy predefined hard constraints of the form g (()()))
can formulate this constrained optimization problem as follows:
arg min eee, a(Yald 5o, ((Yo()())== (d00) E (Yo(o(0),
(I)
where denotes a proper scoring rule.
One widely-used (strictly) proper scoring rule for continuous distributions is the continuous ranked probability score (CRPS) (Gneiting & Raftery, 2007). The CRPS simultaneously evaluates sharpness (how concentrated or "narrow" the distribution is) and calibration (how well the distributional coverage "aligns" with actual observations). More formally, for an observed scalar outcome y and a corresponding probabilistic distributional estimate, Y. the CRPS is defined as:
CRPS(Y. Ey yl - Ey Y'I.
where Y' denotes an i..dd copy of Y. Compared to other scoring rules, e.g., the log probability scoring rules, which require assumptions on the outcome variable, the CRPS is robust to probabilistic model mis-specification. Because of this unique property, the CRPS is widely used the evaluation metric in many applications, c.g., probabilistic time series forecasting (Gasthaus et al., 2019; Rangapuram al., 2021; Park et al., 2022: Olivares et al., 2024b), quantile regression (Fakoor et al., 2023), precipitation nowcasting (Ravuri et al., 2021; Gao et al., 2023) and weather forecasting (Rasp & Lerch, 2018; Kochkov et al., 2024; Price et al., 2025).
We align our training objective with the proposed evaluation metric above. by directly optimizing the CRPS Eq. (2) in Problem We define the loss as the sum of the univariate CRPS:
n
f( (Yo (Y))
Σ CRPS (Ye)(( (.)).
(3)
3

---

Published as a conference paper at ICLR 2026
The CRPS naturally aligns with the goal of producing feasible and well-calibrated predictions, as the metric rewards distributions that closely match observed outcomes. Enforcing our constraints in the distribution space guarantees that every sample from the predicted distribution is physically or operationally valid. Consequently, modeling the loss through the CRPS provides a principled way to reconcile domain constraints with distributional accuracy.
3.2 DIFFERENTIABLE PROBABILISTIC PROJECTION LAYER (DPPL)
We transform the constrained Problem into the unconstrained optimization problem:
arg min (((t)) l (Yo( u,
dc
where θ C θ denotes the feasible parameter space that ensures constraint satisfaction, and É denotes the loss function in Eq. (3). We solve this using two-step procedure: first define a predictive output distribution, then project it onto the constraint manifold using differentiable probabilistic projection layer (DPPL) for end-to-end optimization.
Our framework begins with established probabilistic backbone model. This can be Gaussian Process (Rasmussen & Williams, 2006), neural process (Kim et al., 2019), Dee (Salinas et al., 2019; Rangapuram et al., 2021), ensembles of neural networks operators (Mouli et al., 2024). The base model fa 4 - F predicts the distribution parameters (mean μe (o()) and covariance Ne (o()),
), for 0 E )) without constraint awareness. We then use reparameterization function r: R& x I + IR' define the distribution one of two ways: either as an identity map, where 1= k, that returns f (d6))), Ee(()() for our efficient sample-free paradigm during training; or as map, where n, that combines the distribution parameters with noise & ~ p()) R", where p denotes tractable sampling distribution, and gives sample ze =o(()))) Zo ()())
R=" from the predicted distribution to generate constrained samples at inference. This dual-mode design balances training efficiency with strict constraint feasibility at inference.
The reparameterization function induces the base (unconstrained) distribution parameters or predictive random variable as:
(Ho( Ee (d()))) (Training) (()(1))
r (fo(
(Inference)
Following this Predictor Step above, we use the DPPL in the Corrector Step to restrict the parameter space é S θ, such that for all ()())))) ()())) = satisfied. The DPPL is our core architecture innovation for leveraging the base model to learn predictions that satisfy the given constraints. We define the projected distribution
(o))),) the constraints g(ûg ((()) I and
parameters or projected predictive random variable as:
r(fo((())) (Mo(o())) L. (o(0)), (Training) -( (Yo(o()),
(f)((((((()))))))))))))
(6)
(Inference)
for (of {) in Eq. (5), where fo Φ + y CRC R* denotes the probabilistic model that outputs the constrained distribution parameters (o()), (o((i)) Our DPPL yields a constraint-satisfying realization u² Y Yo(g(()) as the final predictive random variable.
This two-step approach mirrors predictor-corrector methods (Boyd & Vandenberghe, 2004; Bertsekas, 1997). with the DPPL serving as our key architectural innovation for ensuring constraint satisfaction. Equivalently, the DPPL can be formulated as a constrained least squares problem on the samples of Zo((()). ). (See Appendix C for details.) Prior works on imposing hard constraints in time series and solving PDEs (Rangapuram al., 2021; Hansen et al., 2023) reduce special cases of our method with linear constraints. (See Appendix D for details.) We draw (n((()) (o( and solve the following constrained optimization problem:
u (ze (c((()))))
ay(0())ER* (x )) (OM(((()))))) (م) zo (0
arg min
where u² (20(<)()) denotes predicted sample of Y)() and where ||~l|Q v2 Qz for some symmetric positive semi--efinite matrix Q. (See Appendix E for details on the flexibility of learning various forms of Q.)
4

---

Published as a conference paper at ICLR 2026
3.3 DPPL ON THE DISTRIBUTION PARAMETERS FOR LOCATION- -SALLE DISTRIBUTIONS
In this subsection, we detail how to directly compute the parameters for the constrained distribution by applying our DPPL on the base distribution parameters for an efficient, sampling- free during training. To do so, we can assume that the prior distribution JF belongs to multivariate, location- sccllioon-caleeee family, i..., a distribution such that any affine transformation Y of a random variable Z + 1//2 ~ F(u, 2)) and & ~ F(0, 1), remains within the same distribution family F. This an example of how to compute the random variable in Eq. (5) for a multivariate location-scale distribution. A familiar case of this when Z (p, Σ) and Y AZ B is an affine transformation; in which case Y ~ (Ap + B, ). Alternatively, we can show that when Y is nonlinear transformation of Z, it has approximately (to first-order) the same distribution Z, with an appropriately-chone of parameters (given in Eq. (8) below). We state this result more formally in Theorem 3.1. The proof, given in Appendix F, uses a first-order Taylor expansion to linearize the nonlinear function transformation, and is similar to the Multivariate Delta Method (Casella & Berger, 2001).
set
Theorem 3.1. Let Z ~ F(p, Σ) be random variable, where the underlying distribution F belongs to a multivariate location-scae scale family of distributions, with mean p and covariance 2; and let T be function with continuous first derivatives, such that JT (T))EiJT(1) is symmetric positive semi definiee Then, the transformed distribution Y = (Z)) converges in distribution with first- order accuracy to F(p, Σ) with mean ğa = T(u) and covariance (T)))))T( where (M) = (T( denotes the Jacobian of with respect to evaluated at μ.
Let Z ~ F(u, Σ) denote the prior distribution and z ~ Z. We apply Theorem 3.1 with =(=)=u" = ==²
where (z) denotes the solution of the constrained least squares problem in Problem (7). In this case, the projected random variable satisfies Y ~ F(t,, Σ) with updated parameters:
j = T());
E = JT() E
(8)
3.4 DPPL FOR VARIOUS CONSTRAINT TYPES
In this subsection, we discuss how to compute the DPPL for various constraint types (linear equality, nonlinear equality, and convex inequality) for both train and inference modes. Table shows these constraints types require different treatments: linear equality have closed-fore projections, nonlinear equality can be solved with iterative methods, and convex inequality require optimization solvers.
Table I: Summary of DPPL in ProbHardd2E for various constraint types. For linear equality constraints, the oblique projection Po-1 = 1 -Q-1AT 1Q--ATAAQ--AAAA (AQ-IAT)- A; for nonlinear equality constraints, R denotes the first-order optimality conditions.
Constraint Type
Solution u' u2)
Solver Type
Jacobian JT
Linear Equality
Pg 12+ (1 Po 1)A'b
closed-form
Po
Nonlinear Equality
w. s.t. R(u" x=;2) =
nonlinear
implicit differentiation sensitivity analysis; argmin differentiation
Convex Inequality
argmin || - 1/12
convex opt.
g(ū)
3.4.1 LINEAR EQUALITY CONSTRAINTS
For linear equality constraints, we have an underdetermined linear system h(ū) Aá b 0, where A E R . < and has full row rank g. In this case, we can derive a closed-form solution to the constrained least squares Problem (7). this case, both training and inference modes are equivalent since the DPPL projection is exact. (See Appendix C.
3.4.2 NONLINEAR EQUALITY CONSTRAINTS
For nonlinear equality constraints, h(ū) 0, we longer derive the exact closed-form expres- sion for the solution. Instead, we can provide an expression which is satisfied by the optimal solution. In particular, approximate the parameter-level projection training time. This can then be
5

---

Published as a conference paper at ICLR 2026
solved for the posterior mean ====(() (p) in Eq. (8) with the nonlinear transformation T(₁) u (1) with iterative optimization methods, e.g., Newton's 's Method. (We can then compute the posterior covariance Σ in Eq. (8) by estimating the Jacobian (r) by differentiating the nonlinear equations via the implicit function theorem (Blondel et al., 2022), and evaluating it at μ.. (See Appendix C.2.) At inference, we project each sample exactly with our custom, batched optimization solver to ensure strict constraint feasibility.
(²()) = - Q Q-1@@i Vh(h" (2)) (2)) X, h=)" = O with respect
3.4.3 (NONLINEAR) CONVEX INEQUALITY CONSTRAINTS
For convex incquality constraints, á in Problem (7) is in a convex set, C C R". Closed-for formm expressions (such as those in previous subsections for linear and nonlinear equality constraints) do not exist (Boyd & Vandenberghe, 2004). Instead, we rely on convex optimization solvers ensure computational efficiency and scalability to compute the solution u' in training. The gradients of the convex program can be calculated efficiently using sensitivity analysis (Bertsekas, 1997; Bonnans & Shapiro, 2013), argmin differentiation (Sun et al., 2022: Agrawal et al., 2019; Amos & Kolter, 2017; Gould et al., 2016), and/or variational analysis (Rockafellar & Wets, 2009). These techniques provide a means to compute the Jacobian J (1), which represents the sensitivity of the optimal solution u² to changes in the input vector 1u, whose projection we are essentially computing to the convex constraints space. During inference, we solve the convex program per sample. (See Appendix C.3.)
3.5 SAMPLE-FREE WITH CLOSED-FORM CRPS
We use a closed-form expression for the CRPS to enable a computationally efficient and sample-free approach for evaluating the CRPS in the loss function in Eq. (3). Calculating the CRPS for arbi-
trary distribution requires generating samples (Rangapuram et al., 2021: Gneiting & Raftery, 2007), but closed-forn expressions for the CRPS exist for several location-scale distributions (Gaussian, lo- gistic, student's t, beta, gamma, uniform). Most notably, for the univariate Gaussian, the closed-form
CRPS is given as: (R) = 2 (2P(2) - 1) + 2p(2) } where p(2) 1 exp ( (-2² 22)
v2x
denotes the standard normal probability density function (PDF), and P(z) - dyy)dy denotes the standard normal cumulative distribution function (CDF) for 2 ~ (0, 1) (Gneiting et al., 2005;
Taillardat et al., 2016). This sample-free formulation is especially beneficial when the DPPL is computationally intensive, e.g., in the presence of nonlinear constraints.
4 EMPIRICAL RESULTS
In our empirical evaluations, we aim to answer the following five questions about ProbHardE25
(Q1) Does training end-to-end with the imposed hard constraints improve upon the performance of imposing them only at inference time?
(Q2) Is using a general oblique projection more beneficial than using the commonly-used orthog- onal projection, and if so when?
(Q3) Does training with the catributiio-agnosste proper scoring rule, CRPS instead of NLL, improve performance?
(Q4) What are the computational savings of projecting directly on the distribution parameters and using the closed form CRPS vs. projecting the samples?
(Q5) How does P I rdE2E perform when extended to more general constraints, e.g., non- linear equality and convex inequality constraints?
See Appendix G for details on the test datasets, Appendix H for implementation details, and Ap- pendix for additional empirical results.
Test Cases. We demonstrate the efficacy of ProbHardE2E in two constrained optimization applications: PDEs; and hierarchical forecasting. We show that our methodology with DPPL is model-agnostic, as demonstrated through its high-performance integration with different base models across applications. We first consider a series of PDE problems with varying levels of difficulty in learning their solutions, following the empirical evaluation from Hansen et al. (2023). These PDEs are categorized as "easy." "medium," and "hard," with the difficulty level determined by the smoothness or sharpness of the solution. (See Appendix G. for details.) In addition to PDEs,

---

Published as a conference paper at ICLR 2026
also evaluate ProbHardE2E five hierarchical time-series forecasting benchmark datasets from Alexandrov et al. (2019), where the goal is to generate probabilistic predictions that are coherent with known aggregation constraints across cross-sectional hierarchies (Rangapuram et al., 2021). (See Appendix G.2 for details.)
Baselines. We compare two variants of ProbHardE2E, i.c., ProbHardE2E-OOa which uses a general oblique projection (Q= ") projection and our default unless otherwise specified, and ProbHa NobHardEEE-OP which uses an orthogonal projection (Q 7), against several probabilistic deep learning baselines commonly used for uncertainty quantification in constrained PDEs and probabilistic time series forecasting. For PDEs, ProbHardE2 uses VarianceNO (Mouli et al., 2024), which is a probabilisitic extension of the Fourier Neural Operator (FNO) (Li et al., 2021) as the unconstrained model. We compare ProbHardE2E with: (i) HardC, which is based on Négiar et al. (2023); Hansen et al. (2023), and which imposes the orthogonal projection only on the mean, but does not update the covariance; (ii) ProbConserv (Hansen al., 2023), which applies the oblique projection only at inference time, and works only with linear constraints (in the nonlinear constraint case, we compare with ProbHa rdInf, which is variant of P ProCConserv that imposes the nonlinear constraint at inference time only); (iii) SoftC (Hansen et al., 2023), which introduces soft penalty on constraint violation la PINNs (Raissi al., 2019; Li al., 2024) during training but does not guarantee constraint satisfaction at inference; and (iv) the unconstrained model backbone Vari anceNo. For hierarchical time-series forecasting, ProbHardE25 uses DeepVAR (Salinas et al., 2019) as the probabilistic base model. We compare P ProbHa rdE2 with: (i) ProbConserv; (ii) HierE2E (Rangapuram et al., 2021), which enforces linear constraints via an end-to-end orthogonal projection; classical statistical approaches including (iii) ARIMA-NO - veBU, (iv) ETS Pers-Naii NaiveBU (Hyndman et al., 201 I; 2025). (v) PERMBU-MINT (Taieb et al., 2017; Olivares et al., 2022); and (vi) the unconstrained model backbone DeepVAR.
Evaluation. We evaluate ProbHa adE22 using the following metrics: Mean Squared Error (MSE), which measures the mean prediction accuracy; Constraint Error (CF), which measures the constraint errors on the samples (conservation law for PDEs and coherency for hierarchical time series forecast- ing); and Continuous Ranked Probability Score (CRPS), which measures performance in uncertainty quantification (UQ). (See Appendix H.3 for details on the metrics.) For each model, we report these metrics when trained with either CRPS or Negative Log-Likelihood (NLL) as the loss. Although originally optimized with NLL, we also train a CRPS-based variant of P ProbConserv to ensure fair comparison. The experiments are conducted on single NVIDIA V100 GPU in the PDEs case, and on an Intel(R) XR) CPU E5-2603 v4 @ 1.70GHz the time series forecasting case. To ensure scalability, use diagonal covariance matrix Q in Problem 7, following prior work (Hansen et al., 2023; Mouli et al., 2024). (See Appendix E for low-rank and full covariance structures.)
4.1 LINEAR CONSERVATION AND HIERARCHICAL CONSTRAINTS
In this subsection, we test P ProbHardE2E on linear constraints. Table 2 presents our comparative evaluation results across multiple PDE datasets under linear constraints, and Table presents our evaluation results across multiple time series forecasting datasets. We use these results to answer questions (Q1)-(Q4) raised above.
Q1. Regarding (QI) on the benefits of training end-to-end with a hard constraint, the results demon- strate that our method achieves superior performance compared to existing approaches. Specifically, when measured against two accuracy metrics across four PDE datasets in Table 2, our method with either oblique [ProbHardE2E-Ob or orthogonal (ProbHardE2E-Or projection consis- tently outperforms both ProbConserv, which applies constraints only during post-processing, and soffc, which implements constraints as soft penalties, as measured by the desired CRPS metric. This performance advantage directly addresses research question (Q1), showing that our end-to-end approach is more effective than methods that treat constraints as either post-processing steps or soft penalties. In addition, across five diverse hierarchical time-series datasets, our method achieves -fate-of-the-art CRPS on the LABOUR, TOURISM-L. and WIKI datasets. On the TOURISM and TRAFFIC datasets, it remains highly competitive, outperforming traditional approaches, e.g., ARIMA and ETS, and offering comparable performance to specialized methods, e.g., PERMBU-MINT and Hi er.2..
7

---

Published as a conference paper at ICLR 2026
Table 2: Test metries on constrained PDEs across four datasets, which are ordered top to bottom in their learning difficulty. Metrics include MSE x x10-6, constraint (conservation) error (CE) x 110 and CRPS x10-3. -3 Each algorithm is trained with either CRPS or NLL. Best values per row are highlighted in bold.
Datast
Motric
#robitardu2E-00 NLL
Probite.ddz2E-00 Hardc
ProbConsezv CRPS
SoftC
VarianceNu (basc)
CRPS
CRPS
NLL
CRPS
NLL
NLL
CRPS
NLL
CRPS
NLL
0.047
0.301
0.027
1.26
0.156
0.029
2,01
Heat
CE CRPS
0
0
0.852
4.806
1.76
34,3
0,304
0.37
0.271
0.713
0.275
0.452
0.392
4.27
0.354
1.129
0.396
4.39
MSE
9.59
6.16
9.01
11.08
8,870
10.55
8.501
10.5
8.187
7.362
7,945
8.132
PME
17.091
29.31
20.19
27.2
CRPS
2.01
2.65
1.798
1.80
1.785
1.667
2.03
2.49
2.068
2,444
2.02
2,43
MSE
131
262
88.09
$10.82
103.78
458.38
134
277
148.11
599.11
149
603
Advection
182
CRPS
4.19
7.03
2.94
8.669
3,236
11.23
3.88
7,00
3,963
9.702
3.98
10.1
MSB
186
207
194, 84
432.92
394.29
433.28
903
273
431.89 89
429.06
425
425
Stefan
CE
0
0
166.93
168.75
180
169
CRPS
7.52
14.147
14.67
14:432
14.539
9.878
Q2. Regarding (Q2) on the effectiveness of the oblique vs. orthogonal projections, Table 2 shows while both oblique (BrobHarddEEE-O Ob) and orthogonal (ProbHardB2E--O0 variants of ProbHardE2E enforce zero constraint error, their impact on predictive fidelity varies signifi- cantly, depending on the difficulty of the PDE problem. For the "easy" (smooth) Heat equation and "medium" tasks (PME and Advection), orthogonal projection reduces CRPS by 10 - 30% relative to oblique projection and improves MSE by up to 33%. However, in the "hard" (sharp) nonlinear Stefan problem, the oblique projection-based method improves CRPS by more than 50% compared the orthogonal projection. Table 3 shows that Pr ProbHardE2E-- series datasets with pross-sectional hierarchies, as it improves CRPS on LABOUR and TOURISM datasets, compared to P rdE2E Ob. This addresses (Q2) that correcting predictions along covariance- weighted (oblique) directions better preserves the true uncertainty structure around shocks and spikes, performing more effectively on problems with heteroscedastic data.
generally performs better on the time
Table 3: CRPS x x10-3
3 for hierarchical time-series datasets across various probabilistic forecasting algorithms. Constraint (coherency) error (CE) is given in parenthesis and is equal to O for all methods except the base unconstrained DeepVAR. PERMBU CERMBU-MINT is not available for TOURISM-L, because the dataset has a nested hicrarchical structure, and PERMBU U -INNT is not well-defined on this type of dataset (Rangapuram et al., 2021). ANIMA-Na i BU, ETS-Na vBBU and PERMBU-MINT are deterministic models with no model uncertainty.
Durset
Pedbdaa oz
nisceto
An120 Suioweu
E2S 275-suiveeu
Diupyvs (boac)
LAROUR
M611 (0)
18.6666 (0)
458+555(0)
50.5+20.600) (0)
45
43.2
(0)
38.2+4.5 (0.215)
TOURISM
989:130 (W)
$2.4866 (0)
100.7:7.7700
108.1:16.3300 (0)
113.8 (0)
100.5 (0)
77.1
92.5=2.2 (2818.01)
TOURISM-1.
155.243.6 (W)
1564400.400)
1769-21.5
161.361109.00
174.3 (0)
169.0
158 1.8.14100 (70000)
TRaNFC
55.0:1066
00.5278
7100339900
41.817.8800
(0.8
(6,5
63.7
40,0:206 (0.)92)
WiKI
282,1 29/4
21588 16.9
264.7
30.7 (0)
216,5 30)
377.3.00
4673 (0)
281.2
209 15 (8396 6)
Q3. Regarding (Q3) on the training objective, Table 2 shows that training with the proper scoring rule, such as CRPS, improves UQ (measured by CRPS) across nearly all PDE datasets compared to training with the commonly-used NLL in previous SciML works. The only exception is HardC in the PME. The CRPS training objective also improves mean accuracy (measured by MSE) in approximately three-quarters of the dataset-model experiments. In addition, Table 3 shows that on four out of five time series datasets, we improve upon the results of HierE2E. which uses the DeepVAR base model with orthogonal projection on the samples, and which optimizes the sampling- bassd quantile loss by projecting directly on the distribution parameters.
Q4. Regarding (Q4) on the computational efficiency of our sampling-free approach, Fig. ((a) shows the training time per epoch for the hierarchical time series datasets. Models trained for time series and PDEs (see Appendix I.1) with 100 posterior samples per training step incur a 3.3-4.6x increase in epoch time relative to P FrobHardE2E, which avoids sampling altogether by using closed-fore foomm CRPS loss. Note that the computational overhead of azobHardE22 is approximately 2x that of the unconstrained model. However, compared to the sampling-basee probabilistic baselines, our
8

---

Published as a conference paper at ICLR 2026
approach with the DPPL layer that directly projects distribution parameters and a closed-form CRPS objective offers significant training- time speed-ups across all forecasting and PDE testbeds.
4.2 GENERAL NONLINEAR EQUALITY AND CONVEX INEQUALITY CONSTRAINTS (Q5)
In this subsection, test ProbHardE2 to address question (QS) on PDE datasets, as (current) time series forecasting applications usually need predictions satisfy only linear (e.g., hierarchical) constraints.
nonlinear equality and convex inequality constraints
4.2.1 NONLINEAR EQUALITY CONSTRAINTS
We test ProbHa rdE2E with general nonlinear constraints using nonlinear conservation laws from PDEs. (See Appendix 1.2 for details.) Importantly, Table shows that even in this challenging case of nonlinear constraints, the constraint error (CE) the samples is 0 that ensure strict feasibility on the samples. addition, we see superior performance enforcing nonlinear constraints with Pr ProbHardE2E. We see even larger MSE performance improvement of at most 2 15 17x when trained on CRPS, and CRPS performance improvement of at most = 2.5 for m E [2, 3] over the various baselines that apply the nonlinear constraint just inference time rProbHardIrs as a reduced linear constraint (P robConserv) at inference time, and the unconstrained model (VarianceNO), These results highlight the benefit of training end-to-end with the constraint in the nonlinear case. In particular, this validates our dual mode training and inference approach, where we obtain the computational benefits of projecting on the parameters and approximately satisfying the constraint during training, and exact constraint sanctification at inference time while achieving state-of-the-ard performance measured in CRPS. In addition, Fig. ((b) shows that P ProbilardE2) is able to significantly better capture the shock and has tighter uncertainty estimates, leading to lower CRPS values than the baselines.
83-
FootbatE2E
0.6
Tue
15-
FreaCaneel =30 =3c
Prosteewn :le Prot/eidEEE
Eime per terrattonn
0.4
FropHa eãc
33-
00x
0.5
labour
Tourio
Trxfx
Barlens
0.46
0.50
0.32
0.36
0.38
0.55
0.60
0.65
0.70
0.75
0.80
Time Series Dalanet
(a) Linear Equality: Timing
(b) Nonlinear Equality: PME
(c) Convex Inequality: TVD
Figure 1: ProbHardE2E on the various constraint types. (a) Linear Equality: Average time per iteration (in seconds) for ProbHa compared to the HierE2E
five hierarchical time- series datasets; (b) Nonlinear Equality: Mean +3 standard deviation for the PME with conservation constraint time 0.51, with PDE parameter mimin [3, 4] and Mest 3.88; (c) Convex Inequality: Mean #3 standard deviation for linear advection with TVD constraint at time = 0.51, with PDE parameter Btrinn E [1, 2] and Beest 1.5. The horizontal axes in (b)-(c) are zoomed in to highlight the uncertainty near the propagating front.
4.2.2 (NONLINEAR) CONVEX INEQUALITY CONSTRAINTS
We impose convex relaxation of the total variation diminishing (TVD) constraint that has been commonly used in numerical methods for PDEs to minimize artificial oscillations VLeve 2002).
In particular, we impose TVD ΣN= Σ,", (u(tz,z44)) - u(t.,; x:)) as a regularization term. (See Appendix 1.3 for details.) Note that this discrete form is analogous to total variation denois- ing in signal processing (Rudin et al., 1992; Boyd & Vandenberghe, 2004). Fig. I(c) illustrates that imposing this TVD relaxation improves the shock location prediction, compared to the uncon- strained model VarianceNO. We see that RobHardE225 has smaller variance, compared to both ProbConserv, which only enforces the conservation law, and Variances inceen. Most importantly, see that ProbHa F2EE is less likely to predict non-physical negative samples, which violates

---

Published as a conference paper at ICLR 2026
Table 4: Test metrics on the nonlinear PME with PDE coefficient k(u) um, which controls the sharpness of the solution (larger values are "harder"), for NLL and CRPS training. The training and test parameters are sampled from this range of m. Metrics include MSE x10- calibration error (CE) x 10 3 and CRPS 10 Best values per row and metric are bolded.
PME Dataset
Metric
ProbEardi EdE2E-BB
ProbifardE2S Or CRPS
ProbffaddInE
ProBConserv CRPS
Varlancetio (hase)
CRPS
NLL
NLL
CRPS
NLL
NLL
CRPS
NLL
MSE:
5.04
106.638
9.38
43.5
78.185
86.147
88.539
94.467
80.342
89.212
m E [2,3]
CE
0
0
0
0.020
0.028
CRPS
8.648
18.867
11.34
14.8
19.005
32,977
20.672
36.724
20.779
37.140
MSE
296.4
471.3
3.19
134.7
157.8
200.2
184.5
276.4
162
201.5
[3
CRPS
11.23
16.9
7.10
11.18
22.60
30.4
24.7
35.1
23.7
48.5
424.8
716.8
1.59
206.49
280.4
332.3
276.7
619.9
199.2
341.7
m E [4,5]
CF
0
0
0
0
0
22.8
$9.7
CRPS
19,9
9.36
27.2
the positivity of the true solution. In addition, the variance of the ProbHardE22 solution also has smaller peak above the shock, and hence it is less prone to oscillations than the other baselines.
5 CONCLUSION
In this work, we propose ProbHardE2E, which seamlessly incorporates constraints ranging from commonly used linear constraints challenging nonlinear constraints into black-box probabilistic deep learning models using differentiable probabilistic projection layer (DPPL). We show that ProbHarsEE2E applicable in a wide range of scientific and operational domains. ranging from linear coherency constraints in time series forecasting to nonlinear conservation laws in solving PDEs. Contrary to prior works (Hansen al., 2023), which only support linear global conservation constraints, Pr ProbHardEa dEE2 supports nonlinear conservation constraints. This paves the way for future work to enforce conservation locally over sub-regions or over control volumes la finite volume methods (Le Veque, 2002). Future work also includes extending method to handle general non-convex, nonlinear inequality constraints using advanced optimization techniques or relaxations, to richer covariance parameterizations,, e.g., low-rank or dense, and to empirical distributions other than location-scale families by sample projection.
ACKNOWLEDGMENTS
The authors thank Jiayao Zhang, Pedro Eduardo Mercado Lopez and Gaurav Gupta for discussions on this work and related early investigations. The authors would additionally like to thank Chris Rackauckas and Alan Edelman for their insights and support of the project.
REPRODUCIRILITY STATEMENT
Code used to produce the experiments in this paper is available amazon-
10

---

Published as a conference paper at ICLR 2026
REFERENCES
Akshay Agrawal, Brandon Amos, Shane Barratt, Stephen Boyd, Steven Diamond, and J Zico Kolter. Differentiable convex optimization layers. In Advances in Neural Information Processing Systems, volume 32, 2019.
Alexander Alexandrov, Konstantinos Benidis, Michael Bohlke-Schneider, Valentin Flunkert, Jan Gasthaus, Tim Januschowski, Danielle C. Maddix, Syama Rangapuram, David Salinas, Jasper Schulz, et al. Gluonts: Probabilistic and neural time series modeling in python. Journal of Machine Learning Research, 21(116):1-6 2019.
Brandon Amos and J Zico Kolter. Optnet: Differentiable optimization as layer in neural networks. In Proceedings of the 34th International Conference on Machine Learning, volume 70, pp. 136-145. PMLR. 2017.
Oren Anava, Vitaly Kuznetsov, and (Google Inc. Sponsorship). Web traffic time series forecasting, forecast future traffic to wikipedia pages. Kaggle Competition, 2018. URL https: /www.
wwreeeeeeeeeeeeeeeeeeeeeeeeeeeee forecasting/.
Abdul Fatir Ansari, Lorenzo Stella, Caner Turkmen, Xiyuan Zhang, Pedro Mercado, Huibin Shen, Oleksandr Shchur, Syama Sundar Rangapuram, Sebastian Pineda Arango, Shubham Kapoor, et al. Chronos: Learning the language of time series. Transactions on Machine Learning Research, 2024.
George Athanasopoulos, RobJ Hyndman, Nikolaos Kourentzes, and Anastasios Panagiotelis. Forecast reconciliation: review. International Journal of Forecasting, 40(2):430-456, 2024.
Australian Bureau of Statistics. Labour Force, Australia, Accessed Online, 2019, URL
https: / www. abs gov. www.A///TT///aaaa....a nst/DetailsPage 6202. 0Decs 202019922 2o11900ppDDDuueet.t.....
Peter Battaglia, Jessica B Hamrick, Victor Bapst, Alvaro Sanzhez- Vinicius Zambaldi, Mateusz Malinowski, Andrea Tacchetti, David Raposo, Adam Santoro, Ryan Faulkner, al. Relational inductive biases, deep learning, and graph networks. arXiv preprint arXiv:1806.01261, 2018.
Souhaib Ben Taieb and Bonsoo Koo. Regularized regression for hierarchical forecasting without unbiasedness conditions. Proceedings of the 25th ACM SIGKDD international conference on knowledge discovery & data mining, pp. 1337-1347, 2019.
Konstantinos Benidis, Syama Sundar Rangapuram, Valentin Flunkert, Yuyang Wang, Danielle Maddix, Caner Turkmen, Jan Gasthaus, Michael
David Salinas, Lorenzo Stella, François-Xave Xavier Aubet, Laurent Callot, and Tim Januschowski. Deep learning for time series forecasting: Tutorial and literature survey. ACM Comput. Surv., 55(6), 2022.
Dimitri P Bertsekas. Nonlinear programming. Journal of the Operational Research Society, 48(3):
334-334, 1997.
Tom Beucler, Michael Pritchard, Stephan Rasp, Jordan Ott, Pierre Baldi, and Pierre Gentine. Enforc- ing analytic constraints in neural networks emulating physical systems. Physical Review Letters, 126(9):098302, 2021.
Lorenz T Biegler, Omar Ghattas, Matthias Heinkenschloss, and Bart van Bloemen Waanders. Large- scale pde-constrained optimization: an introduction. In Large-scale PDE-constrained optimization, pp. 3-13. Springer, 2003.
Åke Björck. Numerical Methods for Least Squares Problems. Society for Industrial and Applied Mathematics, 1996.
Mathieu Blondel and Vincent Roulet. The elements of differentiable programming, arXiv preprint arXiv:2403. 14606, 2024.
Mathieu Blondel, Quentin Berthet, Marco Cuturi, Roy Frostig, Stephan Hoyer, Felipe Llinares-López Fabian Pedregosa, and JeanlPhilippe Vert. Efficient and modular implicit differentiation. In Advances in Neural Information Processing Systems, volume 35, pp. 5230-5242, 2022.

---

Published as a conference paper at ICLR 2026
Thomas Bolton and Laure Zanna. Applications of Deep Learning to Ocean Data Inference and Subgrid Parameterization. Journal of Advances in Modeling Earth Systems, 11(1):376 399, 2019.
J Frédéric Bonnans and Alexander Shapiro. Perturbation analysis of optimization problems. Springer Science & Business Media, 2013.
Nicholas Boumal. An Introduction Optimization Smooth Manifolds. Cambridge University Press, 2024.
George EP Box, Gwilym M Jenkins, Gregory C Reinsel, and Greta M Ljung. Time series analysis: forecasting and control. John Wiley & Sons, 2015.
Stephen Boyd and Lieven Vandenberghe. Convex optimization. Cambridge university press, 2004.
George Casella and Roger Berger. Statistical Inference. Duxbury Resource Center, June 2001.
Nithin Chalapathi, Yiheng Du, and Aditi Krishnapriyan. Scaling physics-informed hard constraints with mixture-of-experee In International Conference on Learning Representations, 2024.
Hao Chen, Gonzalo E. Constante Flores, and Can Li. Physics-informed neural networks with hard linear equality constraints. Computers and Chemical Engineering. 189:108764, 2024.
Chaoran Cheng, Boran Han, Danielle C. Maddix, Abdul Fatir Ansari, Andrew Stuart, Michael W. Ma- honey, and Yuyang Wang. Gradient-free generation for hard-constrained systems. In International Conference on Learning Representations, 2025.
Abhimanyu Das, Weihao Kong, Rajat Sen, and Yichen Zhou. decoder-only foundation model for time-series forecasting. arXiv preprint arXiv:23100 10. 10688, 2023.
Priya L Donti, David Rolnick, and J Zico Kolter. DC3: A learning method for optimization with hard constraints. In International Conference on Learning Representations, 2021.
Chris Edwards. Neural networks learn to speed up simulations. Communications of the ACM, 65(5): 27-29, 2022.
Carson Eisenach, Yagna Patel, and Dhruv Madeka. Mgtransformer: Multi-horizon forecasts with context dependent and feedback-aware attention. arXiv preprint arXiv:2009. 14799, 2022.
Rasool Fakoor, Taesup Kim, Jonas Mueller, Alexander J. Smola, and Ryan J. Tibshirani. Flexible model aggregation for quantile regression. Journal of Machine Learning Research, 24(162):1-45 2023.
Yuchen Fang, Sen Na, Michael W Mahoney, and Mladen Kolar. Fully stochastic trust-region sequential quadratic programming for equality-onnstrained optimization problems. SIAM Journal on Optimization, 34(2):2007-2037, 2024.
Meire Fortunato, Tobias Pfaff, Peter Wirnsberger, Alexander Pritzel, and Peter Battaglia. Multiscale meshgraphnets. 2nd Al4Science Workshop the 39th International Conference on Machine Learning (ICML, 2022.
Zhihan Gao, Xingjian Shi, Boran Han, Hao Wang, Xiaoyong Jin, Danielle C. Maddix, Yi Zhu, Mu Li, and Yuyang Wang. PreDiff: Precipitation nowcasting with latent diffusion models. In Advances in Neural Information Processing Systems, volume 36, 2023.
Jan Gasthaus, Konstantinos Benidis, Yuyang Wang, Syama Sundar Rangapuram, David Salinas, Valentin Flunkert, and Tim Januschowski. Probabilistic forecasting with spline quantile function rns.. In Proceedings of the Twenty-Secon International Conference Artificial Intelligence and Statistics, volume 89, pp. 1901-1910. PMLR, 2019.
Philip E. Gill and Elizabeth Wong. Sequential quadratic programming methods. In Mixed Integer Nonlinear Programming, pp. 147-224, New York, NY, 2012. Springer New York.
Mark Girolami and Ben Calderhead. Riemann manifold langevin and hamiltonian monte carlo methods. Journal of the Royal Statistical Society Series Statistical Methodology, 13(2):12-214,, 2011.
12

---

Published as a conference paper at ICLR 2026
Tilmann Gneiting and Adrian E Raftery. Strictly proper scoring rules, prediction, and estimation. Journal of the American statistical Association. 102(477):359-3777 0000599-77,077 2007.
Tilmann Gneiting, Adrian E Raftery, Anton Westveld, and Tom Goldman. Calibrated probabilistic forecasting using ensemble model output statistics and minimum erps estimation. Monthly Weather Review, (335))1098-111 2005.
Gene H. Golub and Chen Greif. On solving block-structured indefinite linear systems. SIAM Journal on Scientific Computing, 24(6):276--292, 2003.
Stephen Gould, Basura Fernando, Anoop Cherian, Peter Anderson, Rodrigo Santa Cruz, and Edison Guo. On differentiating parameterized argmin and argmax problems with application to bi-level optimization. arXiv preprint arXiv: 045477 2016.
Stephen Gould, Richard Hartley, and Dylan Campbell. Deep declarative networks. IEEE Trans. Pattern Anal. Mach. Intell., (00e))3388840004 2022.
Andreas Griewank and Andrea Walther. Evaluating derivatives: principles and techniques of algorithmic differentiation. SIAM, 2008
Gaurav Gupta, Xiongye Xiao, and Paul Bogdan. Multiwavelet-basee operator learning for differential equations. In Advances in Neural Information Processing Systems, volume 34, pp. (4048-24062, 2021.
Derek Hansen, Danielle C. Maddix, Shima Alizadeh, Gaurav Gupta, and Michael W. Mahoney. Learn- ing physical models that respect conservation laws. In Proceedings of the 40th International Conference on Machine Learning, volume 202, pp. 12469-12510. PMLR, 2023.
Ami Harten. High resolution schemes for hyperbolic conservation laws. Journal of computational physics, 135(2):260-22) 278, 1997.
llgee Hong, Sen Na, Michael W Mahoney, and Mladen Kolar. Constrained optimization via exact augmented lagrangian and randomized iterative sketching. In Proceedings of the 40th International Conference Machine Learning, volume 202, pp. 13174-13198. PMLR, 2023.
Shi Bin Hoo, Samuel Müller, David Salinas, and Frank Hutter. The tabular foundation model tabpfn outperforms specialized time series forecasting models based on simple features. arXiv preprint arXiv:2501.000 02945, 2025.
Thomas JR Hughes. The finite element method: linear static and dynamic finite element analysis. Courier Corporation, 2003.
Rob J Hyndman and George Athanasopoulos. Forecasting: principles and practice. OTexts, 2018.
Rob J. Hyndman, Roman A. Ahmed, George Athanasopoulos, and Han Lin Shang. Optimal com- bination forecasts for hierarchical time series. Computational statistics data analysis, 55(9): 2579-2589, 2011.
Rob Hyndman, George Athanasopoulos, Azul Garza, Cristian Challu, Max Mergenthaler, and Kin G. Olivares. Forecasting: Principles and Practice, the Pythonic Way, OTexts, Melbourne, Australia, 2025. available at https://oeexts.com//ppp///
Harry Joc. Multivariate models and multivariate dependence concepts. CRC press, 1997.
A. Kadambi, C. de Melo, CJ. Hsieh, M. Srivastava, and S. Soatto. Incorporating physics into data-driven computer vision. Nat Mach Intell, 5:572-580, 2023.
Rudolph Emil Kalman. A new approach to linear filtering and prediction problems. Journal of Basic Engineering, 82(1):35-44 45, 1960.
George Em Karniadakis, Ioannis G Kevrekidis, Lu Lu, Paris Perdikaris, Sifan Wang, and Liu Yang. Physics-informed machine learning. Nature Reviews Physics, ((6):422-440 440, 2021.
Muhammad Firmansyah Kasim and Yi Heng Lim. Constants motion network. In Advances in Neural Information Processing Systems, volume 35, pp. 25295-25305 2022.
13

---

Published as a conference paper at ICLR 2026
Alex Kendall and Yarin Gal. What uncertainties do need in bayesian deep learning for computer vision? In Proceedings of the 31st International Conference on Neural Information Processing Systems, pp. 55580-5590. Curran Associates Inc., 2017.
Alex Kendall, Yarin Gal, and Roberto Cipolla. Multi-task learning using uncertainty to weigh losses for scene geometry and semantics. In 2018 IEEE Conference Computer Vision and Pattern Recognition (CVPR), pp. 7482-7491. Computer Vision Foundation IEEE Computer Society, 2018.
Hyunjik Kim, Andriy Mnih, Jonathan Schwarz, Marta Gamelo,, Ali Eslami, Dan Rosenbaum, Oriol Vinyals, and Yee Whye Teh. Attentive neural processes. arXiv preprint arXiv: 1901 000161,61, 2019.
D. Kingma and J. Ba. A method for stochastic optimization. In International Conference on Learning Representations, 2015.
Dmitrii Kochkov, Janni Yuval, Ian Langmore, Peter Norgaard, Jamie Smith, Griffin Mooers, Milan Klöwer, James Lottes, Stephan Rasp, Peter Düben, al. Neural general circulation models for 2024.
weather and climate. Nature, (32(8027):1000-10
Steven George Krantz and Harold R. Parks. The implicit function theorem: history, theory, and applications. Springer Science & Business Media, 2002.
Aditi Krishnapriyan, Amir Gholami, Shandian Zhe, Robert Kirby, and Michael W Mahoney. Char- acterizing possible failure modes in physics-informed neural networks. In Advances in Neural Information Processing Systems, volume 34, pp. 26548-26560, 2021.
Quoc Le, Alex J Smola, and Stéphane Canu. Heteroscedastic gaussian process regression. In Proceedings of the 22nd International Conference on Machine Learning, pp. 489-496. PMLR, 2005.
Randall Veque. Numerical Methods for Conservation Laws. Lectures in mathematics ETH Zürich. Birkhäuser Verlag, 1990.
Randall J. Veque. Finite volume methods for hyperbolic problems, volume 31. Cambridge university press, 2002.
Randall J. Vequee Finite difference methods for ordinary and partial differential equations: steady-state and rime-dependent problems. SIAM, 2007.
Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Neural operator: Graph kernel network for partial differential equations. arXiv preprint arXiv:2003.03443 2020.
Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Fourier Neural Operator for Parametric Partial Differential Equations. In International Conference on Learning Representations, 2021.
Zongyi Li, Hongkai Zheng, Nikola B. Kovachki, David Jin, Haoxuan Chen, Burigede Liu, Kamyar Azizzadenesheli, and Anima Anandkumar. Physics-informet neural operator for learning partial differential equations. ACM / IMS Data Sci., 13). 2024.
Bryan Lim, Sercan Ô Ank, Nicolas Loeff, and Tomas Pfister. Temporal fusion transformers for interpretable multi-horizon time series forecasting. International Journal of Forecasting, 37(4): 2021.
Konstantin Lipnikov, Gianmarco Manzini, J. David Moulton, and Mikhail Shashkov. The mimetic finite difference method for elliptic and parabolic problems with staggered discretization of diffusion coefficient. Journal of Computational Physics, 305: 105:111-126 2016.
Lu Lu, Pengzhan Jin, Guofei Pang, Zhongqiang Zhang, and George Em Karniadakis. Learning nonlinear operators via DeepONet based on the universal approximation theorem of operators. Nature Machine Intelligence, 3(3):218-229 2021.
14

---

Published as a conference paper at ICLR 2026
Danielle C. Maddix, Luiz Sampaio, and Margot Gerritsen. Numerical artifacts in the generalized porous medium equation: Why harmonic averaging itself is not to blame. Journal of Computational Physies, 361:280-220 298, 2018a.
Danielle C. Maddix, Luiz Sampaio, and Margot Gerritsen. Numerical artifacts in the discontinuous generalized porous medium equation: How to avoid spurious temporal oscillations. Journal of Computational Physics, 368:277-298, 2018b.
James E. Matheson and Robert Winkler. Scoring rules for continuous probability distributions. Management science, 22(10)) 1976.
Takashi Matsubara and Takaharu Yaguchi. Finde: Neural differential equations for finding and preserving invariant quantities. In International Conference on Learning Representations, 2023.
Youngjae Min, Anoopkumar Sonar, and Navid Azizan. Hard-constrained neural networks with universal approximation guarantees. arXiv preprint arXiv:2410. 10807, 2024.
S. Chandra Mouli, Danielle C. Maddix, Shima Alizadeh, Gaurav Gupta, Andrew Stuart, Michael W. Mahoney, and Yuyang Wang. Using Uncertainty Quantification to Characterize and Improve Out-of-Domain Learning for PDEs. In Proceedings of the 40th International Conference on Machine Learning, volume 235, pp. 36372-36418 PMLR, 2024.
Eike Hermann Müller. Exact conservation laws for neural network integrators of dynamical systems. arXiv preprint arXiv:2209. 11661, 2022.
Jorge Nocedal and J. Wright, Steven. Numerical Optimization. Springer, 2006.
Geoffrey Négiar, Michael W. Mahoney, and Aditi S. Krishnapriyan. Learning differentiable solvers for systems with hard constraints. In International Conference on Learning Representations, 2023.
Kin G. Olivares, Federico Garza, David Luo, Cristian Challú, Souhaib Ben Taieb, Shanika I.. Wickra- masuriya, and Artur Dubrawski. HierarchicalForecastt A reference framework for hicrarchical forecasting in python. arXiv preprint arXiv:2207.03517, 2022.
Kin Olivares, O Nganba Meetei, Ruijun Ma. Rohan Reddy, Mengfei Cao, and Lee Dicker. Proba- bilistic hicrarchical forecasting with deep poisson mixtures. International Journal of Forecasting, 40(2):470-49)) 2024a.
Kin G. Olivares, Geoffrey Négiar, Ruijun Ma. Oinam Nganba Meetei, Mengfei Cao, and Michael W Mahoney. CLOVER Probabilistic forecasting with coherent learning objective reparameteri- zation. Transactions Machine Learning Research, 2024b.
Youngsuk Park, Danielle C. Maddix, François-Xaves Xavier Aubet, Kelvin Kan, Jan Gasthaus, and Yuyang Wang. Learning quantile functions without quantile crossing for distribution-fre time series forecasting. Proceedings of The 25th International Conference on Artificial Intelligence and Statistics, volume 151, pp. 8127-8150. PMLR, 2022.
Fotios Petropoulos, Daniele Apiletti, Vassilios Assimakopoulos, Mohamed Zied Babai, Devon Barrow, Souhaib Ben Taieb, Christoph Bergmeir, Ricardo J Bessa, Jakub Bijak, John E Boylan, et al. Forecasting: theory and practice. International Journal of Forecasting, 38(3):705-871, 2022.
Tobias Pfaff, Meire Fortunato, Alvaro Sanchez-Gonzalez, and Peter W. Battaglia. Learning mesh- based simulation with graph networks. International Conference on Learning Representations, 2021.
Iln Price, Alvaro Sanchez- Sanchez-Gonzalez, Ferran Alet, Tom R. Andersson, Andrew El-Kadi, Dominic Masters, Timo Ewalds, Jacklynn Stott, Shakir Mohamed, Peter Battaglia, Remi Lam, and Matthew Willson. Probabilistic weather forecasting with machine learning. Nature, (8744)4)) 2025.
Christopher Rackauckas, Yingbo Ma, Julius Martensen, Collin Warner, Kirill Zubov, Rohit Supekar, Dominic Skinner, Ali Ramadhan, and Alan Edelman. Universal differential equations for scientific machine learning. arXiv preprint arXiv:2001.0483 2020.
15

---

Published as a conference paper at ICLR 2026
Maziar Raissi, Paris Perdikaris, and George E Kariiaaaiis Physics- inssrcs-ifforme neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational physics, 578:686-707, 2019.
Syama Sundar Rangapuram, Lucien D. Werner, Konstantinos Benidis, Pedro Mercado, Jan Gasthaus, and Tim Januschowski. End-to-end learning of coherent probabilistic forecasts for hierarchical time series. In International Conference on Machine Learning, pp. 8832-8843. PMLR, 2021.
Syama Sundar Rangapuram, Shubham Kapoor, Rajbir Singh Nirwan, Pedro Mercado, Tim Januschowski, Yuyang Wang, and Michael SchneieS Coherent probabilistic forecasting of temporal hierarchies. In Proceedings of The 26th International Conference on Artificial Intelligence and Statistics, volume 206, pp. 9362-9376. PMLR, 2023.
Carl Edward Rasmussen and Chrisotpher K.I. Williams. Gaussian Processes for Machine Learning. MIT Press, 2006.
Stephan Rasp and Sebastian Lerch. Neural networks for postprocessing ensemble weather forecasts. Monthly Weather Review, 146(11):3))) 1):3885 3900, 2018
Suman Ravuri, Karel Lenc, Matthew Willson, Dmitry Kangin, Remi Lam, Piotr Mirowski, Megan Fitzsimons, Maria Athanassiadou, Sheleem Kashem, Sam Madge, et al. Skilful precipitation nowcasting using deep generative models of radar. Nature, 597(7878)662-6777, 2021.
Jack Richter-Powell, Yaron Lipman, and Ricky T. Q Chen. Neural conservation laws: A divergence- free perspective. In Advances in Neural Information Processing Systems, volume 35, 2022.
Christian P Robert, George Casella, and George Casella. Monte Carlo statistical methods, volume 2.
Springer, 1999.
R. Tyrrell Rockafellar and Roger J.B. Wets. Variational analysis, volume 317. Springer Science Business Media, 2009,
Jo Bo Rosen. The gradient projection method for nonlinear programming, part i. linear constraints. Journal of the society for industrial and applied mathematics, 8(1):181-217, 1960.
Lconid Rudin, Stanley Osher, and Emad Fatemi. Nonlinear total variation based noise removal algorithms. Physica D: nonlinear phenomena, 60(1-4):259 -2622268,-,, 1992.
Nadim Saad, Gaurav Gupta, Shima Alizadeh, and Danielle C. Maddix. Guiding continuous operator learning through physics-based boundary constraints. In International Conference on Learning Representations, 2023.
David Salinas, Michael Seblke-Schneider, Laurent Callot, Roberto Medico, and Jan Gasthaus. High- dimensional multivariate forecasting with low-rank gaussian copula processes. In Advances in Neural Information Processing Systems, volume 32, 2019.
David Salinas, Valentin Flunkert, Jan Gasthaus, and Tim Januschowski. Deepar: Probabilistic forecasting with autoregressive recurrent networks. International Journal of Forecasting, 36(3): 1181-1191, 2020.
Alexander Schein, Kevin T. Carlberg. and Matthew J. Zahr. Preserving general physical properties in model reduction dynamical systems via Journal for Numerical Methods in Engineering, 122(14):3368-339) 2021.
projection. International
James A. Sethian and John Strain. Crystal growth and dendritic solidification. Journal of Computa- tional Physics, 98(2):231-25) 1992.
Nicki Skafte, Martin Jørgensen, and Søren Hauberg. Reliable training and estimation of variance networks. In Advances in Neural Information Processing Systems, volume 32, 2019.
Andrew Stim, Harm Wessels, Megan Schertzer, Laura Pereira, Neville Sanjana, and David Knowles. Faithful heteroscedastic regression with neural networks. In International Conference on Artificial Intelligence and Statistics, pp. 5593-5613. PMLR, 2023.
16

---

Published as a conference paper at ICLR 2026
Patrick Obin Sturm and Anthony S. Wexler. Conservation laws in neural network architecture: enforcing the atom balance a julia- based photochemical model (v0.2.0). Geosci. Model Dev., 15:3417-3431,
2022. 2022.
Haixiang Sun, Ye Shi, Jingya Wang, Hoang Duong Tuan, H Vincent Poor, and Dacheng Tao. Alternating differentiation for optimization layers. arXiv preprint arXiv:2210.01802, 2022.
Roberto Szechtman and Peter Glynn. Constrained monte carlo and the method control variates. In Proceeding of the 2001 Winter Simulation Conference (Cat. No. 01 CH373044, volume 1, pp. 394-400. IEFE, 2001.
Souhaib Ben Taieb, James W. Taylor, and Rob J. Hyndman. Coherent probabilistic forecasts for hierarchical time series. In International Conference on Machine Learning, pp. 3348-3357. PMIR, 2017.
Maxime Taillardat, Olivier Mestre, Michaël Zamo, and Philippe Naveau. Calibrated ensemble forecasts using quantile regression forests and ensemble model output statisties. Monthly Weather Review, 144(6):2375-2399 2016.
Irina Kalashnikova Tezaur, Jeffrey A. Fike. Kevin Thomas Carlberg, Matthew F. Barone, Danielle Maddix, Erin E. Mussoni, and Maciej Balajewicz. Advanced fluid reduced order models for compressible flow. Sandia National Laboratories Report, Sand No. 2017-10335, 2017.
Tourism Australia, Canberra. Tourism Research Australia (2005), Travel by Australians. https: // //www, www.aa.gga...........kkk....... www.co/m.....oo sblanche/ quarterly
wtw.aaaaaaa..aaaaarrereeeeeeeeee torrytturiii sm- in- australia/, itt.iiiminnauustaaaaaaaaaaaaaaaa 2005.
Mare Toussaint, Kelsey R. Allen, Kevin A. Smith, and Joshua B. Tenenbaum. Differentiable physics and stable modes for tool-use and manipulation planning extended abtract. Proceedings the Eienty-Eighht International Joint Conference Artificial Intelligence, IJCAI-19, pp. 6231-6235. International Joint Conferences on Artificial Intelligence Organization, 2019.
Utkarsh Utkarsh, Valentin Churavy, Yingbo Ma, Tim Besard, Prakitr Srisuma, Tim Gymnich, Adam R Gerlach, Alan Edelman, George Barbastathis, Richard D Braatz, et al. Automated translation and accelerated solving of differential equations on multiple gpu platforms. Computer Methods in Applied Mechanics and Engineering. 119:116591, 2024.
Utkarsh Utkarsh, Pengfei Cai, Alan Edelman, Rafael Gomez-Bombarelli, and Christopher Vincent Rackauckas. Physies-constrained flow matching: Sampling generative models with hard constraints. In Advances in Neural Information Processing Systems, volume 38, 2025.
J.M. van der Meer, J.F.B.M. Kraaijevanger, M. Möller, and J.D. Jansen. Temporal oscillations in the simulation of foam enhanced oil recovery. ECMOR XV- 15th European Conference the Mathematics of Oil Recovery, pp. 1-20. 2016.
Aad W Van der Vaart. Asymptotic statistics, volume 3. Cambridge university press, 2000.
Tim Van Erven and Jairo Cugliari. Game-theoretically optimal reconciliation contemporaneous hierarchical time series forecasts. In Modeling and stochastic learning for forecasting in high dimensions, pp. 297-317. Springer, 2015.
Juan Luis Vázquez. The Porous Medium Equation: Mathematical Theory. The Clarendon Press, Oxford University Press, Oxford, 2007.
Kilian Q. Weinberger and Lawrence K. Saul. Distance metric learning for large margin nearest neighbor classification. Journal of machine learning research, 10(2), 2009.
Ruofeng Wen, Kari Torkkola, Balakrishnan Narayanaswamy, and Dhruv Madeka. A quantile recurrent forecaster. arXiv preprint arXiv:1711.1105: 2018.
Alistair White, Anna Büttner, Maximilian Gelbrecht, Valentin Duruisseaux. Niki Kilbertus, Frank Hellmann, and Niklas Boers. Projected neural differential equations for learning constrained dynamics. arXiv preprint 2rXiv:2410.23667 2024.
17

---

Published as a conference paper at ICLR 2026
Shanika L. Wickramasuriya, George Athanasopoulos, and Rob J. Hyndman and. Optimal forecast reconciliation for hierarchical and grouped time series through trace minimization. Journal of the American Statistical Association, 114(526):804 819, 2019.
Jeffrey Willette, Hae Beom Lee, Juho Lee, and Sung Ju Hwang. Meta learning low rank covariance factors for energy-based deterministic uncertainty. arXiv preprint arXiv:21 2021.
Robert B. Wilson. A Simplicial Method for Convex Programming. PhD thesis, Harvard University, 1963.
Gerald Woo, Chenghao Liu, Akshat Kumar, Caiming Xiong, Silvio Savarese, and Doyen Sahoo. Unified training of universal time series forecasting transformers. arXiv preprint 22X2v2240..02592. 2024.
Matthew J. Zahr and Per-Olof Persson. An adjoint method for a high-order discretization of deforming domain conservation laws for optimization of flow problems. Journal of Computational Physics, 326:516-543, 2016.
Laure Zanna and Thomas Bolton. Data-Driven Equation Discovery of Ocean Mesoscale Closures. Geophysical Research Letters, 47(17), 2020.
18

---

Published as a conference paper at ICLR 2026
A RELATED WORK
In this section, we review works that have been motivated by dealing with hard constraints in various domains including imposing constraints in neural networks (Min al., 2024; Donti et al., 2021), probabilistic time series forecasting (Rangapuram et al., 2023; 2021; Olivares et al., 2024b) and scientific machine learning (Négiar et al., 2023; Hansen et al., 2023). Table 5 summarizes several of these methods. We see that existing general methododologies, e.g., HardNet (Min et al., 2024) and DC3 (Donti et al., 2021), work across various domains and different types of constraints- HardNet handles convex constraints, and DC3 tackles nonlinear ones. The biggest limitation of these methods is that they provide point estimates only. Despite having the point forecast satisfying the constraints, they unsuitable for PDEs and forecasting applications, which generally now require variance estimates. Hier-E2E (Rangapuram et al., 2021) and CLOVER (Olivares et al., 2024b) specialized solutions for forecasting problems, which both deal with probabilistic forecasting under linear constraints. Linear constraints are common in the time series forecasting domain. Both methods require sampling during training, which can be computationally intensive. Within the PDE- focused methods, ProbConserv (Hansen et al., 2023) and HardC (Hansen et al., 2023) handle linear constraints and include variance estimates in their probabilistic models. The training procedure with the constraint is not end-to-end since the constraint is only applied at inference time. PDE-CL (Négiar al., 2023) handles nonlinear constraints and supports end-to-end training, but the cost of not supporting variance estimation.
Table 5: Summary of methods motivated by dealing with hard constraints in various domains: imposing constraints in neural networks (Min et al., 2024; Donti et al., 2021), probabilistic time series forecasting (Rangapuram et al., 2023; 2021; Olivares et al., 2024b) and scientific machine learning (Négiar et al., 2023; Hansen et al., 2023). For models that only provide point estimates, we evaluate their capabilities on sampling-free training and satisfying constraints distributions, while treating the point estimate as a degenerate probabilistic distribution.
Method
Domain
Constraint Type
End-to- End
Prob. Model w/ Variance Estimate
Sampling- free Truining
Constraint on Dstba.
HardNet (Min al., 2024)
General
Convex
X
DC3 (Donti al., 2021)
General
Nonlinear
X
Hicr-E2E (Rangapuram al., 2021; 2023)
Forccasting
Linear
CLOVER (Olivares al., 2024h)
Forecasting
Linear
/
PDE-CL (Négiar et al., 2023)
PDEs
Nonlinear
ProbConserv (Hansen ct al., 2023)
PDEs
Lincar
HardC (Hansen et al., 2023)
PDEE
Linear
ProbilardE2 (Ours)
Gencral
Nonlincar
Our proposed method P ProbHard bridges the gaps left by its predecessors. It combines the flexibility of general domain application with the ability to handle nonlinear constraints, and maintains end-to-end training capability. Perhaps most notably, it achieves this while incorporating probabilistic modeling with variance estimates, supporting sampling-free training, and maintaining constraints distributions.
A.1 IMPOSING DETERMINISTIC CONSTRAINTS NEURAL NETWORKS
Enforcing constraints in neural networks (NNs) has been explored in various forms. In fact, activation functions, e.g., sigmoid, ReLU, and softmax, inherently impose implicit constraints by restricting out- puts to specific intervals. Another well-established method for enforcing constraints in NNs involves differentiating the Karush-Kuhn Tucker (KKT) conditions, which enables backpropagation through optimization problems. This technique has led to the development of differentiable optimization layers (Amos & Kolter, 2017; Agrawal et al., 2019) and projected gradient descent methods (Rosen, 1960).
Most commonly, soft constraint methods, e.g., Lagrange duality based methods, are often employed in ML to balance minimizing the primary objective with satisfying the constraints. These methods typically do so by adding the constraint as a penalty term to the loss function (Battaglia et al., 2018). For example, Lagrange dual methods and relaxed formulations are frequently used allow flexibility
19

---

Published as a conference paper at ICLR 2026
in the optimization process, while still guiding the model toward feasible solutions. These methods enforce-adhereece to the constraints during training; and this lack of
encourage- but do not strictly enforce- adherence strict enforcement can be undesirable in some scientific disciplines, where known constraints must be satisfied exactly (Hansen et al., 2023; Rangapuram et al., 2021).
More recently, there have been approaches that have been motivated by satisfying hard constraints. DC3 (Donti et al., 2021) is a general method for learning a family of constrained optimization problems using a correction and variable completion procedure. The variable completion approach has a strong theoretical and practical foundation. A limitation is that it does require knowledge of the structure of the matrix identify these corresponding predicted and completed variables, which hinders its generalizability. In addition for inequality constraints. it only achieves hard constraint satisfaction asymptotically; that is, the "correction" procedure to enforce inequality constraints is carried out through gradient-deseen optimization algorithms (Min et al., 2024; Donti et al., 2021).
Projection-basse based methods are an alternate method for enforcing hard constraints in NNs. Min et al. (2024) identify cases where the aforementioned DC3 framework (Donti et al., 2021) is outperformed by their proposed HardNet projection layer approach. Additionally, Min et al. (2024) investigate the expressiveness of projection layers, which builds on the foundational work in Agrawal et al. (2019); Amos & Kolter (2017), to further advance the understanding of constraint enforcement in NNs. Projection- basedd methods have also been used to enforce constraints on specific architectures, e.g., neural ordinary differential equations (Neural ODEs) (Kasim & Lim, 2022; Matsubara & Yaguchi, 2023). In particular, White al. (2024) use closed-form projection operator to enforce nonlinear constraint g(u) = 0 in a Neural ODE, using techniques from Boumal (2024). A common limitation of these works is that they impose the constraint deterministically, on point estimates rather than on an entire probability distribution.
A.2 PROBABILISTIC TIME SERIES FORECASTING
Probabilistic time series forecasting extends beyond predicting point estimates, e.g., the mean or median, by providing framework to capture uncertainty, with practical application in estimating high quantiles, e.g., P99. Classical statistical models, e.g., autoregressive integrated moving average (ARIMA) models (Box al., 2015), state-space models (Kalman, 1960), and copula-based models (Joe, 1997) are prominent examples. More recently, deep learning models. e.g., DeepAR (Salinas al.. 2020) and its multivariate extension DeepVAR (Salinas al., 2019), multivariate quantile regression-based models (Wen et al., 2018; Eisenach et al., 2022; Park et al., 2022), temporal fusion transformers (TFT) (Lim et al., 2021), and foundational models based on large language models (LLMs) (Ansari et al., 2024; Hoo et al., 2025; Das et al., 2023; Woo et al., 2024) have shown success. See Benidis et al. (2022) for an overview.
Linear constraints are important in hierarchical time series forecasting, where coherent aggregation constraints are required over regions (Rangapuram et al., 2021; Olivares et al., 2024b) and over temporal hierarchies (Rangapuram et al., 2023). This constraint is critical in scenarios where higher- level forecasts must be aggregates of lower-level ones, which is a common requirement in time-series forecasting. Early works in hierarchical forecasting focus mean forecasts under linear/hierarchical constraints, starting from the naive bottom-up and top-down approachess (Hyndman & Athanasopoulos, 2018). More recently, Hyndman et al. (2011) show that the Middle-Ouu Out projection-baser method yields better forecast accuracy. Since then, projection-baser
reconciliation methods, e.g., GTOP (Van Erven & Cugliari, 2015), MinT (Wickramasuriya et al., 2019), and ERM (Ben Taieb & Koo, 2019) have been developed. These methods leverage generic time series models, e.g., ARIMA and exponential smoothing (ETS), derive the unconstrained mean forecast, and then they use linear projection to map these forecasts to the consistent space. Taieb et al. (2017) further extend the reconciliation method (MinT) to probabilistic forecasting by developing a method called PERMBU that constructs cross-sectional dependence through sequence of permutations. A more thorough review of forecasting reconciliation is provided in Athanasopoulos et al. (2024).
To better handle the trade-off between forecast accuracy and coherence within the model, several works have proposed end-to-end methodologies. For example, Rangapuram et al. (2021) propose Hier-E2E, which is an end-to-end learning approach that imposes constraints via orthogonal projection on samples from the distribution. Hier-EEE produces coherent probabilistic forecasts without requiring explicit post-processing reconciliation. One limitation is that Hier-E2E relies on
20

---

Published as a conference paper at ICLR 2026
expensive sampling techniques achieve this coherence, by projecting directly on the samples rather than on the distribution itself, which has a closed-form expression.
Separately, DPMN (Olivares et al., 2024a) adopts equality constraint completion approach similar to that in DC3 (Donti et al., 2021), rather than a projection-based approach, for satisfying the co- herency constraint. DPMN assumes that the bottom-level series follow a Poisson mixture distribution, with a multivariate discrete distribution on the Poisson rates across bottom-level series. Compared to Hier-E2E, DPMN uses a CNN-based encoder rather than Dee, and it shows improved forecast accuracy over Hier-E2E. As follow-up work to DPMN, Olivares et al. (2024b) propose CLOVER, framework which enforces coherency as a hard constraint in probabilistic hierarchical time series forecasting models using a CNN encoder. In particular, CLOVER only predicts the base forecasts in the first step, and it solves for the aggregate forecasts by leveraging the constraint relation. Finally, CLOVER models the joint distribution of all the forecasts in the scoring function calculation. Similar to Hier-E2E, CLOVER also relies on sampling to enforce hierarchical coherency and to generate uncertainty estimates. This affects the training time, and it requires tuning of the number of samples for an accurate approximation of this scoring function. Although CLOVER admits constraint satisfaction, the exact provably convergent procedure only exists for linear equality constraints (Donti et al., 2021), and it has not been applied to nonlinear equality or convex inequality constraints.
A.3 SCIENTIFIC MACHINE LEARNING (SCIML)
Partial differential equations (PDEs) are ubiquitous in science and engineering applications, and have been used to model various physical phenomena, ranging from nonlinear fluid flows with the Navier- Stokes equations to nonlinear heat transfer. Classical numerical methods to solve PDEs include finite difference (Le Veque, 2007), finite element (Hughes, 2003), and finite volume methods (Le Veque, 2002). These numerical methods discretize the solution on a spatio-temporal mesh, and the accuracy increases the mesh becomes finer. For this reason, numerical methods can be computationally expensive on real-world, time dependent, 3D spatial problems that require fine meshes for high accuracy.
Recently, Scientific Machine Learning (SciML) methods aim to alleviate the high computational requirement of numerical methods. State-of-the-art data-driven methods include operator-based methods, which aim to learn a mapping from PDE parameters or initial/boundary conditions to the PDE solution, e.g., Neural Operators (NOs) (Li al., 2020; 2021; Gupta et al., 2021) and DeepONet (Lu et al., 2021), and bassaag-paasiigg Graph Neural Networks (GNNs))-baeed MeshGraphNets (Pfaff et al., 2021; Fortunato al., 2022). These data-driven methods are not guaranteed to satisfy the PDE or known physical laws exactly, e.g., conservation laws (Hansen et al., 2023; Mouli et al., 2024) or boundary conditions (Saad et al., 2023; Cheng et al., 2025) since they only implicitly encode the physics through the supervised training simulation data (Kadambi et al., 2023).
Similar to imposing constraints on NNs, most work on imposing constraints in SciML has been focused on soft constraints. One well-known type of approach is Physics- Informed Neural Networks (PINNs) (Raissi et al., 2019; Kamiaaakii et al., 2021), which approximates the solution of a PDE as NN. PINNs and similarly Physics-Informed Neural Operators (PINOs) (Li et al., 2024) impose the PDE as an additional term in the loss function, akin to the aforementioned soft constraint regularization. Krishnapriyan et al. (2021); Edwards (2022) identify limitations of this approach on problems with large PDE parameter values, where adding this regularization term can actually make the loss landscape sharp, non-smooth and more challenging to optimize. In addition, Hansen et al. (2023) show that adding the constraint to the loss function does not guarantee exact constraint enforcement, which can be critical in the case of conservation and other physical laws. This constraint violation primarily happens because the Lagrangian duals of the constrained optimization problem are typically non-zero, i.e., the physical constraint is not exactly satisfied.
Recent work has studied imposing physical knowledge as hard constraints on various SciML methods. Négiar et al. (2023) propose PDE-CL, which uses differentiable programming and the implicit function theorem (Krantz& & Parks, 2002) to impose nonlinear PDE constraints directly. Chalapathi et al. (2024), extend this work by leveraging a mixture-of -feeeerrsss (MoE) framework to better scale the method. Similarly, Beucler et al. (2021) impose analytical constraints in NNs with applications to climate modeling. In addition, Universal Differential Equations (UDEs) (Rackauckas et al., 2020; Utkarsh et al., 2024) provide a GPU-compatible and end-to-end differentiable way to learn PDEs
21

---

Published as a conference paper at ICLR 2026
while also enforcing implicit constraints. Chen et al. (2024) propose KKT-hPINN enforce linear equality constraints by using a projection layer that is derived from the KKT conditions. These works show the benefit of imposing the PDE hard constraint rather than soft constraint. A limitation of these methods is that they impose the constraints deterninistically. and they do not provide estimates of the underlying variance uncertainties. To address this, Hansen et al. (2023) propose ProbConserv, which incorporates linear conservation laws as hard constraints on probabilistic models by performing an oblique projection to update the unconstrained mean and variance estimates, Limitations are that this projection is only applied as post post-processing step at inference time, and the method only supports linear constraints.
B UNIVERSAL APPROXIMATION GUARANTEES
In this section, we prove a universal approximation result for our differentiable probabilistic projection layer (DPPL) in the case of convex constraints. As a consequence of this result, our ProbHardE2E in Algorithm is a universal approximator, and can approximate any continuous target function that satisfies the given constraints. Our proof of this result generalizes the analysis of Min et al. (2024) to our broader framework with arbitrary Q. Since Q is symmetric positive
from the case Q definite, we compute its Cholesky factorization Q = LLL", where L denotes a lower triangular matrix. We then show that if fe is universal approximator. i.c., sufficiently wide and deep neural network, then our DPPL preserves this universal approximation capability. Hence, ProbHardED retains the expressiveness of neural networks, both in its probabilistic formulation and in enforcing hard constraints.
Our DPPL in Problem (7) is formulated in terms of projecting the samples 20 20(6()) (o( where Zo(a(t) (o) Eo()()) for some multivariate location-scat distribution F, and where (6l0), (()) D denotes training data from distribution D. The mcan (o(م(() R and covariance Eg (g())) R are the output from deep neural network, fo o - R The value of k depends on the approximation for Ee("")), e.g., k + n² for dense Eo(o()), k = 2n for assume in
De ()())) diag(o? ....,,,....................... ..ff....,?)) ,a)), k for (o(()) 1. For notational simplicity, this section that Se + I2" corresponds to the components that output the mean ())())) By setting ))))) us (x))) fo()()) in Problem (7), our DPPL can also be formulated in terms of projecting the mean u (1) ) as the following constrained least squares problem:
()()()) = arg min ()(a(((0)) - fo()())))
Do (o()
((
where Q denotes symmetric positive definite matrix and 9(-) < 0, 0 denote the convex constraints. In particular, we show that the projected mean ())) is a universal approximator of the true solution u E R". We now state the theorem and provide its proof below.
Theorem B.1. Consider Problem 9 with the projection step defined using symmetric positive definite (SPD) matrix Q iR= a deep neural network that is a universal approximator, fa - "R"" where & C I.mh denotes a compact set, convex constraints g(-) < 0, h, 0, and training data (0(0) () (1())
~ D from distribution D. For any continuous target function that satisfies the constraints, i.e., the true solution u + C R²
C(4),, where denotes the convex set of feasible points defined by the convex constraints and C($) denotes the space of continuous functions on 4, there exists a choice of network parameters for (o( = (o(())) E RR", such that the projected mean, which is composition of fo with the projection step, i.e., 11 (fo(o(1)) R", approximates the target function arbitrarily well, where fo : $ + C Hence, under these conditions, ProbHa rdE2E is universal approximator for constrained mappings.
()) ( E C
Proof. Let C C IR" denote the convex set of feasible points defined by the convex constraints. Consider the projection operator onto with respect to the Q-norm:
T& (v) argmin ao (o(e)))) -112-
(10)
22

---

Published as a conference paper at ICLR 2026
Since Q is symmetric positive definite (SPD). it has the following Cholesky factorization,
LL
where L E R=** denotes a lower triangular matrix with strietly positive diagonal entries, and hence invertible. By the definition of the Q-norm, and then substituting in its Cholesky factorization, have
( v12 (Mo(d!)) ( )
LI. (No()(1) - T..LL (Mo((5(0) -v)
((oo()(1) - L)L (Mo()(()) -)) (L (10((110) - u)) (L (Mo(o))) v)) -I ))) -)
(11)
This shows that the Q-norm is equivalent to the standard Euclidean norm after the linear transforma- tion L
We define the invertible linear mapping R + R by (v) = L Then using Eq. (11), the Q-norm in Eq. (10) can be written the Euclidean norm as follows:
n (v) argmin, bo (ct9))c [fa (f 8) 1
argmin (66)))c V ((a (T))
(12)
-ー (
argmin
where w = (d Hence, the projection can be expressed as the Euclidean projection onto transformed set (C). well known that the Euclidean projection onto closed convex set nonexpansive and is Lipschitz continuous. (See, e.g., Min et al. (2024).)
Now, suppose that (((()) is a deep neural network that is universal approximator, i.e., for any continuous function u > IR", C(Φ), and for any > 0, there exists parameters such that
[lu( (0) - fo fe( <
sup 5
where T C R. denotes compact set and C( denotes the space of continuous functions O. Let u + S R", u E C(Φ), be any continuous target function whose outputs satisfy the constraints. Since I is continuous (as the compasition of the continuous mapping V the Euclidean projection onto (C), and - - it follows by the universal approximation theorem and properties of continuous functions that the projected mean Io (J)()))) (a(((()) can uniformly approximate u((()) arbitrarily well on In other words, for every O, there exists choice of network parameters 0 such that
sup ( - - (Jo((())) <E g0
Thus, the composition of the neural network fo with the Q-norm projection retains the universal approximation property for any continuous target function satisfying the constraints.
C COMPUTATION OF POSTERIOR DISTRIBUTION FOR VARIOUS CONSTRAINT TYPES
In this section, we discuss how to compute the differentiable probabilistic projection layer (DPPL) that projects the distribution parameters (Eq. (6)) in FrobHardE25 for various constraint types, which are summarized Table
C.I LINEAR EQUALITY CONSTRAINTS
In this subsection, provide the closed-form expressions for the constrained posterior distribution parameters, i.e., the mean fi and covariance E in Eq. (8), from the DPPL in ProbHa rdE2F. for linear
23

---

Published as a conference paper at ICLR 2026
equality constraints. Linear equality constraints occur in wide range of applications, including coherency constraints hierarchical time series forecasting (Hyndman et al., 2011; Rangapuram et al., 2021; Petropoulos et al., 2022; Olivares et al., 2024b), divergence- -feee conditions in incompressible fluid flows (Raissi et al., 2019; Richter-Powell et al., 2022), boundary conditions in PDEs (Saad et al., 2023), and global linear conservation law constraints (Hansen et al., 2023; Mouli et al., 2024).
Proposition C.1. For linear equality constraints, h(ūi) = Aú b O, with A E Rex", with full row rank q, where q < n1 and b IRy, the optimal solution u* to Problem (7) is given as u (2) = P (I- PQ-1) where Pa-1 -----1 (AQ-1AT)- ) denotes oblique A,
projection operator, and At denotes the Moore-Penrose inverse. In addition, if F(u, Σ) and L for multivariate, location-scca distribution F. then u² Y, where Y F(p. Σ) and n. E are given in Eq. (8) with T(z) = u u²()), which simplifies to the closed-form expressions. ja Pa- 11 + - Po- A) and Σ = Pq .EPQ
Proof. Using the Lagrange multiplier À R we can form the Lagrangian of Problem (7) with linear constraints obtain:
1 T Qu - =Ta + X (A b). L(ú, 2) = -
2
The sufficient optimality conditions to obtain (u", are the first-order gradient conditions:
V L(, (;=), Qu' Q : + A = 0,
(13a)
V (0, A; 2) Au" - -b=0 0.
(13b)
Since Q is SPD, QT and Q- exists. Then from Eq. (13a), we obtain:
Q(u =) + A X = 0,
which simplifies to the following expression for u":
u' ='=2-Q--A A
(14)
We solve Eq. (13b) for u* using the Moore-Penrose inverse, i.e., w* = Alb, where A AT (AAT)-I Note that AA E R289 is invertible with full rank q since E Raxn has full row rank S n. Substituting this expression into Eq. (14) for u² gives:
A (AAT) > = A
Rearranging for the optimal Lagrange multiplier \", and multiplying both sides by A gives:
(AQ-*AT)X* A= (AA (AA
Now, AQ-IAT E RRx9 is invertible since A has full row rank q. Then obtain:
X" x==AQQ (AQQ-AAT-- -AA b).
Substituting in the expression for X" into Eq. (14) gives the following expression for the optimal
solution:
u = = A (AQ- A - (Az- b),
(15)
=( (T-Q-1 A A (AQ (AQ-'AT)- - + - ATA---AAT
Let
PQ- A -- be an oblique projection. To see that this is a projection, observe that
(16)
p2 = (I -I- (Q="AATAA00-----AANA---41144444
-I-2Q- =1-0""1AA40000000AA0000"""111110 A ATAQ-AAAA-AAAAAAAAATAAT -T)-AA
- -I-QQ-1 "AA0000"""""====="""""00000"1010
= -I-Q- -I-Q-IA(AQQ-1TT)>
Pq
24

---

Published as a conference paper at ICLR 2026
Then, the expression for u" in Eq. (15) simplifies to:
u" (z) PQ-1z Q (AQ " (A)-),))),, PQ-1z+ (A
(17)
PQ-12++ + +1 Pa-1) )Abb
since AA' = AAT (
Since the expression for u' in Eq. (17) is a linear transformation T of F(u, 2), can use Theorem 3.1 with T()) (²) write the expression for w F(ū, E), where:
A (w)) w² (n) (1
(18a)
Σ = JT (1)) = PQ-
(18b)
It can easily be verified that (u) = Pa-i by differentiating Eq. (17) with respect to z. We note that Eq. (18) holds exactly in the case of linear constraints since T is linear transformation of D
We note that our probabilistic method applies to underdetermined linear systems when q < where there is existence of many solutions. When q = n and has full row rank, the solution unique. In this case, the projection Po has the following deterministic solution,
e = A
Σ ==0. 0.
In this case, it reduces non-probabilistic point prediction methods, similar to the HardC baseline (Hansen et al., 2023), where only the mean is updated.
C.2 NONLINEAR EQUALITY CONSTRAINTS
In this subsection, describe how compute the DPPL in ProbHardE22 for general nonlinear equality constraints. Nonlinear equality constraints naturally arise in applications that involve structural, physical, or geometrie consistency. These include closed-loop kinematics robotics (Toussaint et al., 2019), nonlinear conservation laws (LeVeque, 1990) in PDE-constrained surrogate modeling (Biegler et al., 2003; Zahr Persson, 2016; Négiar et al., 2023) with applications in climate modeling (Bolton & Zanna, 2019; Zanna & Bolton, 2020; Beucler et al., 2021), compressible flows in aerodynamics (Tezaur et al., 2017) and atomic modeling (Müller, 2022; Sturm & Wexler, 2022).
Proposition C.2. For nonlinear equality constraints, h(ā) = 0 E 12, where h IR" + R3, the optimal solution (z) to Problem (7) forms pair (w" (2),²*) which satisfies u u²(2)
V-'Vh(u' (=)) X* and hiu" h=)"(=)) O. In addition, if Z ~ F(u, Σ) and ~ Z for multivariate, location-scale distribution F. then u² Y, where Y F(ù Σ) and p. E are given in Eq. (8) with T(z) = u² (2).
Proof. Using the Lagrange multiplier À R", we can form the Lagrangian of Problem (7) with nonlinear equality constraints to obtain:
(( 2td Qa-z" "uQuuee + X h(ū).
The sufficient optimality conditions to obtain (u², (²) are the first-order gradient conditions:
R(u", (+;)) (a = =(""" -) + Vh(u">* 0, (=)(ū,;; = h(u²)= =0
(19)
We solve Eq. (19) via root-finding methods, e.g., Newton's 's method for (u", (²) to obtain u² u²(z) =
arg( à : R(ū, >=; = 0}. where the root-finding solution u² implicitly dependent on z. Since
the expression for u² is a nonlinear transformation T of F(u, E). we can use Theorem 3.1 with T(2) (2) to write the expression for u" F(i, E), where:
j = T(x) = u² (n), Σ J (1)
(20a)
(20b)
hold to first-order accuracy. In the following Proposition C.3, we detail the iterative algorithm to compute the terms u² (p) and Jr (n) in Eq. (20).
25

---

Published as a conference paper at ICLR 2026
Proposition C.3. Let h(ū) e Et be smooth nonlinear equality constraint, where h R² + R Consider the constrained projection problem from Problem (7) with 2 = Ju:
u (1) arg min f(ū),
(21)
&ER" h(a)=0
where Q > and [(ū) = | 112 denotes our quadratic objective.
1. At each iteration, we solve the linearized Karush-Kin Kuhn- (KKT) system using the Schur complement to obtain:
A(+1)= = 11110-1007
( 1:000110 1)))
(22a)
a(++1) = p- 1-Q-111 Q(i- ((+11)
(22b)
where J(0 = Vi(ū(()T E R2xxn At the first iteration with ū( u. Eg. (22b) simplifies
to:
i = - Q-1 (JQ- h(y))
(23)
where J = V
2. At convergence, the Jacobian JT (1) of the projection map T(p) u² (n)) is given by:
(u)(k) 1
JT(N)=
(24)
0p
where J Vh(w*)T.
Proof. We begin with the matrix form of the KKT system derived in Eq. (19):
VaL(u xt) V.L(u. X*]
V((())+ +. X*1 0,
R(u", pp
(25)
h(w²)
with quadratic objective f defined in Problem 21.
1. Iteration Update. We use Newton's Method to linearize the KKT system in Eq. (25) evaluated at (a(++1), ((++) about the past iterate (ū(9) ( For the stationarity condition, which the first com- ponent of R(u", p)) in Eq. (25), we use the first-order Taylor expansion of Ro(ū(++1) 16+1); p) about the past iterate (i(), (()) to obtain:
Au(iti
Ro( ((((++)),x0+1),,,, (1) = Ro( A00 + Vo. A00 :1:
An(i+1)
Val(n(), (1) RL + XO)
A
An(+t17
VaL(W(), ( (V + (2)(2()1110)
11
An A
VaL(m ( V 20) _OTJ
= =0,
(26)
where Aa(++)) (++)) 0 A = (=-1) 10 and (2h(a() Rexaxq denotes the
Hessian of the constraints. Solving for the increments we obtain:
A
[V au L() (())
L
AA(+t)) = -VaL(a
(27)
For the feasibility condition, i.c., the second component of R(u", x; in Eq. (25), we also linearize the constraint as:
R, (0(+1), 10+1);;) h((((+1) h(a()) + ( (g(++1) - -((() 0.
(28)
26

---

Published as a conference paper at ICLR 2026
using first-order Taylor expansion (Newton's Method). Then,
1 A1 = ).
(29)
We can then combine Eq. (27) and Eq. (29) to form the following linearized system of KKT conditions:
V2 (( A(1)
JOST
A
V VaL( 1(0)
(30)
J(
A
h((
Note that the system of equations in Eq. (30) used in Sequential Quadratic Programming (SQP) (Wil- son, 1963; Nocedal & Wright, 2006; Gill & Wong, 2012) when there are inequality constraints. In addition, since our objective f is quadratic, we do not need to compute its second-order Taylor expansion, and only need to linearize the constraints. SQP reduces to Newton's Method when there are no constraints. In particular, Eq. (30) gives the standard unconstrained Newton step V A(((t))) = (((())) when h
Now, we use our quadratic objective f in Problem 21 to compute:
VaL(a 1 - Q(a 1)) + J 1
(31)
V L( (0)) Q+ ()++))(()) @
Note that Q is symmetric positive definite, but ()) is not guaranteed to be positive definite in the general case, especially at every iterate, which could make the Newton step undefined. Regularization may be needed ensure that (² is positive semi-definite. In addition, since (²)(2() Raxxxq is three-dimensional tensor, it is computationally expensive to compute this matrix of second derivatives, especially on our large-scale problem and through auto-differentiation (Griewank & Walther, 2008; Blondel & Roulet, 2024). Similar to the Gauss-Newtor method (Björck, 1996; Nocedal & Wright, 2006) for nonlinear least squares problems, we assume that the constraint is approximately affine near its optimal point u", and use only first-order constraint information. Hence, we set ()) x O. We note that even with these approximations for efficiency on large-scale problems, we still show strong performance the nonlinear constraint results in Table 4. An alternate approach could be to use low-rank approximation to the Hessian done in Quasi-Newton, e.g.. BEGS methods (Nocedal & Wright, 2006).
Using Eq. (31) with setting (²)((()) = O, Eq. (30) simplifies to:
[
Q
JO9T
A
[Q(()) 4) + -OOT h(ā())
(32)
AA(+1)
Then,
Q1 1) + BOOT (1(ttt)t ( (+(++)) = - h((
((i) = - - ()) -(()
(33a)
(33b)
We see that the only terms involving (0) cancel from both sides of the equation. Note that the method does not require tracking the dual variable, so it could also be equivalently reset to X = at each iteration, and we compute (++1) only for computing the primal update in Eq. (22b).
Since Q > 0, it is invertible, we can multiply Eq. (33a) by Q-1 to obtain:
(i+1 (++1) ()) = - (a) 11) - Q- O ((-I)
(34)
Multiplying both sides of Eq. (34) by (i) and using Eq. (33b), we can eliminate (a to obtain: -h((()) - -
p)) Q f(o ( x(++)
(35)
Since Q 0, Q- > O and then --1 > 0. and hence it is invertible. We can then solve
Eq. (35) for ((+1) to obtain:
( = () -1 (
(36)
which gives the desired Eq. (22a). Then, solving Eq. (34) for 16--1) gives:
u (41) 1 - x
(37)
27

---

Published as a conference paper at ICLR 2026
which is the desired Eq. (22b). Lastly, for the first iterate, substituting Eq. (36) into Eq. (37), setting i 0 and ()) = 1e gives the desired Eq. (23).
We can also solve Eq. (32) efficiently using block Gaussian elimination and the Schur complement (
IR (Golub & Greif, 2003). In particular, the block matrix in Eq. (32) can be factored into product of elementary matrices as:
O j
JOO 0
-oo- 9 [
0 -
1
Q J 1
(38)
Since the matrix factorization in Eq. (38) is product of elementary matrices and a diagonal matrix, we can easily compute its inverse as:
- oT
1
- -Q -Q-1
[000
Q-1
0 -J
/ --Q
1
(39)
I
Then multiplying by the right-hand side in Eq. (32) gives the solution:
T00+1)
Q-1JØT T
Q Q-
J Q-1 o Q((00) 11)
16+11
0
I
0
-(J
h(a
Q01
Q-1
Q(a(1) 14)
0
I
-(J0)
-(JQQ-1JOT)
((())-1
(((()) - ( (1(0) 11)
(00
I
1(()
j( - -ni)
0
I
(-1
Q
(J J(T (I(((())))
(M((())))))))))))))))))))))))))) - (( - m))
Using the Schur complement reduces the Newton system from an indefinite (n + q) x (n + solve to n SPD solve with Q-1 and q x q SPD solve with the Schur complement, where q S n. Similarly, the Jacobian expression in Eq. (24), which we will show next, is obtained by implicitly differentiating the linearized KKT conditions and eliminating the dual block, which also avoids the need to invert full (n + 4) x (n + saddle-point or indefinite matrix.
Ou" (ts)/pp of the transformation T(u) = u² (μ)) using implicit differentiation. At convergence, the optimal pair (u", X*) satisfies Eq. (19). Differentiating both sides of the first stationarity equation in Eq. (19) w.r.t j gives:
Jacobian (p)).. Here, we compute the Jacobian ( =
Ru (1:p)== 0u
(Q(u" p1) + Vh(u*)^") =0,
0p
(40)
Ou"
Ox = Q.
(Q+ (Q+v>)(u*)N))
+ Vh(u')
0p
0p
Similar to Eq. (31), we assume h(u") is approximately affine near the optimal point, and we approximate (²)(u() Similarly differentiating both sides of the second feasibility equation in Eq. (19) w.r. gives:
0 h(u²) = Vh(u*) Ou"
Ru X::
0μ = 0.
(41)
Dp
0p
Combining Eq. (40) and Eq. (41) leads the following block linear system:
Q J
J
Ou @8p Q 0
0
0x 0x000 08p
Similar Eq. (38), we can use the Schur complement eliminate the dual term via block substitution. Using the block inverse in Eq. (39) with j(i = J". have
28

---

Published as a conference paper at ICLR 2026
Ou" /0p]
I
q-
I
1
/
0x* 18/
0
I
0
-(
-rQ-1
I
0
0
1
0
-(J'Q-1
---
QQ-1
0 I
(
I 1-Q--JrToN//00 (~'Q-11'))~1
Hence, the Jacobian is given by the first component as:
Ou" (n)
JT(
= Q- (FQ- J
0u
which is the desired Eq. (24).
C.3 (NONLINEAR) CONVEX INEQUALITY CONSTRAINTS
In this subsection, we describe how to compute the DPPL ProbHardE2E for nonlinear convex inequality constraints. Convex inequality constraints arise naturally in many scientific and engineering applications. For example, total variation (TV) regularization widely used to promote smoothness piccewise-constant structure in spatial fields, e.g., image denoising (Rudin et al., 1992; Boyd & Vandenberghe. 2004) and total variation diminishing (TVD) constraints to avoid spurious artificial oscillations in numerical solutions to PDEs (Harten, 1997; Tezaur et al., 2017; Schein et al., 2021). Other common convex constraints include box constraints, which enforce boundedness physical or operational quantities (Bertsekas, 1997).
We consider the constrained projection Problem (7), where Q > 0, and h: R + R R" R' denote smooth functions representing equality and convex inequality constraints, respectively. This is a convex optimization problem due to the strictly convex quadratic objective and the assumption that g(u) is convex. The associated Lagrangian
L(u, v; 2) =(u=)) h()) v. g(((),,
- =)
with Lagrange multipliers À E TR2 for the equality constraints and E IR. for the inequality constraints, where 0. The KKT optimality conditions are given as:
(Stationarity)
Q(u 2) + Vh(u ) + Vg(u")e =
(Primal feasibility)
M(u") = 0, g(u") S 0,
(42)
(Dual feasibility)
v² 2 0,
1s g (u" = for all j = I,....s. Note the first conditions are the same the ones for nonlinear equality constraints with in Eq. (19).
(Complementary slackness)
The KKT conditions in Eq. (42) are necessary and sufficient for optimality, under standard constraint qualifications, e.g., Slater's condition (Boyd & Vandenberghe, 2004). Eq. (42) can be solved by various optimization methods, e.g., stochastic trust-region methods with sequential quadratic pro- gramming (SQP) (Boyd & Vandenberghe, 2004; Hong et al., 2023) and exact augmented Lagrangian (Boyd& & Vandenberghe, 2004; Fang al., 2024). The augmented Lagrangian balances the need for both constraint satisfaction and computational efficiency, which makes it particularly effective in large-scale optimization problems. While the inequality constraints g(u) S 0 are convex by assumption, the equality constraints h(u)) 0 are typically required to be affine to ensure that the feasible set remains convex (Boyd & Vandenberghe, 2004). Nonlinear equalitics generally yield non-convex level sets, which can violate problem convexity even when the objective and inequalities are convex. Although exceptions exist where nonlinear equalities define convex sets, these cases are rare and must be verified explicitly (Bertsekas, 1997; Boyd & Vandenberghe, 2004).
To compute the Jacobian JT == du" (w)/Op of the projection map with respect to the input μ., could, in principle, apply implicit differentiation the KKT conditions in Eq. (42). For general
29

---

Published as a conference paper at ICLR 2026
constrained problems with nonlinear equality and convex inequality constraints, the derivation becomes analytically complex, particularly due to active set variability and non-affine structure. the special case of quadratic programs with affine constraints, OptNet (Amos & Kolter, 2017) provides an explicit expression for the derivatives via KKT conditions. In addition, CVXPYLayers (Agrawal al., 2019) enables gradient-based learning for general convex cone programs by canonicalizing them into a standard conic form. In our implementation, we use CVXPYLayers to enforce the constraints during the projection step. Since CVXPYLayers does currently support full Jacobian extraction or higher-order derivatives, we estimate the variance of the projection map using Monte Carlo methods by applying random perturbations the inputs and computing empirical statistics over repeated forward passes.
D SPECIAL CASES OF PROSHARDE2E
this section, we show applications of ProbHardE22 in two seemingly unrelated but technically re- lated domains: (1) hierarchical time series forecasting with coherency constraints (Rangapuram et al., 2021; Olivares et al., 2024a); (2) solving partial differential equations (PDEs) with global conservation constraints (Hansen et al., 2023; Mouli et al., 2024). Both are special cases of ProbHardE2E with linear equality constraints, and orthogonal (Q 7) and oblique (Q = Σー1) projections, respectively. Fig. 2 illustrates the wide variety of cases that framework covers.
Homo sondact
Predistions
Constraint Type
Linea Commrrana
Pioteti Idd
Deterren
ProbHardEZE
End -----o-ndd
Inferende
Time-sste
Constraint Application
Domain
Figure 2: ProbHardE2E serves as a probabilistic umified framework for learning with hard con- straints.
D.I ENFORCING COMERENCY HIERARCHICAL TIME SERIES FORECASTING
Hierarchical time series forecasting is abundant in several applications, e.g., retail demand forecasting and electricity forecasting. retail demand forecasting, the sales are tracked at various granularities, including item, store, and region levels. In electricity forecasting, the consumption demand is tracked at individual and regional levels. Each time series at time t can be separated into bottom and aggregate levels. Bottom-levels aggregate into higher-level series each time point through known
30

---

Published as a conference paper at ICLR 2026
Total Ot
Region Utz
Region Atg
Region Ota
City
City 4
City 5 bts
City 2 bi
City 6 bss
City 3 bts
be
Figure 3:
Example hierarchical time series structure with at E R. be R and Samm
[
1
0
0
1
1
0
0
0
1
0
0
0
relationships, which can be represented as dependency graphs. Let 2 = [a b,]T E IR", where a, IR9 denotes the aggregate entries, d, IR9 denotes the bottom-level entries, and n q + q. Let Ssum E {0, 1 denote the summation matrix, which defines the relationship between the bottom and aggregate levels a, Stmb This coherency constraint can be equivalently expressed as:
[az bt
[Lg
Ssum]
= Azt = O, Vt.
(43)
where Ig denotes the q x q identity matrix. See Hyndman et al. (2011); Rangapuram et (2021); Olivares et al. (2024b) and the references therein for details, and Fig. 3 for illustration.
HierE2E (Rangapuram et al., 2021) enforces the coherency constraint in Eq. (43) by projecting the multivariate samples 2, onto the null space of the constraint, i.e., Az, 0. It uses the following projection:
u (2+) = (I
A ATA (I A A)
(44)
Pi
where At = (A denotes the right psuedoinverse, and Pi P7 = PJ denotes an orthogonal projector.
We show that HierE2E can be formulated in our ProbHarddEE framework with the following posterior mean and covariance:
PHieEEZ (I (I-A3 A A)μ,
(45a)
c HierEE = ALAE - EALA+ AEAAAA A,
(45b)
where is defined in Eq. (44). In particular, show in Proposition D.1 that the FierEZE posterior update in Eq. (45) is special linear constraint case of our ProbHardE2E method, which uses an orthogonal projection with Q I and b 0.
Proposition D.1. The projected mean and covariance for HierE2E in Eq. (45) is given by the solution to Problem (7) with linear constraints in Proposition C.I, i.e. h(u) = Au = 0, b 0, where Q = for an orthogonal projection and Z (1,, Σ) is a multivariate Gaussian.
Proof. The oblique projection Eq. (16) used in P ProbHardE2E for linear constraints is given as
Pg = - Q A (AQQ 4 A. Setting the expression simplifies to Po P,
T- A (AA A AtA.
31

---

Published as a conference paper at ICLR 2026
The posterior mean j for ProbHardE2E with linear constraints is given in Eq. (18a) with 6 as:
j P
PTH
(46)
AAAA)ı
AiHicEEB:
which is the desired expression in Eq. (45a). Similarly, the posterior covariance for ProbHardE2E in Eq. (18b) is given as:
I = Po- DP-
PIEPT
= PIEPI
= (I - AIA))(II- ATA)
(47)
(I - (-AA))( -AATA)
Σ -'AEE -EA'++++ A'AAAA
Σ HierEE
which is the desired expression in Eq. (45b).
Note that HierE2E does not directly project the distribution parameters, even though a closed form exists, as shown in Eq. (46) and Eq. (47). Instead, it directly projects the samples in Eq. (44). An improvement to HierE2E (that we do in FrobHardE2E) is to eliminate the computationally expensive sampling in the training loop. (See Section 3.5.) HierE2E samples from the parametric distribution generated by DeepVAR (Salinas et al., 2019; 2020; Alexandrov et al., 2019), reconciles these samples, and computes the loss over time using the Continuous Ranked Probability Score (CRPS). Generally, for unknown distributions, the CRPS evaluation requires sampling, which may explain its necessity Gaussian distribution in Hier HierE2E, the CRPS can be computed analytically (Matheson & Winkler, 1976; Taillardat al., 2016) using the mean and covariance of the output distribution.
their framework. For many standard distributions, e.g., the multivariate
D.2 ENFORCING CONSERVATION LAWS IN PDES
In addition to hierarchical forecasting, another (at first -urreltte) application of ProbHardE2E is enforcing conservation laws in solutions to partial differential equations (PDEs). conservation law is given as u, + V F(u) F(u) Veque, 1990). Hansen et al. (2023) propose the ProbConserv method to enforce the integral form of conservation laws from finite volume methods (Le Veque, 2002) as a linear con- straint Au = for specific problems that satisfy a boundary flux linearity assumption. In particular, ProbConserv proposes the following update equations for the posterior mean and covariance matrix:
for unknown unt,, a:) and nonlinear flux function
jppnhConsser F ΣA EA EAEA
b),
(48a)
ErrobConserv = Σ - ES--AT(S (AEAA))
(48b)
given the mean a and the covariance matrix E estimated from black-box probabilistic model, e.g.. Gaussian Process, probabilistic Neural Operators (Mouli et al., 2024) Attentive Neural Process (ANP) (Hansen et al., 2023) or DeeRVAR (Salinas et al., 2019) used in the hierarchical forecasting case.
In Pr ProbConserv, the posterior mean j in Eg. (48a) is shown to be the solution to the constrained least squares problem:
@rProbConserv arg min
2 |1
1-13--1
PER
A
We formulate this optimization problem more generally, and show that by assuming that Z
Mi, Σ) multivariate Gaussian, constrained sample u" (=) Y ~N(i,, $.) in ProbConserv
32

---

Published as a conference paper at ICLR 2026
solution to Problem (7) with Σ-! and linear constraints. In particular, show in Proposition D.2 that the ProbConsery posterior update in Eq. (48) is a special linear constraint case of our ProbHardE2E method, which uses oblique projection with Q E-
Proposition D.2. The projected mean and covariance for ProbConserv in Eq. (48) is given by the solution to Problem (7) with linear constraints in Proposition C.I, i.e., h(u) = =Au = 0, where Q Σ - for an oblique projection and Z (p, 2) is multivariate Gaussian.
Proof. The oblique projection in Eq. (16) used in ProbHardE2E for linear constraints is given as Pa-1 = P3 = I ------A- T (AQ-1AT) - Setting Q Σ-1, we have that Po-1
ΣAT (AEAT)- == Ps. Ps.
The posterior mean ji for ProbHard with linear constraints is given in Eq. (18a)
A d=Pa Pe 11+ (1 - Po 1)AT6
= =I--EA
A)μ+ ( A')
= V1-EA"(AEA -1 A) (ACAT)- AA
= 1 - EAT (Ap (Ap b)
PipotConser::
which is equal to the desired expression in Eq. (48a).
Similarly, the posterior covariance for ProbHardE2E in Eq. (18b) is given as:
Σ = Pq-IEP--
-I- =(1
ΣAT (AEA AAX
=
EA "ATA (I EAT TAA
- 1E-2EAT (AEAT)-AAE) ""AA(""))""AAAAAA
EAT(AA
ProohConsrr
which is equal to the desired expression in Eq. (48b).
Note that the projected distribution parameters in Eq. (48) are applied only at inference time in ProbConserv. In ProbHardE2E, we show the benefits of imposing the constraints at training time as well in an end-to-end manner.
E FLEXIBILITY IN THE CHOICE OF Q AND ITS STRUCTURE
In this section, we discuss the modeling choices for the projection matrix Q in our DPPL, which defines the energy norm in the objective the constrained least squares Problem (7). Its specification significantly influences both the learning dynamics and the inductive biases of the model. Selecting or learning Q offers a principled mechanism to reflect the statistical structure of the data, particularly in settings involving multivariate regression or heteroscedastic noise (Kendall et al., 2018; Stirn et al., 2023). Table summarizes common structure choices for Q and their trade-offs. Of course, many applications, there is single goal for the choice of Q-to to optimize accuracy.
In practice, the space of symmetric positive definite (SPD) matrices too large to be explored (and "learned")) without additional structure, especially in high-dimensionall settings. To address this, structural constraints are often imposed on Q, reducing the number of parameters, and acting as form of regularization (Willette al., 2021). These structures encode modeling assumptions, e.g., output independence, sparsity, or low-rank correlations, and they trade off statistical expressivity against computational efficiency.
In many cases, the choice of Q (or the form of Q) should ideally reflect (knowledge or assumptions hope about) the structure of the underlying data distribution. The simplest choice, Q T,
33

---

Published as a conference paper at ICLR 2026
Structure of Q
Example Form
Merits and Demerits
+ Simplest choice, no parameters
Identity
Q =
Strong regularization
Ignores uncertainty and correlations
Captures heteroscedasticity
Diagonal (learned)
Q = diag( .l.....9... 9n)
Efficient to compute and invert
Ignores correlations
Captures dominant correlations
Low-rank (learned L)
Q = LL L Rexd
+ Fewer parameters than full
Still computationally involved
Fully expressive
Full (learned T.)
Q = T.L., I. Raxa
High memory and compute cost
Prone overfitting
Table 6: Several structure choices for the matrix Q and their associated trade-offs.
assumes isotropy across output dimensions, and is often used for its regularization benefits and ease of implementation. This choice neglects any correlation structure in the data, and it tends to perform poorly in the presence of strong heteroscedassicctyy A diagonal matrix Q = diag .9.............................. !n) introduces per-dimension weighting, and is well-suited to heteroscedastic tasks where the variance differs across outputs (Kendall & Gal, 2017: Skafte et al., 2019). Low-rank approximations provide a compromise between model complexity and expressivity, by capturing dominant correlation directions (Willette et al., 2021). Full-rank matrices allow flexibility and often require strong priors or large datasets avoid overfitting (Weinberger & Saul, 2009).
We focus on two concrete realizations of Q: the identity matrix Q = / that is used in the HierE22 (Rangapuram et al., 2021) (see Appendix D.1), and a diagonal matrix defined as the inverse of predicted diagonal covariance, == Σー1 that is used in ProbConserv (Hansen et al., 2023) (see Appendix D.2), where the model. This latter choice corresponds to heteroscedastic formulation that scales residuals based on their predicted precision, which emphasizes more confident predictions, and down- weinhws less certain ones (Stirn et al., 2023; Le et al., 2005; Hansen et al., 2023).
di . denote the empirical variances output by
F PROOF OF THEOREM 3.1
In this section, we begin by first restating Theorem 3.1, which provides closed-form update for our DPPL in Eq. (8) for a prior distribution that belongs to multivariate local-scale family of distributions; and then we provide its proof.
Theorem 3.1. Let Z ~ F(p, Σ) be random variable, where the underlying distribution F belongs multivariate location-scale family of distributions, with mean p. and covariance Σ; and let T be function with continuous first derivatives, such that JT(H)EJT(H)) is symmetric positive semi-definite. Then, the transformed distribution
T (Z) converges distribution with first- order accuracy to F(ū, E) with mean ğ = T(u) and covariance JT(U)EJJT(p where JT(U) = (T())" denotes the Jacobian of with respect to z evaluated at
Proof. Recall that a family of probability distributions is said to be location-scale family if for any random variable whose distribution belongs to the family ~ F(u, E), then there exists a transformation (c-paaaameterization) of the form
Y d
AZ B,
where denotes a scale transformation matrix, B denotes the location parameter, and denotes equality in distribution.
34

---

Published as a conference paper at ICLR 2026
Let Y TZ) be nonlinear transformation. We calculate the first-order Taylor series expansion to linearize the function about the mean μ us:
Y (Z) T(p) + (m) 1u) JT (1) (T(u)- - JT(()))
(49)
Then, since Z belongs to the location-scan scale family distributions, the linearization of Y F(ù, E) also belongs to the family with mean jú and covariance Σ, which we compute below. Taking the expectation of both sides of Eq. (49) we get:
A = (T) (1) + (Z p4]
E[T(p)) E[J (Z [4)] (by linearity of expectation)
= =(w)) + (T() [E - μL) (since μi is not a random variable)
= T(u)-
(50)
Then, the covariance E. is given as:
Σ B[( - E[T()))) -[T(Z)))
=E((TZ)) (T()) [(x)) (by Eq. (50))
= =R((T(n) + JT
+ (T) (T()))
1) -T())) (by Eq. (49))
- =[(1-())2--))--0)))1 -
=JT(W)E((Z - (Z (((==))]JT(w)
=Ir()EIr(()).
Importantly, the approximation error between the nonlinear transformation and its linearization converges zero in probability (Van der Vaart, 2000), which ensures the validity of this approach asymptotically. We note that this result closely related to the Multivariate Delta Method (Casella Berger, 2001). which shows that for nonlinear function T, the sample mean of T(=1. 2n,) also converges distribution, under mild conditions. Specifically, if the sample mean of n i.dd draws from Z converges multivariate Gaussian (by the CLTT), then the same linearization argument and Slutsky's theorem imply that the sample mean of the projected samples converges to multivariate Gaussian, with parameters given in Eq. (8). Second-ordes approximations (via quadratic expansion of T) yield higher-order corrections, and can lead to non-Gaussian outcomes (e.g., chi-squared) (Casella & Berger, 2001).
G BENCHMARKING DATASETS
In this section, we detail the benchmarking datasets in both applications domains, i.e., PDEs and probabilistic time series forecasting.
G.1 PDES
We consider series of conservative PDEs with varying levels of difficulties, where the goal is to learn an approximation of the solution that satisfies known conservation laws. We follow the empirical evaluation protocol from Hansen et al. (2023). The PDEs we study are conservation laws, which take the following differential form:
U + V F(u)= O,
(51)
for some nonlinear flux function F(u). These equations can be written in their conservative form as:
dt u(t, x)d = F(u(t,ao)) (anpaa)))FF(((((( F(u(t,
(52)
35

---

Published as a conference paper at ICLR 2026
by applying the divergence term in ID over the domain Q = [zo. (LeVeque, 1990; Hansen et al., 2023). This global conservation law states that the rate of change of total mass or energy in this system is given by the difference of the flux into the domain and the flux out of the domain. Note that in higher dimensions, the flux difference on the right-hand side of Eq. (52) can be written as a surface integral along the boundary of the domain. This conservative form is at the heart of numerical finite volume methods (LeVeque, 2002), which discretize the domain into control volumes and solve this equation locally in each control volume, to enforce local conservation, i.c., so that the flux into control volume is equal to the flux out of it. In the following, we summarize the PDE test cases with their initial and boundary conditions, exact solutions, and derived linear conservation constraints from Hansen et al. (2023).
G.1.1 GENERALIZED POROUS MEDIUM EQUATION (GPME)
The Generalized Porous Medium Equation (GPME) is given by the following degenerate parabolic PDE:
Ue V (k(u) Vu) 0,
(53)
where the flux in Eq. (51) is given as F(u) = kku and k(u) denotes the diffusivity parameter. This diffusivity parameter k(u) may depend nonlinearly and/or discontinuously on the solution u. We consider three representative cases within the GPME family, by changing this parameter k(u). Each instance of the GPME increases in difficulty based on the regularity of the solution and the presence of shocks or discontinuities.
Heat Equation ("Easy"") The classical parabolic heat equation arises when the diffusivity is constant, i.e., k(u) = k in Eq. (53). We use the heat equation with the following sinusoidal initial condition and periodic boundary conditions from Krishnapriyan et al. (2021); Hansen et al. (2023):
dt kAu, u(0, x)) = sin(x),
Vr Ω = [0, 2T], Vt E [0,
Vx E [0, 2T),
(54)
(t, 0) u(2x,t) t),
Vt [0, 1].
respectively. The exact solution, which can be solved using the Fourier Transform, is given as:
Mexact (t, x) e sin(x).
The solution is smooth sinusoidal curve that exponentially decays dissipates over time, and has an infinite speed of propagation. With these specific initial and boundary conditions in Eq. (54), the global conservation law in Eq. (52) reduces the following linear equation:
2x
u(t, x) dr = O, Vt E [0. 1],
(55)
since the net flux on the boundaries is O.
Porous Medium Equation (PME) ("Medium")). The PME is a nonlincar degenerate subclass of the GPME, where the diffusivity is a nonlinear, monomial the solution, i.e., k(()) u in Eq. (53). It has been using in modeling nonlinear heat transfer (Vázquez, 2007; Maddix et al., 201 8a). We use the PME with the following initial condition and growing in time left Dirichlet boundary condition from Lipnikov et al. (2016); Maddix et al. (2018a); Hansen et al. (2023):
Ut V (u" Vu) = O,
Væ E Ω = [0, 1], Vt E 1],
u(0, x) O,
Vz [0, 1].
(56)
u(t, 0) = (mt) 1/m
Vt E [0, 1].
The exact solution is given as:
Mexaa (t,x) (m - RoLU((=2))1// 1/m
For small values of k(u), this degenerate parabolic equation behaves hyperbolic in nature. The solution exhibits a sharp front at the degeneracy point = x with a finite speed of propagation. With these specific initial and boundary conditions in Eq. (56), the global conservation law in Eq. (52) reduces to the following linear equation:
S u(t,xx) dr (ml))-1/m
Vt E (0. 1].
(57)
36

---

Published as a conference paper at ICLR 2026
Stefan Equation ("Hard"). The Stefan equation has been used in foam modeling (van der Meer et al., 2016) and crystallization (Sethian & Strain, 1992), and models phase transitions with the following discontinuous diffusivity:
a 2 u² k(u)= = 0, a < ur
1,
2 0,
in Eq. (53). We use the Stefan equation with the following initial condition and Dirichlet boundary conditions from Maddix et al. (2018b); Hansen et (2023):
V (R(u) Vu)
Va f [0, 1]. e [0,
u(0,5 x) = O, Vo E [0, 1].
(58)
ut, 0) V [0, 1].
The exact solution is given
- "²* erf 2VE
(t, = 1 1uzw"
erf(a/2)
where denotes the indicator function, erf(z) = (2/v lo exp(-1²)er dyy denotes the error function, and 2á and á satisfies the following nonlinear equation:
- n² u² erf(ā)āe**
The solution is a rightward moving shock. With these specific initial and boundary conditions in Eq. (58), the global conservation law in Eq. (52) reduces the following linear equation:
L. wlt, x) 2(1 w²) -(1=4) V E [0,
(59)
crf(a/2)
G.1.2 HYPERBOLIC LINEAR ADVECTION EQUATION
The hyperbolic linear advection equation models fluids transported at a constant velocity, and is given by Eq. (51) with linear flux Fuu) = Su. We use the linear advection problem with the following step-function initial condition and inflow Dirichlet boundary conditions from Hansen et al. (2023):
+ fu. 0,
Vr E Ω = [0, 1], Vt E [).
00, x) 1.50.5.
Ve E [0, 1],
(60)
u(t, 0) 1,
Vt E [0, 1].
The exact solution is given as:
u(x = - Bt),
BBt)
where hxx) 1, s.s denotes the initial condition. The solution remains a shock, which travels to the right with finite speed of propagation B. With these specific initial and boundary conditions in Eq. (60), the global conservation law in Eq. (52) reduces to the following linear equation:
L. t((
2 + St,
(61)
which shows that the total mass increases linearly with time due to the fixed inflow.
G.2 PROBABILISTIC TIME SERIES FORECASTING
In addition to PDEs, we also evaluate ProbHardE2E on five hierarchical time series forecasting benchmark datasets, where the goal is to generate probabilistic predictions that coherent with known aggregation constraints across cross-sectional hierarchies (Rangapuram et al., 2021).
Table 7 provides overview of the time series datasets used in our empirical evaluation. For each benchmarking dataset, details the total number of series, the number bottom level series (i.c., the leaf nodes in the hierarchy), the number of series aggregated from the bottom-level series, the depth
37

---

Published as a conference paper at ICLR 2026
of the hierarchy in terms of the number of levels, the number of time series observations, and the prediction horizon т.
We adopt the same dataset configurations in Rangapuram et al. (2021), from which hicrarchical forecasting benchmarks and pre- pre-processing pipeline. These datasets are available in Glu- onT'S package (Alexandrov et al., 2019). The LABOUR dataset (Australian Bureau of Statistics, 2019) contains monthly Australian employment statistics from 1978 to 2020, organized into 57-series hi- erarchy. The TRAFFIC dataset (Ben Taieb & Koo, 2019) includes sub-hourly freeway lane occupancy data, aggregated into daily observations forming a 207-series structure. TOURISM (Tourism Australia, Canberra, 2005) consists of quarterly tourism counts across 89 Australian regions (1998-2006) and the extended TOURISM-L dataset (Wickramasuriya et al., 2019) comprises 555 grouped series based on both geography and travel purpose. Lastly, WIKT contains daily page view counts from 199 Wikipedia pages collected over two years (Anava et al., 2018).
the
Table 7: A summary of the time-series datasets. TOURISM-L has two hierarchies, defined by geography and travel purpose; consequently, it has different numbers of bottom series and different depths in each hierarchy.
Dataset
Total
Bottom
Aggregated
Levels
Obs.
Horizon т
Frequency
TOURISM
89
56
33
36
8
Quarterly
TOURISM-L
555
76; 304
175
4;
5
228
12
Monthly
LABOUR
57
514
8
Monthly
TRAFFIC
207
200
4
366
1
Daily
WIKI
199
150
49
366
1
Daily
H IMPLEMENTATION DETAILS
In this section, we provide the implementation details ProbHardE2E. Fig. illustrates the overall pipeline of P FroHardEEE which integrates probabilistic modeling, constraint enforcement, and loss-based calibration into a unified differentiable architecture. The core contribution lies in the DPPL, which acts as a "corrector" to the "predictor," which is the unconstrained distribution predicted by wide class of models. Conceptually, this layer parallels classical predictor-corrector and primal-dual methods from numerical optimization (Boyd & Vandenberghe, 2004; Bertsekas, 1997), where candidate solution is refined to satisfy known constraints before evaluation.
We evaluate ProbHardE25 on two scientific domains: (1) PDEs, where structured physical con- straints, e.g., conservation laws and boundary conditions, must be enforced (see Appendix G.I), and (2) probabilistic hierarchical time series forecasting, where aggregation coherency is required (see Appendix G.2). We show that ProbHardE2E is aggosti by using a base probabilistic model (predictor) from each application domain, i.e., Var ainneO (Mouli al., 2024) for PDEs and DeepVAR (Salinas et al., 2019) for forecasting. We then enforce the corresponding constraint with our DPPL (corrector). We provide the experimental details for each application in the following subsections.
H.1 PDES
All the experiments are performed on single NVIDIA V100 GPU. We use probabilistic Fourier Neural Operator (FNO) (Li al., 2021), i.e., VarianceNO (Mouli al., 2024) to learn mapping from PDE parameters to solutions, e.g., the diffusivity mapping k(u) + (f,,x) in the (degenerate) parabolic Generalized Porous Medium Equation (GPME), or the velocity mapping 3 + uu(, x)) in the hyperbolic linear advection equation. (See Appendix G.1 for details on the datasets.)
H.1.1 DATASET GENERATION
Table provides an overview of the PDE data generation. For each PDE in Appendix G.I, generate dataset of N 200 parameter-solution pairs ( input PDE parameters, e.g., ht m, u", 8, and (0)) denotes the corresponding spatiotemporal solution
w() Jez D, where の(0) denotes the
38

---

Published as a conference paper at ICLR 2026
User-Specified Constraints
h(w)) 0,
s(u) S 0
(010,u10) D
Differentiable Probabilistic Projection Layer (DPPL)
Base Probabilistic
Zo((
Yo(
CRPS Loss Computation ((Volol0),
(d19)
fs((.E)
DPPL(Z.(0 b())
Unconstrained Random Variable
Constrained Random Variable
(Predictor)
Projection Matrix
(Corrector)
End-to-end training
Figure 4: Schematic representation of ProbHardE25 (see Algorithm 1). Here, known pathwise- differentiable probabilistic model is chosen to predict a (unconstrained) prior distribution. (Optionally, the projection matrix be specified as a part of the prediction from the probabilistic model or modeled separately.) Next, we transform the distribution with our DPPL to obtain the transformed distribution, done empirically via the Delta Method (see Section 3.3), which enforces the constraints. Lastly, we choose an appropriate loss function, e.g., CRPS, to calibrate the transformed distribution with the target variable.
field. Each solution (+) (t, x is simulated over a grid of 100 equidistant points both space and time, yielding total of 100 x 100 observations per instance. During evaluation, predict the final 20 equidistant time slices while conditioning on the earlier time steps.
Table 8: Overview of PDE dataset generation. Each dataset contains 200 samples with a fixed 160/40 train-test split.
PDE
Parameter range
Spatial domain
Time domain
Train/ (%)
Heat PME
h [1,
0, 2T)
0. 1
80/20
m E [2, 3]
0, 1
0,
80/20
Stefan
u* E [0.6, 0.65]
1
0,
80/20
Linear Advection
B [1
0,
0,
80/20
H.1.2 ARCHITECTURAL DETAILS
We VarianceNO (Mouli al., 2024) base unconstrained probabilistic model. VarianceNO an augmented Fourier Neural Operator (FNO) (Li et al., 2021) that updates the last layer output two prediction heads instead of one, i.e one for the mean and the other for the variance of the multivariate Gaussian distribution. Table 9 details the model hyperparameters.
H.1.3 TRAINING AND TESTING SETUP
We follow the standard training procedure for FNO-based models as proposed by Li et al. ()021). Specifically, we use the Adam optimizer (Kingma & Ba, 2015) with weight decay and train using mini- batches of fixed size B 20. A step-based learning rate scheduler is applied, which reduces the learning rate by half every 50 epochs. During evaluation, we uniformly sample parameters from the specified parameter ranges in Table 8 to construct test and compute the evaluation metrics.
39

---

Published as a conference paper at ICLR 2026
Table 9: Hyperparameters for the VarianceNO model.
Hyperparameter
Values
VarianceNO
Number of Fourier layers
4
Channel width Number of Fourier modes
{32, 64)
12
Batch size
20
Learning rate
{10-4, 10-2.
10-2)
H.2 PROBABILISTIC TIME SERIES FORECASTING
The experiments are performed on an Intel(R) Xeon(R) CPU B5-2603 v4 @ 1.70GHz.
H.2.1 DATASET GENERATION
We adopt the hierarchical forecasting benchmarks and preprocessing pipeline introduced in Ranga- puram al. (2021), using five standard datasets: LABOUR, TRAFFIC, TOURISM, TOURISM-L, and WIKI. Each dataset contains a hierarchy of time series with varying depth and number of aggregation levels (see Table 7). The train/test splits, seasonal resolutions, and prediction horizons follow the standardized setup provided in Rangapuram et al. (2021).
H.2.2 ARCHITECTURAL DETAILS
We use DeepVAR (Salinas al., 2019) base unconstrained probabilistic model, which is aligned with Hier-E2E. DeepVAR is a probabilistic autoregressive LSTM-based model that leverages multivariate Gaussian distribution assumption on the multivariate target. DeepVAR models the joint dynamics of all the time series in the hierarchy through latent temporal dependencies, and outputs both the mean and scale of the predictive distribution, by optimizing the negative log likelihood (NLL). Our ProbHardE22 model in the time series application is developed based on Hier-E2E in Gluon'SS (Alexandrov al., 2019). We use default base model architecture DeepVAR, and make further modifications to H errE2... Specifically, we tune the hyperparameters Table 10, and adjust the loss CRPS for structured probabilistic evaluation. We disable sampling- based projection during training because Pr PEE2 optimizes the closed-fiee from CRPS for Gaussian distributions, and our projection methodology ensures that linear constraints are met probabilistically. During inference, report CRPS through samples, in order align the evaluation definition with the various hierarchical forecasting baselines.
Table 10: Key hyperparameters for DeepVAR across hierarchical forecasting datasets.
Dataset
Epochs
Batch Size
Learning Rate
Context Length
No. Prediction Samples
LABOUR
32
0.01
24
400
TRAFFIC
10
32
0.001
40
400
TOURISM
32
200
TOURISM-L
10
4
0.001
36
200
WIKT
25
32
0.001
15
200
H.2.3 TRAINING AND TESTING SETUP
We follow the standard GluonTS (Alexandrov al., 2019) training setup using the Adam opti- mizer (Kingma & Ba, 2015) and mini-batch updates. Each epoch consists 50 batches, with batch size set according Table 10. We run evaluation five times and report the mean and variance of the CRPS values in Table 3.
Unlike Hier-E2E (Rangapuram et al., 2021), which samples forecast trajectories during training and projects them to ensure structural coherence on samples, our method operates entirely in the parameter space during training. We avoid sampling and instead minimize the closed-form CRPS
40

---

Published as a conference paper at ICLR 2026
loss (Gneiting et al., 2005) directly the predieted mean and variance. This makes the training process sampling- free and reduces training time, similar to the PDE case discussed later in Figure 5. This key distinction avoids the use of coherent train samples, described in the Appendix of Rangapuram et al. (2021).
At inference time, because the reported CRPS is computed on the samples in the hierarchical baselines, we enable structured projection by drawing predicted samples from the learned distribution, and apply our DPPL ensure that they satisfy the hierarchical aggregation constraints. This setup parallels the coherent-pree ptr.dttped_aaaplll.. with this approach for experimentation simplicity. Table 10 shows the number of prediction samples in evaluation to compute the CRPS and calibration metrics over the projected outputs. Alternatively, we can also evaluate the CRPS using samples from the projected distribution.
mode in HierE2E, and we implement the inference step
H.3 METRICS
We evaluate P robHardEEE, and the various baselines using the following metrics. We denote the exact solution or ground truth observations u, and report the metrics the mean fi, covariance Σ, and samples {u, J.N. drawn from the constrained multivariate Gaussian distribution ((, Σ).
Mean Squared Error (MSE). The MSE measures the mean prediction accuracy and given as:
MSE(A)
||u Aille
where the Frobenius norm is taken over all the datapoints n in p.
Constraint Error (CE), The CE measures the error in the various equality constraints h(u") = 0. i.e., conservation laws for PDEs and coherency for hierarchical time series forecasting, on the samples, and is given as:
CE(u*)
Σ [
where compute the average error over N 100 samples {u] fi=1
Continuous Ranked Probability Score (CRPS). The CRPS (Gneiting Raftery, 2007) measures the quality of uncertainty quantification by comparing a predictive distribution to a ground-truth observation. For a multivariate Gaussian distribution with independent components (j,, diag(o?)). where o denotes the i-th diagonal entry the predictive covariance Σ, the CRPS is given in closed-form as:
CRPSN(U,G;U)
0 24 (2P(z)) 22(zz)
Vi
where z; (w: - pu)/on. p(z)) (L/V2TT) (xp(=)xxpppppppp) exp(-22/22
denotes the standard normal PDF, and P(z) - ply)dy denotes the standard normal CDF (Gneiting al., 2005; Taillardat al., 2016).
I ADDITIONAL EMPIRICAL RESULTS AND DETAILS
In this section, we include additional empirical results and details for ProbHardE2E with various constraint types, i.e., linear equality, nonlinear equality and convex inequality.
1.1 LINEAR EQUALITY CONSTRAINTS
In this subsection, we show additional results and details for ProbHardE22 with linear equality constraints in both PDEs and hierarchical time series forecasting.
41

---

Published as a conference paper at ICLR 2026
1.1.1 PDES WITH CONSERVATION LAW CONSTRAINTS
Here, we impose the discretized form of the simplified linear global conservation laws given in Appendix G.1 for the heat equation (Eq. (55)), PME (Eq. (57)), Stefan (Eq. (59)) and linear advection equation (Eq. (61)). We use the trapezoidal discretizations of the integrals from Hansen et al. (2023).
Fig. 5 shows the analogous training time per epoch to Fig. ((a) for PDE datasets. Models trained with 100 posterior samples per training step incur a 3.5-3.6x increase in epoch time relative to our ProbHard rdE2E, which avoids sampling altogether by using a closed-form CRPS loss. See Table 2 for the accuracy results.
No Sampling
With Sampling (100)
1.50
19
1.25
0.75
0.2
0.00
PME
Advection
Figure 5: ProbHardE225 PDE timing comparisons for our sampling-free approach.
1.1.2 HIERARCHICAT. TIME SERIES FORECASTING WITH COHERENCY CONSTRAINTS
Here, we test ProbHa rdE2E on probabilistic hierarchical forecasting with coherency constraints. (See Appendix G.2 for details and Table for the results.) We compare the two variants of I-1 (ProbiardE2E-OOs and with orthogonal Q
ProbHardE2E, i... with oblique (ProbiaarEE2EOo)) projection to the following baselines:
DeepVAR (Salinas et al., 2019) is the base unconstrained probabilistic model, which assumes multivariate Gaussian distribution for Z ~ (u,, Σ) with mean μ and diagonal covariance Σ.
Hier-E2E (Rangapuram et al., 2021) uses DeepVAR as the base model, and enforces the exact coherency constraint by applying the orthogonal projection directly the samples in an end-to-end manner. Another difference from their approach is that we use the closed-form
CRPS expression rather than the approximate weighted quantile loss.
Pr ProbConserv (Hansen et al., 2023) enforces the coherency constraint as an oblique projection at inference time only.
ARIMA-Naive NeBB and CTS-Nal veBU are two simple baseline models that use ARIMA and exponential smoothing (ETS), respectively. These methods use a naive bottom-up approach of deriving aggregated level forecasts (Hyndman & Athanasopoulos, 2018).
PERMBU-MINE (Taieb et al., 2017) a hierarchical probabilistic forecasting model that is based on a linear projection method MINT (Wickramasuriya et al., 2019). It generates probabilistic forecasts for aggregated series using permuted bottom-level forecasts.
We do not include DPMN (Olivares al., 2024a) CLOVER (Olivares et al., 2024b) in our experi- ments because the implementations are proprietary. Given that Hier-E2E is the best open-access hierarchical forecasting model, through GluonTS (Alexandrov et al., 2019), to the best of our knowl- edge, we use the same base model to Hier-E2 E2E (i.e., DeepVAR), and we evaluate forecast accuracy compared Hier-E2E assess the added value of our P ProbHardE25
accessc
I.2 NONLINEAR EQUALITY CONSTRAINTS
In this subsection, we impose the discretized form of the general nonlinear global linear conservation laws from Eq. (52) in Appendix G.1 for the PME with various ranges for the parameter m. (See Table 4 for the results and Fig. 1(b) for the solution profile.) For the PME, the flux in Eq. (52) is
42

---

Published as a conference paper at ICLR 2026
given F(u) VV, where k(u)= = um Substituting this flux into Eq. (52) and integrating in time gives the general conservation law for the PME as:
J ut,, ()tt
/
[u" ro (t, u" EN (N
vt E [0, 1].
Similar to the linear equality constraint case, we discretize the integral using the trapezoidal rule. Unlike P ProbConserv (Hansen et al., 2023), which requires an analytical flux expression to evaluate the right-hand side, our ProbHardE28 can enforce arbitrary (nonlinear) conservation laws directly. In addition, P ProbHardE2E with nonlinear constraints can be applied arbitrary PDEs with any initial or boundary conditions. We impose the initial and boundary conditions as additional linear constraints and enforce positivity the parameter m, i.c., m E [2, 3/. [3, 4] and [4, 5]. As the exponent m is increased, the degeneracy increases, and as a result the solution becomes sharper and more challenging to solve (Maddix et al., 2018a; 8a; Hansen et al., 2023). We see in Table 4 that across all values of m, either our oblique ProbHar rew. or orthogonal projection ProbHardE2E-O0 variants of our method perform better than all the baselines.
the solution. We test on various training and testing ranges for
1.3 (NONLINEAR) CONVEX INEQUALITY CONSTRAINTS
In this subsection, we impose a convex total variation diminishing (TVD) inequality constraint. (See Fig. 1(c) for the solution profile.) TVD numerical schemes have been commonly using in solving hyperbolic conservation laws with shocks to minimize numerical oscillations from dispersion (Harten, 1997; Le Veque, 1990; Tezaur et al., 2017). The total variation (TV) is defined in its continuous form as:
Ou
TV(u(t,))
dl2 dz
This integral can be approximated as the discrete form of the total variation (TV) used in image processing
TV(u(t)) TV(u(t,;))) Ju(t, x++1) utt, mx)1,
(62)
where we discretize the spatial domain f [Z1 .N. into N gridpoints. A numerical scheme is called TVD if:
TV(u(tn+1)) S (V(((t))), Vn = 1n...... N₁ where we discretize the temporal domain [0, 7] [,,,...;AA, AN.] into N, gridpoints, and TV denotes the discretized form defined in Eq. (62).
(63)
The TVD constraint in Eq. (63) is nonlinear inequality constraint, and enforcing it hard constraint is challenging with current frameworks, e.g., DCL (Agrawal et al., 2019). To address this, we perform a convex relaxation of the constraint by imposing:
N
N₂ ju(tn, T++1) - u()n,
TVD
n=1
as regularization term. This approach is analogous to total variation denoising in signal processing (Rudin et al., 1992; Boyd & Vandenberghe, 2004).
Fig. 1(c) demonstrates the application of the modified TVD constraint, resulting in more physically- meaningful solutions by decreasing both the artificial oscillations and probability of negative samples, which violate the monotonicity and positivity properties of the true solution, respectively. addition, ProbHardE2E leads to improved (tighter) uncertainty estimates.
J THE USE OF LARGE LANGUAGE MODELS (LLMs)
LLMs are used for grammatical corrections and minor formatting of the paper. They not used in the conceptualization and implementation research.
43


Reviews:
[{"review_id": "UPtTz3V807", "reviewer": "Reviewer_LXk5", "summary": "This paper proposes an end-to-end probabilistic framework that strictly enforces hard constraints during training while producing full predictive distributions (mean + covariance) for uncertainty quantification (UQ). The technical core is the Differentiable Probabilistic Projection Layer (DPPL), which projects an unconstrained probabilistic prediction onto the constraint manifold. For equality/linear constraints, the projection/transformation has closed forms. For nonlinear/convex constraints, the DPPL layer uses Newton/KKT iterations to compute the projection and implicit differentiation (linearized KKT systems) to get the Jacobian for covariance propogation. Then, using CRPS and its closed forms for common local-scale families, it further enables sample-free and end-to-end training, which is faster compared to the sampling-based UQ baselines. The method is shown to imporve forecasting/UQ metrics on several benchmarks in hierarchical time series forecasting and PDE solving problems.", "questions": "1. Is it possible to provide a synthetic analysis of the effect of the first-order approximations on the results, or compute the approximation errors? Have you considered other methods than the linearized KKT, such as Gauss-Newton augmented Lagrangian for nonlinear constraints? The experimental results are not strong enough to support the superior performance of the proposed method (details in the 2nd question below), so maybe relaxing some of the approximations would improve the peformance.\n \n2. The paper overclaims its performance based on the experimental results in some places. When analyzing these results,  sometimes it is hard to relate the numbers in the tables to the conclusions stated in the paper. For example, \n   - When addressing Q1, the paper stated \"Specifically, when measured against two accuracy metrics across four PDE datasets in Table 2, our method with either oblique (ProbHardE2E-Ob) or orthogonal (ProbHardE2E-Or) projection consistently outperforms both ProbConserv…\". However, in Table 2, when comparing the results on Heat, the MSEs for ProbHardE2E-Ob/ProbHardE2E-Or trained with CRPS are higher than ProbConserv trained with CRPS (i.e., 0.036 and 0.031 compared to 0.027); when comparing the results on PME, the MSEs for ProbHardE2E-Ob/ProbHardE2E-Or trained with CRPS are higher than the other methods trained with CRPS (i.e., 9.59 and 9.01 compared to 8.801, 8.187 and 7.945). These results conflict with your statement that your method \"consistently outperforms\" the other methods. \n   - In Q3, the target is to compare performance using CRPS as training objective to that using NLL, but Table 3 does not provide such experiments or information. I do not see the improvement upon HierE2E addressed Q3, either.\n   - In Section 4.2.1, the paper stated \"We see an even larger MSE performance improvement of ≈ 15 − 17× when trained on CRPS, and CRPS performance improvement of ≈ 2.5× over the various baselines\". This statement only holds true for the case when $m\\in[2,3]$, but not for the other two cases. \n   - In Figure 1(b), it seems to me that the uncertainty estimates from the three methods are roughly the same (maybe the blue region is a bit wider than the others; it is hard to tell the difference between the red and purple regions). In Figure 2\\(c\\), it seems that all the three curves perform roughly the same, none of them captured the ground truth curve.\n \n3. Typos/suggested corrections in the manuscript:\n   - In Algorithm 1, it would be more clear if you indicate what $\\mathbf{Y}_\\theta(\\phi)$ and $\\mathbf{Z}_\\theta(\\phi)$ represent (constrained/unconstrained random variable), even though you have explained them in the main text. In addition, I think at the end of step 7, \"where $z_\\theta(\\phi)\\sim\\mathbf{Z}_\\theta(\\phi)$\" should be moved to the end of step 8.\n   - The paper frequently mentioned \"Problem 7\" while I think it would be better to use \"Problem (7)\" as it appears. \n   - In line 322-323 on page 6, the classical statistical approaches should be listed separately instead of being listed under \"(ii) HierE2E\".\n   - Be consistent with using \"HierE2E\" and \"Hier-E2E\" (the paper is using both).\n   - In line 424, \"approximating\" should be \"approximately\".\n   - In lines 1327-1332, $\\nabla L_{\\hat{u}}$ should be $\\nabla_{\\hat{u}}L$.\n   - In line 1448, in the last vector, $-J^{(*)}$ should be $J^{*}$.\n   - In Eq. (44), add definition of $P_I=I-A^{\\top}(AA^{\\top})A$ (as you mentioned $P_I$ later but it was undefined).", "limitations": "", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "1. The main nolvety of this paper is to propose an end-to-end enforcement of hard constraints, where constraints are optimized as part of learning, not only enforced at post-processing or inference as in convention.\n\n2. It also enables joint UQ and constraint satisfaction, thus suitable for safety/physics/constrained engineering tasks.\n\n3. The sample-free & closed-form CRPS, combined with analytic projection propagations gives substantial training speedups.\n\n4. Generality: the framework covers linear equalities, nonlinear equalities and convex inequalities (with iterative solvers for the nonlinear/convex cases). In addition, it clearly demonstrates how the method works in each scenario with examples and application areas provided.", "weaknesses": "1. The major concern is that the method relies on first-order approximations to linearize the nonlinear function transformation, the KKT system, and the constraints, which together may bring too much estimation errors. Specifically, the covariance propagation relies on a first-order approximation of the Jacobian, thus the UQ under strong nonlinearity constraint can be misestimated. Also the linearized KKT system is only locally valid and can fail to converge or even diverge if the constraint is highly nonlinear or the initial point is far from feasible. \n\n2. The closed-form sample-free CRPS loss depends on using tractable distribution families (e.g., Gaussian). For multimodal/complex posteriors, it still requires sampling or approximations, thus the advantage of computational efficiency is not guaranteed.\n\n3. In some places, the experimental results are not consistent with the claims made in the paper.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "FFpKDAAk2y", "reviewer": "Reviewer_8rdZ", "summary": "The paper proposes ProbHardE2E, a unified probabilistic framework for learning under hard constraints with uncertainty quantification. The core module DPPL turns an unconstrained probabilistic predictor (e.g., DeepVAR for time series or VarianceNO for PDEs) into a constraint-satisfying probabilistic model in an end-to-end pipeline. Training is sample-free: the model predicts location–scale parameters $(\\mu, \\Sigma)$ and DPPL projects these to constrained parameters $(\\hat{\\mu}, \\hat{\\Sigma})$ via a delta-method linearization of a projection map $T$, and the paper optimizes CRPS on the constrained distribution. At inference, they sample from the prior and exactly project each sample by solving a constrained least-squares problem to enforce feasibility. The approach aims to handle linear equalities, nonlinear equalities, and convex inequalities within one recipe.", "questions": "1. Can you provide error bounds (or at least empirical error studies) for the delta-method approximation vs exact sampling across varying nonlinearity/variance.\n2. Can you report multivariate calibration checks (e.g., energy score/logS), not just diagonal-CRPS, given the fact that you are doing multivariate forecasting, where CRPS is not strictly proper. \n3. Include ablations on $Q$ (identity vs diagonal vs low-rank) with runtime/accuracy trade-offs, since $Q$ defines the projection geometry. \n4. The claim \"We demonstrate the importance of using a strictly proper scoring rule for evaluating probabilistic predictions, e.g., the CRPS, rather than negative log-likelihood (NLL).\" is completely wrong, NLL is strictly proper.", "limitations": "", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "- The method is clearly formulated in a principled “predictor–corrector” view.\n- The authors provide an exact handling for linear constraints.\n- The training objective is close-formed and sampling free.", "weaknesses": "- First-order DPPL approximation with no error bounds (risk under strong nonlinearity/variance). \n- Loss/evaluation assume independence while the projected posterior is generally correlated. \n- Practical experiments lock $Q$ to diagonal/identity, limiting the benefits of oblique projections and potentially biasing outcomes. \n- Nonlinear equality projection claims “strict feasibility,” but non-convexity and active-set issues are largely handled numerically without theoretical or robustness analysis. \n- Inequality-constraint metrics are thin (CE defined for equalities), so support for those claims is more qualitative than quantitative.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "0SBGP0SgIY", "reviewer": "Reviewer_971x", "summary": "The authors present ProbHardE2E, a probabilistic forecasting framework that incorporates hard operational and physical constraints while providing uncertainty quantification. Their methodology introduces a novel differentiable probabilistic projection layer (DPPL) that can be integrated with a wide range of neural network architectures. The DPPL enables the model to learn the system in an end-to-end manner, in contrast to other approaches where constraints are enforced only through post-processing or during inference.\n\nProbHardE2E optimizes a strictly proper scoring rule without making any distributional assumptions about the target, allowing it to produce robust distributional estimates. This stands in contrast to existing likelihood-based methods, which are often biased by their distributional assumptions and model choices. Moreover, the framework can incorporate a variety of non-linear constraints, thereby enhancing modeling power and flexibility.\n\nThe authors demonstrate the effectiveness of ProbHardE2E in learning partial differential equations with uncertainty estimates and in probabilistic time-series forecasting, establishing it as a broadly applicable general framework that bridges these seemingly disparate domains.", "questions": "please refer to the weakness", "limitations": "", "rating": 8, "confidence": 3, "soundness": 4, "presentation": 3, "contribution": 4, "strengths": "1. The proposed method uses strictly proper scoring rules (e.g., CRPS) instead of log-likelihood objectives, reducing the learning bias caused by incorrect distributional assumptions.\n2. The training process can be sample-free, offering potential efficiency advantages.\n3. In both time-series forecasting and PDE-solving tasks, the proposed method demonstrates strong empirical performance.", "weaknesses": "1. During inference, the authors indicate that a projection must be computed at every step. Does this imply that the computational cost during inference could be substantial? Could the authors provide an analysis of the time complexity for both training and inference? Similarly, although the training process is sampling-free, repeatedly computing the Jacobian matrix can also increase computational time, especially when dealing with high-dimensional data or scenarios involving multiple constraints.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "uO2x1QlzzE", "reviewer": "Reviewer_MxEb", "summary": "The paper presents *ProbHardE2E*, an end-to-end probabilistic prediction framework that incorporates hard constraints. The approach first trains an unconstrained probabilistic backbone using the CRPS, chosen for its robustness and existence properties compared to the NLL. The main contribution is the *differentiable probabilistic projection layer (DDPL)*, which projects unconstrained predictions onto constrained probabilistic outputs (mean and covariance) during training and onto corresponding random samples during inference.  \n\nFor a location–scale predictive family, the authors derive closed-form solutions of the push-forward through the optimal constrained objective under linear, nonlinear, and convex constraints, allowing sample-free training. The DDPL thus produces uncertainty-aware predictions that satisfy hard constraints. Numerical experiments across time-series and PDE tasks demonstrate improved accuracy and uncertainty metrics over several baselines.\n\nOverall, this is a well-written and valuable paper with a clear contribution, sound theoretical derivations, and strong empirical validation. I have a few main comments that, if clarified, could further strengthen the paper. Given satisfactory responses, I would be willing to raise my score.", "questions": "## CRPS notation\nEquation (2) uses the CRPS from Gneiting & Raftery (2007):  \n$$\n\\mathrm{CRPS}(Y,y)=\\tfrac12\\mathbb{E}|Y-Y'|-\\mathbb{E}|Y-y|.\n$$\nHowever, this corresponds to a proper scoring rules that is defined such that $S(Q,Q)\\ge S(P,Q)$, i.e., the true distribution maximizes the score. This conflicts with the minimization objectives in (3) and (4).  \nSince the closed-form in §3.5 is derived for minimization, as opposed to (2), the sign in (2) should be flipped for consistency. This is likely a notational issue and easily corrected.\n\n---\n\n## Notation in §3.1\nThe notation in §3.1 is somewhat inconsistent.\nYou define $\\mathcal{Y} \\subset \\mathbb{R}^{k}$ as the space of predicted distribution parameters,(l. 104), but later (l. 111) use $\\mathcal{Y}$ for the space of constrained distribution samples. \nMoreover, ${\\mathbf{Y}\\_{\\theta}} \\in \\mathbb{R}^{n}$ is introduced as a realization in $\\mathbb{R}^{n}$, which is a notation conflict with the parameter space $\\mathcal{Y}$. The way I understood it, $\\mathbf{Y}\\_{\\theta}$ is supposed to be a random variable. \n\nHowever, with the formulation $\\mathbf{Y}_{\\theta}(\\phi^{(i)})$ it appears like a functional, although it is probably meant as the conditional predictive distribution given $\\phi^{(i)}$. Clarity would improve by reserving distinct symbols: e.g., use $\\Theta$ for the parameter space, $\\Phi$ for network parameters, and define all variables upon first use. This would make it explicit that the network outputs distribution parameters, which then define a predictive measure from which realizations can be sampled.\n\n\n---\n\n## Framing of log-likelihood and CRPS\nLine 19 states that likelihood-based objectives are “heavily biased by their distributional assumptions.” This phrasing feels too strong since the proposed approach also assumes a distributional form (e.g., location-scale or Gaussian), for which the log-likelihood would be an equally valid objective. The CRPS is indeed advantageous, especially under misspecification, but I suggest softening the claim in the abstract and §3.1 or framing it in terms of the observed empirical performance benefits.  \nAdditionally, it could be interesting to evaluate *ProbHardE2E* with alternative proper scoring rules—such as the energy or kernel scores (Gneiting & Raftery 2007)—as you did for NLL in Q3. Showing that the approach generalizes to various scores would reinforce its flexibility.\n\n---\n\n## Discussion of unique solutions\nA brief discussion on the *existence and uniqueness* of the constrained optimization problem (eqs 1 and 7) would enhance completeness.  \nIf the feasible set is empty, the algorithm naturally terminates; but if the solution is unique, uncertainty vanishes and the probabilistic framework offers no benefit.  \nFor instance, with linear equalities where $q=n$ and $A$ is full rank, the projection $P_{Q^{-1}}=0$ yields the deterministic solution  \n$$\n\\hat{\\mu}=A^{\\dagger}b,\\quad \\Sigma=0,\n$$\nwhich is a degenerate distribution and technically not part of the underlying distributional family anymore (maybe depending a bit on the chosen definition). Certainly, in this case uncertainty quantification is not possible and the approach offers no advantage against non-probabilistic solvers. While this is not a limitation of the method, clarifying that *ProbHardE2E* mainly targets under-determined $q<n$ or non-unique settings would improve conceptual precision.\nMore generally, for any constraint, the transformation $\\mathcal{T}$ projects the distribution onto the subspace induced by the constraints and $\\Sigma$ does in general not have full rank, i.e. the underlying distribution is degenerate (by design of course). Maybe it would be worth to mention and/or analyze some implication of this, for example that uncertainty cannot be (meaningfully) analyzed across the whole space of the initially learned distribution, but only across the subspace of the transformed distribution.\n\n\n---\n\n## Minor comments\n- $\\theta$ appears in l. 106 before being defined.  \n- The domains of $\\phi^{(i)}$ and $u^{(i)}$ (l. 106) are missing.  \n- The spaces of $\\xi$ and $z_\\theta$ (l. 168–169) should be specified.\n\n---\n\n## Reference\nGneiting, T., & Raftery, A. E. (2007). *Strictly Proper Scoring Rules, Prediction, and Estimation.* **JASA**, 102 (477), 359–378.  \n<https://doi.org/10.1198/016214506000001437>", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "The paper provides an interesting combination of probabilistic learning and optimization under hard constraints. To my best judgement, the author have a good grasp of the theory and provide detailed proofs for their theoretical claims. I particularly enjoyed the idea of analyzing the push-forward distribution of the constrained optimization problem and showing ways to derive the anyltical solution for a location-scale family.\nThe numerical experiments support the author's claims, showing improved performance while fulfilling the hard constraints.", "weaknesses": "I believe, the presentation and readability could be improved, and some details of the method might need further clarification, see questions.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Poster)