End-to-End Probabilistic Framework for Learning with Hard Constraints

## Abstract

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 are heavily 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.

## 1 Introduction

Recently, machine learning (ML) models have been applied to a variety of engineering and scientific tasks, including probabilistic time series forecasting (Rangapuram et al., 2021; Hyndman et al., 2011; Taieb et al., 2017; Olivares et al., 2024b) and scientific applications (Krishnapriyan et al., 2021; Hansen et al., 2023; Négiar et 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 ML 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 are 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 a 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 and Kolter, 2017; Agrawal et al., 2019; Min et al., 2024), or 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 probabilistic forecasts.
Generating output distribution statistics under hard constraints is often computationally expensive or yields only approximate solutions (Robert et al., 1999; Szechtman and Glynn, 2001; Girolami and Calderhead, 2011).
There have been domain-specific 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 a 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 a post-processing 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 a 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, ProbHardE2E, that integrates a broad class of hard constraints (including non-linear constraints) in an end-to-end fashion, while incorporating UQ.
By leveraging key results from statistics and optimization in a novel way, we predict both the mean and covariance of the output data, moving beyond point estimate predictions.
ProbHardE2E enforces nonlinear constraints with an efficient sampling-free method to generate distribution statistics. 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 ProbHardE2E can incorporate a broad class of deep learning backbone models.
- •
The key technical novelty of ProbHardE2E is a *differentiable probabilistic projection layer* (DPPL) that extends standard projection methods to accommodate UQ while enforcing hard constraints. ProbHardE2E can handle constraints ranging from linear equality to general nonlinear equality to convex inequality constraints.
- •
We use the DPPL to impose constraints directly on the marginals of the multivariate distribution for an efficient sampling-free approach for posterior distribution estimation, which reduces the computational overhead by up to $3\text{--}5\times$ during training.
- •
We show that ProbHardE2E is effective in two (seemingly-unrelated, but technically-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 $15\times$ lower mean-squared error (MSE) in mean forecast and $2.5\times$ improved uncertainty estimates, measured by the Continuous Ranked Probability Score (CRPS), compared to the baseline methods.
- •
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).
While the need for this is well-known in, e.g., time series forecasting, previous PDE learning works commonly use NLL-based metrics for UQ.

