HEAP: Hyper Extended A-PDHG Operator for Constrained High-dim PDEs
Authors: Mingquan Feng
Weixin Liao
Yixin Huang
Yifan Fu
Qifu Zheng
Junchi Yan
Venue: ICML 2025 poster (2025)
Keywords: PDE
Neural Operator
Constraints
High-dimension
Primary area: applications->chemistry_physics_and_earth_sciences

Abstract: Neural operators have emerged as a promising approach for solving high-dimensional partial differential equations (PDEs). However, existing neural operators often have difficulty in dealing with constrained PDEs, where the solution must satisfy additional equality or inequality constraints beyond the governing equations. To close this gap, we propose a novel neural operator, Hyper Extended Adaptive PDHG (HEAP) for constrained high-dim PDEs, where the learned operator evolves in the parameter space of PDEs. We first show that the evolution operator learning can be formulated as a quadratic programming (QP) problem, then unroll the adaptive primal-dual hybrid gradient (APDHG) algorithm as the QP-solver into the neural operator architecture. This allows us to improve efficiency while retaining theoretical guarantees of the constrained optimization. Empirical results on a variety of high-dim PDEs show that HEAP outperforms the state-of-the-art neural operator model.

HEAP: Hyper Extended A-PDHG Operator for Constrained High-dim PDEs
Mingquan Feng" Weixin Liao
Yixin Huang
Yifan Fu' Qifu Zheng' Junchi Yan
Abstract
To tackle the dimensionality challenge, recent advances in learning-based solvers often try to leverage the specific structure of PDE (Han et al., 2017; Wang et al., 2022; Yu et al., 2018) and develop memory memory--fficien training tech- niques (Hu et al., 2024; Shang et al., 2023). These methods have shown empirical success on a variety of PDEs. How- ever, they are designed to handle specific instances, meaning that they require expensive retraining from scratch when- ever there changes in the initial or boundary conditions. Neural operators (Gaby et al., 2024; Gaby & Ye, 2024) over- come these issues by learning the evolution operator of the solution function. The training is solely based on the PDE residual loss, no needing spatial discretization or solving the PDE directly. Once trained, the model (ideally) generalizes to new initial conditions with a single inference pass.
Neural operators have emerged as promising approach for high-dimensional partial differen- tial equations (PDEs). However, existing neural operators often have difficulty in dealing with constrained PDEs, which is practical setting where the solution must satisfy additional equal- ity or inequality constraints beyond the govern- ing equations. To close this gap, we propose a novel neural operator, Hyper Extended Adaptive PDHG (HEAP) for constrained high-dim PDEs, where the learned operator evolves in the param- eter space of PDEs. We first show that the evo- lution operator learning can be formulated as a quadratic programming (QP) problem, then un- roll the adaptive primal-dual hybrid gradient (A- PDHG) algorithm as QP-solver into the neural operator architecture. It allows improve effi- ciency while retaining theoretical guarantees of the constrained optimization. Empirical results on a variety of high-dim PDEs show that HEAP outperforms the state-of-the-an neural operator model.
Despite such progress, notable bottleneck (we also note the early efforts date back to (Lagaris et al., 1998)) to neural PDE solvers is the (in)equality constraints, which are com- mon in both real-world scenarios and numerical modeling, e.g. heat diffusion with temperature bounds and the Euler equation with entropy conditions. Simple methods incor- porating these constraints can result in unstable training, slow convergence, and violation of the constraints (Wang et al., 2021). While more advanced techniques have been proposed to handle constraints in either soft penalty (Wang et al., 2024) or hard encoding (Lu et al., 2021b), these methods come with either significant computational costs or network architecture restrictions and are limited to low- dimensional PDEs as discussed in (Feng et al., 2025).
1. Introduction
Partial differential equations (PDEs) are widely used across various disciplines in science, engineering, economics, and finance (Strauss, 2007; Folland, 2020). Since many PDEs lack analytical solutions, numerical techniques like finite element methods (Huebner et al., 2001) and recently neural network-based solvers (Blechschmidt & Ernst, 2021) have been developed. Although these approaches have witnessed significant progress, they face the curse of dimensionality, where the number of required data points grows exponen- tially the dimension d of the independent variables in- creases (Bellman, 1966).
In this work, propose novel quadratic programming- based neural operator, HEAP, for constrained high- dimensional PDEs, which can efficiently handle constraints in high-dimensional PDEs beyond existing SOTA neural operators. We first show that the evolution operator in the parameter space of constrained PDEs can be formulated as QP problem, where the QP constraints correspond the quadratic linear constraints in the original space. We then unroll the adaptive primal-dual hybrid gradient (A-PDHG) algorithm into the neural operator architecture to improve computational efficiency. Moreover, we theoretically prove the approximation capability of the proposed method and provide empirical results on variety of high-dim PDEs to demonstrate that HEAP outperforms existing SOTA neural operator model. The contributions are follows:
'Equal contribution 'Sch. of Computer Science & Sch. of Artificial Intelligence, Shanghai Jiao Tong University, China. Cor- respondence to: Junchi Yan <yanjunchi @ @et. edu..edd.cuu...... This work was partly supported by NSFC (62222607, 623B1009). Code avail- able at GitHub.
P Proceedings of the 42nd International Conference Machine L Learning. Vancouver, Canada. PMLR 267, 2025. Copyright 2025 by the author(s).

---

HEAP: Hyper Extended A A-PDHG Operator for Constrained High-dim PDEs
First time identify and formulate the evolution operator in the parameter space of constrained high-dimensional PDEs as QP optimization problem.
PDE. For general form PDBs, Zang et al. (2020) proposed adversarial network that solves PDEs. Hu et al. (2024) designed stochastic dimension gradient descent to reduce memory. Shi et al. designs an efficient amortized algorithm for differential operations in solving high-dim PDEs.
Propose a hyper-cxtended A-PDHG operator architecture for QP problem end-to-end learning and solving, with theoretical justification.
Another approach is to learn the solution operator, e.g. the mapping from initial conditions to the solution functions. The pioneering work is (Gaby et al., 2024), which learns the high-dimensional PDE solution operator by reducing the solution operator to the evolution of the parameters of a reduced-orde model (a neural network) and introduces another network to learn the evolution in discrete time steps by minimizing PDE residual. The second work (Gaby & Ye, 2024) extends the first work to continuous evolution by adapting Neural ODE (NODE) (Chen al., 2018) framework, improving the sample efficiency. We follow this line of work and propose novel neural operator that can handle constraints in high-dimensional PDEs.
Solve various types of PDEs on 5-20 dimensions with various constraints accurately and stably, generalizing across initial conditions.
2. Related Work and Preliminaries
Recent neural PDE solvers have shown the potential to overcome the limitations of traditional numerical solvers, especially for high-dim and constrained problems.
2.1. From Low-dim to High-dim Neural PDEs Solvers
The neural-networkk-based approaches for low dimensional PDEs can be mainly divided into two categories: physics- informed network and data-driven operator(Karmiadkki et al., 2021). Specifically, Innfore-In Neural Net- works (PINNs) (Raissi al., 2019) leverage neural networks to approximate the solution of PDE and enforce boundary loss and PDE residual loss at few selected points. This ap- proach has various extensions such as weak form (De Ryck et al., 2024), adaptive sample weight (McClenny & Braga- Neto, 2020), and integral loss (Salch et al., 2023). These methods may suffer from slow convergence and re-training requirements for new initial conditions. Alternatively, Deep- ONet (Chen & Chen, 1995; Lu et al., 2021a) a) leverages the universal approximation theorem of infiite-dimensional operators and directly fits the PDE solution operators in a data-driven manner, solving a family of PDEs with a single model. However, this approach requires a large amount of data for training and may not generalize well to new ini- tial conditions. itysics-informed neural operators (PINO) (Li et al., 2024) is a hybrid approach incorporating data and PDE residuals in the loss function. There exist many neural operators with other architectures, such as Fourier transform layer (Li et al., 2020), and graph neural networks (Ales)(Alet et al., 2019). The existing GNN-based models describe irregular spacial grids and their relative position as graph (Lötzzch et al., 2022: Horie & Mitsume, 2022; Bryutkin et al., 2024).
2.2. Neural PDE Solvers with Constraints
When PDE admits multiple weak or classical solutions, additional constraints- either equalities or inequalities act as regularisers that single out the physically (or financially) meaningful branch of the solution manifold.
Initial and boundary conditions are the most common equal- ity constraints, yet they may be partially unknown or insuf- ficient to guarantee uniqueness. In such cases inequality constraints can still be derived from first principles. A cel- Davier-Stokes
ebrated example is the three dieee-iinniional equation, whose uniqueness remains a Millennium Prize Problem; imposing an upper bound on the total kinetic energy is one standard way to discard non-physical blow- up solutions (Fefferman, 2006). Similarly, the multi-asset Black-Scholes equation lacks a unique solution in prac- tice; no-arbitrage bounds on the option's 's theta restore well- posedness and are routinely used in quantitative finance (Hull & Basu, 2016). We visualise this effect in Fig. 3, where distinct constraint sets lead to markedly different Black-Scholes solution surfaces.
A few recent A14PDE studies have begun encode in- equality information, e.g. derivative bounds convex con- straints, into physics-informed losses (Hoshisashi et al., 2024; Moro & Chamon, 2025). However, these methods rely on dense spatial discretisation and have only been demon- strated on low-dimensional (d < 3) problems, Scalability to the high-dim regime (d 2 10) targeted in this work re- mains open; HEAP addresses this gap by reformulating the constraint enforcement differentiable QP layer and un- rolling an adaptive PDHG solver inside the neural operator.
For high-dimensional PDEs, the curse of dimensionality poses challenge, and specialized methods are designed which mainly fall in the physics-informed category. (Han et al., 2018; 2017) proposed the DeepBSDE, reforming a special class of hyperbolic PDEs as backward stochastic dif- ferential equations (BSDEs) by the Feynman- Kac formula. The deep Ritz method (Yu al., 2018) solves high-dim variation PDEs by minimizing the energy functional of the
Previous works incorporate constraints into neural solvers in two ways: soft and hard. Soft methods introduce penalty
2

---

HEAP: Hyper Extended A-PDHG Operator for Constrained High-dim PDEs
terms into the objective. For instance, physics-infommn -nnfroeed neural networks (PINNs) (Raissi et al., 2019) apply initial and boundary conditions penalty terms. Other examples include enforcing turbulence models (List et al., 2022), com- plex geometries (Wang et al., 2024), and spectral represen- tation (Du et al., 2023). However, in many real-world cases, the constraints need to be strictly enforced which also the focus of this paper. Hard constraints can be implemented in various ways, e.g. encoding structure of constraints into network architecture (Lu et al., 2021b; lb; Richter-Powell et al., 2022), using an implicit constraint layer (Négiar et al., 2022), incorporating numerical solvers into the network (Saad et al., 2022; Chalapathi et al., 2024). However, these methods are limited to low-dim PDEs and do not scale well to high-dim problems. In our work, we propose the first soft-constrained end-to-end neural solver for high-dimensional PDEs.
models for solving:
/'" Voue: 0, 10-u.ull dt
min
E OUND
Jo
1
(2(
(2)
s.t.
H Muoo.1(x 2 0, x SH te (0, (),
0, = Vo.). t (0,T).
In this formulation, the PDE operator learning is converted to the parameter evolution operator learning, which is a linearly constrained nonlinear optimization problem.
Approximated QP Operator. To design the architecture of Var in principled way, we exploit the quadratic structure of Eq. 2 and approximately reformulate Vw as a QP problem w.r.t. O,. A general QP problem can be formulated as:
min xE
2 Txxx+t cxx, s.t. Ax b,
(3)
2.3. Genreal Formulation of Constrained PDEs
where Q is a symmetrie, positive semi-definite matrix. V can be implemented as a QP solver, after three approxima- tion steps. Firstly, the ( norm in the objective function and 92, in the constraint are approximated by Monte Carlo sampling, i.e. evaluation over sets of randomly sampled points X, Xx from 22, QH respectively. Secondly, ue, in constraint A is replaced by its first-order Taylor expansion with small time step At, i.e. no, + Vouo, At. Finally, we add lower and upper bounds I, u of 6g to the constraints to enhance the numerical stability. The resulting operator, termed as Vor, leads to a new optimization problem:
We consider nonlinear PDE initial value problem with linear constraints in its general form:
0,u(x, t) = (x,,
xe (0,
"(x, 0) = g(x), H[u[(x, t) 20 0
x XEΩ,
(1)
x QH t (0,T)
where u(x, t) is the solution function, F is a possible non- linear operator, g(x) is the initial condition, Q C Rd is the spatial domain, and T is the time horizon. In addi- tion to the PDE, have inequality constraints H[u] 2 b, where H is a linear operator. The constraints are defined on the subdomain Uzz The constraints might be rooted in the multi-solution nature of some PDEs, associated with physical laws, or numerical stability requirements.
Veue, (X) - Flus, ||
VOP (0,) = argmin 6.E(Luu
(4)
H[u, + Voue, -..:t 2 0.
s.t.
The problem above is a QP instance, by keeping the bound (1, u) and assembling the components (Q, A, b, c) as:
3. Methods
Q V Veue, (X) Voue, (X),
A = H Ar (X,,,
3.1. Problem Formulation
(5)
- Flue,] Vous, (X),
Evolution Operator in the Parameter Space. Under the high-dimensional settings, it is usually infeasible to numeri- cally calculate the PDE solution directly. In this paper, we approximately surrogate the solution function u(;,t) with a neural network ansatz ug (-), and learn the evolution of the parameters 0 starting from an initial parameter 0g. The evolution governed by a hyper-operator Vu, which takes the current parameter 0r and outputs its temporal derivative, i.e. 0, (O.). In the definition of the constraint PDEs, we use the formulation in Equation. 1.
b = - -H(u,, (XHH.
Adaptive primal-dual hybrid gradient (PDHG) Algo- rithm for QP. The QP problem Eq. 3 can be solved by various classic algorithms, e.g. matrix-faatorizaaion--asss methods including interior-point and simplex, and first- order methods based matrix-vector production including PDHG and ADMM. Considering the high-dimension na- ture, we choose PDHG for its simplicity and scalability. It converts the QP problem into a saddle-poin problem:
The overall target is to fit the hyper-operator Vú minimizing the expected PDE residual over a distribution D of f and under the constraints. Formally, we drive the following
min max xE[La] y20 x 1x c x -Ax (Ax-b), b), (6)
3

---

HEAP: Hyper Extended A-PDHG Operator for Constrained High-dim PDEs
Algorithm 1 A-PDHG embodiment for solving QP in Eq. 3
Algorithm 2 Hyper-Extended A-PDHG (i.e. HEAP)
1: Input: Step-size {(Tk, 0k, 9k, nh),. initial points
1: Input: ansatz parameter O,, hypernet Na., iteration K
(x", y²),, QP instance (Q,
Assemble QP instance (Q, A, b, c, u) from Eq. 5
2: Initialize x" = x²
3: Estimate initial points (X", YO), step-sizes and expan-
3: for k 0, 11,,.... do
sion weights (W) by (0,)
4: md = (1 - 3x)x** + AEx**
4: Encode X x (1 + W Y y W
5:
x x Projj[.ui fx* Th (Qxmd c- A )
5: Initialize X = X
6:
7:
y yktl -Projk, [y
y k+1 = y² + ok A(m (x xx*) + x(tl))
md
7:
6: for k 0, 1, ...,.... K do
X md = (1 - Bx + 3x
8:
ght1_ = -1- 3x)x* + 3
8:
X md X TTk (Qx md Wh
md tc- A AY
X X
Clamp [1.1] (X
9:
if(x" y") converges then break
10: Output: Solution x";
10:
Y rk+1 = Y + 0x A(Ta(xHt1 - x²) +
X
W
11:
Yk+1 ReLL (+++1)
12:
X (1 3r
add X
+ 3
and iteratively updates the primal and dual variables
13: Decode x - -
W
y,y from initial guesses, leveraging gradient information. Among the variants of PDHG, the momentum momentum-acceleratte accelerated PDHG (A-PDHG) (Lu & Yang. 2023) is chosen for its fast convergence, as summarized in Algorithm 1.
14: Output: Parameter evolution estimation 0 xx
intermediate variables.
3) Estimate the initial points yx), step-sizes and ex- pansion weights {W} by a hypernetork (w. The hypernetwork weights w are trainable.
The Advantages of A-PDHG. Prior high-dim neural op- erators cast parameter evolution as unconstrained least- squares fit or treat it as black-box regression, making it impossible to encode hard safety or no-arbitrage rules. By rewriting the evolution step as a quadratic programming (QP) we can natively embed linear inequality information and thus extend a broader spectrum of constrained PDEs.
4) Convert the projection operators to differentiable func- tions, e.g. ReLU, Clamp, to enable end-to-ent training.
The HEAP is shown in Algorithm 2, with the modifications marked as blue, and the workflow is illustrated in Fig. 1.
A vanilla A-PDHG solver typically needs 100+ primal-dual iterations our largest QP instances--rrohibitive prohibitively expen- sive during back propaaprpaa.tion.. HEAP truncates the algo- rithm to K 3 steps and learns step sizes, momentum and expansion weights from data, retaining algorithmic priors while reducing compute time by two orders of magnitude. Table 6 confirms that K <3 already matches, or surpasses, the accuracy of the full A-PDHG loop.
Training and Testing, The training loss is the accumulated residual r along the time horizon, whose increment é is the weighted sum of the PDE residual fr.p the physical constraint violation pe. and the numerical constraints ir,
Qú. + c TT.1"² |
fr,pde
(min(AAd b, b.0) /
it.po
3.2. Unroll A-PDHG into Neural Architecture
|min((, - 1,u Úr 0) 2 /1.11' 7. X Apetie.p ge + Aact t,me
Despite the fast theoretical convergence rate of A-PDHG, it still suffers long-tail convergence in practice, partially due to the difficulty in tuning the hyperparameters. In our PDE evolution operator settings, the long iterations lead to slow and unstable backpropagation. To address this issue, we truncate and extend A-PDHG into a hyper-network archi- tecture, termed as Hyper-Extended A-PDHG (HEAP), with the following modifications:
where the {A} are tunable scalar weights. The forward pass is the NeuralODE (Chen et al., 2018) integration of both the residual r; and the ansatz parameter 0,;
TT = ['" Ti dt
1) Truncate the iterations to fixed number K. This is motivated by the observation that the A-PDHG algorithm converges quickly in the first few iterations.
Or 80 + 1- T
(8)
(HHEA((0)) dt
2) Extend the computation to a learnable latent space, com- pensating for the inaccuracy due to truncation. The latent space is parameterized by linear expansion weights (W)
All operations in the HEAP model are differentiable, and the gradients are computed by backpropagation via the ad- joint method (Chen et al., 2018). The adjoint method
4

---

HEAP: Hyper Extended A-PDHG Operator for Constrained High-dim PDEs
Extended APDHG iteration x
PDE & Constraints
Q, c, W
HEAP Evolution Operator
Assemble
QP Instance Q.A.b,c,l,u
oo
x
xk=1
Encode
W
Project
Decode
d0
1t
dt
Extend Params
yo
A, W
Hypemet NN'
yk
yk+1
W, W,x",
V
Ansatz
U
PDE
u(x,
u(x = =t)
u(xt = T)
Residual
dr, R(O, t)
Constraints Violation
r(t =0
r(t=
Backpropagate
Figure The PDE evolution embedded into the ansatz: parameter space, approximated as QP operator and solved by the HEAP Its operator unrolled from the A-PDHG algorithm, with the iterations truncated and extended to learnable latent space. The hypernetwork Nr., estimates the initial points, step sizes, and expansion weights. The projection operators are replaced by differentiable functions.
an efficient technique for computing gradients by solving backward equations, which propagate gradient information from the output to the input.
Theorem 3.1. [Existence of HEAP] Given any QP instance M (Q, A,b, l,,,,u),, u), HyperNet N. and its corresponding primal-dual sequence (Iy, yk Wk)ssK generated by the PDQP algorithum within K iterations, there exists a K-layer HEAP with parameter assignment θHEAP that can output the same ilerative solution sequence.
In inference, first the initial parameters 0g are obtained by fitting the ansatz with the initial function. Then 0g is fed to the trained HEAP model to get the solution function. The test metric is the PDE relative residual at the final time T:
Approximation Capability. Secondly, verify the ap- proximation capacity of HEAP. The proof of this theorem is given in Appendix A.
L2R
||u0r Upeeel//|ueerll
(9)
Theorem 3.2. [Linear Convergence] Given an approxima- tion error bound e, there exists a K-layer HEAP with O neurons that exhibits linear convergence, i.e. for any QP instance M (Q, A, c, u) and (z,y) )) R R
Intuitive view, A-PDHG algorithm. Adaptive PDHG alter- nates primal descent 0 (reducing mismatch) with a dual ascent on the Lagrange multipliers y (penalising violations). Each update is projected onto the feasible set.
satisfying < u, it follows that
HEAP network. HEAP unrolls a fixed number K of A-PDHG iterations into a feed-forward network. Pri- mal and dual variables are lifted from vectors to matri- ces-effectively solving K QP steps in parallel-and the lifting weights are produced by hyper-network that per- ceives the current PDE state. The entire construct is end-to- end differentiable.
L(X My - L(x, Y M) ; <e
(10)
4. Experiments
All experiments were conducted a machine with 1TB memory, 144 cores Intel Xeon Platinum 8352V CPU, and 8 GPUs (NV GeForce RTX 4090).
3.3. Theoretical Analysis
We provide our theoretical study follows.
4.1. PDEs Selection
Alignment with A-PDHG Algorithm. Firstly, we prove HEAP's alignment with A-PDHG, i.c. there exists a HEAP that could replicate the primal-dual sequence generated by A-PDHG. A stricter expression is provided as follows, and the proof of this theorem given in Appendix
We benchmark HEAP on four prototype equations, each tested in both an unconstrained form (suffix -UC) and constrained form (suffix -C). For -C cases we en- force inequality constraints that must hold everywhere in
5

---

HEAP: Hyper Extended A-PDHG Operator for Constrained High-dim PDEs
space-time; for -CC we drop those constraints but keep the same governing PDE. What's more, the numerical constraints which are applied in both -C and -CC cases are the bound 1, u of parameters 8,, which is chosen as
Table I. Hyper-parameters for CSO on Burgers; HEAP uses iden- tical values. T represents the time range of the equation, steps denotes the number iterations of the equation in time.
dim
steps
input dim
dropout
learning rate
1 1= 1 u = 1. Unless noted otherwise, all experiments use spatial dimension de (5, 10, 15, 20).
5
0.2
10
1921
0.30
5 x 10
10
0.2
10
2121
0.30
x 10-5
15
0.2
20
2321
0.30
2 x 10-4
Heat Equation.
20
0.2
20
2521
0.30
3 10-4
xE [-1,1]], LE (0,7), 0,u= u,
d
u(1,0 0) = Theat (I) = a II sin (xx,),
Parameters a, denote drift, volatility and the risk-free rate; 6, Uh U; The 7I are calibration constants. k E [0.5, 1.5] is the the adjustable coefficient. Constraints (T1)-(2) en- code standard no-arbitrage bounds on the option's theta.
0,u= S Aheas
(only for Heat-c).
(17)
where a E [0, 1] is a tunable parameter.
Burgers Equation.
4.2. Experiments Implementation Details
d.u 2 + (u
2+d 2d
Oz U,
Setup. The extended dimension of HEAP is 20; the itera- tion is 3, and the hypernetwork Ni,, is a 5-layer ResNet with 1000 hidden units. The HEAP is trained with Adam opti- mizer, default batch size 64, and 100 epochs. E.g., Table I summarises the settings used for CSO our Burgers type equation experiments; HEAP employs the same configura- tion. The ODE is integrated with the Runge- Kutta (Butcher, 1996) with a time step 1 1/10 of the time horizon.
1
u(x,0) 0)- = 9burg (2) =
+ exp(k: k: t 2 Ex x4))
Oeu S Aburg (only for Burgers-C).
(18)
where & E [0.8, 1.3] is the the adjustable coefficient.
Reaction-Diffusinn (RD) Equation.
Sampling of initial parameters. The training set consists of 150K randomly generated 0 from normal distribution, and the test set is 250 randomly generated 0 from another normal distribution with slightly different mean and vari- ance. Unless otherwise noted, the initial ansatz parameters Gg are drawn i.i.d. from a Gaussian distribution with mean O and standard deviation v0.5 for training, and from a wider ((, 1) for test; this encourages robustness to distribution shift while keeping the training residual numerically stable.
0,u = 1 + min [1, (t,
ut) = KK-- 1-
Fi
(1-T)
u(x,0) = (x)(x) = + x
+ sin
0,u S Azd (only for RD-C).
(19)
Composite loss weights. The accumulated residual r's (see Eq. (7)) is weighted sum of three terms: (i) PDE residual, (ii) physical constraint violation, and (iii) numerical box- constraint violation. To balance their magnitudes during optimisation we set Apde = 10-2, Ape = 1, and A = 10, respectively, for all experiments. The same coefficients are applied the baseline CSO to ensure fair comparison.
where A 1/Vd is constant representing the wave fre-
quency, K E [0.5, 1.5] is the tunable parameters.
Black-Scholes (BS) Equation.
d
1
d
0,u = Σ μ R. 0z u 20-0"
02
x
2
i=1
- [[11- (1-5)QQ4 Q1) +rlu,
4.3. Peer Methods and Protocols
Quu) ReLU ReLU(u- vn)
Yh 7
To our best knowledge, only a handful of neural-operator models can handle high gigh--imensional PDEs, and none of them explicitly target strong inequality constraints. We thus adopt the recently proposed Control-based barrol Solution Opera- for (CSO) (Gaby & Ye, 2024) as the SOTA baseline. CSO parametrises the parameter trajectory O with a black-box network; for fair comparison we reuse HEAP's hypernet- work backbone, identical optimiser settings, batch size and training horizon across all tasks.
7h
7
+ n,
Uh 11
u(z,0) 0) = (x) lo min Zi.
0,u, 2 ru
(TI, only for BS-C),
0,u, 2 r Abs
8 (T2, only for BS-C).
(20)
6

---

HEAP: Hyper Extended A-PDHG Operator for Constrained High-dim PDEs
Table PDE residual Fl on unconstrained equations com- parison between the strong and recently proposed SOTA baseline CSO (Gaby Ye, 2024) and our proposed HEAP on various PDEs.
Table 3. PDE residual 1.2 (Only for Heat-CC) and accumu- lated PDE residual r for (Others) on constrained equations com- parison between the strong and recently proposed SOTA baseline CSO (Gaby & Ye, 2024) and our proposed HEAP on various PDEs.
Eqn
Heat
Eqn
Constrained Heat
Dim
d=5
d=10
d=15
d=20
Dim
d=5
d=10
d=15
d=20
CSO
0.0222
0.02188
0.0218
0.0273
HEAP
0.0041
0.0068
0.0106
0.0141
CSO
0.0030
0.0187
0.0136
0.0407
HEAP
0.0025
0.0044
0.0092
0.0169
Eqn
Burgers type
Eqn
Constrained Burgers Type
CSO
0,0028
0.0058
0.0093
0.0095
HEAP
0.0020
0.0032
0.0079
0.0089
CSO
0.9798
0.9125
0.9764
0.9529
HEAP
0.7377
0.6281
0.6926 0.8778
Eqn
Reaction Diffusion
Eqn
Constrained Reaction Diffusion
CSO
0.0028
0.0024
0.0076
0.0021
HEAP
0.0021
0.0016
0.0024
0.0020
CSO
0.3560
0.4746
0.5687
0.5102
Eqn
Equation
HEAP
0.3560
0.3768
0.4730
0.3900
Eqn
Constrained
Type 1
CSO
0.0329
0.0439
0.0391
0.0494
HEAP
0.0270
0.0308
0.0340
0.0433
CSO
0,9462
0.8649
0,9145
1.000
HEAP
0.7022
0.5667
0.5232
0.4756
Eqn
Constrained Black-Schole Type 2
Metrics with ground-truth reference. For those PDEs whose exact or high-fidelity reference solu- tions are available (Heat-UC/C, Burgers- UC, RD-UC, BS- UC), we report the 1.2 relative error of the PDE residual defined in Eq. (9). References are analytical for Heat, and produced by high-resolution DeepBSDE solver for Burgers, RD and BS.
CSO
0,9570
0.8625
0,9260
1.0000
HEAP
0.6046
0.4241
0.4617
0.4670
Table 2, it can be found that under constrained conditions, the improvement of the HEAP to CSO is more significant, demonstrating the effectiveness of HEAP in solving high- dimensional PDEs with constraints.
Metrics without ground-truth reference. For Burgers-C, C, RD-C and BS-C the reference solution is unknown; we therefore adopt the accumulated residual A (see Eq. (7)), which jointly measures the physics error and constraint mismatch along the trajectory.
The normalized terminal time physical constraint violation .t.pe on constrained equations are also compared in Table 4. Experimental results show that the Heap performs excel- lent in satisfying the physical constraints of the equation. Especially in the Black-Scholes equation with more com- plex constraints, Heap significantly reduces the physical constraint violation compared with the CSO, showing the potential of Heap to solve complex linear constraints.
Constraint satisfaction.
In every constrained task we additionally report the nor- malised terminal-time physical constraint violation hr.ee (Eq. (7)), irrespective of the availability of ground truth.
To intuitively show the advantages of the Heap over the CSO baseline, we visualized the solution and resisual of the con- strained 10 10-dimensional Black-Scholes equation on Fig.2. We changed two of the dimensions and set the other dimen- sions to O, and performed grid solution on a [0. 1] plane. The image indicate that the accumulated PDE residual of the Heap is significantly smaller than that of CSO.
4.4. Results and Discussion
Results for Unconstrained Equations. The PDE residual on unconstrained equations (Heat-UC, Burgers-UC, RD- UC, BS-UC) are compared in Table 2. The L2RE of PDE of HEAP is consistently lower than that of CSO, indicating better accuracy and generalization.
Results for Constrained Equations. The PDE residual on constrained equations (Heat-C, Burgers-C, RD-C, BS-C-TI, BS-C-T2) are compared in Table 3. As mentioned in 4.1.1, we set some reasonable constraints for each type of equation, where the Black-Scholee equation contains two relatively more complex no-arbitrage constraints. The L2RE of PDE of HEAP is consistently lower than that of CSO, indicating better accuracy and generalization. By comparison with
Time and Space Usage. The table 5 shows the performance indicators of Heap and CSO methods for Dim=5, 10, 15, 20, 50, 100, include time consumption, CPU memory us- age(CMem), and GPU memory usage(GMem). Without any optimization applied, the HEAP costs 1.6x GPU memory and about 2x training/inference time compared to the base- line CSO. The time and memory usage are nearly linearly increasing with the dimension, thus HEAP is scalable to
7

---

HEAP: Hyper Extended A-PDHG Operator for Constrained High-dim PDEs
Table physical constraint violation Fx2 comparison he- tween CSO and HEAP constrained PDEs. The best is marked.
Table 5. Comparison of runtime (sec.) and memory usage (MB) for CSO and HEAP on Burgers-type equations. The test time is evaluated with batch, and training time calculated with and 156 batches (156 batches is approximately correspond to
Eqn
Constrained Heat
Dim
d=5
d=10
d=15
d=20
epoch). The batchsize is set as 64. The symbol * indicates that
CSO
0.2944
0.2843
0.0080
1.58 x 10-4
training not completed due to GPU memory overflow.
HEAP
0.2540
0.2234
0.0168
3.21
10
Dim
Method
CMmem
GMem
T
Tirsin
Eqn
Constrained Burgers Type
5
CSO
4.47
12336
0.90
24.12
596
CSO
0.4544
0.5606
0,4335
0,4098
5
HEAP
4.50
15865
2.94
67.48
1071
HEAP
0.4224
0.4603
0.5948
0.3345
10
CSO
4.21
14157
1.21
32.06
699
Eqn
Constrained Reaction Diffusion
10
HEAP
4.15
18696
3.52
79.63
1136
CSO
0.2556
0.7593
0.6725
0.6936
CSO
4.44
16075
42.02
902
HEAP
0.2555
0.2585
0.2359
0.1794
HEAP
4.57
21927
4.31
97.87
1668
Eqn
Constrained Black-Schole Type 1
20
CSO
4.15
18085
2.31
53.44
1119
HEAP
4.22
25405
112.83
1896
CSO
0.6368
0.9386
0.6837
0.694
HEAP
0.2843
0.1932
0.3310
0,2231
50
CSO
4.39
32166
4.85
109.90
50
HEAP
4.50
51654
11.19
198.46
Eqn
Constrained Black-Schole Scholes Type
100
CSO
4.05
71958
11.54
CSO
0.7466
0.7919
0.6546
0.6936
100
HEAP
4.29
80000
25.46
HEAP
0.5595
0.3865
0.0719
0.0165
Pred u (CSO)
Pred (HEAP)
Table 6. Ablation study on the number of unrolling layers and extened dims in HEAP on BS equation with constraint type 1. The PDE residuals are normalized by constant scale.
Unrolling Layers
224
Extended Dim
Residual
2
3
5
PDE
0.9630
0.5205
1.7511
4.5298
Con
0.5084
0.1163
0,4339
0,3460
Res (CSO)
Res (HEAP)
10
PDE
0.8557
0.5246
0.5857
1.5917
Jaar
Con
0.5869
0,3310
0,7472
0.6787
20
PDE
0.9181
0.4298
0.5283
1.5857
4230
Con
0.3455
0.1203
0.2708
0.6474
AM
ever, K = is relatively robust choice across different PDE families.
Figure 2. The predict solution (Pred.u) and accumulated PDE residual (Res) Heap and CSO Black-Scholes Equation (dim- 10) with constraint type 2.
5. Conclusion
We propose novel neural operator for Constrained High- dim PDEs and devise the HEAP architecture for end-to-end training. We theoretically prove the approximation capabil- ity and provide empirical results on a variety of high-dim PDEs to demonstrate that it outperforms existing state-of- the-art neural operator models and shows better constraint satisfaction and generalization.
higher dimensions in terms of computation.
Ablation Experiment for the hypernetwork. The number of unrolling layers and extended dimension in the hyper- network are the core hyperparameters of the Heap method. We show the results of the ablation experiment on these two hyperparameters in Table 6 (taking the 15 15-dimensional Black-Scholee equation as example), whose metrics value is also the accumulated residual i; (see Eq. (7)). The results show that both unrolling layers = 1,2.3,4 and extended di- mension 5,10.20 in the hypernetwork have a significant impact on final accuracy and constraint satisfaction. How-
Impact Statement
This paper presents the work that aims to advance the field of Machine Learning for Science (ML4Sci). There are many potential societal consequences of our work, none which we feel must be specifically highlighted here.
8

---

HEAP: Hyper Extended A-PDHG Operator for Constrained High-dim PDEs
References
Gaby, N. and Ye, X. Approximation of solution operators for high-dimensional pdes. arXiv preprint arX1v:2001.10385, 2024.
Alet, F., Jeewajee, A. K., Villalonga, M. B., Rodriguez, A., Lozano-Perez, T., and Kaelbling, L. Graph element networks: adaptive, structured computation and memory. In International Conference on Machine Learning, pp- 212-222. PMLR, 2019.
Gaby, N., Ye, X., and Zhou, H. Neural control of para- metric solutions for high-dimensional evolution pdes. SIAM Journal on Scientific Computing, 46(2):C155- C185, 2024.
Bellman, R. Dynamic programming. science, 153(3731): 34-37, 1966.
Han, J., Jentzen, A., et al. Deep toarinng-baser numerical methods for figh-dimensional parabolic partial differen- tial equations and backward stochastic differential equa- tions. Communications mathematics and statistics, (4):349-380, 2017.
Blechschmidt, J. and Ernst, O. G. Three ways to solve partial differential equations with neural networks review. Muttllnitillnngen, 24(2)=2021000066 2021.
Bryutkin, A., Huang, J., Deng, Z., Yang. G., Schönlieb, C.-B., and Aviles-Rivero, A. Hamlet: Graph transformer neural operator for partial differential equations. arXiv preprint xxXv:2402..03411, 2024.
Han, J., Jentzen, A., and E, W. Solving high-dimensional partial differential equations using deep learning. Pro- ceedings of the National Academy of Sciences, 115(34): 8505-8510, 2018.
Butcher, J. C. A history of runge-kutta methods. Applied namerical mathematics, 20(3):247-260, 1996.
Horie, M. and Mitsume, N. Physics-embedded neural net- works: Graph neural pde solvers with mixed boundary conditions. Advances in Neural Information Processing Systems, 35:23218 2222299, 2022.
Chalapathi, N., Du, Y., and Krishnapriyan, A. Scal- ing physics-informed hard constraints with mixture-of- experts. arXiv preprint arXiv:2402.13411 2024.
Chambolle, A. and Pock, T. On the ergodic convergence rates of first-order primal-dual algorithm. Mathematical Programming, (59(1))233-238) 2016.
Hoshisashi, K., Phelan, C. E., and Barucca, P. Physics- informed neural networks for derivative-oonstranedd pdes. Proceedings of the ICML 2024 Workshop on Al for Science, 2024. to appear.
Chen, R. T., Rubanova, Y., Bettencourt, J., and Duvenaud, D. K. Neural ordinary differential equations. Advances in neural information processing systems, 31, 2018.
Hu, Z., Shukla, K., Karniadakis, G. E.,, and Kawaguchi, K. Tackling the curse of dimensionality with physics- informed neural networks. Neural Networks, 176: 106369, 2024.
Chen, T. and Chen, H. Universal approximation nonlinear operators by neural networks with arbitrary activation functions and its application to dynamical systems. IEEE Transactions Neural Networks, 6(4):911 917-777, 1995.
Huebner, K. H., Dewhirst, D. L., Smith, D. E., and Byrom, T. G. The finite element method engineers. John Wiley & Sons, 2001.
De Ryck, T., Mishra, S., and Molinaro, R. wpinns: Weak physics informed neural networks for approximating en- tropy solutions of hyperbolic conservation laws. SIAM Journal on Numerical Analysis, 62(2):811-88) 1-841, 2024.
Hull, J. C. and Basu, S. Options, Futures, and Other Deriva- tives. Pearson Education, Boston, 9 edition, 2016.
Du, Y., Chalapathi, N., and Krishnapriyan, A. Neural spec- tral methods: Self-supervised learning in the spectral domain. arXiv preprint arXiv:23/2.05225, 2023.
Karniadakis, G. E., Kevrekidis, I. G., Lu, L., Perdikaris, P., Wang, S., and Yang, L. Physics-informed machine learning. Nature Reviews Physics, 3(6):422-440, 2021.
Fefferman, C. L. Existence and smoothness of the navier- stokes equation. In Carlson, J., Jaffe, A., and Wiles, A. (eds.), The Millennium Prize Problems, pp. 57-67. Clay Mathematics Institute, 2006.
Lagaris, I. E., Likas, A., and Fotiadis, D. I. Artificial neu- ral networks for solving ordinary and partial differential equations. IEEE transactions on neural networks, 9(5): 987-1000, 1998.
Feng. M., Huang, Y., Liu, Y., Liao, W., Liu, Y., and Yan, Singer: Stochastic network graph evolving operator for high dimensional pdes. ICLR, 2025.
Li, 7.,, Kovachki, N. B., Arizzadenesheli, K., Bhattacharya, K., Stuart, A., Anandkumar, A., et al. Fourier neural operator for parametric partial differential equations. In International Conference on Learning Representations, 2020.
Folland, G. B. Introduction to partial differential equations. Princeton university press, 2020.

---

HEAP: Hyper Extended A-PDHG Operator for Constrained High-dim PDEs
Li, Z., Zheng, H., Kovachki, N., Jin, D., Chen, H., Liu, B., Azizzadenesheli, K., and Anandkumar, A. Physies- informed neural operator for learning partial differential equations. ACM/JMS Journal of Data Science, 1(3):1-27, 2024.
Saleh, E., Ghaffari, S., Bretl, T., Olson, L.,, and West, M. Learning from integral losses in physics informed neural networks. arXiv preprint 2rxiv:2305.17338 2023.
Shang, Y., Wang, F., and Sun, Randomized neural net- work with petrov-galerkin methods for solving linear and nonlinear partial differential equations. Communications in Nonlinear Science and Numerical Simulation, 127: 107518, 2023.
List, B., Chen, L.-W., and Thuerey, N. Learned turbulence modelling with differentiable fluid solvers: physics-basee -asses loss functions and optimisation horizons. Journal Fluid Mechanics, 949:A25, 2022.
Shi, Z., Hu, Z., Lin, M., and Kawaguchi, K. Stochastic taylor derivative estimator: Efficient amortization for ar- bitrary differential operators. In The Thirty-eighth Annual Conference on Neural Information Processing Systems.
Lötzsch, W., Ohler, S., and Otterbach, J. S. Learning the solution operator of boundary value problems using graph neural networks. arXiv preprint arXiv:2206. 2022.
Lu, H. and Yang. A practical and optimal first-order method for large-scale convex quadratic programming. arXiv preprint arXiv:2311.077100 2023.
Strauss, W. A. Partial differential equations: An introduc- tion. John Wiley & Sons, 2007.
Wang, H., Li, Dwivedi, A., Hara, K., and Wu, Beno: neural operators for elliptic pdes. arXiv preprint arXiv:2401.0932. 2024.
Lu, L., Jin, P., Pang, G., Zhang, Z., and Karniadakis, G. E. Learning nonlinear operators via deeponet based on the universal approximation theorem of operators. Nature machine intelligence, 3(3):218-22) 2021a.
Wang. S., Teng, Y., and Perdikaris, P. Understanding and mitigating gradient flow pathologies in physics-informed neural networks. SIAM Journal on Scientific Computing, (((55:33055AA300 2021.
Lu, L., Pestourie, R., Yao, W., Wang. Z., Verdugo, F., and Johnson, S. G. Physics-informed neural networks with hard constraints for inverse design. SIAM Journal Scientific Computing, 43(6):B1100 B1132, 2021b.
Wang, Y., Jin, P., and Xie, H. Tensor neural network and its numerical integration. arXiv preprint arXiv:2207.02754 2022.
McClenny, L. and Braga-Neto, U. Self-adaptive physics- informed neural networks using a soft attention mecha- nism. arXiv preprint arXiv:2009.0455
Yu, B. et al. The deep ritz method: deep learning- caaningghhasdd nu- merical algorithm for solving variational problems. Com- munications in Mathematics and Statistics, 6(1):1-12, 2018.
2020.
Moro, V. and Chamon, L. F. O. Solving differential equa- tions with constrained learning. In Proceedings of the 13th International Conference Learning Representa- tions (ICLR), 2025. to appear.
Zang, Y., Bao, G., Ye, X., and Zhou, H. Weak adversarial networks for high-ddimensionl partial differential equa- tions. Journal of Computational Physics, 411::09409,, 2020.
Négiar, G., Mahoney, M. W., and Krishnapriyan, A. S. Learning differentiable solvers for systems with hard con- straints. arXiv preprint arXiv:2207.08675. 2022.
Raissi, M., Perdikaris, P., and Karniadakis, G. E. 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.
J., Lipman, Y., and Chen, R. Neural conservation laws: A divergence-free perspective. Ad- vances in Neural Information Processing Systems, 35: 38075-38088, 2022.
Saad, N., Gupta, G., Alizadeh, S., and Maddix, D. C. Guid- ing continuous operator learning through physics- based boundary constraints, arXiv preprint ar20v.2212.07477, 2022.
10

---

HEAP: Hyper Extended A-PDHG Operator for Constrained High-dim PDEs
A. Theorem Proof
HEAP exhibit linear convergence, i.e.
Proof of Theorem 3.1. We use mathematical induction.
L(X y M) - L(x, Y :M) [ M) [1 ; ;M)
Base Step. At iteration & = 0, both APDQP and HEAP are initialized by hypernet Ne., and the base case holds trivially.
||==="|| + 11v-0"1l"" (y ))) A(xx xx))
Inductive Step. Assume that for some k 2 0, X = x² and y = =y We now show that X = x*+1 and y
(E(
(15)
gx for a suitable choice of parameters in HHEAP.
| |a=="012
11y y y0 (u- - T A(x - x20) is
Since 26
For the APDQP update, we know:
constant given pair of (x,, v). we can take K(c) = 1 +1 to obtain:
ktt Projji.u (x² Th 1 A-AA))) 1²))
+c-
L(X*, y; -) Lxx, M) S 6
2 <E
(16)
g = Proj" + ox(6b Alrx(1ht1 - x²)) + x ()))
(11)
Because the network width is a fixed number, the number of ncurons at each layer is bounded by a constant Cram Thus, the number of neurons in HEAP is bounded by Crum K(). which of O(
The corresponding HEAP updates are:
x Clamp (x* Tx(2 W + -c- ATx*w()))
y = RReLU(Y* + (b - A( x²*
B. PDE Solutions under Different Constraints
+ +
(12)
As discussed in Section 2.3, the solution to the multi-asset Black-Scholer equation may not always be unique in prac- tice. Constraints can be employed to obtain the desired solution. In Figure 3, we compare the convergence of the solution under different constraint conditions (using the Black-Schole equation from Equation 20, with con- straints and with constraint T2 at varying values of An.).. The comparison reveals that, even after excluding the in- fluence of solution residuals, the equation yields different solutions in certain regions.
By setting / = 3 = I, and ensuring that Tki 0k and nn in HEAP match the corresponding step sizes in APDQP, it can be verified that:
x k+1 and ykh+1 = +++
(13)
Thus, by induction, HEAP can exactly replicate the primal- dual sequence generated by the APDQP algorithm, indicat- ing HEAP's expressivity.
D
Proof of Theorem 3.2. Firstly, we introduce a proposition that describes A-PDHG's convergence rate.
C. The Converge of Learnable Step Sizes
As outlined in Algorithm 2, the A-PDHG method adaptively learns the update step sizes for the X and Y layers (r and n) to ensure effective training. Through experiments in Figture 6, we have verified that the update step sizes are indeed adjusted during training and converge to a fixed value after approximately 60 batches.
Proposition A.1. (Chambolle & Pock, 2016) Let (of.yy)) be the primal-dual variables generated by the PDHG algorithum for the QP problem M (Q, A, e. I, u). If the step sizes a,a satisfy To || 1, then for any (x. v) E Ra x R satisfying 1 I x u, the primal-dual gup satisfies
D. Training Trajectory
L(**, M) - L(x, :M)
The training trajectory is visualized in Figures 4 and 5, where we record the PDE residuals at training checkpoints. These figures show that HEAP converges faster and more stably than CSO. Our method's advantage in these bench- marks arises from its principled design and theoretical soundness, particularly through the elegant incorporation of constraints into the neural networks. In contrast, CSO incorporates constraints as a penalty term in the objective, which may lead to constraint violations on the test set.
1
xx- x01 212 +
1 - 1012
(y 1 10) A(x - =0) TATT
<
2k
T
(14)
" k and L(a, yr
where
t
gb =
=
1
Qx + x - y TAx b).
We can conclude from the proposition above that A-PDHG converges linearly under appropriate conditions. Hence, if we set T, a to satisfy ro| [Al < 1 in both APDQP and HEAP, according to Thm. 3. the dual series of
Although our method demonstrates improvements, we ac- knowledge that constrained high- nigh-dimensional PDEs remain

---

HEAP: Hyper Extended A-PDHG Operator for Constrained High-dim PDEs
Sel co
Sel. c1
Sel
128F Compasioon on "waiiing
C50
Res, cD
Res
Res c2
22-
een
Sol.co c0
Sol c1
Sol. c2
Figure 5. Burgers L2REpte comparison between CSO and our HEAP during the training. The CSO model oscillates and even collapses the test set.
Ros co
Bes c1
Res
Figure 3. Visualization of predict solution (Pred_u) and accumu- lated PDE residual (Res) of Heap and CSO in Black- Scholes Equa- tion (dim- 10) with no constraints (cO) and two different instances of constraints T2.
0.004
0.042
L2RE Comparison on Training
0.034
0.10
C50
16:3P
0.034
10.5
Yaining
M
0.00
300
8
Figure 4. Constrained Heat T.2 comparison hetween CSO and our HEAP during the training. Our methods quickly achieve stable (relative) loss over the test set.
74
60
100
320
challenging due their inherent nonlinearity. Learning methods still face the risk of overfitting to the training set, which highlights the need for further refinement. It is also worth noting that, to date, we have not extensively fine-tuned the hyperparameters for each specific case.
Figure 6. During the training of the 10-dimensional Reaction Dif- fusion equation, the iteration of the adaptive step size parameter (r and m) finally converges to a certain value.
12


Reviews:
[{"review_id": "6TAjW7rFop", "reviewer": "Reviewer_yj2S", "summary": "The paper introduces a novel method called HEAP, designed to solve high-dimensional PDEs that include additional constraints. The PDE solution is approximated with the evolution of the neural network parameters, which is formulated as a quadratic programming problem. To solve the QP efficiently, the method unrolls a fixed number of iterations of the APDHG algorithm. This makes the iterative solver a differentiable module that can be trained end-to-end. A hypernetwork is incorporated to estimate initial values, step sizes, and latent weights needed for the APDHG iterations. \n\nThe paper presents theoretical results demonstrating that HEAP can replicate the iterative sequence of the APDHG algorithm and achieves linear convergence under suitable conditions. Experiments on several PDEs (including constrained and unconstrained heat equations, Burgers equation, and reaction-diffusion equations) demonstrate that HEAP achieves lower PDE residuals and better satisfaction of constraints compared to CSO.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thanks for your suggestions. The additional experiment results and visualizations are reported in the rebuttal supplementary material (RSM), available at <https://anonymous.4open.science/r/HEAP/RSM.pdf>, and will be added to the draft once allowed. Here we provide responses to your questions point by point:\n\n> Q1: A more detailed explanation of the experiment problems and possible visualization.\n\nA1: Here, we briefly explain two of the existing experiment problems and a newly added problem.\n\n1) Heat equation: The heat equation serves as a fundamental prototype PDEs in physics and mathematics. As a paradigmatic parabolic PDE, it governs ubiquitous physical phenomena ranging from thermal diffusion in materials to probability distribution evolution in stochastic processes. However, solving high-dim instances presents formidable challenges: both traditional mesh-based discretization techniques and neural network methods suffer from the \"curse of dimensionality\", with storage and computational costs growing exponentially as dimension increases.\n\n2) Heat equation with constraints: The incorporation of monotonically decreasing temperature fields introduces a physically grounded regularization to the classical heat equation, elevating both its modeling fidelity and computational complexity. This constraint encodes the thermodynamic irreversibility of cooling processes in materials with latent heat barriers. Traditional numerical time-stepping schemes may violate the monotonicity condition unless rigorously coupled with projection operators or barrier methods.\n\n3) Black-Sholes (newly added, with visualization): The Black-Scholes equation (Nobel prize 1997) is the cornerstone of modern quantitative finance, providing a rigorous framework for pricing financial derivatives under idealized market conditions. In practice, the no-arbitrage principle ensures the absence of trivial risk-free profit opportunities in the financial market. The no-arbitrage constraint arises naturally when modeling structured products with embedded downside protections or regulatory circuit breakers. The visualizations are Fig. 1 and 2 of RSM.\n\n---\n\n> Q2: A clearer explanation of how QP formulation and HEAP algorithm overcome limitations in existing approaches.\n\nA2: \n\n1) QP formulation: the existing approaches of high-dim PDEs formulate the parameter evolution either as a least square problem or a black-box optimization, without explicitly considering the potential constraints. The QP formulation allows the incorporation of linear constraints into the parameter evolution, extending the applicability of NN-based solvers to a broader range of high-dim PDEs.\n\n2) HEAP algorithm: the existing methods for solving large-scale QP problem are either iterative solvers like APDHG or vanilla neural operators. The former requires a large number (100+) of iterations to converge, while the latter suffers from the generalization issues. HEAP combines the best of both: it unrolls the APDHG algorithm into a neural network, which achieves better accuracy than vanilla neural operator due to algorithm priors and solves within 3 iterations, much fewer iterations than the APDHG due to learned knowledge from data.\n\n---\n\n> Q3: Intuitive explanations for concepts.\n\nA3: Intuitions of key components:\n1) Evolving operator: the solution of a high-dim PDE is surrogated by a subnetwork, whose parameters evolve over time to approximate PDE. When constraints are added, the parameter evolution is formulated as a QP problem, where the objective is to minimize the discrepancy between the left and right-hand side of PDE, and the linear constraints represent the given constraints.\n\n2) APDHG algorithm: the APDHG algorithm is a primal-dual iteration algorithm for QP, where the primal variable x is the solution variable and the dual variable y is the Lagrange multiplier or the penalty of constraints. The primal step updates x by a gradient descent of the objective function and penalty term, while the dual step updates y by a gradient ascent of the penalty term. Both updates are projected to the feasible set of the constraints.\n\n3) HEAP network: HEAP unrolls a fixed number of iterations of the APDHG algorithm into a neural network, where the primal and dual variables are extended from a vector to a matrix, i.e. parallelized over a new dimension. The extension is parameterized by the output of the hypernetwork, which takes the information of the problem state as input. The HEAP network can be trained end-to-end by backpropagation.\n\n4) Theorems: Theorem 3.1 shows HEAP network is truly an extended version of APDHG algorithm, since the HEAP network can be reduced to the APDHG algorithm when the parameters are set to ones. Theorem 3.2 shows that HEAP convergence of solving QP problems in terms of the required number of parameters. The proof naturally inherits from the convergence rate of the APDHG algorithm.", "claims_and_evidence": "1. Theoretical proof is provided for claims: HEAP aligns with the APDHG algorithm, the approximation capacity of HEAP. I have not fully checked the proof.\n\n2. Empirical results are provided to demonstrate that HEAP outperforms the baseline CSO in terms of lower PDE residuals and improved constraint satisfaction. The experimental evidence is clear, showing nearly consistent improvements over the baseline.", "methods_and_evaluation_criteria": "The proposed methods and the evaluation criteria are suited to the problem of solving high-dimensional constrained PDEs.", "theoretical_claims": "I have not fully checked the proof.", "experimental_designs_or_analyses": "The experimental design is sound for the problem at hand. However, a more detailed explanation of these problems and possible visualization can help the readers to understand the importance and difficulty of these problems.", "relation_to_broader_scientific_literature": "The paper is related to AI for PDE research. No apparent relevance to broader scientific literature.", "essential_references_not_discussed": "I can't think of apparent references that are missing.", "other_strengths_and_weaknesses": "The method’s motivation of reformulating the parameter evolution as a QP problem and unrolling the APDHG algorithm is interesting, but the reasoning behind these choices is not strongly conveyed. A clearer explanation of how these design choices overcome limitations in existing approaches would help readers, especially those from the general machine learning community.\n\nThe presentation is technical and dense, which makes it difficult for readers who are not already familiar with these methods to understand. The paper does not fully utilize its available space to offer intuitive explanations or visual aids that could help clarify the concepts.", "other_comments_or_suggestions": "I have no additional comments.", "questions_for_authors": "I don't have any further questions.", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "GFBmUfsh2y", "reviewer": "Reviewer_xiA3", "summary": "This work provides a principled way of handling constraints in the framework of Control-based solution operator (CSO) for learning PDE solution operators under constraints, called HEAP. In the original CSO, the evolution of the network parameter $\\theta$ is governed by a neural network $V$ called the neural control fields that directly map network parameter $\\theta$ to its time derivative $\\dot{\\theta}$. Constraints can be incorporated via soft penalties.\n\nThis paper instead handles the problem by performing constrained quadratic programming (QP). Assuming the constraints are linear, the CSO constraint optimization objective can be identified as a QP problem. This paper modifies an existing QP solver, adaptive primal-dual hybrid gradient (APDHG) for the CSO QP problem, where a matrices $\\{W\\}$ corresponding to the linear mapping to and from a linear latent space, the initial point of QP solver and the step sizes are predicted from a NN $V$. The output of $V$ is then used in the modified ADPHG, which runs for a fixed number of iterations and produces the solution, which is the time derivative $\\dot{\\theta}$. By identifying the modified ADPHG iterations as layers, $V$ can be seen as a hyper-network that predicts $\\dot{\\theta}$ from $\\theta$. The paper demonstrates the superior performance of HEAP over CSO in terms of accuracy, measured by L2 relative error, and constraint satisfaction on a set of constrained PDEs of dimension up to 20.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for the insightful questions and suggestions. The additional experiment results are reported in the rebuttal supplementary material (RSM), available at <https://anonymous.4open.science/r/HEAP/RSM.pdf>. Here we provide brief responses to your questions point by point:\n\n\n> Q1: Computational efficiency: memory usage and training/inference speed for d<=20 and higher dimensions scalability.\n\nA1: We provide results for d=5, 10, 15, 20, 50, 100 in the Tab. 4 of RSM. The HEAP costs 1.6x GPU memory and about 2x training/inference time compared to the baseline CSO. The time and memory usage are nearly linearly increasing with the dimension, thus HEAP is scalable to higher dimensions in terms of computation.\n\n---\n\n> Q2: Experimental design: parameters in the Gaussian distribution, hyperparameters for CSO, penalty coefficient in CSO\n\nA2: The parameters in the Gaussian distribution are mean=0, std=$\\sqrt{0.5}$ for training and mean=0 std=1 for testing. The penalty coefficients are PDE residual loss (with weight $w_1$), the constraint violation loss ($w_2$), and the numerical constraint loss ($w_3$). To balance their magnitudes during training, the weights are set as $w_1 = 10^{-2}$, $w_2 = 1.0$, and $w_3 = 10$. The remaining hyperparameters are displayed in Tab. 3 of RSM.\n\n---\n\n> Q3: Hyperparameter Ablation: iteration count K, extended dimension in the hypernetwork.\n\nA3: As shown in Tab. 5 of RSM, both iteration count K = 1,2,3,4 and extended dimension = 5,10,20 in the hypernetwork have a significant impact on the final PDE accuracy and constraint satisfaction. However, K=3 is a relatively robust choice across different PDE families.\n\n---\n\n\n> Q4: Handling more general, non-quadratic, or nonlinear constraints\n\nA4: The effectiveness of HEAP largely depends on its underlying APDHG algorithm and the structure of constraints. For convex but non-linear constraints, modifications for constraint projections may still preserve convergence guarantees, thus HEAP can be easily extended. For non-convex constraints, however, theoretical convergence becomes challenging, though heuristic adaptations like penalty methods or relaxed Lagrangian formulations might work in practice. \n\n---\n\n\n> Q5: Convergences of learnable step sizes (tau, sigma, etc.)\n\nA5: As shown in Fig. 3 of RSM, the learnable step sizes converge to stable values during training after 60 batches. \n\n---\n\n\n> Q6: Rationale for using ResNet\n\nA6: We follow the backbone choice of previous work, CSO, where ResNet is chosen for its simplicity and effectiveness in training deep networks. The skip connections in ResNet help to alleviate the vanishing gradient problem. Other architectures like Transformer or LSTM are also possible for HEAP, but we leave them for future work.\n\n---\n\n\n> Q7: Add related work STDE.\n\nA7: STDE is a state-of-the-art method for efficiently estimating extra-high-dimensional gradients, which can be integrated into HEAP as a gradient estimator. We will add a discussion and cite to STDE in the final version.\n\n---\n\n\n> Q8: Validity of the approximation by truncating the Taylor expansion of the constraints (Eqn. 4).\n\nA8: The original constraint is about the function (infinite dimension), which must be approximated by some finite-dimensional numerical scheme for computation. The truncating Taylor expansion is actually an implicit Euler method, in the sense that the constraint is enforced at the next time step. The implicit Euler method is a widely used stable scheme for PDEs, and the error can be controlled by the square of the step size.", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "The paper tested their algorithm HEAP again the baseline method CSO on (1) constrained heat equations, (2) unconstrained heat, (3) Burgers, (4) reaction–diffusion PDEs in 5, 10, 15, and 20 spatial dimensions, where training set is generated by sampling the initial condition $\\theta_{0}$ from a Gaussian distribution, and the test set is generated by sampling $\\theta_{0}$ from another Gaussian distribution.\n\nThe main evaluation criteria are:\n- PDE residual error (L2REpde) after a learned final solution,\n- Constraint violation (L2REcon) for physically constrained PDEs.\n\nI think the PDEs chosen and the evaluation criteria are valid, however, the dimension is only up to $20$, which still seems relatively low, and the authors did not report any metric on the computational efficiency like memory usage and training/inference speed. This is crucial since the improvements in L2RE for some cases are marginal, and if the proposed method incurs heavy computational cost, then the usefulness of the proposed method would be discounted.", "theoretical_claims": "Yes, I went through all theorems in section 3.3. except for proposition 3.3. which comes from another paper. They all seem to be coherent.", "experimental_designs_or_analyses": "- I find that there is too little information provided on the experimental design. What were the parameters in the Gaussian distribution used for generating the training and test set? What were the hyperparameters used for implementing the baseline method CSO? The author should, at the very least, report the penalty coefficient used in CSO, which could greatly impact the performance.\n- As mentioned before, the authors did not report any metric on the computational efficiency, like memory usage and training/inference speed.", "relation_to_broader_scientific_literature": "The paper references neural PDE-operator methods (DeepONet, PINOs, etc.) and high-dimensional PDE solvers (Han et al., Yu et al., among others).\nThey specifically build upon the recent line of \"parametric PDE operator\" approaches (like the CSO or Neural Control of PDE solutions).\nFor constrained PDE training, they mention relevant works on soft/hard constraints in PDE/inverse design but highlight how those typically remain in lower dimensions.", "essential_references_not_discussed": "None that I can think of.", "other_strengths_and_weaknesses": "Strengths:\n\n- The QP formulation + APDHG unrolling is a more principled way to enforce constraints than naive penalty-based methods used in CSO, and the experiment shows that the proposed method does outperform the baseline method in terms of accuracy.\n\nWeaknesses / Questions:\n\n- As mentioned, there is too little information provided on the experimental design, which greatly reduces the validity of the authors' claim.\n- As mentioned, the authors do not provide explicit complexity or memory scaling analysis.\n- The authors mention that a few iterations (K=3) suffice, but it remains unclear how robust this is across PDE families with more complicated constraints or stiffer PDE dynamics.\n- The method only works for linear constraints. It is unclear how to extend the QP formulation to nonlinear constraints.\n- Even for linear constraint, the proposed method uses an approximation by truncating the Taylor expansion of $u_{\\theta_{t}}$ (Eqn. 4), which modifies the original constraint. There was no discussion on how this approximation is valid throughout the paper.", "other_comments_or_suggestions": "1. The paper STDE (https://openreview.net/forum?id=J2wI2rCG2u) from NeurIPS 2024 seems to be the current state-of-the-art method for solving high-dimensional with PINN, which you could include in section 2.2. for a complete reference list.", "questions_for_authors": "1. Nonlinear Constraints: How easily can HEAP be extended to handle more general, non-quadratic, or nonlinear constraints, which might no longer yield a strict QP form in parameter space?\n2. Scalability for d>20: Have you tested or estimated memory/time usage for truly large dimensions, e.g., d=50+ or 100, to see if the cost of forming Q, A becomes limiting?\n3. Hyperparameter Ablation: How does the choice of iteration count K or extended dimension in the hypernetwork impact final PDE accuracy or constraint satisfaction?\n4. Adaptive Step Sizing: The approach includes some learnable step sizes ($\\tau$, $\\sigma$, etc.). Do you observe they converge to stable values, or do they exhibit large variance across training?\n5. What is the rationale for using ResNet?", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "TE4rNfTLUr", "reviewer": "Reviewer_fun5", "summary": "The paper introduces HEAP (Hyper Extended Adaptive PDHG), a new neural operator designed to solve constrained high-dimensional PDEs, where solutions must meet additional constraints beyond the governing equations. HEAP learns the evolution of PDE parameters and formulates this process as a quadratic programming (QP) problem. To efficiently solve this, the method unrolls the adaptive primal-dual hybrid gradient (APDHG) algorithm into the neural network. This approach improves efficiency while ensuring theoretical guarantees for constrained optimization. Experiments on various high-dimensional PDEs show that HEAP outperforms existing neural operators in accuracy and efficiency.", "questions": "", "limitations": "", "rating": 2, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for the insightful question. The additional experiment results are reported in the rebuttal supplementary material (RSM), available at <https://anonymous.4open.science/r/HEAP/RSM.pdf>. Regarding your question, we provide the following explanations:\n\n> Q1: What does constraint mean here, if PDE has a fixed solution?\n\nA1: The purpose of constraints, both equality and inequality, is to regularize the solution space and select the meaningful PDE solutions when multiple solutions exist.\n\nFirstly, in a general sense, the IC/BCs are also constraints of PDE, albeit equality constraints. \n\nSecondly, in practice, proper IC/BCs are sometimes not easy to obtain, while inequality is usually more accessible. For example, under which IC/BCs the Navier-Stokes equation has a unique solution is still an open problem (the Millennium Prize Problems) [1]. In this case, the inequality constraints are added to ensure a physically reasonable solution, e.g., the upper-bounded total energy (Eq. 7 in [1]).  Constraints do affect the solution of PDEs, as we show in Fig. 2 of RSM, where the solution of the Black-Scholes (BS) equation with different financial constraints is compared. The BS equation is a fundamental PDE in quantitative finance (Nobel prize 1997) [2], which does not have a unique solution in practice.\n\nThirdly, the constrained PDE problems have been formulated and studied in previous AI4PDE works [3, 4]. However, the existing methods are not scalable to high-dimensional problems, which is the focus of our work.\n\n---\n\n> Q2: Discuss the Evolutional deep neural network (EDNN) paper\n\nA2:\nThe EDNN also formulates the evolution of surrogate model parameters as an optimization problem, but it 1) solves the problem by a numerical solver instead of a neural network, 2) considers only equality constraints, and 3) is not directly scalable to high-dimensional problems. We will add a comparison with EDNN in the final version.\n\n---\n\n> Q3: Is there any more complex problem? Are all the testing cases with exact or trusted numerical solutions?\n\nA3:\nWe provide four additional experiments, two of which are Black-Scholes with no-arbitrage constraints (Eq. 1-5 of RSM), and the other two are burgers and reaction-diffusion with rate-of-change constraints. The results are shown in Tab. 1,2,6,7 of RSM. All the additional testing cases are without exact solutions or trusted numerical solutions. \n\n\nWe have added all the above discussions and RSM to the draft, but we are not able to provide the updated draft here due to the rebuttal rules this year. We hope this response addresses your concerns, and please let us know if you have any further questions.\n\n\nReference:\n\n[1] Fefferman, Charles L. \"Existence and smoothness of the Navier-Stokes equation.\" The millennium prize problems 57.67 (2006): 22.\n\n[2] Hull, John C., and Sankarshan Basu. Options, futures, and other derivatives. Pearson Education India, 2016.\n\n[3] Hoshisashi, Kentaro, Carolyn E. Phelan, and Paolo Barucca. \"Physics-Informed Neural Networks for Derivative-Constrained PDEs.\" ICML 2024 AI for Science Workshop.\n\n[4] Moro, Viggo, and Luiz FO Chamon. \"Solving Differential Equations with Constrained Learning.\" The Thirteenth International Conference on Learning Representations.", "claims_and_evidence": "Why constraint is so important for solving PDE? If the PDE and proper initial and boundary conditions are given, the solution is fixed, what does constraint mean here? Are you doing optimization? Could you make the motivation and problem setting more clear?", "methods_and_evaluation_criteria": "\"The reference solution is obtained by either the explicit solution or the PINN-based numerical solver\", since the method solves equations, we should only compare with high-order numerical solutions or exact solutions. Are the authors doing so?", "theoretical_claims": "3.1 to 3.3 look fine to me.", "experimental_designs_or_analyses": "Is there any more complex problem? Are all the testing cases with exact or trusted numerical solutions?", "relation_to_broader_scientific_literature": "The contribution to me is unclear, it needs to be emphasized that what do the constraints mean? Are they residuals? If so, aren't we solving them?", "essential_references_not_discussed": "'Evolutional deep neural network', need to be cited and discussed.", "other_strengths_and_weaknesses": "I am not sure the motivation of the paper, what does the constraint mean here? Why we need it for solving PDE instead of optimization?", "other_comments_or_suggestions": "No such.", "questions_for_authors": "What does the constraint mean here? Why we need it for solving PDE instead of optimization?", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "tAPXHAKGQ2", "reviewer": "Reviewer_m3Db", "summary": "This paper focuses on solving high-dimensional time dependent PDE with constraints. Traditional method suffers from curse of dimensionality. The proposed method effectively solves the problem by combining quadratic programing (QP) with NeuralODE. The approach can be summarized as:\n1. Reformulating the PDE with constraints into a QP problem;\n2. Modifying existing QP solver APDHG with learnable parameters to improve efficiency, named the method as HEAP;\n3. Training solution ansatz parameters and HEAP parameters with NeuralODE method in a differentiable way by minimizing residuals.\n\n## update after rebuttal\n\nI had a specific question on how to measure performance on BS equations without explicit solution, but the authors did not reply. This is my main concern. Hence, I change my score to 3.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We appreciate your insightful comments and suggestions. Here are our responses to your questions:\n\n> Q1: It would strengthen the effectiveness of HEAP if some other constraints were demonstrated\n\nA1: The additional experiment results are reported in the rebuttal supplementary material (RSM), available at <https://anonymous.4open.science/r/HEAP/RSM.pdf>. We provide four additional experiments, two of which are Black-Scholes with no-arbitrage constraints (Eq. 1-5 of RSM), and the other two are burgers and reaction-diffusion with rate-of-change constraints. The results are shown in Tables 1, 2, 6, 7 of RSM.\n\n---\n\n> Q2: Does HEAP have some boundaries for applications? Are there any constraints that HEAP cannot handle?\n\nA2: HEAP is based on a PDHG algorithm variant, designed for convex constraints. Therefore, HEAP can handle any convex constraints in principle. For more general non-convex constraints, the convergence of HEAP is not guaranteed. \n\nWe will add this discussion and fix the typos in the final version once the paper update is allowed.", "claims_and_evidence": "Yes. The experiments show their method being effective in a clear and convincing way.", "methods_and_evaluation_criteria": "Yes. The experiments are benchmarked on heat equations", "theoretical_claims": "I checked both Theorem 3.1 & 3.2 and found no issue.", "experimental_designs_or_analyses": "Yes, I find no issue with the experimental design. Since there is little work on solving high-dim PDE with constraints, the authors only chose one baseline for demonstration, and I understand their choice.", "relation_to_broader_scientific_literature": "N.A.", "essential_references_not_discussed": "N.A.", "other_strengths_and_weaknesses": "**Originality** is high. There is little work for high-dim PDE with constraints for several reasons. Curse of dimensionality rules out traditional numerical methods. Meanwhile, constraints are not well handled by methods such as PINN or neural operators. Therefore, this work innovatively tackles the problem by combining traditional QP method with NeuralODE, which fills the blank to some degree.\n\n**Weakness** in my opinion is about the experiment. The authors chose the constraints as upper and lower bounds for $\\theta_t$. However, it seems the proposed method HEAP should not be limited to such simple constraints. It would strengthen the effectiveness of HEAP if some other constraints are demonstrated, for example, some typical linear constraints $H[u](x, t)\\geq 0$ as indicated in equation (1). Also, does HEAP have some boundary for application? Is there any constraints that HEAP cannot handle?\n\nTypo:\n1. Line 164, $\\Omega_A$ -> $\\Omega_H$?\n2. Line 127 (right), \"Constraint $A$\" -> Constraint $H$?", "other_comments_or_suggestions": "N.A.", "questions_for_authors": "1. As discuss in weakness, is there any linear constraints that HEAP cannot handle?", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Meta-review: The paper proposes HEAP, designed to solve high-dimensional PDEs in the presence of additional constraints. The proposed method proposes considerable performance along with a theoretical guarantee on the soundness of the method.
Decision: Accept (poster)