## 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](https://arxiv.org/html/2506.07003v2#A1.T5)
in [Appendix A](https://arxiv.org/html/2506.07003v2#A1)
summarizes some advantages and disadvantages of these methods that are motivated by enforcing hard constraints in these domains.
(See [Appendix A](https://arxiv.org/html/2506.07003v2#A1) for additional details.)

## 3 ProbHardE2E : A Unified Probabilistic Optimization Framework

In this section, we introduce ProbHardE2E.
See [Algorithm 1](https://arxiv.org/html/2506.07003v2#alg1) for a summary.
(See also [Appendix B](https://arxiv.org/html/2506.07003v2#A2) for a universal approximation guarantee.)
In [Section 3.1](https://arxiv.org/html/2506.07003v2#S3.SS1), we discuss the proper evaluation metric for a constrained probabilistic learner, and we define our objective function that corresponds to that evaluation metric.
In [Section 3.2](https://arxiv.org/html/2506.07003v2#S3.SS2), we propose our differentiable probabilistic projection layer (DPPL) that enforces the hard constraints.
In [Section 3.3](https://arxiv.org/html/2506.07003v2#S3.SS3), we describe how to compute the parameters of the resulting constrained posterior distribution.
In [Section 3.4](https://arxiv.org/html/2506.07003v2#S3.SS4), we discuss update rules for various types of constraints (linear equality, nonlinear equality, and convex inequality constraints).
In [Section 3.5](https://arxiv.org/html/2506.07003v2#S3.SS5), we propose a sample-free formulation for satisfying the constraints while optimizing for the objective.

Figure: Algorithm 1 ProbHardE2E: Training and Inference

### 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
$\hat{f}_{\theta}:\Phi\rightarrow\mathcal{Y},$
where $\Phi$ $\subset\mathbb{R}^{m}$ denotes the input space, and $\mathcal{Y}$ $\subset\mathbb{R}^{k}$
denotes the space of predicted distribution parameters that meet the constraints.
Given a multivariate distribution class, these learned parameters induce a predictive multivariate $\mathbf{Y}_{\theta}(\phi^{(i)})\in\mathbb{R}^{n}$, where $({\phi}^{(i)},u^{(i)})\sim\mathcal{D}$ denotes training data from a distribution $\mathcal{D}$.
Each realization of $\hat{u}(\phi^{(i)})\sim\mathbf{Y}_{\theta}(\phi^{(i)})$ is required to satisfy predefined hard constraints of the form $g(\hat{u}(\phi^{(i)}))\leq 0$ and $h(\hat{u}(\phi^{(i)}))=0$.
We can formulate this constrained optimization problem as follows:

$$ $\displaystyle\operatorname*{argmin}_{\begin{subarray}{c}\theta\in{\Theta},\hskip 2.84544ptg(\mathbf{Y}_{\theta}(\phi^{(i)}))\leq 0,\hskip 2.84544pth(\mathbf{Y}_{\theta}(\phi^{(i)}))=0\end{subarray}}\mathbb{E}_{({\phi}^{(i)},u^{(i)})\sim\mathcal{D}}\ \ell\big(\mathbf{Y}_{\theta}({\phi}^{(i)}),u^{(i)}\big),$ (1) $$

where ${\Theta}$
denotes the parameter space, and
$\ell$ denotes a proper scoring rule.

One widely-used (strictly) proper scoring rule for continuous distributions is the continuous ranked probability score (CRPS) (Gneiting and 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:

$$ $\text{CRPS}(Y,y)=\tfrac{1}{2}\mathbb{E}_{Y}|Y-Y^{\prime}|-\mathbb{E}_{Y}|Y-y|,$ (2) $$

where $Y^{\prime}$ denotes an i.i.d. copy of $Y$.
Compared to other scoring rules, e.g., the log probability scoring rules, which require strong assumption on the outcome variable, the CRPS is robust to probabilistic model mis-specification.
Because of this unique property, the CRPS is widely used as the evaluation metric in many applications, e.g., probabilistic time series forecasting (Gasthaus et al., 2019; Rangapuram et 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 and 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 in [Equation 2](https://arxiv.org/html/2506.07003v2#S3.E2) in Problem [1](https://arxiv.org/html/2506.07003v2#S3.E1). We define the loss as the sum of the univariate CRPS:

$$ $\ell\big(\mathbf{Y}_{\theta}({\phi}^{(i)}),u^{(i)}\big)=\sum_{j=1}^{n}\text{CRPS}((\mathbf{Y}_{\theta}({\phi}^{(i)}))_{j},u_{j}^{(i)}).$ (3) $$

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 [1](https://arxiv.org/html/2506.07003v2#S3.E1) into the unconstrained optimization problem:

$$ $\underset{\theta\in\bar{\Theta}}{\operatorname*{argmin}}\ \mathbb{E}_{({\phi}^{(i)},u^{(i)})\sim\mathcal{D}}\ \ell\big(\mathbf{Y}_{\theta}({\phi}^{(i)}),u^{(i)}\big),$ (4) $$

where $\bar{\Theta}\subseteq\Theta$ denotes the feasible parameter space that ensures constraint satisfaction,
and $\ell$ denotes the loss function in [Equation 3](https://arxiv.org/html/2506.07003v2#S3.E3). We solve this using a two-step procedure: first define a predictive output distribution, then project it onto the constraint manifold using a *differentiable probabilistic projection layer (DPPL)* for end-to-end optimization.

Our framework begins with an established probabilistic backbone model. This can be a Gaussian Process (Rasmussen and Williams, 2006), neural process (Kim et al., 2019), DeepVAR (Salinas et al., 2019; Rangapuram et al., 2021), or ensembles of neural networks or operators (Mouli et al., 2024). The base model $f_{\theta}:\Phi\rightarrow\mathbb{R}^{k}$ predicts the distribution parameters (mean $\mu_{\theta}(\phi^{(i)})$ and covariance $\Sigma_{\theta}(\phi^{(i)})$, for $\theta\in\Theta$) – without constraint awareness. We then use a reparameterization function $r:\mathbb{R}^{k}\times\mathbb{R}^{n}\rightarrow\mathbb{R}^{l}$ to define the distribution in one of two ways: either as an identity map, where $l=k$, that returns $f_{\theta}(\phi^{(i)})=\big(\mu_{\theta}(\phi^{(i)}),\Sigma_{\theta}(\phi^{(i)})\big)$ for our efficient sample-free paradigm during training; or as a map, where $l=n$, that combines the distribution parameters with noise $\xi\sim p(\xi)$, where $p$ denotes a tractable sampling distribution, and gives a sample $z_{\theta}(\phi^{(i)})\sim\mathbf{Z}_{\theta}(\phi^{(i)})$ 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:

$$ $r(f_{\theta}(\phi^{(i)}),\xi)=\begin{cases}\big(\mu_{\theta}(\phi^{(i)}),\Sigma_{\theta}(\phi^{(i)})\big),\hfill\text{ (Training)}\\ \mathbf{Z}_{\theta}(\phi^{(i)}),\hfill\text{ (Inference)}\end{cases}$ (5) $$

Following this Predictor Step above, we use the DPPL in the Corrector Step to restrict the parameter space to $\bar{\Theta}\subseteq\Theta$, such that for all $\hat{u}_{\theta}(\phi^{(i)})\sim\mathbf{Y}_{\theta}(\phi^{(i)})$, the constraints $g(\hat{u}_{\theta}(\phi^{(i)}))\leq 0$ and $h(\hat{u}_{\theta}(\phi^{(i)}))=0$ are 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 parameters or projected predictive random variable as:

$$ $\text{DPPL}(r(f_{\theta}(\phi^{(i)}),\xi),g(\cdot),h(\cdot))=r(\hat{f}_{\theta}(\phi^{(i)}),\xi)=\begin{cases}\big(\hat{\mu}_{\theta}(\phi^{(i)}),\hat{\Sigma}_{\theta}(\phi^{(i)})\big),\hfill\text{ (Training)}\\ \mathbf{Y}_{\theta}(\phi^{(i)}),\hfill\text{ (Inference)}\end{cases}$ (6) $$

for $r(f_{\theta}(\phi^{(i)}),\xi)$ in [Equation 5](https://arxiv.org/html/2506.07003v2#S3.E5),
where $\hat{f}_{\theta}:\Phi\rightarrow\mathcal{Y}\subset\mathbb{R}^{k}$ denotes the probabilistic model that outputs the constrained distribution parameters $(\hat{\mu}_{\theta}(\phi^{(i)}),\hat{\Sigma}_{\theta}(\phi^{(i)}))$.
Our DPPL yields a constraint-satisfying realization $u^{*}\sim\mathbf{Y}_{\theta}(\phi^{(i)})$
as the final predictive random variable.

This two-step approach mirrors predictor-corrector methods (Boyd and 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 $\mathbf{Z_{\theta}}(\phi^{(i)})$.
(See [Appendix C](https://arxiv.org/html/2506.07003v2#A3) for details.)
Prior works on imposing hard constraints in time series and solving PDEs (Rangapuram et al., 2021; Hansen et al., 2023) reduce to special cases of our method with linear constraints.
(See [Appendix D](https://arxiv.org/html/2506.07003v2#A4) for details.)
We draw $z_{\theta}(\phi^{(i)})\sim\mathbf{Z}_{\theta}(\phi^{(i)})$, and we solve the following constrained optimization problem:

$$ $\displaystyle u^{*}(z_{\theta}(\phi^{(i)})):=\operatorname*{argmin}_{\begin{subarray}{c}\hat{u}_{\theta}(\phi^{(i)})\in\mathbb{R}^{n},g(\hat{u}_{\theta}(\phi^{(i)}))\leq 0,h(\hat{u}_{\theta}(\phi^{(i)}))=0\end{subarray}}\lVert\hat{u}_{\theta}(\phi^{(i)})-z_{\theta}(\phi^{(i)})\rVert^{2}_{Q},$ (7) $$

where $u^{*}(z_{\theta}(\phi^{(i)}))$ denotes a predicted sample of $\mathbf{Y}_{\theta}(\phi^{(i)})$, and where $\lVert x\rVert_{Q}=\sqrt{x^{\top}Qx}$ for some symmetric positive semi-definite matrix $Q$.
(See [Appendix E](https://arxiv.org/html/2506.07003v2#A5) for details on the flexibility of learning various forms of $Q$.)

### 3.3 DPPL on the Distribution Parameters for Location-Scale 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 $\mathcal{F}$ belongs to a multivariate, location-scale family, i.e.,
a distribution such that
any affine transformation $\mathbf{Y}$ of a random variable $\mathbf{Z}=\mu+\Sigma^{1/2}\xi\sim\mathcal{F}(\mu,\Sigma)$ and $\xi\sim\mathcal{F}(0,1)$,
remains within the same distribution family $\mathcal{F}$. This is an example of how to compute the random variable in [Equation 5](https://arxiv.org/html/2506.07003v2#S3.E5) for a multivariate location-scale distribution.
A familiar case of this is when $\mathbf{Z}\sim\mathcal{N}(\mu,\Sigma)$ and $\mathbf{Y}=A\mathbf{Z}+B$ is an affine transformation; in which case $\mathbf{Y}\sim\mathcal{N}(A\mu+B,A\Sigma A^{\top})$.
Alternatively,
we can show that when $\mathbf{Y}$ is a nonlinear transformation of $\mathbf{Z}$, it has approximately (to first-order) the same distribution $\mathbf{Z}$, with an appropriately-chosen set of parameters (given in [Equation 8](https://arxiv.org/html/2506.07003v2#S3.E8) below).
We state this result more formally in [Theorem 3.1](https://arxiv.org/html/2506.07003v2#S3.Thmtheorem1).
The proof, given in [Appendix F](https://arxiv.org/html/2506.07003v2#A6), uses a first-order Taylor expansion to linearize the nonlinear function transformation, and is similar to the Multivariate Delta Method (Casella and Berger, 2001).

###### Theorem 3.1 .

Let $\mathbf{Z}\sim\mathcal{F}(\mu,\Sigma)$ be a random variable, where the underlying distribution $\mathcal{F}$ belongs to a multivariate location-scale family of distributions, with mean $\mu$ and covariance $\Sigma$; and let $\mathcal{T}$ be a function with continuous first derivatives, such that ${J_{\mathcal{T}}(\mu)}\Sigma J_{\mathcal{T}}(\mu)^{\top}$ is symmetric positive semi-definite.
Then, the transformed distribution $\mathbf{Y}=\mathcal{T}(\mathbf{Z})$ converges in distribution with first-order accuracy to $\mathcal{F}(\hat{\mu},\hat{\Sigma})$ with mean $\hat{\mu}=\mathcal{T}(\mu)$ and covariance $\hat{\Sigma}=J_{\mathcal{T}}(\mu)\Sigma J_{\mathcal{T}}(\mu)^{\top}$, where $J_{\mathcal{T}}(\mu)={\nabla\mathcal{T}(\mu)}^{\top}$ denotes the Jacobian of $\mathcal{T}$ with respect to $z$ evaluated at $\mu$.

Let $\mathbf{Z}\sim\mathcal{F}(\mu,\Sigma)$ denote the prior distribution and $z\sim\mathbf{Z}$.
We apply [Theorem 3.1](https://arxiv.org/html/2506.07003v2#S3.Thmtheorem1) with $\mathcal{T}(z)=u^{*}(z)$, where $u^{*}(z)$ denotes the solution of the constrained least squares problem in Problem [7](https://arxiv.org/html/2506.07003v2#S3.E7).
In this case, the projected random variable satisfies
$\mathbf{Y}\sim\mathcal{F}(\hat{\mu},\hat{\Sigma})$
with updated parameters:

$$ $\hat{\mu}=\mathcal{T}(\mu),\qquad\hat{\Sigma}=J_{\mathcal{T}}(\mu)\,\Sigma\,J_{\mathcal{T}}(\mu)^{\top}.$ (8) $$

###### Theorem 3.1 .

### 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 1](https://arxiv.org/html/2506.07003v2#S3.T1) shows these constraints types require different treatments: linear equality have closed-form projections, nonlinear equality can be solved with iterative methods, and convex inequality require optimization solvers.

**Table 1: Summary of DPPL in ProbHardE2E for various constraint types. For linear equality constraints, the oblique projection $P_{Q^{-1}}=I-Q^{-1}A^{\top}(AQ^{-1}A^{\top})^{-1}A$; for nonlinear equality constraints, $R$ denotes the first-order optimality conditions.**
| Constraint Type | Solution $u^{*}(z)$ | Solver Type | Jacobian $J_{\mathcal{T}}$ |
| --- | --- | --- | --- |
| Linear Equality | $P_{Q^{-1}}z+(I-P_{Q^{-1}})A^{{\dagger}}b$ | closed-form | $P_{Q^{-1}}$ |
| Nonlinear Equality | $(u^{*},\lambda^{*})$ s.t. $R(u^{*},\lambda^{*};z)=0$ | nonlinear | implicit differentiation |
| Convex Inequality | $\underset{h(\hat{u})=0,\,g(\hat{u})\leq 0}{\text{argmin}}\|\hat{u}-z\|_{Q}^{2}$ | convex opt. | sensitivity analysis;<br>argmin differentiation |

#### 3.4.1 Linear Equality Constraints

For linear equality constraints, $h(\hat{u})=A\hat{u}-b=0$, where $A\in\mathbb{R}^{q\times n},q\leq n$, has full row rank $q$, we can derive a closed-form solution to the constrained least squares Problem [7](https://arxiv.org/html/2506.07003v2#S3.E7). In this case, both training and inference modes are equivalent since the DPPL projection is exact. (See [Section C.1](https://arxiv.org/html/2506.07003v2#A3.SS1).)

#### 3.4.2 Nonlinear Equality Constraints

For nonlinear equality constraints, $h(\hat{u})=0$, we can no longer derive the exact closed-form expression for the solution.
Instead, we can provide an expression which is satisfied by the optimal solution. In particular, we approximate the parameter-level projection at training time.
This can then be solved for the posterior mean $\hat{\mu}=u^{*}(\mu)$ in [Equation 8](https://arxiv.org/html/2506.07003v2#S3.E8) with the nonlinear transformation $\mathcal{T}(\mu)=u^{*}(\mu)$ with iterative optimization methods, e.g., Newton’s Method.
(We can then compute the posterior covariance $\hat{\Sigma}$ in [Equation 8](https://arxiv.org/html/2506.07003v2#S3.E8) by estimating the Jacobian $J_{\mathcal{T}}(\mu)$ by differentiating the nonlinear equations $u^{*}(z)=z-Q^{-1}\nabla h(u^{*}(z))^{\top}\lambda\,,\ h(u^{*}(z))=0$ with respect to $z$ via the implicit function theorem (Blondel et al., 2022), and evaluating it at $\mu$.
(See [Section C.2](https://arxiv.org/html/2506.07003v2#A3.SS2).) At inference, we project each sample exactly with our custom, batched optimization solver to ensure strict constraint feasibility.

#### 3.4.3 (Nonlinear) Convex Inequality Constraints

For convex inequality constraints, $\hat{u}$ in Problem [7](https://arxiv.org/html/2506.07003v2#S3.E7) is in a convex set, $\mathcal{C}\subset\mathbb{R}^{n}$.
Closed-form expressions (such as those in previous subsections for linear and nonlinear equality constraints) do not exist (Boyd and Vandenberghe, 2004). Instead, we rely on convex optimization solvers to 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 and Shapiro, 2013), argmin differentiation (Sun et al., 2022; Agrawal et al., 2019; Amos and Kolter, 2017; Gould et al., 2016), and/or variational analysis (Rockafellar and Wets, 2009).
These techniques provide a means to compute the Jacobian $J_{\mathcal{T}}(\mu)$, which represents the sensitivity of the optimal solution $u^{*}$ to changes in the input vector $\mu$, whose projection we are essentially computing to the convex constraints space. During inference, we solve the convex program per sample. (See [Section C.3](https://arxiv.org/html/2506.07003v2#A3.SS3).)

### 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 $\ell$ in [Equation 3](https://arxiv.org/html/2506.07003v2#S3.E3). Calculating the CRPS for an arbitrary distribution requires generating samples (Rangapuram et al., 2021; Gneiting and Raftery, 2007), but closed-form expressions for the CRPS exist for several location-scale distributions (Gaussian, logistic, student’s t, beta, gamma, uniform).
Most notably, for the univariate Gaussian, the closed-form CRPS is given as:
$\text{CRPS}_{\mathcal{N}}(z)=\left[z\cdot\left(2\mathrm{P}(z)-1\right)+2p(z)-\frac{1}{\sqrt{\pi}}\right],$
where $p(z)=\frac{1}{\sqrt{2\pi}}\exp{(-z^{2}/2)}$ denotes the standard normal probability density function (PDF), and $\mathrm{P}(z)=\int_{-\infty}^{z}p(y)dy$ denotes the standard normal cumulative distribution function (CDF) for $z\sim\mathcal{N}(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 ProbHardE2E:

- (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 orthogonal projection, and if so when?
- (Q3)
Does training with the distribution-agnostic 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 on the samples?
- (Q5)
How does ProbHardE2E perform when extended to more general constraints, e.g., nonlinear equality and convex inequality constraints?

See [Appendix G](https://arxiv.org/html/2506.07003v2#A7)
for details on the test datasets, [Appendix H](https://arxiv.org/html/2506.07003v2#A8) for implementation details, and [Appendix I](https://arxiv.org/html/2506.07003v2#A9) 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 [Section G.1](https://arxiv.org/html/2506.07003v2#A7.SS1) for details.) In addition to PDEs, we also evaluate ProbHardE2E on 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 [Section G.2](https://arxiv.org/html/2506.07003v2#A7.SS2) for details.)

##### Baselines.

We compare two variants of ProbHardE2E, i.e., ProbHardE2E-Ob, which uses a general oblique projection ($Q=\Sigma^{-1}$) projection and is our default unless otherwise specified, and ProbHardE2E-Or, which uses an orthogonal projection ($Q=I$), against several probabilistic deep learning baselines commonly used for uncertainty quantification in constrained PDEs and probabilistic time series forecasting. For PDEs, ProbHardE2E 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 et 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 ProbHardInf, which is a variant of ProbConserv that imposes the nonlinear constraint at inference time only);
(iii) SoftC (Hansen et al., 2023), which introduces a soft penalty on constraint violation
à la PINNs (Raissi et al., 2019; Li et al., 2024) during training but does not guarantee constraint satisfaction at inference; and (iv)
the unconstrained model backbone VarianceNO. For hierarchical time-series forecasting, ProbHardE2E uses DeepVAR (Salinas et al., 2019) as the probabilistic base model.
We compare ProbHardE2E with:
(i) ProbConserv;
(ii) HierE2E (Rangapuram et al., 2021), which enforces linear constraints via an end-to-end orthogonal projection;
classical statistical approaches including ARIMA-NaiveBU, ETS-NaiveBU  (Hyndman et al., 2011, 2025), (iii) PERMBU-MINT (Taieb et al., 2017; Olivares et al., 2022); and (iv)
the unconstrained model backbone DeepVAR.

##### Evaluation.

We evaluate ProbHardE2E using the following metrics:
Mean Squared Error (MSE), which measures the mean prediction accuracy;
Constraint Error (CE), which measures the constraint errors on the samples (conservation law for PDEs and coherency for hierarchical time series forecasting); and
Continuous Ranked Probability Score (CRPS), which measures performance in uncertainty quantification (UQ). (See [Section H.3](https://arxiv.org/html/2506.07003v2#A8.SS3) 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 ProbConserv to ensure a fair comparison. The experiments are conducted on a single NVIDIA V100 GPU in the PDEs case, and on an Intel(R) Xeon(R) CPU E5-2603 v4 @ 1.70GHz in the time series forecasting case. To ensure scalability, we use a diagonal covariance matrix $Q$ in Problem [7](https://arxiv.org/html/2506.07003v2#S3.E7), following prior work (Hansen et al., 2023; Mouli et al., 2024). (See [Appendix E](https://arxiv.org/html/2506.07003v2#A5) for low-rank and full covariance structures.)

### 4.1 Linear Conservation and Hierarchical Constraints

In this subsection, we test ProbHardE2E on linear constraints. Table [2](https://arxiv.org/html/2506.07003v2#S4.T2) presents our comparative evaluation results across multiple PDE datasets under linear constraints, and Table [3](https://arxiv.org/html/2506.07003v2#S4.T3) presents our evaluation results across multiple time series forecasting datasets. We use these results to answer questions [(Q1)](https://arxiv.org/html/2506.07003v2#S4.I1.i1)-[(Q4)](https://arxiv.org/html/2506.07003v2#S4.I1.i4) raised above.

**Table 2: Test metrics on constrained PDEs across four datasets, which are ordered top to bottom in their learning difficulty. Metrics include MSE $\times 10^{-5}$, constraint (conservation) error (CE) $\times 10^{-3}$, and CRPS $\times 10^{-3}$. Each algorithm is trained with either CRPS or NLL. Best values per row are highlighted in bold.**
| Dataset | Metric | ProbHardE2E-Ob | ProbHardE2E-Or | HardC | ProbConserv | SoftC | VarianceNO (base) |  |  |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
|  |  | CRPS | NLL | CRPS | NLL | CRPS | NLL | CRPS | NLL | CRPS | NLL | CRPS | NLL |
| Heat | MSE | 0.036 | 0.047 | 0.031 | 0.301 | 0.031 | 0.090 | 0.027 | 1.26 | 0.051 | 0.156 | 0.029 | 2.01 |
| CE | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0.852 | 4.806 | 1.76 | 34.3 |  |
| CRPS | 0.304 | 0.37 | 0.271 | 0.713 | 0.275 | 0.452 | 0.392 | 4.27 | 0.354 | 1.129 | 0.396 | 4.39 |  |
| PME | MSE | 9.59 | 6.16 | 9.01 | 11.08 | 8.870 | 10.55 | 8.801 | 10.5 | 8.187 | 7.362 | 7.945 | 8.132 |
| CE | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 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.065 | 2.444 | 2.02 | 2.43 |  |
| Advection | MSE | 131 | 262 | 88.09 | 310.82 | 103.78 | 458.38 | 134 | 277 | 148.11 | 599.11 | 149 | 605 |
| CE | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 19.334 | 182.99 | 18.9 | 182 |  |
| CRPS | 4.19 | 7.03 | 2.94 | 8.669 | 3.236 | 11.23 | 3.88 | 7.90 | 3.963 | 9.702 | 3.98 | 10.1 |  |
| Stefan | MSE | 186 | 207 | 394.84 | 432.92 | 394.29 | 433.28 | 303 | 273 | 431.89 | 429.06 | 425 | 425 |
| CE | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 166.93 | 168.75 | 180 | 169 |  |
| CRPS | 7.52 | 7.85 | 14.147 | 14.67 | 14.432 | 14.539 | 7.85 | 8.33 | 9.878 | 10.062 | 9.51 | 10.2 |  |

##### Q1.

Regarding [(Q1)](https://arxiv.org/html/2506.07003v2#S4.I1.i1) on the benefits of training end-to-end with a hard constraint, the results demonstrate that our method achieves superior performance compared to existing approaches.
Specifically, when measured against two accuracy metrics across four PDE datasets in Table [2](https://arxiv.org/html/2506.07003v2#S4.T2), our method
with either oblique (ProbHardE2E-Ob) or orthogonal (ProbHardE2E-Or) projection
consistently outperforms both ProbConserv, which applies constraints only during post-processing, and SoftC, which implements constraints as soft penalties. This performance advantage directly addresses research question [(Q1)](https://arxiv.org/html/2506.07003v2#S4.I1.i1), 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 state-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 HierE2E.

##### Q2.

Regarding [(Q2)](https://arxiv.org/html/2506.07003v2#S4.I1.i2) on the effectiveness of the oblique vs. orthogonal projections, [Table 2](https://arxiv.org/html/2506.07003v2#S4.T2) shows while both oblique (ProbHardE2E-Ob) and orthogonal (ProbHardE2E-Or) variants of ProbHardE2E enforce zero constraint error, their impact on predictive fidelity varies significantly, 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 to the orthogonal projection. [Table 3](https://arxiv.org/html/2506.07003v2#S4.T3) shows that ProbHardE2E-Or generally performs better on the time series datasets with cross-sectional hierarchies, as it improves CRPS on Labour and Tourism datasets, compared to ProbHardE2E-Ob. This addresses [(Q2)](https://arxiv.org/html/2506.07003v2#S4.I1.i2) 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.

**Table 3: CRPS $\times 10^{-3}$ for hierarchical time‑series datasets across various probabilistic forecasting algorithms. Constraint (coherency) error (CE) is given in parenthesis and is equal to 0 for all methods except the base unconstrained DeepVAR. PERMBU-MINT is not available for Tourism-L, because the dataset has a nested hierarchical structure, and PERMBU-MINT is not well-defined on this type of dataset (Rangapuram et al., 2021). ARIMA-NaiveBU, ETS-NaiveBU and PERMBU-MINT are deterministic models with no model uncertainty.**
| Dataset | ProbHardE2E-Ob | ProbHardE2E-Or | ProbConserv | HierE2E | ARIMA-NaiveBU | ETS-NaiveBU | PERMBU-MINT | DeepVAR (base) |
| --- | --- | --- | --- | --- | --- | --- | --- | --- |
| Labour | 36.1$\pm$ 2.7 (0) | 28.6$\pm$6.5 (0) | 45.8$\pm$6.5 (0) | 50.5$\pm$20.6 (0) | 45.3 (0) | 43.2 (0) | 39.3 (0) | 38.2$\pm$4.5 (0.215) |
| Tourism | 98.9$\pm$13.0 (0) | 82.4$\pm$6.6 (0) | 100.7$\pm$7.7 (0) | 103.1$\pm$16.3 (0) | 113.8 (0) | 100.8 (0) | 77.1 (0) | 92.5$\pm$2.2 (2818.01) |
| Tourism-L | 155.2$\pm$3.6 (0) | 156.4$\pm$9.4 (0) | 176.9$\pm$21.5 (0) | 161.3$\pm$10.9 (0) | 174.1 (0) | 169.0 (0) | – | 158.1$\pm$10.2 ($70000$) |
| Traffic | 55.0$\pm$10.6 (0) | 60.6$\pm$7.8 (0) | 71.0$\pm$3.9 (0) | 41.8$\pm$7.8 (0) | 80.8 (0) | 66.5 (0) | 67.7 (0) | 40.0$\pm$2.6 (0.192) |
| Wiki | 212.1$\pm$29.4 (0) | 215.8$\pm$16.9 (0) | 264.7$\pm$ 30.7 (0) | 216.5$\pm$26.7 (0) | 377.2 (0) | 467.3 (0) | 281.2 (0) | 229.4$\pm$15.8 (8398.6) |

##### Q3.

Regarding [(Q3)](https://arxiv.org/html/2506.07003v2#S4.I1.i3) on the training objective, [Table 2](https://arxiv.org/html/2506.07003v2#S4.T2) 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](https://arxiv.org/html/2506.07003v2#S4.T3) shows that on four out of five time series datasets, we improve upon the results of Hier-E2E, which uses the DeepVAR base model with an orthogonal projection on the samples, and which optimizes the sampling-based quantile loss by projecting directly on the distribution parameters.

##### Q4.

Regarding [(Q4)](https://arxiv.org/html/2506.07003v2#S4.I1.i4) on the computational efficiency of our sampling-free approach, LABEL:fig:pde_timing shows the training time per epoch for the hierarchical time series datasets. Models trained for time series and PDEs (see [Section I.1](https://arxiv.org/html/2506.07003v2#A9.SS1)) with 100 posterior samples per training step incur a $3.3\text{--}4.6\times$ increase in epoch time relative to ProbHardE2E, which avoids sampling altogether by using a closed-form CRPS loss. Note that the computational overhead of ProbHardE2E is approximately $2\times$ that of the unconstrained model. However, compared to the sampling-based probabilistic baselines, our 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, we test ProbHardE2E on nonlinear equality and convex inequality constraints to address question [(Q5)](https://arxiv.org/html/2506.07003v2#S4.I1.i5) on PDE datasets, as (current) time series forecasting applications usually need predictions to satisfy only linear (e.g., hierarchical) constraints.

#### 4.2.1 Nonlinear Equality Constraints

We test ProbHardE2E with general nonlinear constraints using nonlinear conservation laws from PDEs. (See [Section I.2](https://arxiv.org/html/2506.07003v2#A9.SS2) for details.) Importantly, [Table 4](https://arxiv.org/html/2506.07003v2#S4.T4) shows that even in this challenging case of nonlinear constraints, the constraint error (CE) on the samples is 0 so that we ensure strict feasibility on the samples. In addition, we see superior performance of enforcing nonlinear constraints with ProbHardE2E. We see an even larger MSE performance improvement of $\approx 15-17\times$ when trained on CRPS, and CRPS performance improvement of $\approx 2.5\times$ over the various baselines that apply the nonlinear constraint just at inference time (ProbHardInf), as a reduced linear constraint (ProbConserv) 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 approximating satisfying the constraint during training, and exact constraint sanctification at inference time while achieving state-of-the-art performance measured in CRPS.
In addition, LABEL:fig:nonlin_pme shows that ProbHardE2E is able to significantly better capture the shock and has tighter uncertainty estimates, leading to lower CRPS values than the baselines.

Figure: (a) ProbHardE2E on the various constraint types. (a) Linear Equality: Average time per iteration (in seconds) for ProbHardE2E, compared to the HierE2E on five hierarchical time-series datasets; (b) Nonlinear Equality: Mean $\pm$3 standard deviation for the PME with conservation constraint at time $t=0.51$, with PDE parameter $m_{\text{train}}\in[3,4]$ and $m_{\text{test}}=3.88$; (c) Convex Inequality: Mean $\pm$3 standard deviation for linear advection with TVD constraint at time $t=0.51$, with PDE parameter $\beta_{\text{train}}\in[1,2]$ and $\beta_{\text{test}}=1.5$. The horizontal axes in (b)-(c) are zoomed in to highlight the uncertainty near the propagating front.
Refer to caption: x1.png

#### 4.2.2 (Nonlinear) Convex Inequality Constraints

We impose a convex relaxation of the total variation diminishing (TVD) constraint that has been commonly used in numerical methods for PDEs to minimize artificial oscillations (LeVeque, 2002). In particular, we impose $\text{TVD}=\sum_{n=1}^{N_{t}}\sum_{i=1}^{N_{x}}\left|u(t_{n},x_{i+1})-u(t_{n},x_{i})\right|$ as a regularization term. (See [Section I.3](https://arxiv.org/html/2506.07003v2#A9.SS3) for details.) Note that this discrete form is analogous to total variation denoising in signal processing (Rudin et al., 1992; Boyd and Vandenberghe, 2004).
LABEL:fig:lin_adv_tvd illustrates that imposing this TVD relaxation improves the shock location prediction, compared to the unconstrained model VarianceNO. We see that ProbHardE2E has smaller variance, compared to both ProbConserv, which only enforces the conservation law, and VarianceNO. Most importantly, we see that ProbHardE2E is less likely to predict non-physical negative samples, which violates the positivity of the true solution. In addition, the variance of the ProbHardE2E solution also has a 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 to challenging nonlinear constraints into black-box probabilistic deep learning models using a differentiable probabilistic projection layer (DPPL). We show that ProbHardE2E is 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 et al., 2023), which only support linear global conservation constraints, ProbHardE2E can support 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 (LeVeque, 2002).
Future work also includes extending our 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.

##### Acknowledgements.

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.

## References

- Agrawal et al. (2019)
A. Agrawal, B. Amos, S. Barratt, S. Boyd, S. Diamond, and J. Z. Kolter.
Differentiable convex optimization layers.
In *Advances in neural information processing systems*, volume 32, 2019.
- Alexandrov et al. (2019)
A. Alexandrov, K. Benidis, M. Bohlke-Schneider, V. Flunkert, J. Gasthaus, T. Januschowski, D. C. Maddix, S. Rangapuram, D. Salinas, J. Schulz, et al.
Gluonts: Probabilistic and neural time series modeling in python.
*Journal of Machine Learning Research*, 21(116):1–6, 2019.
- Amos and Kolter (2017)
B. Amos and J. Z. Kolter.
Optnet: Differentiable optimization as a layer in neural networks.
In *Proceedings of the 34th International Conference on Machine Learning*, volume 70, pages 136–145. PMLR, 2017.
- Anava et al. (2018)
O. Anava, V. Kuznetsov, and (Google Inc. Sponsorship).
Web traffic time series forecasting, forecast future traffic to wikipedia pages.
Kaggle Competition, 2018.
URL [https://www.kaggle.com/c/web-traffic-time-series-forecasting/](https://www.kaggle.com/c/web-traffic-time-series-forecasting/).
- Ansari et al. (2024)
A. F. Ansari, L. Stella, C. Turkmen, X. Zhang, P. Mercado, H. Shen, O. Shchur, S. S. Rangapuram, S. P. Arango, S. Kapoor, et al.
Chronos: Learning the language of time series.
*Transactions on Machine Learning Research*, 2024.
- Athanasopoulos et al. (2024)
G. Athanasopoulos, R. J. Hyndman, N. Kourentzes, and A. Panagiotelis.
Forecast reconciliation: A review.
*International Journal of Forecasting*, 40(2):430–456, 2024.
- Australian Bureau of Statistics (2019)
Australian Bureau of Statistics.
Labour Force, Australia.
Accessed Online, 2019.
URL [https://www.abs.gov.au/AUSSTATS/abs@.nsf/DetailsPage/6202.0Dec%202019?OpenDocument](https://www.abs.gov.au/AUSSTATS/abs@.nsf/DetailsPage/6202.0Dec%202019?OpenDocument).
- Battaglia et al. (2018)
P. W. Battaglia, J. B. Hamrick, V. Bapst, A. Sanchez-Gonzalez, V. Zambaldi, M. Malinowski, A. Tacchetti, D. Raposo, A. Santoro, R. Faulkner, et al.
Relational inductive biases, deep learning, and graph networks.
*arXiv preprint arXiv:1806.01261*, 2018.
- Ben Taieb and Koo (2019)
S. Ben Taieb and B. Koo.
Regularized regression for hierarchical forecasting without unbiasedness conditions.
In *Proceedings of the 25th ACM SIGKDD international conference on knowledge discovery & data mining*, pages 1337–1347, 2019.
- Benidis et al. (2022)
K. Benidis, S. S. Rangapuram, V. Flunkert, Y. Wang, D. Maddix, C. Turkmen, J. Gasthaus, M. Bohlke-Schneider, D. Salinas, L. Stella, F.-X. Aubet, L. Callot, and T. Januschowski.
Deep learning for time series forecasting: Tutorial and literature survey.
*ACM Comput. Surv.*, 55(6), 2022.
- Bertsekas (1997)
D. P. Bertsekas.
Nonlinear programming.
*Journal of the Operational Research Society*, 48(3):334–334, 1997.
- Beucler et al. (2021)
T. Beucler, M. Pritchard, S. Rasp, J. Ott, P. Baldi, and P. Gentine.
Enforcing analytic constraints in neural networks emulating physical systems.
*Physical Review Letters*, 126(9):098302, 2021.
- Biegler et al. (2003)
L. T. Biegler, O. Ghattas, M. Heinkenschloss, and B. van Bloemen Waanders.
Large-scale pde-constrained optimization: an introduction.
In *Large-scale PDE-constrained optimization*, pages 3–13. Springer, 2003.
- Björck (1996)
d. Björck.
*Numerical Methods for Least Squares Problems*.
Society for Industrial and Applied Mathematics, 1996.
- Blondel and Roulet (2024)
M. Blondel and V. Roulet.
The elements of differentiable programming.
*arXiv preprint arXiv:2403.14606*, 2024.
- Blondel et al. (2022)
M. Blondel, Q. Berthet, M. Cuturi, R. Frostig, S. Hoyer, F. Llinares-López, F. Pedregosa, and J.-P. Vert.
Efficient and modular implicit differentiation.
In *Advances in Neural Information Processing Systems*, volume 35, pages 5230–5242, 2022.
- Bolton and Zanna (2019)
T. Bolton and L. Zanna.
Applications of Deep Learning to Ocean Data Inference and Subgrid Parameterization.
*Journal of Advances in Modeling Earth Systems*, 11(1):376–399, 2019.
- Bonnans and Shapiro (2013)
J. F. Bonnans and A. Shapiro.
*Perturbation analysis of optimization problems*.
Springer Science & Business Media, 2013.
- Boumal (2024)
N. Boumal.
*An Introduction to Optimization on Smooth Manifolds*.
Cambridge University Press, 2024.
- Box et al. (2015)
G. E. Box, G. M. Jenkins, G. C. Reinsel, and G. M. Ljung.
*Time series analysis: forecasting and control*.
John Wiley & Sons, 2015.
- Boyd and Vandenberghe (2004)
S. P. Boyd and L. Vandenberghe.
*Convex optimization*.
Cambridge university press, 2004.
- Casella and Berger (2001)
G. Casella and R. Berger.
*Statistical Inference*.
Duxbury Resource Center, June 2001.
- Chalapathi et al. (2024)
N. Chalapathi, Y. Du, and A. Krishnapriyan.
Scaling physics-informed hard constraints with mixture-of-experts.
In *International Conference on Learning Representations*, 2024.
- Chen et al. (2024)
H. Chen, G. E. C. Flores, and C. Li.
Physics-informed neural networks with hard linear equality constraints.
*Computers and Chemical Engineering*, 189:108764, 2024.
- Cheng et al. (2025)
C. Cheng, B. Han, D. C. Maddix, A. F. Ansari, A. Stuart, M. W. Mahoney, and Y. Wang.
Gradient-free generation for hard-constrained systems.
In *International Conference on Learning Representations*, 2025.
- Das et al. (2023)
A. Das, W. Kong, R. Sen, and Y. Zhou.
A decoder-only foundation model for time-series forecasting.
*arXiv preprint arXiv:2310.10688*, 2023.
- Donti et al. (2021)
P. L. Donti, D. Rolnick, and J. Z. Kolter.
DC3: A learning method for optimization with hard constraints.
In *International Conference on Learning Representations*, 2021.
- Edwards (2022)
C. Edwards.
Neural networks learn to speed up simulations.
*Communications of the ACM*, 65(5):27–29, 2022.
- Eisenach et al. (2022)
C. Eisenach, Y. Patel, and D. Madeka.
Mqtransformer: Multi-horizon forecasts with context dependent and feedback-aware attention.
*arXiv preprint arXiv:2009.14799*, 2022.
- Fakoor et al. (2023)
R. Fakoor, T. Kim, J. Mueller, A. J. Smola, and R. J. Tibshirani.
Flexible model aggregation for quantile regression.
*Journal of Machine Learning Research*, 24(162):1–45, 2023.
- Fang et al. (2024)
Y. Fang, S. Na, M. W. Mahoney, and M. Kolar.
Fully stochastic trust-region sequential quadratic programming for equality-constrained optimization problems.
*SIAM Journal on Optimization*, 34(2):2007–2037, 2024.
- Fortunato et al. (2022)
M. Fortunato, T. Pfaff, P. Wirnsberger, A. Pritzel, and P. Battaglia.
Multiscale meshgraphnets.
*2nd AI4Science Workshop at the 39th International Conference on Machine Learning (ICML*, 2022.
- Gao et al. (2023)
Z. Gao, X. Shi, B. Han, H. Wang, X. Jin, D. C. Maddix, Y. Zhu, M. Li, and Y. Wang.
PreDiff: Precipitation nowcasting with latent diffusion models.
In *Advances in Neural Information Processing Systems*, 2023.
- Gasthaus et al. (2019)
J. Gasthaus, K. Benidis, Y. Wang, S. S. Rangapuram, D. Salinas, V. Flunkert, and T. Januschowski.
Probabilistic forecasting with spline quantile function rnns.
In *Proceedings of the Twenty-Second International Conference on Artificial Intelligence and Statistics*, volume 89, pages 1901–1910. PMLR, 2019.
- Gill and Wong (2012)
P. E. Gill and E. Wong.
Sequential quadratic programming methods.
In *Mixed Integer Nonlinear Programming*, pages 147–224, New York, NY, 2012. Springer New York.
- Girolami and Calderhead (2011)
M. Girolami and B. Calderhead.
Riemann manifold langevin and hamiltonian monte carlo methods.
*Journal of the Royal Statistical Society Series B: Statistical Methodology*, 73(2):123–214, 2011.
- Gneiting and Raftery (2007)
T. Gneiting and A. E. Raftery.
Strictly proper scoring rules, prediction, and estimation.
*Journal of the American statistical Association*, 102(477):359–378, 2007.
- Gneiting et al. (2005)
T. Gneiting, A. E. Raftery, A. H. Westveld, and T. Goldman.
Calibrated probabilistic forecasting using ensemble model output statistics and minimum crps estimation.
*Monthly Weather Review*, 133(5):1098–1118, 2005.
- Golub and Greif (2003)
G. H. Golub and C. Greif.
On solving block-structured indefinite linear systems.
*SIAM Journal on Scientific Computing*, 24(6):2076–2092, 2003.
- Gould et al. (2016)
S. Gould, B. Fernando, A. Cherian, P. Anderson, R. S. Cruz, and E. Guo.
On differentiating parameterized argmin and argmax problems with application to bi-level optimization.
*arXiv preprint arXiv:1607.05447*, 2016.
- Gould et al. (2022)
S. Gould, R. Hartley, and D. Campbell.
Deep declarative networks.
*IEEE Trans. Pattern Anal. Mach. Intell.*, 44(8):3988–4004, 2022.
- Griewank and Walther (2008)
A. Griewank and A. Walther.
*Evaluating derivatives: principles and techniques of algorithmic differentiation*.
SIAM, 2008.
- Gupta et al. (2021)
G. Gupta, X. Xiao, and P. Bogdan.
Multiwavelet-based operator learning for differential equations.
In *Advances in Neural Information Processing Systems*, volume 34, pages 24048–24062, 2021.
- Hansen et al. (2023)
D. Hansen, D. C. Maddix, S. Alizadeh, G. Gupta, and M. W. Mahoney.
Learning physical models that can respect conservation laws.
In *Proceedings of the 40th International Conference on Machine Learning*, volume 202, pages 12469–12510. PMLR, 2023.
- Harten (1997)
A. Harten.
High resolution schemes for hyperbolic conservation laws.
*Journal of computational physics*, 135(2):260–278, 1997.
- Hong et al. (2023)
I. Hong, S. Na, M. W. Mahoney, and M. Kolar.
Constrained optimization via exact augmented lagrangian and randomized iterative sketching.
In *Proceedings of the 40th International Conference on Machine Learning*, volume 202, pages 13174–13198. PMLR, 2023.
- Hoo et al. (2025)
S. B. Hoo, S. Müller, D. Salinas, and F. Hutter.
The tabular foundation model tabpfn outperforms specialized time series forecasting models based on simple features.
*arXiv preprint arXiv:2501.02945*, 2025.
- Hughes (2003)
T. J. Hughes.
*The finite element method: linear static and dynamic finite element analysis*.
Courier Corporation, 2003.
- Hyndman and Athanasopoulos (2018)
R. J. Hyndman and G. Athanasopoulos.
*Forecasting: principles and practice*.
OTexts, 2018.
- Hyndman et al. (2011)
R. J. Hyndman, R. A. Ahmed, G. Athanasopoulos, and H. L. Shang.
Optimal combination forecasts for hierarchical time series.
*Computational statistics & data analysis*, 55(9):2579–2589, 2011.
- Hyndman et al. (2025)
R. J. Hyndman, G. Athanasopoulos, A. Garza, C. Challu, M. Mergenthaler, and K. G. Olivares.
*Forecasting: Principles and Practice, the Pythonic Way*.
OTexts, Melbourne, Australia, 2025.
available at https://otexts.com/fpppy/.
- Joe (1997)
H. Joe.
*Multivariate models and multivariate dependence concepts*.
CRC press, 1997.
- Kadambi et al. (2023)
A. Kadambi, C. de Melo, C. Hsieh, M. Srivastava, and S. Soatto.
Incorporating physics into data-driven computer vision.
*Nat Mach Intell*, 5:572–580, 2023.
- Kalman (1960)
R. E. Kalman.
A new approach to linear filtering and prediction problems.
*Journal of Basic Engineering*, 82(1):35–45, 1960.
- Karniadakis et al. (2021)
G. E. Karniadakis, I. G. Kevrekidis, L. Lu, P. Perdikaris, S. Wang, and L. Yang.
Physics-informed machine learning.
*Nature Reviews Physics*, 3(6):422–440, 2021.
- Kasim and Lim (2022)
M. F. Kasim and Y. H. Lim.
Constants of motion network.
In *Advances in Neural Information Processing Systems*, volume 35, pages 25295–25305, 2022.
- Kendall and Gal (2017)
A. Kendall and Y. Gal.
What uncertainties do we need in bayesian deep learning for computer vision?
In *Proceedings of the 31st International Conference on Neural Information Processing Systems*, page 5580–5590. Curran Associates Inc., 2017.
- Kendall et al. (2018)
A. Kendall, Y. Gal, and R. Cipolla.
Multi-task learning using uncertainty to weigh losses for scene geometry and semantics.
In *2018 IEEE Conference on Computer Vision and Pattern Recognition (CVPR)*, pages 7482–7491. Computer Vision Foundation / IEEE Computer Society, 2018.
- Kim et al. (2019)
H. Kim, A. Mnih, J. Schwarz, M. Garnelo, A. Eslami, D. Rosenbaum, O. Vinyals, and Y. W. Teh.
Attentive neural processes.
*arXiv preprint arXiv:1901.05761*, 2019.
- Kingma and Ba (2015)
D. Kingma and J. Ba.
A method for stochastic optimization.
In *International Conference on Learning Representations*, 2015.
- Kochkov et al. (2024)
D. Kochkov, J. Yuval, I. Langmore, P. Norgaard, J. Smith, G. Mooers, M. Klöwer, J. Lottes, S. Rasp, P. Düben, et al.
Neural general circulation models for weather and climate.
*Nature*, 632(8027):1060–1066, 2024.
- Krantz and Parks (2002)
S. G. Krantz and H. R. Parks.
*The implicit function theorem: history, theory, and applications*.
Springer Science $\&$ Business Media, 2002.
- Krishnapriyan et al. (2021)
A. Krishnapriyan, A. Gholami, S. Zhe, R. Kirby, and M. W. Mahoney.
Characterizing possible failure modes in physics-informed neural networks.
In *Advances in neural information processing systems*, volume 34, pages 26548–26560, 2021.
- Le et al. (2005)
Q. V. Le, A. J. Smola, and S. Canu.
Heteroscedastic gaussian process regression.
In *Proceedings of the 22nd International Conference on Machine Learning*, pages 489–496. PMLR, 2005.
- LeVeque (1990)
R. J. LeVeque.
*Numerical Methods for Conservation Laws*.
Lectures in mathematics ETH Zürich. Birkhäuser Verlag, 1990.
- LeVeque (2002)
R. J. LeVeque.
*Finite volume methods for hyperbolic problems*, volume 31.
Cambridge university press, 2002.
- LeVeque (2007)
R. J. LeVeque.
*Finite difference methods for ordinary and partial differential equations: steady-state and time-dependent problems*.
SIAM, 2007.
- Li et al. (2020)
Z. Li, N. Kovachki, K. Azizzadenesheli, B. Liu, K. Bhattacharya, A. Stuart, and A. Anandkumar.
Neural operator: Graph kernel network for partial differential equations.
*arXiv preprint arXiv:2003.03485*, 2020.
- Li et al. (2021)
Z. Li, N. Kovachki, K. Azizzadenesheli, B. Liu, K. Bhattacharya, A. Stuart, and A. Anandkumar.
Fourier Neural Operator for Parametric Partial Differential Equations.
In *International Conference on Learning Representations*, 2021.
- Li et al. (2024)
Z. Li, H. Zheng, N. B. Kovachki, D. Jin, H. Chen, B. Liu, K. Azizzadenesheli, and A. Anandkumar.
Physics-informed neural operator for learning partial differential equations.
*ACM / IMS J. Data Sci.*, 1(3), 2024.
- Lim et al. (2021)
B. Lim, S. Ö. Arık, N. Loeff, and T. Pfister.
Temporal fusion transformers for interpretable multi-horizon time series forecasting.
*International Journal of Forecasting*, 37(4):1748–1764, 2021.
- Lipnikov et al. (2016)
K. Lipnikov, G. Manzini, J. D. Moulton, and M. Shashkov.
The mimetic finite difference method for elliptic and parabolic problems with a staggered discretization of diffusion coefficient.
*Journal of Computational Physics*, 305:111–126, 2016.
- Lu et al. (2021)
L. Lu, P. Jin, G. Pang, Z. Zhang, and G. E. Karniadakis.
Learning nonlinear operators via DeepONet based on the universal approximation theorem of operators.
*Nature Machine Intelligence*, 3(3):218–229, 2021.
- Maddix et al. (2018a)
D. C. Maddix, L. Sampaio, and M. Gerritsen.
Numerical artifacts in the generalized porous medium equation: Why harmonic averaging itself is not to blame.
*Journal of Computational Physics*, 361:280–298, 2018a.
- Maddix et al. (2018b)
D. C. Maddix, L. Sampaio, and M. Gerritsen.
Numerical artifacts in the discontinuous generalized porous medium equation: How to avoid spurious temporal oscillations.
*Journal of Computational Physics*, 368:277–298, 2018b.
- Matheson and Winkler (1976)
J. E. Matheson and R. L. Winkler.
Scoring rules for continuous probability distributions.
*Management science*, 22(10):1087–1096, 1976.
- Matsubara and Yaguchi (2023)
T. Matsubara and T. Yaguchi.
Finde: Neural differential equations for finding and preserving invariant quantities.
In *International Conference on Learning Representations*, 2023.
- Min et al. (2024)
Y. Min, A. Sonar, and N. Azizan.
Hard-constrained neural networks with universal approximation guarantees.
*arXiv preprint arXiv:2410.10807*, 2024.
- Mouli et al. (2024)
S. C. Mouli, D. C. Maddix, S. Alizadeh, G. Gupta, A. Stuart, M. W. Mahoney, and Y. 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, pages 36372–36418. PMLR, 2024.
- Müller (2022)
E. H. Müller.
Exact conservation laws for neural network integrators of dynamical systems.
*arXiv preprint arXiv:2209.11661*, 2022.
- Nocedal and Wright (2006)
J. Nocedal and J. Wright, Steven.
*Numerical Optimization*.
Springer, 2006.
- Négiar et al. (2023)
G. Négiar, M. W. Mahoney, and A. S. Krishnapriyan.
Learning differentiable solvers for systems with hard constraints.
In *International Conference on Learning Representations*, 2023.
- Olivares et al. (2022)
K. G. Olivares, F. Garza, D. Luo, C. Challú, M. Mergenthaler, S. Ben Taieb, S. L. Wickramasuriya, and A. Dubrawski.
HierarchicalForecast: A reference framework for hierarchical forecasting in python.
*Work in progress paper, submitted to Journal of Machine Learning Research.*, abs/2207.03517, 2022.
URL [https://arxiv.org/abs/2207.03517](https://arxiv.org/abs/2207.03517).
- Olivares et al. (2024a)
K. G. Olivares, O. N. Meetei, R. Ma, R. Reddy, M. Cao, and L. Dicker.
Probabilistic hierarchical forecasting with deep poisson mixtures.
*International Journal of Forecasting*, 40(2):470–489, 2024a.
- Olivares et al. (2024b)
K. G. Olivares, G. Négiar, R. Ma, O. N. Meetei, M. Cao, and M. W. Mahoney.
$\clubsuit$ CLOVER $\clubsuit$: Probabilistic forecasting with coherent learning objective reparameterization.
*Transactions on Machine Learning Research*, 2024b.
- Park et al. (2022)
Y. Park, D. C. Maddix, F.-X. Aubet, K. Kan, J. Gasthaus, and Y. Wang.
Learning quantile functions without quantile crossing for distribution-free time series forecasting.
In *Proceedings of The 25th International Conference on Artificial Intelligence and Statistics*, volume 151, pages 8127–8150. PMLR, 2022.
- Petropoulos et al. (2022)
F. Petropoulos, D. Apiletti, V. Assimakopoulos, M. Z. Babai, D. K. Barrow, S. B. Taieb, C. Bergmeir, R. J. Bessa, J. Bijak, J. E. Boylan, et al.
Forecasting: theory and practice.
*International Journal of Forecasting*, 38(3):705–871, 2022.
- Pfaff et al. (2021)
T. Pfaff, M. Fortunato, A. Sanchez-Gonzalez, and P. W. Battaglia.
Learning mesh-based simulation with graph networks.
In *International Conference on Learning Representations*, 2021.
- Price et al. (2025)
I. Price, A. Sanchez-Gonzalez, F. Alet, T. Andersson, A. El-Kadi, D. Masters, T. Ewalds, J. Stott, S. Mohamed, P. Battaglia, R. Lam, and M. Willson.
Probabilistic weather forecasting with machine learning.
*Nature*, 637(8044):84–90, 2025.
- Rackauckas et al. (2020)
C. Rackauckas, Y. Ma, J. Martensen, C. Warner, K. Zubov, R. Supekar, D. Skinner, A. Ramadhan, and A. Edelman.
Universal differential equations for scientific machine learning.
*arXiv preprint arXiv:2001.04385*, 2020.
- Raissi et al. (2019)
M. Raissi, P. Perdikaris, and G. E. Karniadakis.
Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations.
*Journal of Computational physics*, 378:686–707, 2019.
- Rangapuram et al. (2021)
S. S. Rangapuram, L. D. Werner, K. Benidis, P. Mercado, J. Gasthaus, and T. Januschowski.
End-to-end learning of coherent probabilistic forecasts for hierarchical time series.
In *International Conference on Machine Learning*, pages 8832–8843. PMLR, 2021.
- Rangapuram et al. (2023)
S. S. Rangapuram, S. Kapoor, R. S. Nirwan, P. Mercado, T. Januschowski, Y. Wang, and M. Bohlke-Schneider.
Coherent probabilistic forecasting of temporal hierarchies.
In *Proceedings of The 26th International Conference on Artificial Intelligence and Statistics*, volume 206, pages 9362–9376. PMLR, 2023.
- Rasmussen and Williams (2006)
C. Rasmussen and C. Williams.
*Gaussian Processes for Machine Learning*.
MIT Press, 2006.
- Rasp and Lerch (2018)
S. Rasp and S. Lerch.
Neural networks for postprocessing ensemble weather forecasts.
*Monthly Weather Review*, 146(11):3885 – 3900, 2018.
- Ravuri et al. (2021)
S. Ravuri, K. Lenc, M. Willson, D. Kangin, R. Lam, P. Mirowski, M. Fitzsimons, M. Athanassiadou, S. Kashem, S. Madge, et al.
Skilful precipitation nowcasting using deep generative models of radar.
*Nature*, 597(7878):672–677, 2021.
- Richter-Powell et al. (2022)
J. Richter-Powell, Y. Lipman, and R. T. Q. Chen.
Neural conservation laws: A divergence-free perspective.
In *Advances in Neural Information Processing Systems*, 2022.
- Robert et al. (1999)
C. P. Robert, G. Casella, and G. Casella.
*Monte Carlo statistical methods*, volume 2.
Springer, 1999.
- Rockafellar and Wets (2009)
R. T. Rockafellar and R. J.-B. Wets.
*Variational analysis*, volume 317.
Springer Science & Business Media, 2009.
- Rosen (1960)
J. B. 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.
- Rudin et al. (1992)
L. I. Rudin, S. Osher, and E. Fatemi.
Nonlinear total variation based noise removal algorithms.
*Physica D: nonlinear phenomena*, 60(1-4):259–268, 1992.
- Saad et al. (2023)
N. Saad, G. Gupta, S. Alizadeh, and D. C. Maddix.
Guiding continuous operator learning through physics-based boundary constraints.
In *International Conference on Learning Representations*, 2023.
- Salinas et al. (2019)
D. Salinas, M. Bohlke-Schneider, L. Callot, R. Medico, and J. Gasthaus.
High-dimensional multivariate forecasting with low-rank gaussian copula processes.
In *Advances in neural information processing systems*, volume 32, 2019.
- Salinas et al. (2020)
D. Salinas, V. Flunkert, J. Gasthaus, and T. Januschowski.
Deepar: Probabilistic forecasting with autoregressive recurrent networks.
*International Journal of Forecasting*, 36(3):1181–1191, 2020.
- Schein et al. (2021)
A. Schein, K. T. Carlberg, and M. J. Zahr.
Preserving general physical properties in model reduction of dynamical systems via constrained-optimization projection.
*International Journal for Numerical Methods in Engineering*, 122(14):3368–3399, 2021.
- Sethian and Strain (1992)
J. A. Sethian and J. Strain.
Crystal growth and dendritic solidification.
*Journal of Computational Physics*, 98(2):231–253, 1992.
- Skafte et al. (2019)
N. Skafte, M. Jørgensen, and S. Hauberg.
Reliable training and estimation of variance networks.
In *Advances in Neural Information Processing Systems*, volume 32, 2019.
- Stirn et al. (2023)
A. Stirn, H. Wessels, M. Schertzer, L. Pereira, N. Sanjana, and D. Knowles.
Faithful heteroscedastic regression with neural networks.
In *International Conference on Artificial Intelligence and Statistics*, pages 5593–5613. PMLR, 2023.
- Sturm and Wexler (2022)
P. O. Sturm and A. S. Wexler.
Conservation laws in a neural network architecture: enforcing the atom balance of a julia-based photochemical model (v0.2.0).
*Geosci. Model Dev.*, 15:3417–3431, 2022.
- Sun et al. (2022)
H. Sun, Y. Shi, J. Wang, H. D. Tuan, H. V. Poor, and D. Tao.
Alternating differentiation for optimization layers.
*arXiv preprint arXiv:2210.01802*, 2022.
- Szechtman and Glynn (2001)
R. Szechtman and P. W. Glynn.
Constrained monte carlo and the method of control variates.
In *Proceeding of the 2001 Winter Simulation Conference (Cat. No. 01CH37304)*, volume 1, pages 394–400. IEEE, 2001.
- Taieb et al. (2017)
S. B. Taieb, J. W. Taylor, and R. J. Hyndman.
Coherent probabilistic forecasts for hierarchical time series.
In *International Conference on Machine Learning*, pages 3348–3357. PMLR, 2017.
- Taillardat et al. (2016)
M. Taillardat, O. Mestre, M. Zamo, and P. Naveau.
Calibrated ensemble forecasts using quantile regression forests and ensemble model output statistics.
*Monthly Weather Review*, 144(6):2375–2393, 2016.
- Tezaur et al. (2017)
I. K. Tezaur, J. A. Fike, K. T. Carlberg, M. F. Barone, D. Maddix, E. E. Mussoni, and M. Balajewicz.
Advanced fluid reduced order models for compressible flow.
*Sandia National Laboratories Report, Sand No. 2017-10335*, 2017.
- Tourism Australia, Canberra (2005)
Tourism Australia, Canberra.
Tourism Research Australia (2005), Travel by Australians.
[https://www.kaggle.com/luisblanche/quarterly-tourism-in-australia/](https://www.kaggle.com/luisblanche/quarterly-tourism-in-australia/), 2005.
- Toussaint et al. (2019)
M. Toussaint, K. R. Allen, K. A. Smith, and J. B. Tenenbaum.
Differentiable physics and stable modes for tool-use and manipulation planning - extended abtract.
In *Proceedings of the Twenty-Eighth International Joint Conference on Artificial Intelligence, IJCAI-19*, pages 6231–6235. International Joint Conferences on Artificial Intelligence Organization, 2019.
- Utkarsh et al. (2024)
U. Utkarsh, V. Churavy, Y. Ma, T. Besard, P. Srisuma, T. Gymnich, A. R. Gerlach, A. Edelman, G. Barbastathis, R. D. Braatz, et al.
Automated translation and accelerated solving of differential equations on multiple gpu platforms.
*Computer Methods in Applied Mechanics and Engineering*, 419:116591, 2024.
- Utkarsh et al. (2025)
U. Utkarsh, P. Cai, A. Edelman, R. Gomez-Bombarelli, and C. V. Rackauckas.
Physics-constrained flow matching: Sampling generative models with hard constraints.
*arXiv preprint arXiv:2506.04171*, 2025.
- van der Meer et al. (2016)
J. van der Meer, J. Kraaijevanger, M. Möller, and J. Jansen.
Temporal oscillations in the simulation of foam enhanced oil recovery.
*ECMOR XV - 15th European Conference on the Mathematics of Oil Recovery*, pages 1–20, 2016.
- Van der Vaart (2000)
A. W. Van der Vaart.
*Asymptotic statistics*, volume 3.
Cambridge university press, 2000.
- Van Erven and Cugliari (2015)
T. Van Erven and J. Cugliari.
Game-theoretically optimal reconciliation of contemporaneous hierarchical time series forecasts.
In *Modeling and stochastic learning for forecasting in high dimensions*, pages 297–317. Springer, 2015.
- Vázquez (2007)
J. Vázquez.
*The Porous Medium Equation: Mathematical Theory*.
The Clarendon Press, Oxford University Press, Oxford, 2007.
- Weinberger and Saul (2009)
K. Q. Weinberger and L. K. Saul.
Distance metric learning for large margin nearest neighbor classification.
*Journal of machine learning research*, 10(2), 2009.
- Wen et al. (2018)
R. Wen, K. Torkkola, B. Narayanaswamy, and D. Madeka.
A multi-horizon quantile recurrent forecaster.
*arXiv preprint arXiv:1711.11053*, 2018.
- White et al. (2024)
A. White, A. Büttner, M. Gelbrecht, V. Duruisseaux, N. Kilbertus, F. Hellmann, and N. Boers.
Projected neural differential equations for learning constrained dynamics.
*arXiv preprint arXiv:2410.23667*, 2024.
- Wickramasuriya et al. (2019)
S. L. Wickramasuriya, G. Athanasopoulos, and R. J. H. and.
Optimal forecast reconciliation for hierarchical and grouped time series through trace minimization.
*Journal of the American Statistical Association*, 114(526):804–819, 2019.
- Willette et al. (2021)
J. Willette, H. B. Lee, J. Lee, and S. J. Hwang.
Meta learning low rank covariance factors for energy-based deterministic uncertainty.
*arXiv preprint arXiv:2110.06381*, 2021.
- Wilson (1963)
R. Wilson.
*A Simplicial Method for Convex Programming*.
PhD thesis, Harvard University, 1963.
- Woo et al. (2024)
G. Woo, C. Liu, A. Kumar, C. Xiong, S. Savarese, and D. Sahoo.
Unified training of universal time series forecasting transformers.
*arXiv preprint arXiv:2402.02592*, 2024.
- Zahr and Persson (2016)
M. Zahr and P.-O. 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.
- Zanna and Bolton (2020)
L. Zanna and T. Bolton.
Data-Driven Equation Discovery of Ocean Mesoscale Closures.
*Geophysical Research Letters*, 47(17), 2020.

## Appendix 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 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].
[Table 5](https://arxiv.org/html/2506.07003v2#A1.T5) 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 are unsuitable for PDEs and forecasting applications, which generally now require variance estimates. Hier-E2E [Rangapuram et al., 2021] and CLOVER [Olivares et al., 2024b] are 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 et al., 2023] handles nonlinear constraints and supports end-to-end training, but at 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 on 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 Training | Constraint on Dstbn. |
| --- | --- | --- | --- | --- | --- | --- |
| HardNet [Min et al., 2024] | General | Convex | ✓ | ✗ | ✓ | ✓ |
| DC3 [Donti et al., 2021] | General | Nonlinear | ✓ | ✗ | ✓ | ✓ |
| Hier-E2E [Rangapuram et al., 2021, 2023] | Forecasting | Linear | ✓ | ✓ | ✗ | ✗ |
| CLOVER [Olivares et al., 2024b] | Forecasting | Linear | ✓ | ✓ | ✗ | ✓ |
| PDE-CL [Négiar et al., 2023] | PDEs | Nonlinear | ✓ | ✗ | ✓ | ✓ |
| ProbConserv [Hansen et al., 2023] | PDEs | Linear | ✗ | ✓ | ✓ | ✓ |
| HardC [Hansen et al., 2023] | PDEs | Linear | ✗ | ✓ | ✓ | ✓ |
| ProbHardE2E (Ours) | General | Nonlinear | ✓ | ✓ | ✓ | ✓ |

Our proposed method ProbHardE2E bridges the gaps left by its predecessors. It combines the flexibility of general domain application with the ability to handle nonlinear constraints, and it 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 on distributions.

### A.1 Imposing Deterministic Constraints on 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 outputs 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 and 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 to allow flexibility in the optimization process, while still guiding the model toward feasible solutions.
These methods encourage—but do not strictly enforce—adherence to the constraints during training; and this lack of 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 $A$ to 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-descent optimization algorithms [Min et al., 2024, Donti et al., 2021].

Projection-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 and Kolter [2017], to further advance the understanding of constraint enforcement in NNs.
Projection-based methods have also been used to enforce constraints on specific architectures, e.g., neural ordinary differential equations (Neural ODEs) [Kasim and Lim, 2022, Matsubara and Yaguchi, 2023]. In particular, White et al. [2024] use a closed-form projection operator to enforce a 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 a 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 et 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 et al., 2020] and its multivariate extension DeepVAR [Salinas et 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 on mean forecasts under linear/hierarchical constraints, starting from the naive bottom-up and top-down approaches [Hyndman and Athanasopoulos, 2018]. More recently, Hyndman et al. [2011] show that the Middle-Out projection-based method yields better forecast accuracy. Since then, projection-based reconciliation methods, e.g., GTOP [Van Erven and Cugliari, 2015], MinT [Wickramasuriya et al., 2019], and ERM [Ben Taieb and Koo, 2019] have been developed. These methods leverage generic time series models, e.g., ARIMA and exponential smoothing (ETS), to derive the unconstrained mean forecast, and then they use a 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 a 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 an orthogonal projection on samples from the distribution. Hier-E2E produces coherent probabilistic forecasts without requiring explicit post-processing reconciliation. One limitation is that Hier-E2E relies on expensive sampling techniques to 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 an equality constraint completion approach similar to that in DC3 [Donti et al., 2021], rather than a projection-based approach, for satisfying the coherency 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 DeepVAR, and it shows improved forecast accuracy over Hier-E2E. As follow-up work to DPMN, Olivares et al. [2024b] propose CLOVER, a 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 [LeVeque, 2007], finite element [Hughes, 2003], and finite volume methods [LeVeque, 2002]. These numerical methods discretize the solution on a spatio-temporal mesh, and the accuracy increases as 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 et al., 2020, 2021, Gupta et al., 2021] and DeepONet [Lu et al., 2021], and message-passing Graph Neural Networks (GNNs)-based MeshGraphNets [Pfaff et al., 2021, Fortunato et 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, Karniadakis et al., 2021], which approximates the solution of a PDE as a 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 and Parks, 2002] to impose nonlinear PDE constraints directly. Chalapathi et al. [2024], extend this work by leveraging a mixture-of-experts (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 while also enforcing implicit constraints. Chen et al. [2024] propose KKT-hPINN to 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 as a hard constraint rather than as a soft constraint. A limitation of these methods is that they impose the constraints deterministically, and they do not provide estimates of the underlying variance or 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 a post-processing step at inference time, and the method only supports linear constraints.

## Appendix 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 1](https://arxiv.org/html/2506.07003v2#alg1) 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] from the case $Q=I$ to our broader framework with arbitrary $Q$. Since $Q$ is symmetric positive definite, we compute its Cholesky factorization $Q=LL^{T}$, where $L$ denotes a lower triangular matrix. We then show that if $f_{\theta}$ is a universal approximator, i.e., a sufficiently wide and deep neural network, then our DPPL preserves this universal approximation capability. Hence, ProbHardE2E retains the expressiveness of neural networks, both in its probabilistic formulation and in enforcing hard constraints.

Our DPPL in Problem [7](https://arxiv.org/html/2506.07003v2#S3.E7) is formulated in terms of projecting the samples $z_{\theta}(\phi^{(i)})\sim\mathbf{Z}_{\theta}(\phi^{(i)})$, where $\mathbf{Z}_{\theta}(\phi^{(i)})\sim\mathcal{F}(\mu_{\theta}(\phi^{(i)}),\Sigma_{\theta}(\phi^{(i)}))$ for some multivariate location-scale distribution $\mathcal{F}$, and where $({\phi}^{(i)},u^{(i)})\sim\mathcal{D}$ denotes training data from a distribution $\mathcal{D}$. The mean $\mu_{\theta}(\phi^{(i)})\in\mathbb{R}^{n}$ and covariance $\Sigma_{\theta}(\phi^{(i)})\in\mathbb{R}^{n\times n}$ are the output from a deep neural network, $f_{\theta}:\Phi\to\mathbb{R}^{k}$. The value of $k$ depends on the approximation for $\Sigma_{\theta}(\phi^{(i)})$, e.g., $k=n+n^{2}$ for dense $\Sigma_{\theta}(\phi^{(i)})$, $k=2n$ for $\Sigma_{\theta}(\phi^{(i)})=\text{diag}(\sigma_{1}^{2},\dots,\sigma_{n}^{2})$, or $k=n$ for $\Sigma_{\theta}(\phi^{(i)})=I$. For notational simplicity, we assume in this section that $f_{\theta}:\Phi\to\mathbb{R}^{n}$ corresponds to the components that output the mean $\mu_{\theta}(\phi^{(i)})$. By setting $z_{\theta}(\phi^{(i)})=\mu_{\theta}(\phi^{(i)})=f_{\theta}(\phi^{(i)})$ in Problem [7](https://arxiv.org/html/2506.07003v2#S3.E7), our DPPL can also be formulated in terms of projecting the mean $\mu_{\theta}(\phi^{(i)})$ as the following constrained least squares problem:

$$ $\displaystyle\hat{\mu}_{\theta}(\phi^{(i)}):=\operatorname*{argmin}_{\begin{subarray}{c}\tilde{\mu}_{\theta}(\phi^{(i)})\in\mathbb{R}^{n}\\ g(\tilde{\mu}_{\theta}(\phi^{(i)}))\leq 0\\ h(\tilde{\mu}_{\theta}(\phi^{(i)}))=0\end{subarray}}\lVert\tilde{\mu}_{\theta}(\phi^{(i)})-f_{\theta}(\phi^{(i)})\rVert^{2}_{Q},$ (9) $$

where $Q$ denotes a symmetric positive definite matrix and $g(\cdot)\leq 0,h(\cdot)=0$ denote the convex constraints. In particular, we show that the projected mean $\hat{\mu}_{\theta}(\phi^{(i)})$ is a universal approximator of the true solution $u\in\mathbb{R}^{n}$. We now state the theorem and provide its proof below.

###### Theorem B.1 .

Consider Problem [9](https://arxiv.org/html/2506.07003v2#A2.E9) with the projection step defined using a symmetric positive definite (SPD) matrix $Q\in\mathbb{R}^{n\times n}$, a deep neural network that is a universal approximator, $f_{\theta}:\Phi\to\mathbb{R}^{n}$, where $\Phi\subset\mathbb{R}^{m}$ denotes a compact set, convex constraints $g(\cdot)\leq 0,h(\cdot)=0$, and training data $({\phi}^{(i)},u^{(i)})\sim\mathcal{D}$ from a distribution $\mathcal{D}$.
For any continuous target function that satisfies the constraints, i.e., the true solution $u:\Phi\to\mathcal{C}\subseteq\mathbb{R}^{n}$, $u\in C(\Phi)$, where $\mathcal{C}$ denotes the convex set of feasible points defined by the convex constraints and $C(\Phi)$ denotes the space of continuous functions on $\Phi$, there exists a choice of network parameters for $f_{\theta}(\phi^{(i)})=\mu_{\theta}(\phi^{(i)})\in\mathbb{R}^{n}$, such that the projected mean, which is composition of $f_{\theta}$ with the projection step, i.e., $\Pi_{\mathcal{C}}^{Q}\bigl(f_{\theta}(\phi^{(i)})\bigr)=\hat{f}_{\theta}(\phi^{(i)})=\hat{\mu}_{\theta}(\phi^{(i)})\in\mathcal{C}\subseteq\mathbb{R}^{n}$, approximates the target function arbitrarily well, where $\hat{f}_{\theta}:\Phi\rightarrow\mathcal{C}\subseteq\mathbb{R}^{n}$. Hence, under these conditions, ProbHardE2E is a universal approximator for constrained mappings.

###### Proof.

Let $\mathcal{C}\subseteq\mathbb{R}^{n}$ denote the convex set of feasible points defined by the convex constraints. Consider the projection operator onto $\mathcal{C}$ with respect to the $Q$-norm:

$$ $\Pi_{\mathcal{C}}^{Q}(v)=\text{argmin}_{\tilde{\mu}_{\theta}(\phi^{(i)})\in\mathcal{C}}\|\tilde{\mu}_{\theta}(\phi^{(i)})-v\|^{2}_{Q}.$ (10) $$

Since $Q$ is symmetric positive definite (SPD), it has the following Cholesky factorization,

$$ $Q=LL^{\top},$ $$

where $L\in\mathbb{R}^{n\times n}$ denotes a lower triangular matrix with strictly positive diagonal entries, and hence is invertible.
By the definition of the $Q$-norm, and then substituting in its Cholesky factorization, we have

$$ $\displaystyle\|\tilde{\mu}_{\theta}(\phi^{(i)})-v\|^{2}_{Q}$ $\displaystyle=(\tilde{\mu}_{\theta}(\phi^{(i)})-v)^{\top}Q(\tilde{\mu}_{\theta}(\phi^{(i)})-v)$ (11) $\displaystyle=(\tilde{\mu}_{\theta}(\phi^{(i)})-v)^{\top}LL^{\top}(\tilde{\mu}_{\theta}(\phi^{(i)})-v)$ $\displaystyle=\big((\tilde{\mu}_{\theta}(\phi^{(i)})-v)^{\top}L)(L^{\top}(\tilde{\mu}_{\theta}(\phi^{(i)})-v)\big)$ $\displaystyle=\big(L^{\top}(\tilde{\mu}_{\theta}(\phi^{(i)})-v))^{\top}(L^{\top}(\tilde{\mu}_{\theta}(\phi^{(i)})-v)\big)$ $\displaystyle=\|L^{\top}(\tilde{\mu}_{\theta}(\phi^{(i)})-v)\|^{2}_{2}.$ $$

This shows that the $Q$-norm is equivalent to the standard Euclidean norm after the linear transformation $L^{\top}$.

We define the invertible linear mapping $\Psi:\mathbb{R}^{n}\to\mathbb{R}^{n}$ by $\Psi(v)=L^{\top}v$.
Then using [Equation 11](https://arxiv.org/html/2506.07003v2#A2.E11), the $Q$-norm in [Equation 10](https://arxiv.org/html/2506.07003v2#A2.E10) can be written as the Euclidean norm as follows:

$$ $\displaystyle\Pi_{\mathcal{C}}^{Q}(v)$ $\displaystyle=\text{argmin}_{\tilde{u}_{\theta}(\phi^{(i)})\in\mathcal{C}}\|L^{T}(\tilde{\mu}_{\theta}(\phi^{(i)})-v)\|^{2}_{2}$ (12) $\displaystyle=\text{argmin}_{\tilde{u}_{\theta}(\phi^{(i)})\in\mathcal{C}}\|\Psi(\tilde{\mu}_{\theta}(\phi^{(i)}))-\Psi(v)\|^{2}_{2}$ $\displaystyle=\Psi^{-1}\Bigl(\text{argmin}_{w\in\Psi(\mathcal{C})}\|w-\Psi(v)\|^{2}_{2}\Bigr),$ $$

where $w=\Psi(\tilde{\mu}_{\theta}(\phi^{(i)}))$.
Hence, the projection can be expressed as the Euclidean projection onto the transformed set $\Psi(\mathcal{C})$. It is well known that the Euclidean projection onto a closed convex set is nonexpansive and is Lipschitz continuous. (See, e.g., Min et al. [2024].)

Now, suppose that $f_{\theta}(\phi^{(i)})$ is a deep neural network that is a universal approximator, i.e., for any continuous function
$u:\Phi\to\mathbb{R}^{n}$, $u\in C(\Phi)$,
and for any $\epsilon>0$, there exists parameters $\theta$ such that

$$ $\sup_{\phi^{(i)}\in\Phi}\|u(\phi^{(i)})-f_{\theta}(\phi^{(i)})\|<\epsilon,$ $$

where $\Phi\subset\mathbb{R}^{m}$ denotes a compact set and $C(\Phi)$ denotes the space of continuous functions on $\Phi$.
Let $u:\Phi\to\mathcal{C}\subseteq\mathbb{R}^{n}$, $u\in C(\Phi)$, be any continuous target function whose outputs satisfy the constraints.
Since $\Pi_{\mathcal{C}}^{Q}$ is continuous (as the composition of the continuous mapping $\Psi$, the Euclidean projection onto $\Psi(\mathcal{C})$, and $\Psi^{-1}$), it follows by the universal approximation theorem and properties of continuous functions that the projected mean $\Pi_{\mathcal{C}}^{Q}\bigl(f_{\theta}(\phi^{(i)})\bigr)=\hat{\mu}(\phi^{(i)})$ can uniformly approximate $u(\phi^{(i)})$ arbitrarily well on
$\Phi$. In other words, for every $\epsilon>0$, there exists a choice of network parameters $\theta$ such that

$$ $\sup_{\phi^{(i)}\in\Phi}\|u(\phi^{(i)})-\Pi_{\mathcal{C}}^{Q}(f_{\theta}(\phi^{(i)}))\|<\epsilon.$ $$

Thus, the composition of the neural network $f_{\theta}$ with the $Q$-norm projection retains the universal approximation property for any continuous target function satisfying the constraints.
∎

###### Theorem B.1 .

###### Proof.

## Appendix 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 ([Equation 6](https://arxiv.org/html/2506.07003v2#S3.E6)) in ProbHardE2E for various constraint types, which are summarized in [Table 1](https://arxiv.org/html/2506.07003v2#S3.T1).

### C.1 Linear Equality Constraints

In this subsection, we provide the closed-form expressions for the constrained posterior distribution parameters, i.e., the mean $\hat{\mu}$ and covariance $\hat{\Sigma}$ in [Equation 8](https://arxiv.org/html/2506.07003v2#S3.E8),
from the DPPL in ProbHardE2E for linear equality constraints. Linear equality constraints occur in a wide range of applications, including coherency constraints in hierarchical time series forecasting [Hyndman et al., 2011, Rangapuram et al., 2021, Petropoulos et al., 2022, Olivares et al., 2024b], divergence-free 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(\hat{u})=A\hat{u}-b=0$, with $A\in\mathbb{R}^{q\times n}$, with full row rank $q$, where $q\leq n$, and $b\in\mathbb{R}^{q}$, the optimal solution $u^{*}$ to Problem [7](https://arxiv.org/html/2506.07003v2#S3.E7) is given as
$u^{*}(z)=P_{Q^{-1}}z+(I-P_{Q^{-1}})A^{{\dagger}}b$,
where
$P_{Q^{-1}}=I-Q^{-1}A^{\top}(AQ^{-1}A^{\top})^{-1}A,$
denotes an oblique projection operator, and $A^{{\dagger}}$ denotes the Moore-Penrose inverse.
In addition, if $\mathbf{Z}\sim\mathcal{F}(\mu,\Sigma)$ and $z\sim\mathbf{Z}$ for multivariate, location-scale distribution $\mathcal{F}$, then $u^{*}\sim\mathbf{Y}$, where $\mathbf{Y}\sim\mathcal{F}(\hat{\mu},\hat{\Sigma})$ and $\hat{\mu},\hat{\Sigma}$ are given in [Equation 8](https://arxiv.org/html/2506.07003v2#S3.E8) with $\mathcal{T}(z)=u^{*}(z)$, which simplifies to the closed-form expressions, $\hat{\mu}=P_{Q^{-1}}\mu+(I-P_{Q^{-1}})A^{{\dagger}}b$ and $\hat{\Sigma}=P_{Q^{-1}}\Sigma P^{\top}_{Q^{-1}}$.

###### Proof.

Using the Lagrange multiplier $\lambda\in\mathbb{R}^{q}$, we can form the Lagrangian of Problem [7](https://arxiv.org/html/2506.07003v2#S3.E7) with linear constraints to obtain:

$$ $L(\hat{u},\lambda;z)=\frac{1}{2}\hat{u}^{\top}Q\hat{u}-z^{\top}Q\hat{u}+\lambda^{\top}(A\hat{u}-b).$ $$

The sufficient optimality conditions to obtain $(u^{*},\lambda^{*})$ are the first-order gradient conditions:

$$ $\displaystyle\nabla_{\hat{u}}L(\hat{u},\lambda;z)|_{u^{*},\lambda^{*}}=Qu^{*}-Q^{\top}z+A^{\top}\lambda^{*}=0,$ (13a) $\displaystyle\nabla_{\lambda}L(\hat{u},\lambda;z)|_{u^{*},\lambda^{*}}=Au^{*}-b=0.$ (13b) $$

Since $Q$ is SPD, $Q=Q^{\top}$ and $Q^{-1}$ exists. Then from [Equation 13a](https://arxiv.org/html/2506.07003v2#A3.E13.1), we obtain:

$$ $\displaystyle Q(u^{*}-z)+A^{\top}\lambda^{*}=0,$ $$

which simplifies to the following expression for $u^{*}$:

$$ $u^{*}=z-Q^{-1}A^{\top}\lambda^{*}.$ (14) $$

We solve [Equation 13b](https://arxiv.org/html/2506.07003v2#A3.E13.2) for $u^{*}$ using the Moore-Penrose inverse, i.e., $u^{*}=A^{{\dagger}}b$, where $A^{{\dagger}}=A^{\top}(AA^{\top})^{-1}$. Note that $AA^{\top}\in\mathbb{R}^{q\times q}$ is invertible with full rank $q$ since $A\in\mathbb{R}^{q\times n}$ has full row rank $q\leq n$. Substituting this expression into [Equation 14](https://arxiv.org/html/2506.07003v2#A3.E14) for $u^{*}$ gives:

$$ $\displaystyle A^{\top}(AA^{\top})^{-1}b=z-Q^{-1}A^{\top}\lambda^{*}.$ $$

Rearranging for the optimal Lagrange multiplier $\lambda^{*}$, and multiplying both sides by $A$ gives:

$$ $\displaystyle(AQ^{-1}A^{\top})\lambda^{*}=Az-\underbrace{(AA^{\top})(AA^{\top})^{-1}}_{I}b.$ $$

Now, $AQ^{-1}A^{\top}\in\mathbb{R}^{q\times q}$ is invertible since $A$ has full row rank $q$.
Then we obtain:

$$ $\displaystyle\lambda^{*}=(AQ^{-1}A^{\top})^{-1}(Az-b).$ $$

Substituting in the expression for $\lambda^{*}$ into [Equation 14](https://arxiv.org/html/2506.07003v2#A3.E14) gives the following expression for the optimal solution:

$$ $\displaystyle u^{*}$ $\displaystyle=z-Q^{-1}A^{\top}(AQ^{-1}A^{\top})^{-1}(Az-b),$ (15) $\displaystyle=(I-Q^{-1}A^{\top}(AQ^{-1}A^{\top})^{-1}A)z+Q^{-1}A^{\top}(AQ^{-1}A^{\top})^{-1}b.$ $$

Let

$$ $P_{Q^{-1}}=I-Q^{-1}A^{\top}(AQ^{-1}A^{\top})^{-1}A,$ (16) $$

be an oblique projection.
To see that this is a projection, observe that

$$ $\displaystyle P_{Q^{-1}}^{2}$ $\displaystyle=(I-Q^{-1}A^{\top}(AQ^{-1}A^{\top})^{-1}A)(I-Q^{-1}A^{\top}(AQ^{-1}A^{\top})^{-1}A)$ $\displaystyle=I-2Q^{-1}A^{\top}(AQ^{-1}A^{\top})^{-1}A+Q^{-1}A^{\top}(AQ^{-1}A^{\top})^{-1}AQ^{-1}A^{\top}(AQ^{-1}A^{\top})^{-1}A$ $\displaystyle=I-2Q^{-1}A^{\top}(AQ^{-1}A^{\top})^{-1}A+Q^{-1}A^{\top}(AQ^{-1}A^{\top})^{-1}A$ $\displaystyle=I-Q^{-1}A^{\top}(AQ^{-1}A^{\top})^{-1}A$ $\displaystyle=P_{Q^{-1}}.$ $$

Then, the expression for $u^{*}$ in [Equation 15](https://arxiv.org/html/2506.07003v2#A3.E15) simplifies to:

$$ $\displaystyle u^{*}(z)$ $\displaystyle=P_{Q^{-1}}z+Q^{-1}A^{\top}(AQ^{-1}A^{\top})^{-1}(AA^{{\dagger}})b,$ (17) $\displaystyle=P_{Q^{-1}}z+(Q^{-1}A^{\top}(AQ^{-1}A^{\top})^{-1}A)A^{{\dagger}}b,$ $\displaystyle=P_{Q^{-1}}z+(I-{P}_{Q^{-1}})A^{{\dagger}}b,$ $$

since $AA^{{\dagger}}=AA^{\top}(AA^{\top})^{-1}=I$.

Since the expression for $u^{*}$ in [Equation 17](https://arxiv.org/html/2506.07003v2#A3.E17) is a linear transformation $\mathcal{T}$ of $z\sim\mathcal{F}(\mu,\Sigma)$, we can use [Theorem 3.1](https://arxiv.org/html/2506.07003v2#S3.Thmtheorem1) with $\mathcal{T}(z)=u^{*}(z)$ to write the expression for $u^{*}\sim\mathcal{F}(\hat{\mu},\hat{\Sigma})$, where:

$$ $\displaystyle\hat{\mu}=\mathcal{T}(\mu)=u^{*}(\mu)=P_{Q^{-1}}\mu+(I-P_{Q^{-1}})A^{{\dagger}}b,$ (18a) $\displaystyle\hat{\Sigma}=J_{\mathcal{T}}(\mu)\Sigma J_{\mathcal{T}}(\mu)^{\top}=P_{Q^{-1}}\Sigma P^{\top}_{Q^{-1}}.$ (18b) $$

It can easily be verified that $J_{\mathcal{T}}(\mu)=P_{Q^{-1}}$ by differentiating [Equation 17](https://arxiv.org/html/2506.07003v2#A3.E17) with respect to $z$. We note that [Equation 18](https://arxiv.org/html/2506.07003v2#A9.EGx9) holds exactly in the case of linear constraints since $\mathcal{T}$ is a linear transformation of $z$.
∎

###### Proposition C.1 .

###### Proof.

### C.2 Nonlinear Equality Constraints

In this subsection, we describe how to compute the DPPL in ProbHardE2E for general nonlinear equality constraints. Nonlinear equality constraints naturally arise in applications that involve structural, physical, or geometric consistency.
These include closed-loop kinematics in robotics [Toussaint et al., 2019], nonlinear conservation laws [LeVeque, 1990] in PDE-constrained surrogate modeling [Biegler et al., 2003, Zahr and Persson, 2016, Négiar et al., 2023] with applications in climate modeling [Bolton and Zanna, 2019, Zanna and Bolton, 2020, Beucler et al., 2021], compressible flows in aerodynamics [Tezaur et al., 2017] and atomic modeling [Müller, 2022, Sturm and Wexler, 2022].

###### Proposition C.2 .

For nonlinear equality constraints, $h(\hat{u})=0\in\mathbb{R}^{q}$, where $h:\mathbb{R}^{n}\to\mathbb{R}^{q}$, the optimal solution $u^{*}(z)$ to Problem [7](https://arxiv.org/html/2506.07003v2#S3.E7) forms a pair $(u^{*}(z),\lambda^{*})$ which satisfies $u^{*}(z)=z-Q^{-1}\nabla h(u^{*}(z))^{\top}\lambda^{*}$ and $h(u^{*}(z))=0$.
In addition, if $\mathbf{Z}\sim\mathcal{F}(\mu,\Sigma)$ and $z\sim\mathbf{Z}$ for multivariate, location-scale distribution $\mathcal{F}$, then $u^{*}\sim\mathbf{Y}$, where $\mathbf{Y}\sim\mathcal{F}(\hat{\mu},\hat{\Sigma})$ and $\hat{\mu},\hat{\Sigma}$ are given in [Equation 8](https://arxiv.org/html/2506.07003v2#S3.E8) with $\mathcal{T}(z)=u^{*}(z)$.

###### Proof.

Using the Lagrange multiplier $\lambda\in\mathbb{R}^{q}$, we can form the Lagrangian of Problem [7](https://arxiv.org/html/2506.07003v2#S3.E7) with nonlinear equality constraints to obtain:

$$ $L(\hat{u},\lambda;z)=\frac{1}{2}\hat{u}^{\top}Q\hat{u}-z^{\top}Q\hat{u}+\lambda^{\top}h(\hat{u}).$ $$

The sufficient optimality conditions to obtain $(u^{*},\lambda^{*})$ are the first-order gradient conditions:

$$ $R(u^{*},\lambda^{*};z)\;=\;\left\{\begin{aligned} &\nabla_{\hat{u}}L(\hat{u},\lambda;z)|_{u^{*},\lambda^{*}}=Q(u^{*}-z)+\nabla h(u^{*})\lambda^{*}=0,\\ &\nabla_{\lambda}L(\hat{u},\lambda;z)|_{u^{*},\lambda^{*}}=h(u^{*})=0.\end{aligned}\right.$ (19) $$

We solve [Equation 19](https://arxiv.org/html/2506.07003v2#A3.E19) via root-finding methods, e.g., Newton’s method for $(u^{*},\lambda^{*})$ to obtain $u^{*}(z)\;=\;\arg\{\,\hat{u}:R(\hat{u},\lambda^{*};z)=0\}\,$, where the root-finding solution $u^{*}$ is implicitly dependent on $z$.
Since the expression for $u^{*}$ is a nonlinear transformation $\mathcal{T}$ of $z\sim\mathcal{F}(\mu,\Sigma)$, we can use [Theorem 3.1](https://arxiv.org/html/2506.07003v2#S3.Thmtheorem1) with $\mathcal{T}(z)=u^{*}(z)$ to write the expression for $u^{*}\sim\mathcal{F}(\hat{\mu},\hat{\Sigma})$, where:

$$ $\displaystyle\hat{\mu}=\mathcal{T}(\mu)=u^{*}(\mu),$ (20a) $\displaystyle\hat{\Sigma}=J_{\mathcal{T}}(\mu)\Sigma J_{\mathcal{T}}(\mu)^{\top},$ (20b) $$

hold to first-order accuracy. In the following [Proposition C.3](https://arxiv.org/html/2506.07003v2#A3.Thmproposition3), we detail the iterative algorithm to compute the terms $u^{*}(\mu)$ and $J_{\mathcal{T}}(\mu)$ in [Equation 20](https://arxiv.org/html/2506.07003v2#A9.EGx10).
∎

###### Proposition C.3 .

Let $h(\hat{u})=0\in\mathbb{R}^{q}$ be a smooth nonlinear equality constraint, where $h:\mathbb{R}^{n}\to\mathbb{R}^{q}$. Consider the constrained projection problem from Problem [7](https://arxiv.org/html/2506.07003v2#S3.E7) with $z=\mu$:

$$ $u^{*}(\mu)=\operatorname*{argmin}_{\begin{subarray}{c}\hat{u}\in\mathbb{R}^{n}\\ h(\hat{u})=0\end{subarray}}\ f(\hat{u}),$ (21) $$

where $Q\succ 0$ and $f(\hat{u})=\frac{1}{2}\ \|\hat{u}-\mu\|_{Q}^{2}$ denotes our quadratic objective.

- 1.
At each iteration, we solve the linearized Karush-Kuhn-Tucker (KKT) system using the Schur complement to obtain:
$\displaystyle\lambda^{(i+1)}$
$\displaystyle=\left(J^{(i)}Q^{-1}J^{(i)\top}\right)^{-1}\left(h(\hat{u}^{(i)})-J^{(i)}(\hat{u}^{(i)}-\mu)\right),$
(22a)
$\displaystyle\hat{u}^{(i+1)}$
$\displaystyle=\mu-Q^{-1}J^{(i)\top}\lambda^{(i+1)},$
(22b)
where $J^{(i)}=\nabla h(\hat{u}^{(i)})^{\top}\in\mathbb{R}^{q\times n}$.
At the first iteration with $\hat{u}^{(0)}=\mu$, [Equation 22b](https://arxiv.org/html/2506.07003v2#A3.E22.2) simplifies to:
$\hat{u}^{(1)}=\mu-Q^{-1}J^{\top}\left(JQ^{-1}J^{\top}\right)^{-1}h(\mu),$
(23)
where $J=\nabla h(\mu)^{\top}$.
- 2.
At convergence, the Jacobian $J_{\mathcal{T}}(\mu)$ of the projection map $\mathcal{T}(\mu):=u^{*}(\mu)$ is given by:
$\displaystyle J_{\mathcal{T}}(\mu):=\frac{\partial u^{*}(\mu)}{\partial\mu}$
$\displaystyle=I-Q^{-1}J^{*\top}(J^{*}Q^{-1}J^{*\top})^{-1}J^{*}\in\mathbb{R}^{n\times n},$
(24)

where $J^{*}=\nabla h(u^{*})^{\top}$.

###### Proof.

We begin with the matrix form of the KKT system derived in [Equation 19](https://arxiv.org/html/2506.07003v2#A3.E19):

$$ $R(u^{*},\lambda^{*};\mu)=\begin{bmatrix}\nabla_{\hat{u}}L(u^{*},\lambda^{*})\\ \nabla_{\lambda}L(u^{*},\lambda^{*})\end{bmatrix}=\begin{bmatrix}\nabla f(u^{*})+J^{*\top}\lambda^{*}\\ h(u^{*})\end{bmatrix}=0,$ (25) $$

with quadratic objective $f$ defined in Problem [21](https://arxiv.org/html/2506.07003v2#A3.E21).

1. Iteration Update. We use Newton’s Method to linearize the KKT system in [Equation 25](https://arxiv.org/html/2506.07003v2#A3.E25) evaluated at $(\hat{u}^{(i+1)},\lambda^{(i+1)})$ about the past iterate $(\hat{u}^{(i)},\lambda^{(i)})$. For the stationarity condition, which is the first component of $R(u^{*},\lambda^{*};\mu)$ in [Equation 25](https://arxiv.org/html/2506.07003v2#A3.E25),
we use the first-order Taylor expansion of $R_{0}(\hat{u}^{(i+1)},\lambda^{(i+1)};\mu)$ about the past iterate $(\hat{u}^{(i)},\lambda^{(i)})$ to obtain:

$$ $\displaystyle R_{0}(\hat{u}^{(i+1)},\lambda^{(i+1)};\mu)$ $\displaystyle=R_{0}(\hat{u}^{(i)},\lambda^{(i)};\mu)+\nabla_{\hat{u},\lambda}R_{0}(\hat{u}^{(i)},\lambda^{(i)};\mu)^{\top}\begin{bmatrix}\Delta\hat{u}^{(i+1)}\\ \Delta\lambda^{(i+1)}\end{bmatrix}$ (26) $\displaystyle=\nabla L_{\hat{u}}(\hat{u}^{(i)},\lambda^{(i)})+\nabla_{\hat{u},\lambda}(\nabla f(\hat{u}^{(i)})+J^{(i)\top}\lambda^{(i)})^{\top}\begin{bmatrix}\Delta\hat{u}^{(i+1)}\\ \Delta\lambda^{(i+1)}\end{bmatrix}$ $\displaystyle=\nabla L_{\hat{u}}(\hat{u}^{(i)},\lambda^{(i)})+\begin{bmatrix}(\nabla^{2}f(\hat{u}^{(i)})+\nabla^{2}h(\hat{u}^{(i)})\lambda^{(i)})&J^{(i)\top}\end{bmatrix}\begin{bmatrix}\Delta\hat{u}^{(i+1)}\\ \Delta\lambda^{(i+1)}\end{bmatrix}$ $\displaystyle=\nabla L_{\hat{u}}(\hat{u}^{(i)},\lambda^{(i)})+\begin{bmatrix}\nabla_{\hat{u}\hat{u}}^{2}L(\hat{u}^{(i)},\lambda^{(i)})&J^{(i)\top}\end{bmatrix}\begin{bmatrix}\Delta\hat{u}^{(i+1)}\\ \Delta\lambda^{(i+1)}\end{bmatrix}=0,$ $$

where $\Delta\hat{u}^{(i+1)}=\hat{u}^{(i+1)}-\hat{u}^{(i)}$, $\Delta\lambda^{(i+1)}=\lambda^{(i+1)}-\lambda^{(i)}$, and $\nabla^{2}h(\hat{u}^{(i)})\in\mathbb{R}^{n\times n\times q}$ denotes the Hessian of the constraints.
Solving for the increments we obtain:

$$ $\displaystyle\begin{bmatrix}\nabla_{\hat{u}\hat{u}}^{2}L(\hat{u}^{(i)},\lambda^{(i)})\quad J^{(i)\top}\end{bmatrix}\begin{bmatrix}\Delta\hat{u}^{(i+1)}\\ \Delta\lambda^{(i+1)}\end{bmatrix}=-\nabla_{\hat{u}}L(\hat{u}^{(i)},\lambda^{(i)}).$ (27) $$

For the feasibility condition, i.e., the second component of $R(u^{*},\lambda^{*};\mu)$ in [Equation 25](https://arxiv.org/html/2506.07003v2#A3.E25), we also linearize the constraint as:

$$ $\displaystyle R_{1}(\hat{u}^{(i+1)},\lambda^{(i+1)};\mu)=h(\hat{u}^{(i+1)})$ $\displaystyle=h(\hat{u}^{(i)})+J^{(i)}(\hat{u}^{(i+1)}-\hat{u}^{(i)})=0,$ (28) $$

using first-order Taylor expansion (Newton’s Method).
Then,

$$ $\displaystyle J^{(i)}\Delta\hat{u}^{(i+1)}=-h(\hat{u}^{(i)}).$ (29) $$

We can then combine [Equation 27](https://arxiv.org/html/2506.07003v2#A3.E27) and [Equation 29](https://arxiv.org/html/2506.07003v2#A3.E29) to form the following linearized system of KKT conditions:

$$ $\begin{bmatrix}\nabla_{\hat{u}\hat{u}}^{2}L(\hat{u}^{(i)},\lambda^{(i)})&J^{(i)\top}\\ J^{(i)}&0\end{bmatrix}\begin{bmatrix}\Delta\hat{u}^{(i+1)}\\ \Delta\lambda^{(i+1)}\end{bmatrix}=-\begin{bmatrix}\nabla_{\hat{u}}L(\hat{u}^{(i)},\lambda^{(i)})\\ h(\hat{u}^{(i)})\end{bmatrix}.$ (30) $$

Note that the system of equations in [Equation 30](https://arxiv.org/html/2506.07003v2#A3.E30) is used in Sequential Quadratic Programming (SQP) [Wilson, 1963, Nocedal and Wright, 2006, Gill and Wong, 2012] when there are no 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, [Equation 30](https://arxiv.org/html/2506.07003v2#A3.E30) gives the standard unconstrained Newton step $\nabla^{2}f(\hat{u}^{(i)})\Delta\hat{u}^{(i+1)}=-\nabla f(\hat{u}^{(i)})$ when $h=0$.

Now, we use our quadratic objective $f$ in Problem [21](https://arxiv.org/html/2506.07003v2#A3.E21) to compute:

$$ $\displaystyle\nabla_{\hat{u}}L(\hat{u}^{(i)},\lambda^{(i)})$ $\displaystyle=Q(\hat{u}^{(i)}-\mu)+J^{(i)\top}\lambda^{(i)},$ (31) $\displaystyle\nabla_{\hat{u}\hat{u}}^{2}L(\hat{u}^{(i)},\lambda^{(i)})$ $\displaystyle=Q+\nabla^{2}h(\hat{u}^{(i)})\lambda^{(i)}\approx Q.$ $$

Note that $Q$ is symmetric positive definite, but $\nabla^{2}h(\hat{u}^{(i)})$ 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 to ensure that $\nabla^{2}h(\hat{u}^{(i)})$ is positive semi-definite. In addition, since $\nabla^{2}h(\hat{u}^{(i)})\in\mathbb{R}^{n\times n\times q}$ is a 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 and Walther, 2008, Blondel and Roulet, 2024]. Similar to the Gauss-Newton method [Björck, 1996, Nocedal and Wright, 2006] for nonlinear least squares problems, we assume that the constraint $h$ is approximately affine near its optimal point $u^{*}$, and use only first-order constraint information. Hence, we set $\nabla^{2}h(\hat{u}^{(i)})\approx 0$. We note that even with these approximations for efficiency on large-scale problems, we still show strong performance in the nonlinear constraint results in [Table 4](https://arxiv.org/html/2506.07003v2#S4.T4). An alternate approach could be to use a low-rank approximation to the Hessian as done in Quasi-Newton, e.g., BFGS methods [Nocedal and Wright, 2006].

Using [Equation 31](https://arxiv.org/html/2506.07003v2#A3.E31) with setting $\nabla^{2}h(\hat{u}^{(i)})=0$, [Equation 30](https://arxiv.org/html/2506.07003v2#A3.E30) simplifies to:

$$ $\begin{bmatrix}Q&J^{(i)\top}\\ J^{(i)}&0\end{bmatrix}\begin{bmatrix}\Delta\hat{u}^{(i+1)}\\ \Delta\lambda^{(i+1)}\end{bmatrix}=-\begin{bmatrix}Q(\hat{u}^{(i)}-\mu)+J^{(i)\top}\lambda^{(i)}\\ h(\hat{u}^{(i)})\end{bmatrix}.$ (32) $$

Then,

$$ $\displaystyle Q\Delta\hat{u}^{(i+1)}+J^{(i)\top}(\lambda^{(i+1)}-\lambda^{(i)})$ $\displaystyle=-Q(\hat{u}^{(i)}-\mu)-J^{(i)\top}\lambda^{(i)},$ (33a) $\displaystyle J^{(i)}\Delta\hat{u}^{(i+1)}$ $\displaystyle=-h(\hat{u}^{(i)}).$ (33b) $$

We see that the only terms involving $\lambda^{(i)}$ 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 $\lambda^{(i)}=0$ at each iteration, and we compute $\lambda^{(i+1)}$ only for computing the primal update in [Equation 22b](https://arxiv.org/html/2506.07003v2#A3.E22.2).

Since $Q\succ 0$, it is invertible, we can multiply [Equation 33a](https://arxiv.org/html/2506.07003v2#A3.E33.1) by $Q^{-1}$ to obtain:

$$ $\displaystyle\Delta\hat{u}^{(i+1)}$ $\displaystyle=\hat{u}^{(i+1)}-\hat{u}^{(i)}=-(\hat{u}^{(i)}-\mu)-Q^{-1}J^{(i)\top}\lambda^{(i+1)}.$ (34) $$

Multiplying both sides of [Equation 34](https://arxiv.org/html/2506.07003v2#A3.E34) by $J^{(i)}$ and using [Equation 33b](https://arxiv.org/html/2506.07003v2#A3.E33.2), we can eliminate $\Delta\hat{u}^{(i+1)}$ to obtain:

$$ $-h(\hat{u}^{(i)})=-J^{(i)}(\hat{u}^{(i)}-\mu)-J^{(i)}Q^{-1}J^{(i)\top}\lambda^{(i+1)}.$ (35) $$

Since $Q\succ 0$, $Q^{-1}\succ 0$ and then $J^{(i)}Q^{-1}J^{(i)\top}\succ 0$, and hence it is invertible. We can then solve [Equation 35](https://arxiv.org/html/2506.07003v2#A3.E35) for $\lambda^{(i+1)}$ to obtain:

$$ $\lambda^{(i+1)}=\left(J^{(i)}Q^{-1}J^{(i)\top}\right)^{-1}\left(h(\hat{u}^{(i)})-J^{(i)}(\hat{u}^{(i)}-\mu)\right),$ (36) $$

which gives the desired [Equation 22a](https://arxiv.org/html/2506.07003v2#A3.E22.1).
Then, solving [Equation 34](https://arxiv.org/html/2506.07003v2#A3.E34) for $\hat{u}^{(i+1)}$ gives:

$$ $u^{(i+1)}=\mu-Q^{-1}J^{(i)\top}\lambda^{(i+1)},$ (37) $$

which is the desired [Equation 22b](https://arxiv.org/html/2506.07003v2#A3.E22.2). Lastly, for the first iterate, substituting [Equation 36](https://arxiv.org/html/2506.07003v2#A3.E36) into [Equation 37](https://arxiv.org/html/2506.07003v2#A3.E37), setting $i=0$ and $u^{(0)}=\mu$ gives the desired [Equation 23](https://arxiv.org/html/2506.07003v2#A3.E23).

We can also solve [Equation 32](https://arxiv.org/html/2506.07003v2#A3.E32) efficiently using block Gaussian elimination and the Schur complement $(J^{(i)}Q^{-1}J^{(i)\top})^{-1}\in\mathbb{R}^{q\times q}$  [Golub and Greif, 2003]. In particular, the block matrix in [Equation 32](https://arxiv.org/html/2506.07003v2#A3.E32) can be factored into a product of elementary matrices as:

$$ $\begin{bmatrix}Q&J^{(i)\top}\\ J^{(i)}&0\end{bmatrix}=\begin{bmatrix}I&0\\ J^{(i)}Q^{-1}&I\end{bmatrix}\begin{bmatrix}Q^{-1}&0\\ 0&-J^{(i)}Q^{-1}J^{(i)\top}\end{bmatrix}\begin{bmatrix}I&Q^{-1}J^{(i)\top}\\ 0&I\end{bmatrix}.$ (38) $$

Since the matrix factorization in [Equation 38](https://arxiv.org/html/2506.07003v2#A3.E38) is a product of elementary matrices and a diagonal matrix, we can easily compute its inverse as:

$$ $\begin{bmatrix}Q&J^{(i)\top}\\ J^{(i)}&0\end{bmatrix}^{-1}=\begin{bmatrix}I&-Q^{-1}J^{(i)\top}\\ 0&I\end{bmatrix}\begin{bmatrix}Q^{-1}&0\\ 0&-(J^{(i)}Q^{-1}J^{(i)\top})^{-1}\end{bmatrix}\begin{bmatrix}I&0\\ -J^{(i)}Q^{-1}&I\end{bmatrix}$ (39) $$

Then multiplying by the right-hand side in [Equation 32](https://arxiv.org/html/2506.07003v2#A3.E32) gives the solution:

$$ $\displaystyle\begin{bmatrix}\hat{u}^{(i+1)}\\ \lambda^{(i+1)}\end{bmatrix}$ $\displaystyle=\begin{bmatrix}\hat{u}^{(i)}\\ 0\end{bmatrix}-\begin{bmatrix}I&-Q^{-1}J^{(i)\top}\\ 0&I\end{bmatrix}\begin{bmatrix}Q^{-1}&0\\ 0&-(J^{(i)}Q^{-1}J^{(i)\top})^{-1}\end{bmatrix}\begin{bmatrix}I&0\\ -J^{(i)}Q^{-1}&I\end{bmatrix}\begin{bmatrix}Q(\hat{u}^{(i)}-\mu)\\ h(\hat{u}^{(i)})\end{bmatrix}$ $\displaystyle=\begin{bmatrix}\hat{u}^{(i)}\\ 0\end{bmatrix}-\begin{bmatrix}I&-Q^{-1}J^{(i)\top}\\ 0&I\end{bmatrix}\begin{bmatrix}Q^{-1}&0\\ 0&-(J^{(i)}Q^{-1}J^{(i)\top})^{-1}\end{bmatrix}\begin{bmatrix}Q(\hat{u}^{(i)}-\mu)\\ h(\hat{u}^{(i)})-J^{(i)}(\hat{u}^{(i)}-\mu)\end{bmatrix}$ $\displaystyle=\begin{bmatrix}\hat{u}^{(i)}\\ 0\end{bmatrix}-\begin{bmatrix}I&-Q^{-1}J^{(i)\top}\\ 0&I\end{bmatrix}\begin{bmatrix}\hat{u}^{(i)}-\mu\\ -(J^{(i)}Q^{-1}J^{(i)\top})^{-1}(h(\hat{u}^{(i)})-J^{(i)}(\hat{u}^{(i)}-\mu))\end{bmatrix}$ $\displaystyle=\begin{bmatrix}\mu-Q^{-1}J^{(i)\top}\lambda^{(i+1)}\\ (J^{(i)}Q^{-1}J^{(i)\top})^{-1}(h(\hat{u}^{(i)})-J^{(i)}(\hat{u}^{(i)}-\mu))\end{bmatrix}.$ $$

Using the Schur complement reduces the Newton system from an indefinite $(n+q)\times(n+q)$ solve to a $n\times n$ SPD solve with $Q^{-1}$ and $q\times q$ SPD solve with the Schur complement, where $q\leq n$. Similarly, the Jacobian expression in [Equation 24](https://arxiv.org/html/2506.07003v2#A3.E24), 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 a full $(n+q)\times(n+q)$ saddle-point or indefinite matrix.

2. Jacobian $J_{\mathcal{T}}(\mu)$. Here, we compute the Jacobian $J_{\mathcal{T}}(\mu):=\partial u^{*}(\mu)/\partial\mu$ of the transformation $\mathcal{T}(\mu)=u^{*}(\mu)$ using implicit differentiation. At convergence, the optimal pair $(u^{*},\lambda^{*})$ satisfies [Equation 19](https://arxiv.org/html/2506.07003v2#A3.E19). Differentiating both sides of the first stationarity equation in [Equation 19](https://arxiv.org/html/2506.07003v2#A3.E19) w.r.t. $\mu$ gives:

$$ $\displaystyle\frac{\partial}{\partial\mu}R_{0}(u^{*},\lambda^{*};\mu)$ $\displaystyle=\frac{\partial}{\partial\mu}(Q(u^{*}-\mu)+\nabla h(u^{*})\lambda^{*})=0,$ (40) $\displaystyle\iff\big(Q+\nabla^{2}h(u^{*})\lambda^{*}\big)\frac{\partial u^{*}}{\partial\mu}+\nabla h(u^{*})\frac{\partial\lambda^{*}}{\partial\mu}=Q.$ $$

Similar to [Equation 31](https://arxiv.org/html/2506.07003v2#A3.E31), we assume $h(u^{*})$ is approximately affine near the optimal point, and we approximate $\nabla^{2}h(u^{*})\approx 0$.

Similarly differentiating both sides of the second feasibility equation in [Equation 19](https://arxiv.org/html/2506.07003v2#A3.E19) w.r.t $\mu$ gives:

$$ $\displaystyle\frac{\partial}{\partial\mu}R_{1}(u^{*},\lambda^{*};\mu)=\frac{\partial}{\partial\mu}h(u^{*})=\nabla h(u^{*})^{\top}\frac{\partial u^{*}}{\partial\mu}=0.$ (41) $$

Combining [Equation 40](https://arxiv.org/html/2506.07003v2#A3.E40) and [Equation 41](https://arxiv.org/html/2506.07003v2#A3.E41) leads to the following block linear system:

$$ $\begin{bmatrix}Q&J^{*\top}\\ J^{*}&0\end{bmatrix}\begin{bmatrix}\partial u^{*}/\partial\mu\\ \partial\lambda^{*}/\partial\mu\end{bmatrix}=\begin{bmatrix}Q\\ 0\end{bmatrix}.$ $$

Similar to [Equation 38](https://arxiv.org/html/2506.07003v2#A3.E38), we can use the Schur complement to eliminate the dual term via block substitution. Using the block inverse in [Equation 39](https://arxiv.org/html/2506.07003v2#A3.E39) with $J^{(i)}=J^{*}$, we have

$$ $\displaystyle\begin{bmatrix}\partial u^{*}/\partial\mu\\ \partial\lambda^{*}/\partial\mu\end{bmatrix}$ $\displaystyle=\begin{bmatrix}I&-Q^{-1}J^{*\top}\\ 0&I\end{bmatrix}\begin{bmatrix}Q^{-1}&0\\ 0&-(J^{*}Q^{-1}J^{*\top})^{-1}\end{bmatrix}\begin{bmatrix}I&0\\ -J^{*}Q^{-1}&I\end{bmatrix}\begin{bmatrix}Q\\ 0\end{bmatrix}$ $\displaystyle=\begin{bmatrix}I&-Q^{-1}J^{*\top}\\ 0&I\end{bmatrix}\begin{bmatrix}Q^{-1}&0\\ 0&-(J^{*}Q^{-1}J^{*\top})^{-1}\end{bmatrix}\begin{bmatrix}Q\\ -J^{(*)}\end{bmatrix}$ $\displaystyle=\begin{bmatrix}I&-Q^{-1}J^{*\top}\\ 0&I\end{bmatrix}\begin{bmatrix}I\\ (J^{*}Q^{-1}J^{*\top})^{-1}J^{*}\end{bmatrix}$ $\displaystyle=\begin{bmatrix}I-Q^{-1}J^{*\top}\partial\lambda^{*}/\partial\mu\\ (J^{*}Q^{-1}J^{*\top})^{-1}J^{*}\end{bmatrix}.$ $$

Hence, the Jacobian is given by the first component as:

$$ $J_{\mathcal{T}}(\mu)=\frac{\partial u^{*}(\mu)}{\partial\mu}=I-Q^{-1}J^{*\top}\left(J^{*}Q^{-1}J^{*\top}\right)^{-1}J^{*},$ $$

which is the desired [Equation 24](https://arxiv.org/html/2506.07003v2#A3.E24).
∎

###### Proposition C.2 .

###### Proof.

###### Proposition C.3 .

###### Proof.

### C.3 (Nonlinear) Convex Inequality Constraints

In this subsection, we describe how to compute the DPPL in ProbHardE2E for nonlinear convex inequality constraints. Convex inequality constraints arise naturally in many scientific and engineering applications. For example, total variation (TV) regularization is widely used to promote smoothness or piecewise-constant structure in spatial fields, e.g., image denoising [Rudin et al., 1992, Boyd and 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 of physical or operational quantities [Bertsekas, 1997].

We consider the constrained projection Problem [7](https://arxiv.org/html/2506.07003v2#S3.E7),
where $Q\succ 0$, and $h:\mathbb{R}^{n}\to\mathbb{R}^{q}$, $g:\mathbb{R}^{n}\to\mathbb{R}^{s}$ 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 is

$$ $L(u,\lambda,\nu;z)=\tfrac{1}{2}(u-z)^{\top}Q(u-z)+\lambda^{\top}h(u)+\nu^{\top}g(u),$ $$

with Lagrange multipliers $\lambda\in\mathbb{R}^{q}$ for the equality constraints and $\nu\in\mathbb{R}^{s}$ for the inequality constraints, where $\nu\geq 0$. The KKT optimality conditions are given as:

$$ $\displaystyle\begin{split}\text{(Stationarity)}\quad&Q(u^{*}-z)+\nabla h(u^{*})\lambda^{*}+\nabla g(u^{*})\nu^{*}=0,\\ \text{(Primal feasibility)}\quad&h(u^{*})=0,\quad g(u^{*})\leq 0,\\ \text{(Dual feasibility)}\quad&\nu^{*}\geq 0,\\ \text{(Complementary slackness)}\quad&\nu_{j}^{*}\cdot g_{j}(u^{*})=0\quad\text{for all }j=1,\dots,s.\end{split}$ (42) $$

Note the first two conditions are the same as the ones for nonlinear equality constraints with $\nu=0$, in [Equation 19](https://arxiv.org/html/2506.07003v2#A3.E19).

The KKT conditions in [Equation 42](https://arxiv.org/html/2506.07003v2#A3.E42) are necessary and sufficient for optimality, under standard constraint qualifications, e.g., Slater’s condition [Boyd and Vandenberghe, 2004]. [Equation 42](https://arxiv.org/html/2506.07003v2#A3.E42) can be solved by various optimization methods, e.g., stochastic trust-region methods with sequential quadratic programming (SQP) [Boyd and Vandenberghe, 2004, Hong et al., 2023] and exact augmented Lagrangian [Boyd and Vandenberghe, 2004, Fang et 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)\leq 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 and Vandenberghe, 2004]. Nonlinear equalities 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 and Vandenberghe, 2004].

To compute the Jacobian $J_{\mathcal{T}}(\mu):=\partial u^{*}(\mu)/\partial\mu$ of the projection map with respect to the input $\mu$, we could, in principle, apply implicit differentiation to the KKT conditions in [Equation 42](https://arxiv.org/html/2506.07003v2#A3.E42). For general constrained problems with nonlinear equality and convex inequality constraints, the derivation becomes analytically complex, particularly due to active set variability and non-affine structure. In the special case of quadratic programs with affine constraints, OptNet [Amos and Kolter, 2017] provides an explicit expression for the derivatives via KKT conditions. In addition, CVXPYLayers [Agrawal et 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 not 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 to the inputs and computing empirical statistics over repeated forward passes.

## Appendix D Special Cases of ProbHardE2E

In this section, we show applications of ProbHardE2E in two seemingly unrelated but technically related 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=I$) and oblique ($Q=\Sigma^{-1}$) projections, respectively.
[Figure 3](https://arxiv.org/html/2506.07003v2#S4.F3) illustrates the wide variety of cases that our framework covers.

Figure: Figure 3: ProbHardE2E serves as a probabilistic unified framework for learning with hard constraints.
Refer to caption: x1.png

### D.1 Enforcing Coherency in Hierarchical Time Series Forecasting

Hierarchical time series forecasting is abundant in several applications, e.g., retail demand forecasting and electricity forecasting. In 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 at each time point through known relationships, which can be represented as dependency graphs. Let $z_{t}=[a_{t}\quad b_{t}]^{\top}\in\mathbb{R}^{n}$, where $a_{t}\in\mathbb{R}^{q}$ denotes the aggregate entries, $b_{t}\in\mathbb{R}^{\tilde{q}}$ denotes the bottom-level entries, and $n=q+\tilde{q}$. Let $S_{\text{sum}}\in\{0,1\}^{q\times\tilde{q}}$ denote the summation matrix, which defines the relationship between the bottom and aggregate levels as $a_{t}=S_{\text{sum}}b_{t}$. This coherency constraint can be equivalently expressed as:

$$ $\begin{bmatrix}I_{q}&-S_{\text{sum}}\end{bmatrix}\begin{bmatrix}a_{t}\\ b_{t}\end{bmatrix}=0\Leftrightarrow Az_{t}=0,\quad\forall t,$ (43) $$

where $I_{q}$ denotes the $q\times q$ identity matrix. See Hyndman et al. [2011], Rangapuram et al. [2021], Olivares et al. [2024b] and the references therein for details, and [Figure 4](https://arxiv.org/html/2506.07003v2#A4.F4) for an illustration.

Figure: Figure 4: Example hierarchical time series structure with $a_{t}\in\mathbb{R}^{4}$, $b_{t}\in\mathbb{R}^{6}$ and $S_{\text{sum}}=\begin{bmatrix}1&1&1&1&1&1\\ 1&0&0&1&1&0\\ 0&0&0&1&1&1\\ 0&1&1&0&0&0\end{bmatrix}$.

HierE2E [Rangapuram et al., 2021] enforces the coherency constraint in [Equation 43](https://arxiv.org/html/2506.07003v2#A4.E43) by projecting the multivariate samples $z_{t}$ onto the null space of the constraint, i.e., $Az_{t}=0$. It uses the following projection:

$$ $\displaystyle\begin{split}&u^{*}(z_{t})=(I-A^{\top}(AA^{\top})^{-1}A)z_{t}=(I-A^{{\dagger}}A)z_{t},\end{split}$ (44) $$

where $A^{{\dagger}}=A^{\top}(AA^{\top})^{-1}$ denotes the right psuedoinverse, and $P_{I}=P_{I}^{2}=P_{I}^{\top}$ denotes an orthogonal projector.

We show that HierE2E can be formulated in our ProbHardE2E framework with the following posterior mean and covariance:

$$ $\displaystyle\hat{\mu}_{\text{HierE2E}}$ $\displaystyle=(I-A^{{\dagger}}A)\mu,$ (45a) $\displaystyle\hat{\Sigma}_{\text{HierE2E}}$ $\displaystyle=\Sigma-A^{{\dagger}}A\Sigma-\Sigma A^{{\dagger}}A+A^{{\dagger}}A\Sigma A^{{\dagger}}A,$ (45b) $$

where $P_{I}$ is defined in [Equation 44](https://arxiv.org/html/2506.07003v2#A4.E44).
In particular, we show in [Proposition D.1](https://arxiv.org/html/2506.07003v2#A4.Thmproposition1) that the HierE2E posterior update in [Equation 45](https://arxiv.org/html/2506.07003v2#A4.E45)
is a 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 [Equation 45](https://arxiv.org/html/2506.07003v2#A4.E45) is given by the solution to Problem [7](https://arxiv.org/html/2506.07003v2#S3.E7) with linear constraints in [Proposition C.1](https://arxiv.org/html/2506.07003v2#A3.Thmproposition1), i.e., $h(u)=Au=0$, $b=0$, where $Q=I$ for an orthogonal projection and $\mathbf{Z}\sim\mathcal{N}(\mu,\Sigma)$ is a multivariate Gaussian.

###### Proof.

The oblique projection in [Equation 16](https://arxiv.org/html/2506.07003v2#A3.E16) used in ProbHardE2E for linear constraints is given as $P_{Q^{-1}}=I-Q^{-1}A^{\top}(AQ^{-1}A^{\top})^{-1}A$. Setting $Q=I$, the expression simplifies to $P_{Q^{-1}}=P_{I}=I-A^{\top}(AA^{\top})^{-1}A=I-A^{{\dagger}}A$.

The posterior mean $\hat{\mu}$ for ProbHardE2E with linear constraints is given in [Equation 18a](https://arxiv.org/html/2506.07003v2#A3.E18.1) with $b=0$ as:

$$ $\displaystyle\hat{\mu}$ $\displaystyle=P_{Q^{-1}}\mu$ (46) $\displaystyle=P_{I}\mu$ $\displaystyle=(I-A^{{\dagger}}A)\mu$ $\displaystyle=\hat{\mu}_{\text{Hier-E2E}},$ $$

which is the desired expression in [Equation 45a](https://arxiv.org/html/2506.07003v2#A4.E45.1).

Similarly, the posterior covariance $\hat{\Sigma}$ for ProbHardE2E in [Equation 18b](https://arxiv.org/html/2506.07003v2#A3.E18.2) is given as:

$$ $\displaystyle\hat{\Sigma}$ $\displaystyle=P_{Q^{-1}}\Sigma P^{\top}_{Q^{-1}}$ (47) $\displaystyle=P_{I}\Sigma P_{I}^{\top}$ $\displaystyle=P_{I}\Sigma P_{I}$ $\displaystyle=(I-A^{{\dagger}}A)\Sigma(I-A^{{\dagger}}A)$ $\displaystyle=(I-A^{{\dagger}}A)(\Sigma-\Sigma A^{{\dagger}}A)$ $\displaystyle=\Sigma-A^{{\dagger}}A\Sigma-\Sigma A^{{\dagger}}A+A^{{\dagger}}A\Sigma A^{{\dagger}}A$ $\displaystyle=\hat{\Sigma}_{\text{Hier-E2E}},$ $$

which is the desired expression in [Equation 45b](https://arxiv.org/html/2506.07003v2#A4.E45.2).
∎

Note that HierE2E does not directly project the distribution parameters, even though a closed form exists, as shown in [Equation 46](https://arxiv.org/html/2506.07003v2#A4.E46) and [Equation 47](https://arxiv.org/html/2506.07003v2#A4.E47). Instead, it directly projects the samples in [Equation 44](https://arxiv.org/html/2506.07003v2#A4.E44). An improvement to HierE2E (that we do in ProbHardE2E) is to eliminate the computationally expensive sampling in the training loop. (See [Section 3.5](https://arxiv.org/html/2506.07003v2#S3.SS5).) 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 in their framework. For many standard distributions, e.g., the multivariate Gaussian distribution in HierE2E, the CRPS can be computed analytically [Matheson and Winkler, 1976, Taillardat et al., 2016] using the mean and covariance of the output distribution.

###### Proposition D.1 .

###### Proof.

### D.2 Enforcing Conservation Laws in PDEs

In addition to hierarchical forecasting, another (at first seemingly-unrelated) application of ProbHardE2E is enforcing conservation laws in solutions to partial differential equations (PDEs). A conservation law is given as $u_{t}+\nabla\cdot F(u)=0$, for unknown $u(t,x)$ and nonlinear flux function $F(u)$ [LeVeque, 1990]. Hansen et al. [2023] propose the ProbConserv method to enforce the integral form of conservation laws from finite volume methods [LeVeque, 2002] as a linear constraint $Au=b$ 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:

$$ $\displaystyle\hat{\mu}_{\text{ProbConserv}}$ $\displaystyle=\mu-\Sigma A^{\top}(A\Sigma A^{\top})^{-1}(A\mu-b),$ (48a) $\displaystyle\hat{\Sigma}_{\text{ProbConserv}}$ $\displaystyle=\Sigma-\Sigma A^{\top}(A\Sigma A^{\top})^{-1}A\Sigma,$ (48b) $$

given the mean $\mu$ and the covariance matrix $\Sigma$ estimated from a black-box probabilistic model, e.g., Gaussian Process, probabilistic Neural Operators [Mouli et al., 2024] or Attentive Neural Process (ANP) [Hansen et al., 2023] or DeepVAR [Salinas et al., 2019] used in the hierarchical forecasting case.

In ProbConserv, the posterior mean $\hat{\mu}$ in [Equation 48a](https://arxiv.org/html/2506.07003v2#A4.E48.1) is shown to be the solution to the constrained least squares problem:

$$ $\displaystyle\hat{\mu}_{\text{ProbConserv}}=\operatorname*{argmin}_{\begin{subarray}{c}\tilde{\mu}\in\mathbb{R}^{n}\\ A\tilde{\mu}=b\end{subarray}}\frac{1}{2}||\tilde{\mu}-\mu||^{2}_{\Sigma^{-1}}.$ $$

We formulate this optimization problem more generally, and show that by assuming that $z\sim\mathbf{Z}\sim\mathcal{N}(\mu,\Sigma)$ is a multivariate Gaussian, a constrained sample $u^{*}(z)\sim\mathbf{Y}\sim\mathcal{N}(\hat{\mu},\hat{\Sigma})$ in ProbConserv is a solution to our Problem [7](https://arxiv.org/html/2506.07003v2#S3.E7) with $Q=\Sigma^{-1}$ and linear constraints. In particular, we show in [Proposition D.2](https://arxiv.org/html/2506.07003v2#A4.Thmproposition2) that the ProbConserv posterior update in [Equation 48](https://arxiv.org/html/2506.07003v2#A4.E48) is a special linear constraint case of our ProbHardE2E method, which uses an oblique projection with $Q=\Sigma^{-1}$.

###### Proposition D.2 .

The projected mean and covariance for ProbConserv in [Equation 48](https://arxiv.org/html/2506.07003v2#A4.E48) is given by the solution to Problem [7](https://arxiv.org/html/2506.07003v2#S3.E7) with linear constraints in [Proposition C.1](https://arxiv.org/html/2506.07003v2#A3.Thmproposition1), i.e., $h(u)=Au-b=0$, where $Q=\Sigma^{-1}$ for an oblique projection and $\mathbf{Z}\sim\mathcal{N}(\mu,\Sigma)$ is a multivariate Gaussian.

###### Proof.

The oblique projection in [Equation 16](https://arxiv.org/html/2506.07003v2#A3.E16) used in ProbHardE2E for linear constraints is given as $P_{Q^{-1}}=P_{\Sigma}=I-Q^{-1}A^{\top}(AQ^{-1}A^{\top})^{-1}A$. Setting $Q=\Sigma^{-1}$, we have that $P_{Q^{-1}}=I-\Sigma A^{\top}(A\Sigma A^{\top})^{-1}A=P_{\Sigma}$.

The posterior mean $\hat{\mu}$ for ProbHardE2E with linear constraints is given in [Equation 18a](https://arxiv.org/html/2506.07003v2#A3.E18.1) as:

$$ $\displaystyle\hat{\mu}=P_{Q^{-1}}\mu+(I-P_{Q^{-1}})A^{{\dagger}}b$ $\displaystyle=(I-\Sigma A^{\top}(A\Sigma A^{\top})^{-1}A)\mu+(I-(I-\Sigma A^{\top}(A\Sigma A^{\top})^{-1}A))A^{{\dagger}}b$ $\displaystyle=(I-\Sigma A^{\top}(A\Sigma A^{\top})^{-1}A)\mu+\Sigma A^{\top}(A\Sigma A^{\top})^{-1}\underbrace{AA^{{\dagger}}}_{I}b$ $\displaystyle=\mu-\Sigma A^{\top}(A\Sigma A^{\top})^{-1}(A\mu-b)$ $\displaystyle=\hat{\mu}_{\text{ProbConserv}},$ $$

which is equal to the desired expression in [Equation 48a](https://arxiv.org/html/2506.07003v2#A4.E48.1).

Similarly, the posterior covariance $\hat{\Sigma}$ for ProbHardE2E in [Equation 18b](https://arxiv.org/html/2506.07003v2#A3.E18.2) is given as:

$$ $\displaystyle\hat{\Sigma}$ $\displaystyle=P_{Q^{-1}}\Sigma P^{\top}_{Q^{-1}}$ $\displaystyle=(I-\Sigma A^{\top}(A\Sigma A^{\top})^{-1}A)\Sigma(I-A^{\top}(A\Sigma A^{\top})^{-1}A\Sigma)$ $\displaystyle=(I-\Sigma A^{\top}(A\Sigma A^{\top})^{-1}A)(\Sigma-\Sigma A^{\top}(A\Sigma A^{\top})^{-1}A\Sigma)$ $\displaystyle=\Sigma-2\Sigma A^{\top}(A\Sigma A^{\top})^{-1}A\Sigma+\Sigma A^{\top}(A\Sigma A^{\top})^{-1}(A\Sigma A^{\top})(A\Sigma A^{\top})^{-1}A\Sigma$ $\displaystyle=\Sigma-\Sigma A^{\top}(A\Sigma A^{\top})^{-1}A\Sigma$ $\displaystyle=\hat{\Sigma}_{\text{ProbConserv}},$ $$

which is equal to the desired expression in [Equation 48b](https://arxiv.org/html/2506.07003v2#A4.E48.2).
∎

Note that the projected distribution parameters in [Equation 48](https://arxiv.org/html/2506.07003v2#A4.E48) 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.

###### Proposition D.2 .

###### Proof.

## Appendix E Flexibility in the choice of Q 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 in the constrained least squares Problem [7](https://arxiv.org/html/2506.07003v2#S3.E7). 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 [6](https://arxiv.org/html/2506.07003v2#A5.T6) summarizes common structure choices for $Q$ and their trade-offs.
Of course, in many applications, there is a single goal for the choice of Q—to optimize accuracy.

**Table 6: Several structure choices for the matrix $Q$ and their associated trade-offs.**
| Structure of $Q$ | Example Form | Merits and Demerits |
| --- | --- | --- |
| Identity | $Q=I$ | + Simplest choice, no parameters<br>+ Strong regularization<br>– Ignores uncertainty and correlations |
| Diagonal (learned) | $Q=\text{diag}(q_{1},\dots,q_{n})$ | + Captures heteroscedasticity<br>+ Efficient to compute and invert<br>– Ignores correlations |
| Low-rank (learned $L$) | $Q=LL^{\top},\,L\in\mathbb{R}^{n\times d}$ | + Captures dominant correlations<br>+ Fewer parameters than full<br>– Still computationally involved |
| Full (learned $L$) | $Q=LL^{\top},\,L\in\mathbb{R}^{n\times n}$ | + Fully expressive<br>– High memory and compute cost<br>– Prone to overfitting |

In practice, the space of symmetric positive definite (SPD) matrices is too large to be explored (and “learned”) without additional structure, especially in high-dimensional settings.
To address this, structural constraints are often imposed on $Q$, reducing the number of parameters, and acting as a form of regularization [Willette et 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 or hope about) the structure of the underlying data distribution.
The simplest choice, $Q=I$, 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 heteroscedasticity. A diagonal matrix $Q=\text{diag}(q_{1},\dots,q_{n})$ introduces per-dimension weighting, and is well-suited to heteroscedastic tasks where the variance differs across outputs [Kendall and 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 to avoid overfitting [Weinberger and Saul, 2009].

We focus on two concrete realizations of $Q$: the identity matrix $Q=I$ that is used in the HierE2E [Rangapuram et al., 2021] (see [Section D.1](https://arxiv.org/html/2506.07003v2#A4.SS1)), and a diagonal matrix defined as the inverse of a predicted diagonal covariance, $Q={\Sigma}^{-1}$ that is used in ProbConserv [Hansen et al., 2023] (see [Section D.2](https://arxiv.org/html/2506.07003v2#A4.SS2)), where ${\Sigma}=\text{diag}(\sigma_{1}^{2},\dots,\sigma_{d}^{2})$ denote the empirical variances output by the model. This latter choice corresponds to a heteroscedastic formulation that scales residuals based on their predicted precision, which emphasizes more confident predictions, and down-weights less certain ones [Stirn et al., 2023, Le et al., 2005, Hansen et al., 2023].

## Appendix F Proof of Theorem 3.1

In this section, we begin by first restating Theorem [3.1](https://arxiv.org/html/2506.07003v2#S3.Thmtheorem1), which provides a closed-form update for our DPPL in [Equation 8](https://arxiv.org/html/2506.07003v2#S3.E8) for a prior distribution that belongs to a multivariate local-scale family of distributions; and then we provide its proof.

See [3.1](https://arxiv.org/html/2506.07003v2#S3.Thmtheorem1)

###### Proof.

Recall that a family of probability distributions is said to be a location-scale family if for any random variable $\mathbf{Z}$ whose distribution belongs to the family $\mathbf{Z}\sim\mathcal{F}(\mu,\Sigma)$, then there exists a transformation (re-parameterization) of the form

$$ $\mathbf{Y}\stackrel{{\scriptstyle d}}{{=}}A\mathbf{Z}+B,$ $$

where $A$ denotes a scale transformation matrix, $B$ denotes the location parameter, and $\stackrel{{\scriptstyle d}}{{=}}$ denotes equality in distribution.

Let $\mathbf{Y}=\mathcal{T}(\mathbf{Z})$ be a nonlinear transformation. We calculate the first-order Taylor series expansion to linearize the function about the mean $\mu$ as:

$$ $\displaystyle\mathbf{Y}=\mathcal{T}(\mathbf{Z})$ $\displaystyle\approx\mathcal{T}(\mu)+J_{\mathcal{T}}(\mu)(\mathbf{Z}-\mu)$ (49) $\displaystyle=\underbrace{J_{\mathcal{T}}}_{A}(\mu)\mathbf{Z}+\underbrace{(\mathcal{T}(\mu)-J_{\mathcal{T}}(\mu)\mu)}_{B}.$ $$

Then, since $\mathbf{Z}$ belongs to the location-scale family of distributions, the linearization of $\mathbf{Y}\sim\mathcal{F}(\hat{\mu},\hat{\Sigma})$ also belongs to the family with mean $\hat{\mu}$ and covariance $\hat{\Sigma}$, which we compute below.

Taking the expectation of both sides of [Equation 49](https://arxiv.org/html/2506.07003v2#A6.E49) we get:

$$ $\displaystyle\hat{\mu}=\mathbb{E}[\mathcal{T}(\mathbf{Z})]$ $\displaystyle\approx\mathbb{E}[\mathcal{T}(\mu)+J_{\mathcal{T}}(\mu)(\mathbf{Z}-\mu)]$ $\displaystyle=\mathbb{E}[\mathcal{T}(\mu)]+\mathbb{E}[J_{\mathcal{T}}(\mu)(\mathbf{Z}-\mu)]\text{ (by linearity of expectation)}$ $\displaystyle=\mathcal{T}(\mu)+J_{\mathcal{T}}(\mu)\underbrace{(\mathbb{E}[\mathbf{Z}]-\mu)}_{0}\text{ (since $\mu$ is not a random variable)}$ $\displaystyle=\mathcal{T}(\mu).$ (50) $$

Then, the covariance $\hat{\Sigma}$ is given as:

$$ $\displaystyle\hat{\Sigma}$ $\displaystyle=\mathbb{E}[(\mathcal{T}(\mathbf{Z})-\mathbb{E}[\mathcal{T}(\mathbf{Z})])(\mathcal{T}(\mathbf{Z})-\mathbb{E}[\mathcal{T}(\mathbf{Z})])^{\top}]$ $\displaystyle=\mathbb{E}[(\mathcal{T}(\mathbf{Z})-\mathcal{T}(\mu))(\mathcal{T}(\mathbf{Z})-\mathcal{T}(\mu))^{\top}]\text{ (by \lx@cref{creftype~refnum}{eqn: delta_mean})}$ $\displaystyle\approx\mathbb{E}[(\mathcal{T}(\mu)+J_{\mathcal{T}}(\mu)(\mathbf{Z}-\mu)-\mathcal{T}(\mu))(\mathcal{T}(\mu)+J_{\mathcal{T}}(\mu)(\mathbf{Z}-\mu)-\mathcal{T}(\mu))^{\top}]\text{ (by \lx@cref{creftype~refnum}{eqn:taylor_theorem})}$ $\displaystyle=\mathbb{E}[(J_{\mathcal{T}}(\mu)(\mathbf{Z}-\mu))(J_{\mathcal{T}}(\mu)(\mathbf{Z}-\mu))^{\top}]$ $\displaystyle=J_{\mathcal{T}}(\mu)\mathbb{E}[(\mathbf{Z}-\mu)(\mathbf{Z}-\mu)^{\top}]J_{\mathcal{T}}(\mu)^{\top}$ $\displaystyle=J_{\mathcal{T}}(\mu)\Sigma J_{\mathcal{T}}(\mu)^{T}.$ $$

∎

Importantly, the approximation error between the nonlinear transformation and its linearization converges to zero in probability [Van der Vaart, 2000], which ensures the validity of this approach asymptotically. We note that this result is closely related to the Multivariate Delta Method [Casella and Berger, 2001], which shows that for a nonlinear function $\mathcal{T}$, the sample mean of $\mathcal{T}(z_{1},\dots,z_{n})$
also converges in distribution, under mild conditions.
Specifically, if the sample mean of $n$ i.i.d. draws from $\mathbf{Z}$ converges to a multivariate Gaussian (by the CLT), then the same linearization argument and Slutsky’s theorem imply that the sample mean of the projected samples converges to a multivariate Gaussian, with parameters given in [Equation 8](https://arxiv.org/html/2506.07003v2#S3.E8).
Second-order approximations (via a quadratic expansion of $\mathcal{T}$) yield higher-order corrections, and can lead to non-Gaussian outcomes (e.g., chi-squared) [Casella and Berger, 2001].

###### Proof.

## Appendix 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 a 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_{t}+\nabla\cdot F(u)=0,$ (51) $$

for some nonlinear flux function $F(u)$. These equations can be written in their conservative form as:

$$ $\frac{d}{dt}\int_{\Omega}u(t,x)d\Omega=F(u(t,x_{0}))-F(u(t,x_{N})),$ (52) $$

by applying the divergence term in 1D over the domain $\Omega=[x_{0},x_{N}]$ [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 [Equation 52](https://arxiv.org/html/2506.07003v2#A7.E52) 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.e., so that the flux into a 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:

$$ $u_{t}-\nabla\cdot(k(u)\nabla u)=0,$ (53) $$

where the flux in [Equation 51](https://arxiv.org/html/2506.07003v2#A7.E51) is given as $F(u)=-k(u)\nabla u$, 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 [Equation 53](https://arxiv.org/html/2506.07003v2#A7.E53). We use the heat equation with the following sinusoidal initial condition and periodic boundary conditions from Krishnapriyan et al. [2021], Hansen et al. [2023]:

$$ $\displaystyle u_{t}$ $\displaystyle=k\Delta u,\hskip 21.33955pt\forall x\in\Omega=[0,2\pi],\ \forall t\in[0,1],$ (54) $\displaystyle u(0,x)$ $\displaystyle=\sin(x),\hskip 15.6491pt\forall x\in[0,2\pi],$ $\displaystyle u(t,0)$ $\displaystyle=u(2\pi,t),\quad\forall t\in[0,1],$ $$

respectively.
The exact solution, which can be solved using the Fourier Transform, is given as:

$$ $u_{\text{exact}}(t,x)=e^{-kt}\sin(x).$ $$

The solution is a smooth sinusoidal curve that exponentially decays or dissipates over time, and has an infinite speed of propagation. With these specific initial and boundary conditions in [Equation 54](https://arxiv.org/html/2506.07003v2#A7.E54), the global conservation law in [Equation 52](https://arxiv.org/html/2506.07003v2#A7.E52) reduces to the following linear equation:

$$ $\int_{0}^{2\pi}u(t,x)\,dx=0,\quad\forall t\in[0,1],$ (55) $$

since the net flux on the boundaries is 0.

##### Porous Medium Equation (PME) (“Medium”).

The PME is a nonlinear degenerate subclass of the GPME, where the diffusivity is a nonlinear, monomial of the solution, i.e., $k(u)=u^{m}$ in [Equation 53](https://arxiv.org/html/2506.07003v2#A7.E53). It has been using in modeling nonlinear heat transfer [Vázquez, 2007, Maddix et al., 2018a]. 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]:

$$ $\displaystyle u_{t}-\nabla\cdot(u^{m}\nabla u)$ $\displaystyle=0,\hskip 55.48277pt\forall x\in\Omega=[0,1],\ \forall t\in[0,1],$ (56) $\displaystyle u(0,x)$ $\displaystyle=0,\hskip 55.48277pt\forall x\in[0,1],$ $\displaystyle u(t,0)$ $\displaystyle=(mt)^{1/m},\hskip 25.6073pt\forall t\in[0,1].$ $$

The exact solution is given as:

$$ $u_{\text{exact}}(t,x)=\left(m\,\mathrm{ReLU}(t-x)\right)^{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 $t=x$ with a finite speed of propagation. With these specific initial and boundary conditions in [Equation 56](https://arxiv.org/html/2506.07003v2#A7.E56), the global conservation law in [Equation 52](https://arxiv.org/html/2506.07003v2#A7.E52) reduces to the following linear equation:

$$ $\int_{0}^{1}u(t,x)\,dx=\frac{(mt)^{1+1/m}}{m+1},\quad\forall t\in[0,1].$ (57) $$

##### Stefan Equation (“Hard”).

The Stefan equation has been used in foam modeling [van der Meer et al., 2016] and crystallization [Sethian and Strain, 1992], and models phase transitions with the following discontinuous diffusivity:

$$ $\displaystyle k(u)$ $\displaystyle=\begin{cases}1,&u\geq u^{*}\\ 0,&u<u^{*}\end{cases},\quad u^{*}\geq 0,$ $$

in [Equation 53](https://arxiv.org/html/2506.07003v2#A7.E53).
We use the Stefan equation with the following initial condition and Dirichlet boundary conditions from Maddix et al. [2018b], Hansen et al. [2023]:

$$ $\displaystyle u_{t}-\nabla\cdot(k(u)\nabla u)$ $\displaystyle=0,\quad\forall x\in\Omega=[0,1],\ t\in[0,1],$ (58) $\displaystyle u(0,x)$ $\displaystyle=0,\quad\forall x\in[0,1],$ $\displaystyle u(t,0)$ $\displaystyle=1,\quad\forall t\in[0,1].$ $$

The exact solution is given as:

$$ $u_{\text{exact}}(t,x)=\mathbf{1}_{u\geq u^{*}}\left[1-\frac{1-u^{*}}{\mathrm{erf}(\alpha/2)}\mathrm{erf}\left(\frac{x}{2\sqrt{t}}\right)\right],$ $$

where $\mathbf{1}$ denotes the indicator function, $\mathrm{erf}(z)=(2/\sqrt{\pi})\int_{0}^{z}\exp(-y^{2})dy$ denotes the error function, and $\alpha=2\tilde{\alpha}$ and $\tilde{\alpha}$ satisfies the following nonlinear equation:

$$ $\frac{1-u^{*}}{\sqrt{\pi}}=u^{*}\,\mathrm{erf}(\tilde{\alpha})\tilde{\alpha}e^{\tilde{\alpha}^{2}}.$ $$

The solution is a rightward moving shock. With these specific initial and boundary conditions in [Equation 58](https://arxiv.org/html/2506.07003v2#A7.E58), the global conservation law in [Equation 52](https://arxiv.org/html/2506.07003v2#A7.E52) reduces to the following linear equation:

$$ $\int_{0}^{1}u(t,x)\,dx=\frac{2(1-u^{*})}{\mathrm{erf}(\alpha/2)}\sqrt{\frac{t}{\pi}},\quad\forall t\in[0,1].$ (59) $$

#### G.1.2 Hyperbolic Linear Advection Equation

The hyperbolic linear advection equation models fluids transported at a constant velocity, and is given by [Equation 51](https://arxiv.org/html/2506.07003v2#A7.E51) with linear flux $F(u)=\beta u$. We use the
1D linear advection problem with the following step-function initial condition and inflow Dirichlet boundary conditions from Hansen et al. [2023]:

$$ $\displaystyle u_{t}+\beta u_{x}$ $\displaystyle=0,\hskip 28.45274pt\forall x\in\Omega=[0,1],\forall t\in[0,1],$ (60) $\displaystyle u(0,x)$ $\displaystyle=\mathbf{1}_{x\leq 0.5},\hskip 5.69046pt\forall x\in[0,1],$ $\displaystyle u(t,0)$ $\displaystyle=1,\hskip 28.45274pt\forall t\in[0,1].$ $$

The exact solution is given as:

$$ $u(x,t)=h(x-\beta t),$ $$

where $h(x)=\mathbf{1}_{x\leq 0.5}$ denotes the initial condition. The solution remains a shock, which travels to the right with a finite speed of propagation $\beta$. With these specific initial and boundary conditions in [Equation 60](https://arxiv.org/html/2506.07003v2#A7.E60), the global conservation law in [Equation 52](https://arxiv.org/html/2506.07003v2#A7.E52) reduces to the following linear equation:

$$ $\int_{0}^{1}u(x,t)\,dx=\frac{1}{2}+\beta t,$ (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 are coherent with known aggregation constraints across cross-sectional hierarchies [Rangapuram et al., 2021].

[Table 7](https://arxiv.org/html/2506.07003v2#A7.T7) provides an overview of the time series datasets used in our empirical evaluation. For each benchmarking dataset, it details the total number of series, the number of bottom level series (i.e., the leaf nodes in the hierarchy), the number of series aggregated from the bottom-level series, the depth of the hierarchy in terms of the number of levels, the number of time series observations, and the prediction horizon $\tau$.

We adopt the same dataset configurations as in Rangapuram et al. [2021], from which we use the hierarchical forecasting benchmarks and pre-processing pipeline. These datasets are available in GluonTS package [Alexandrov et al., 2019]. The Labour dataset [Australian Bureau of Statistics, 2019] contains monthly Australian employment statistics from 1978 to 2020, organized into a 57-series hierarchy. The Traffic dataset [Ben Taieb and 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, Wiki contains daily page view counts from 199 Wikipedia pages collected over two years [Anava et al., 2018].

**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 $\tau$ | Frequency |
| --- | --- | --- | --- | --- | --- | --- | --- |
| Tourism | 89 | 56 | 33 | 4 | 36 | 8 | Quarterly |
| Tourism-L | 555 | 76; 304 | 175 | 4; 5 | 228 | 12 | Monthly |
| Labour | 57 | 32 | 25 | 4 | 514 | 8 | Monthly |
| Traffic | 207 | 200 | 7 | 4 | 366 | 1 | Daily |
| Wiki | 199 | 150 | 49 | 5 | 366 | 1 | Daily |

## Appendix H Implementation Details

In this section, we provide the implementation details of ProbHardE2E. [Figure 5](https://arxiv.org/html/2506.07003v2#A8.F5) illustrates the overall pipeline of ProbHardE2E, 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 a wide class of models. Conceptually, this layer parallels classical predictor-corrector and primal-dual methods from numerical optimization [Boyd and Vandenberghe, 2004, Bertsekas, 1997], where a candidate solution is refined to satisfy known constraints before evaluation.

We evaluate ProbHardE2E on two scientific domains: (1) PDEs, where structured physical constraints, e.g., conservation laws and boundary conditions, must be enforced (see [Section G.1](https://arxiv.org/html/2506.07003v2#A7.SS1)), and (2) probabilistic hierarchical time series forecasting, where aggregation coherency is required (see [Section G.2](https://arxiv.org/html/2506.07003v2#A7.SS2)). We show that ProbHardE2E is model-agnostic by using a base probabilistic model (predictor) from each application domain, i.e., VarianceNO [Mouli et 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.

Figure: Figure 5: Schematic representation of ProbHardE2E (see [Algorithm 1](https://arxiv.org/html/2506.07003v2#alg1)). Here, a known pathwise-differentiable probabilistic model is chosen to predict a (unconstrained) prior distribution. (Optionally, the projection matrix can 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 or via the Delta Method (see [Section 3.3](https://arxiv.org/html/2506.07003v2#S3.SS3)), which enforces the constraints. Lastly, we choose an appropriate loss function, e.g., CRPS, to calibrate the transformed distribution with the target variable.
Refer to caption: x2.png

### H.1 PDEs

All the experiments are performed on a single NVIDIA V100 GPU. We use a probabilistic Fourier Neural Operator (FNO) [Li et al., 2021], i.e., VarianceNO [Mouli et al., 2024] to learn a mapping from PDE parameters to solutions, e.g., the diffusivity mapping $k(u)\mapsto u(t,x)$ in the (degenerate) parabolic Generalized Porous Medium Equation (GPME), or the velocity mapping $\beta\mapsto u(t,x)$ in the hyperbolic linear advection equation. (See [Section G.1](https://arxiv.org/html/2506.07003v2#A7.SS1) for details on the datasets.)

#### H.1.1 Dataset Generation

[Table 8](https://arxiv.org/html/2506.07003v2#A8.T8) provides an overview of the PDE data generation.
For each PDE in [Section G.1](https://arxiv.org/html/2506.07003v2#A7.SS1), we generate a dataset of $N=200$ parameter-solution pairs $\{\phi^{(i)},u^{(i)}\}_{i=1}^{N}\sim\mathcal{D}$, where $\phi^{(i)}$ denotes the input PDE parameters, e.g., $k,m,u^{*},\beta$, and $u^{(i)}$ denotes the corresponding spatiotemporal solution field. Each solution $u^{(i)}(t,x)$ is simulated over a grid of 100 equidistant points in both space and time, yielding a total of $100\times 100$ observations per instance. During evaluation, we 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/Test (%) |
| --- | --- | --- | --- | --- |
| Heat | $k\in[1,5]$ | $[0,2\pi]$ | $[0,1]$ | 80/20 |
| PME | $m\in[2,3]$ | $[0,1]$ | $[0,1]$ | 80/20 |
| Stefan | $u^{*}\in[0.6,0.65]$ | $[0,1]$ | $[0,1]$ | 80/20 |
| Linear Advection | $\beta\in[1,2]$ | $[0,1]$ | $[0,1]$ | 80/20 |

#### H.1.2 Architectural Details

We use VarianceNO [Mouli et al., 2024] as our base unconstrained probabilistic model. VarianceNO is an augmented Fourier Neural Operator (FNO) [Li et al., 2021] that updates the last layer to 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](https://arxiv.org/html/2506.07003v2#A8.T9) details the model hyperparameters.

**Table 9: Hyperparameters for the VarianceNO model.**
| Hyperparameter | Values |
| --- | --- |
| VarianceNO |  |
| Number of Fourier layers | 4 |
| Channel width | {32, 64} |
| Number of Fourier modes | 12 |
| Batch size | 20 |
| Learning rate | {10^-4, 10^-3, 10^-2} |

#### H.1.3 Training and Testing Setup

We follow the standard training procedure for FNO-based models as proposed by Li et al. [2021]. Specifically, we use the Adam optimizer [Kingma and 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](https://arxiv.org/html/2506.07003v2#A8.T8) to construct test sets and compute the evaluation metrics.

### H.2 Probabilistic Time Series Forecasting

The experiments are performed on an Intel(R) Xeon(R) CPU E5-2603 v4 @ 1.70GHz.

#### H.2.1 Dataset Generation

We adopt the hierarchical forecasting benchmarks and preprocessing pipeline introduced in Rangapuram et 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](https://arxiv.org/html/2506.07003v2#A7.T7)). 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 et al., 2019] as our base unconstrained probabilistic model, which is aligned with Hier-E2E.  DeepVAR is a probabilistic autoregressive LSTM-based model that leverages a 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 ProbHardE2E model in the time series application is developed based on Hier-E2E in GluonTS [Alexandrov et al., 2019]. We use the default base model architecture DeepVAR, and make further modifications to Hier-E2E. Specifically, we tune the hyperparameters in [Table 10](https://arxiv.org/html/2506.07003v2#A8.T10), and adjust the loss to CRPS for structured probabilistic evaluation. We disable sampling-based projection during training because ProbHardE2E optimizes the closed-from CRPS for Gaussian distributions, and our projection methodology ensures that linear constraints are met probabilistically.
During inference, we report CRPS through samples, in order to 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. of Prediction Samples |
| --- | --- | --- | --- | --- | --- |
| Labour | 5 | 32 | 0.01 | 24 | 400 |
| Traffic | 10 | 32 | 0.001 | 40 | 400 |
| Tourism | 10 | 32 | 0.01 | 24 | 200 |
| Tourism-L | 10 | 4 | 0.001 | 36 | 200 |
| Wiki | 25 | 32 | 0.001 | 15 | 200 |

#### H.2.3 Training and Testing Setup

We follow the standard GluonTS [Alexandrov et al., 2019] training setup using the Adam optimizer [Kingma and Ba, 2015] and mini-batch updates. Each epoch consists of 50 batches, with batch size set according to [Table 10](https://arxiv.org/html/2506.07003v2#A8.T10). We run our evaluation five times and report the mean and variance of the CRPS values in [Table 3](https://arxiv.org/html/2506.07003v2#S4.T3).

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 loss [Gneiting et al., 2005] directly on the predicted mean and variance. This makes the training process sampling-free and reduces training time, similar to the PDE case discussed later in Figure [6](https://arxiv.org/html/2506.07003v2#A9.F6). This key distinction avoids the use of coherent_train_samples, as 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 we apply our DPPL to ensure that they satisfy the hierarchical aggregation constraints. This setup parallels the coherent_pred_samples mode in HierE2E, and we implement the inference step with this approach for experimentation simplicity. [Table 10](https://arxiv.org/html/2506.07003v2#A8.T10) 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.

### H.3 Metrics

We evaluate ProbHardE2E and the various baselines using the following metrics. We denote the exact solution or ground truth observations as $u$, and we report the metrics on the mean $\hat{\mu}$, covariance $\hat{\Sigma}$, and samples $\{u^{*}_{i}\}_{i=1}^{N}$ drawn from the constrained multivariate Gaussian distribution $\mathcal{N}(\hat{\mu},\hat{\Sigma})$.

##### Mean Squared Error (MSE).

The MSE measures the mean prediction accuracy and is given as:

$$ $\text{MSE}(\hat{\mu})=\frac{1}{n}\|u-\hat{\mu}\|^{2}_{F},$ $$

where the Frobenius norm is taken over all the datapoints $n$ in $\hat{\mu}$.

##### 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:

$$ $\text{CE}(u^{*})=\sum_{i=1}^{N}\|h(u^{*}_{i})\|_{2}^{2},$ $$

where we compute the average error over $N=100$ samples $\{u^{*}_{i}\}_{i=1}^{N}$.

##### Continuous Ranked Probability Score (CRPS).

The CRPS [Gneiting and 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 $\mathcal{N}(\mu,\mathrm{diag}(\hat{\sigma}^{2}))$, where ${\hat{\sigma}}_{ii}^{2}$ denotes the $i$-th diagonal entry of the predictive covariance ${\hat{\Sigma}}$, the CRPS is given in closed-form as:

$$ $\text{CRPS}_{\mathcal{N}}({\hat{\mu}},{\hat{\sigma}};u)=\sum_{i=1}^{n}{\hat{\sigma}}_{ii}\left[z_{i}\left(2P(z_{i})-1\right)+2p(z_{i})-\frac{1}{\sqrt{\pi}}\right],$ $$

where $z_{i}=(u_{i}-{\hat{\mu}}_{i})/{\hat{\sigma}}_{ii}$, $p(z_{i})=(1/\sqrt{2\pi})\exp(-z_{i}^{2}/2)$ denotes the standard normal PDF, and $P(z_{i})=\int^{z_{i}}_{-\infty}p(y)dy$ denotes the standard normal CDF [Gneiting et al., 2005, Taillardat et al., 2016].

## Appendix 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.

### I.1 Linear Equality Constraints

In this subsection, we show additional results and details for ProbHardE2E with linear equality constraints in both PDEs and hierarchical time series forecasting.

#### I.1.1 PDEs with Conservation Law Constraints

Here, we impose the discretized form of the simplified linear global conservation laws given in [Section G.1](https://arxiv.org/html/2506.07003v2#A7.SS1) for the heat equation ([Equation 55](https://arxiv.org/html/2506.07003v2#A7.E55)), PME ([Equation 57](https://arxiv.org/html/2506.07003v2#A7.E57)), Stefan ([Equation 59](https://arxiv.org/html/2506.07003v2#A7.E59)) and linear advection equation ([Equation 61](https://arxiv.org/html/2506.07003v2#A7.E61)). We use the trapezoidal discretizations of the integrals from Hansen et al. [2023].

[Figure 6](https://arxiv.org/html/2506.07003v2#A9.F6) shows the analogous training time per epoch to LABEL:fig:pde_timing for PDE datasets. Models trained with 100 posterior samples per training step incur a $3.5\text{--}3.6\times$ increase in epoch time relative to our ProbHardE2E, which avoids sampling altogether by using a closed-form CRPS loss. See [Table 2](https://arxiv.org/html/2506.07003v2#S4.T2) for the accuracy results.

Figure: Figure 6: ProbHardE2E: PDE timing comparisons for our sampling-free approach.
Refer to caption: x3.png

#### I.1.2 Hierarchical Time Series Forecasting with Coherency Constraints

Here, we test ProbHardE2E on probabilistic hierarchical forecasting with coherency constraints. (See [Section G.2](https://arxiv.org/html/2506.07003v2#A7.SS2) for details and [Table 3](https://arxiv.org/html/2506.07003v2#S4.T3) for the results.) We compare the two variants of ProbHardE2E, i.e., with oblique $Q=\Sigma^{-1}$ (ProbHardE2E-Ob) and with orthogonal $Q=I$ (ProbHardE2E-Or) projection to the following baselines:

- •
DeepVAR [Salinas et al., 2019] is the base unconstrained probabilistic model, which assumes a multivariate Gaussian distribution for $\mathbf{Z}\sim\mathcal{N}(\mu,\Sigma)$ with mean $\mu$ and diagonal covariance $\Sigma$.
- •
Hier-E2E [Rangapuram et al., 2021] uses DeepVAR as the base model, and enforces the exact coherency constraint by applying the orthogonal projection directly on 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.
- •
ProbConserv [Hansen et al., 2023] enforces the coherency constraint as an oblique projection at inference time only.
- •
ARIMA-NaiveBU and ETS-NaiveBU 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 and Athanasopoulos, 2018].
- •
PERMBU-MINT
[Taieb et al., 2017] is 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 et al., 2024a] or CLOVER [Olivares et al., 2024b] in our experiments 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 knowledge, we use the same base model to Hier-E2E (i.e., DeepVAR), and we evaluate forecast accuracy compared to Hier-E2E to assess the added value of our ProbHardE2E.

### I.2 Nonlinear Equality Constraints

In this subsection, we impose the discretized form of the general nonlinear global linear conservation laws from [Equation 52](https://arxiv.org/html/2506.07003v2#A7.E52) in [Section G.1](https://arxiv.org/html/2506.07003v2#A7.SS1) for the PME with various ranges for the parameter $m$. (See [Table 4](https://arxiv.org/html/2506.07003v2#S4.T4) for the results and LABEL:fig:nonlin_pme for the solution profile.)
For the PME, the flux in [Equation 52](https://arxiv.org/html/2506.07003v2#A7.E52) is given as $F(u)=-k\nabla u$, where $k(u)=u^{m}$. Substituting this flux into [Equation 52](https://arxiv.org/html/2506.07003v2#A7.E52) and integrating in time gives the general conservation law for the PME as:

$$ $\int_{\Omega}u(t,x)d\Omega=\int_{0}^{t}[u^{m}(t,x_{0})\nabla u(t,x_{0})-u^{m}(t,x_{N})\nabla u(t,x_{N})]dt,\quad\forall t\in[0,1].$ $$

Similar to the linear equality constraint case, we discretize the integral using the trapezoidal rule. Unlike ProbConserv [Hansen et al., 2023], which requires an analytical flux expression to evaluate the right‐hand side, our ProbHardE2E can enforce arbitrary (nonlinear) conservation laws directly. In addition, ProbHardE2E with nonlinear constraints can be applied to arbitrary PDEs with any initial or boundary conditions. We impose the initial and boundary conditions as additional linear constraints and enforce positivity on the solution. We test on various training and testing ranges for the parameter $m$, i.e., $m\in[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, Hansen et al., 2023]. We see in [Table 4](https://arxiv.org/html/2506.07003v2#S4.T4) that across all values of $m$, either our oblique ProbHardE2E-Ob or orthogonal projection ProbHardE2E-Or variants of our method perform better than all the baselines.

### I.3 (Nonlinear) Convex Inequality Constraints

In this subsection, we impose a convex total variation diminishing (TVD) inequality constraint. (See LABEL:fig:lin_adv_tvd 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, LeVeque, 1990, Tezaur et al., 2017].
The total variation (TV) is defined in its continuous form as:

$$ $\text{TV}(u(t,\cdot))=\int_{\Omega}\bigg|\frac{\partial u}{\partial x}\bigg|\,d\Omega.$ $$

This integral can be approximated as the discrete form of the total variation (TV) used in image processing as:

$$ $\text{TV}(u(t))=\text{TV}(u(t,\cdot,))=\sum_{i=1}^{N_{x}}\left|u(t,x_{i+1})-u(t,x_{i})\right|,$ (62) $$

where we discretize the spatial domain $\Omega=[x_{1},\dots,x_{N_{x}}]$ into $N_{x}$ gridpoints. A numerical scheme is called TVD if:

$$ $\text{TV}(u(t_{n+1}))\leq\text{TV}(u(t_{n})),\quad\forall\ n=1,\dots,N_{t},$ (63) $$

where we discretize the temporal domain $[0,T]=[t_{1},\dots,t_{N_{t}}]$ into $N_{t}$ gridpoints, and TV denotes the discretized form defined in [Equation 62](https://arxiv.org/html/2506.07003v2#A9.E62).

The TVD constraint in [Equation 63](https://arxiv.org/html/2506.07003v2#A9.E63) is a nonlinear inequality constraint, and enforcing it as a 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:

$$ $\text{TVD}=\sum_{n=1}^{N_{t}}\sum_{i=1}^{N_{x}}\left|u(t_{n},x_{i+1})-u(t_{n},x_{i})\right|,$ $$

as a regularization term. This approach is analogous to total variation denoising in signal processing [Rudin et al., 1992, Boyd and Vandenberghe, 2004].

LABEL:fig:lin_adv_tvd 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. In addition, ProbHardE2E leads to improved (tighter) uncertainty estimates.

## Peer Reviews

confidence: 4
contribution: 3
initial_rating: 6
presentation: 4
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.
 
2. 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, 
   - 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. 
   - 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.
   - 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. 
   - 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.
 
3. Typos/suggested corrections in the manuscript:
   - 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.
   - The paper frequently mentioned "Problem 7" while I think it would be better to use "Problem (7)" as it appears. 
   - In line 322-323 on page 6, the classical statistical approaches should be listed separately instead of being listed under "(ii) HierE2E".
   - Be consistent with using "HierE2E" and "Hier-E2E" (the paper is using both).
   - In line 424, "approximating" should be "approximately".
   - In lines 1327-1332, $\nabla L_{\hat{u}}$ should be $\nabla_{\hat{u}}L$.
   - In line 1448, in the last vector, $-J^{(*)}$ should be $J^{*}$.
   - In Eq. (44), add definition of $P_I=I-A^{\top}(AA^{\top})A$ (as you mentioned $P_I$ later but it was undefined).
soundness: 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.

2. It also enables joint UQ and constraint satisfaction, thus suitable for safety/physics/constrained engineering tasks.

3. The sample-free & closed-form CRPS, combined with analytic projection propagations gives substantial training speedups.

4. 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.
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.
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. 

2. 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.

3. In some places, the experimental results are not consistent with the claims made in the paper.

---

confidence: 4
contribution: 3
initial_rating: 6
presentation: 3
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.
2. 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. 
3. Include ablations on $Q$ (identity vs diagonal vs low-rank) with runtime/accuracy trade-offs, since $Q$ defines the projection geometry. 
4. 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.
soundness: 3
strengths: - The method is clearly formulated in a principled “predictor–corrector” view.
- The authors provide an exact handling for linear constraints.
- The training objective is close-formed and sampling free.
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.
weaknesses: - First-order DPPL approximation with no error bounds (risk under strong nonlinearity/variance). 
- Loss/evaluation assume independence while the projected posterior is generally correlated. 
- Practical experiments lock $Q$ to diagonal/identity, limiting the benefits of oblique projections and potentially biasing outcomes. 
- Nonlinear equality projection claims “strict feasibility,” but non-convexity and active-set issues are largely handled numerically without theoretical or robustness analysis. 
- Inequality-constraint metrics are thin (CE defined for equalities), so support for those claims is more qualitative than quantitative.

---

confidence: 3
contribution: 4
initial_rating: 8
presentation: 3
questions: please refer to the weakness
soundness: 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.
2. The training process can be sample-free, offering potential efficiency advantages.
3. In both time-series forecasting and PDE-solving tasks, the proposed method demonstrates strong empirical performance.
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.

ProbHardE2E 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.

The 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.
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.

---

confidence: 3
contribution: 3
initial_rating: 6
presentation: 2
questions: ## CRPS notation
Equation (2) uses the CRPS from Gneiting & Raftery (2007):  
$$
\mathrm{CRPS}(Y,y)=\tfrac12\mathbb{E}|Y-Y'|-\mathbb{E}|Y-y|.
$$
However, 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).  
Since 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.

---

## Notation in §3.1
The notation in §3.1 is somewhat inconsistent.
You 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. 
Moreover, ${\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. 

However, 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.


---

## Framing of log-likelihood and CRPS
Line 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.  
Additionally, 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.

---

## Discussion of unique solutions
A brief discussion on the *existence and uniqueness* of the constrained optimization problem (eqs 1 and 7) would enhance completeness.  
If the feasible set is empty, the algorithm naturally terminates; but if the solution is unique, uncertainty vanishes and the probabilistic framework offers no benefit.  
For instance, with linear equalities where $q=n$ and $A$ is full rank, the projection $P_{Q^{-1}}=0$ yields the deterministic solution  
$$
\hat{\mu}=A^{\dagger}b,\quad \Sigma=0,
$$
which 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.
More 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.


---

## Minor comments
- $\theta$ appears in l. 106 before being defined.  
- The domains of $\phi^{(i)}$ and $u^{(i)}$ (l. 106) are missing.  
- The spaces of $\xi$ and $z_\theta$ (l. 168–169) should be specified.

---

## Reference
Gneiting, T., & Raftery, A. E. (2007). *Strictly Proper Scoring Rules, Prediction, and Estimation.* **JASA**, 102 (477), 359–378.  
<https://doi.org/10.1198/016214506000001437>
soundness: 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.
The numerical experiments support the author's claims, showing improved performance while fulfilling the hard constraints.
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.  

For 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.

Overall, 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.
weaknesses: I believe, the presentation and readability could be improved, and some details of the method might need further clarification, see questions.


Decision: Accept (Poster)