Bi-level Physics-Informed Neural Networks for PDE Constrained Optimization using Broyden's Hypergradients
Authors: Zhongkai Hao
Chengyang Ying
Hang Su
Jun Zhu
Jian Song
Ze Cheng
Venue: ICLR 2023 poster (2023)
Keywords: PINN
machine learning
bi-level optimization

Abstract: Deep learning based approaches like Physics-informed neural networks (PINNs) and DeepONets have shown promise on solving PDE constrained optimization (PDECO) problems. 
However, existing methods are insufficient to handle those PDE constraints that have a complicated or nonlinear dependency on optimization targets. In this paper, we present a novel bi-level optimization framework to resolve the challenge by decoupling the optimization of the targets and constraints. For the inner loop optimization, we adopt PINNs to solve the PDE constraints only. For the outer loop, we design a novel method by using Broyden's method based on the Implicit Function Theorem (IFT), which is efficient and accurate for approximating hypergradients. We further present theoretical explanations and error analysis of the hypergradients computation. Extensive experiments on multiple large-scale and nonlinear PDE constrained optimization problems demonstrate that our method achieves state-of-the-art results compared with strong baselines.

Published as a conference paper at ICLR 2023
BI-LEVEL PHYSICS-INFORMED NEURAL NETWORKS FOR PDE CONSTRAINED OPTIMIZATION USING BROYDEN'S HYPERGRADIENTS
Zhongkai Hao1.2.3, Chengyang Ying', Hang Su².4, Jun Zhu Jian Song", Ze Cheng² l Dept. of Comp. Sei. & Tech., Institute for AI, BNRist Center, THBI Lab, Tsinghua- Banghua-Boschh Joint Center for MI., Tsinghua University Dept. of EE, Tsinghua University.
3 RealAI, 'Pazhou Lab, Guangzhou, 510330, China, Bosch China Investment Ltd (h2 ycy21 @mail. t sighuua. edu. en, [des2], suhangss, isong @t stsiinhu edu. en, cheng@ www.:eccn..sccccoooo
ABSTRACT
Deep learning based approaches like Physics-informed neural networks (PINNs) and DeepONets have shown promise solving PDE constrained optimization (PDECO) problems. However, existing methods are insufficient to handle those PDE constraints that have complicated nonlinear dependency optimization targets. In this paper, we present a novel bi-level optimization framework to resolve the challenge by decoupling the optimization of the targets and constraints. For the inner loop optimization, we adopt PINNs to solve the PDE constraints only. For the outer loop, we design novel method by using Broyden's 's method based the Implicit Function Theorem (IFT), which is efficient and accurate for approximating hypergradients. We further present theoretical explanations and error analysis of the hypergradients computation. Extensive experiments multiple large-scale and nonlinear PDE constrained optimization problems demonstrate that our method achieves state-of- sate--of-theart the-art results compared with strong baselines.
INTRODUCTION
PDE constrained optimization (PDECO) aims at optimizing the performance of a physical system constrained by partial differential equations (PDEs) with desired properties. It is a fundamental task in numerous areas of science (Chakrabarty & Hanson, 2005; Ng & Dubljevic, 2012) and engineering (Hicks & Henne, 1978: Chen et al., 2009), with wide range of important applications including image denoising in computer vision (De los Reyes & Schönlieb, 2013), design of aircraft wings in aerodynamics (Hicks & Henne, 1978), and drug delivery (Chakrabarty Hanson, 2005) in biology etc. These problem have numerous inherent challenges due to the diversity and complexity of physical constraints and practical problems.
Traditional numerical methods like adjoint methods (Herzog & Kunisch, 2010) based on finite element methods (FEMs) (Zienkiewiez al., 2005) have been studied for decades. They could be divided into continuous and discretized adjoint methods (Mitusch et al., 2019). The former one requires complex handcraft derivation of adjoint PDEs and the latter one more flexible and more frequently used. However, the computational cost of FEMs grows quadratically to cubically (Xue et al., 2020) w.r.t mesh sizes. Thus compared with other constrained optimization problems, is much more expensive or even intractable to solve high dimensional PDECO problems with large search space or mesh size.
To mitigate this problem, neural network methods like DeepONet (Lu et al., 2019) have been proposed as surrogate models of FEMs recently. DeepONet learns mapping from control (decision) variables to solutions of PDEs and further replaces PDE constraints with the operator network. But these methods require pretraining a large operator network which is non-trivial and inefficient. Moreover,
"Corresponding author.

---

Published as a conference paper at ICLR 2023
its performance may deteriorate if the optimal solution is out of the training distribution (Lanthaler et al., 2022). Another approach of neural methods (Lu et al., 2021: Mowlavi & Nabi, 2021) proposes to use a single PINN (Raissi et al., 2019) to solve the PDECO problem instead of pretraining an operator network. It uses the method of Lagrangian multipliers to treat the PDE constraints as regularization terms, and thus optimize the objective and PDE loss simultaneously. However, such methods introduce a trade-off' between optimization targets and regularization terms (i.c., PDE losses) which is crucial for the performance (Nandwani et al., 2019). It is generally non-trivial to set proper weights for balancing these terms due to the lack of theoretical guidance. Existing heuristic approaches for selecting the weights may usually yield an unstable training process. Therefore, is imperative to develop effective strategy to handle PDE constraints for solving PDECO problems.
To address the aforementioned challenges, we propose a novel bi-level optimization framework named Bi-level Physics-informed informed Neural networks with Broyden's hypergradients (BPN) for solving PDE constrained optimization problems. Specifically, we first present a bi-level formulation of the PDECO problems, which decouples the optimization of the targets and PDE constraints, thereby naturally addressing the challenge of loss balancing in regularization based methods. To solve the bi-level optimization problem, develop an iterative method that optimizes PINNS with PDE constraints in the inner loop while optimizes the control variables for objective functions in the outer loop using hypergradients. In general, it is nontrivial to compute hypergradients in bi-level optimization for control variables especially if the inner loop optimization is complicated (Lorraine et al., 2020). To address this issue, we further propose novel strategy based on implicit differentiation using Broyden's 's method which is a scalable and efficient quasi-Newton method in practice (Kelley, 1995; Bai et al., 2020). We then theoretically prove an upper bound for the approximation of hypergradients under mild assumptions. Extensive experiments on several benchmark PDE constrained optimization problems show that our method is more effective and efficient compared with the alternative methods.
We summarize our contributions as follows:
To the best of our knowledge, it is the first attempt that solves general PDECO problems based on deep learning using a bi-level optimization framework that enjoys scalability and theoretical guarantee.
We propose a novel and efficient method for hypergradients computation using Broyden's method to solve the bi-level optimization.
We conduct extensive experiments and achieve state-of-the-art results among deep learning methods on several challenging PDECO problems with complex geometry or non- non-lincar Naiver-Stokes equations.
2 RELATED WORK
Neural Networks Approaches for PDE Constrained Optimization. Surrogate modeling is an important class of methods for PDE constrained optimization (Queipo et al., 2005). Physics-informed neural networks (PINNs) are powerful and flexible surrogates to represent the solutions of PDEs (Raissi et al., 2019). hPINN (Lu et al., 2021) treats PDE constraints as regularization terms and optimizes the control variables and states simultaneously. It uses the penalty method and the Lagrangian method to adjust the weights of multipliers. (Mowlavi & Nabi, 2021) also adopts the same formulation but uses line search to find the largest weight when the PDE error is within a certain range. The key limitation of these approaches is that heuristically choosing methods for tuning weights of multipliers might be sub-optimal and unstable. Another class of methods train an operator network from control variables to solutions of PDEs or objective functions. Several works (Xue et al., 2020; Sun et al., 2021: Beatson et al., 2020) use mesh-based methods and predict states on all mesh points from control variables at the same time. PI-DeepONet (Wang et al., 2021a;o) adopts the architecture of DeepONet (L.u et al., 2019) and trains the network using physics-informed losses (PDE losses). However, they produce unsatisfactory results if the optimal solution is out of the distribution (Lanthaler et al., 2022).
Bi-level Optimization in Machine Learning. Bi-level optimization is widely used in various machine learning tasks, e.g., neural architecture search (Liu et al., 2018), meta learning (Rajeswaran et al., 2019) and hyperparameters optimization (Lorraine et al., 2020; Bao al., 2021). One of the key challenges is to compute hypergradients with respect to the inner loop optimization (Liu et al., 2021). Some previous works (Maclaurin et al., 2015; Liu al., 2018) use unrolled optimization or
2

---

Published as a conference paper at ICLR 2023
truncated unrolled optimization which is to differentiate the optimization process. However, this is not scalable if the inner loop optimization is a complicated process. Some other works (Lorraine al., 2020; Clarke et al., 2021; Rajeswaran al., 2019) compute the hypergradient based on implicit function theorem. This requires computing of the inverse hessian-vector product (inverse-HVP), In (Lorraine et al., 2020) it proposes to neumann series approximate hypergradients. Some works also use the conjugated gradient method (Pedregosa, 2016). The approximation for implicit differentiation crucial for the accuracy of hypergradients computation (Grazzi et al., 2020).
3 METHODOLOGY
3.1 PRELIMINARIES
Let V, V be three Banach spaces. The solution fields of PDEs are called state variables, i.e., y E Yaa C Y, and functions or variables we can control are control variables, i.e., u E Uad U where Yad. add called admissible spaces, e.g. subspace parameterized by neural networks or finite element basis. The PDE constrained optimization can be formulated as:
min yc UE
I [y, u), s.t. e(y, u)
where J Y X is the objective function and e: Y x U + V are PDE constraints. Usually, the PDE system of e(y, u) contains multiple equations and boundaryinnitall conditions
((y, 0, Va E u)y, = 0, Va 82,
where (2 Rd) the differential operator representing PDEs and B: Y U (Ω + Rd2) represents boundary/initial conditions. Existing methods based on regularization (e.g. the penalty method) solve the PDÉCO problem by minimizing the following objective (Lu et al., 2021):
min J I [[w uo) + / [^1 F(ywo (x)( 10
/
1A2 B(y us) us) (z)
(3)
where the solutions y and control variables are respectively parameterized by w E Th and Th" with w being the weights of PINNs. A Rd are hyper-parameters balancing these terms of the optimization targets. One main difficulty is that X; are hard to set and the results are sensitive to them due to the complex nature of regularization terms (PDE constraints). In general, large A; makes difficult to optimize the objective J. while small X; can result in a nonphysical solution of yw- Besides, the optimal \; may also vary with the different phases of training.
3.2 REFORMULATING PDECO AS BI-LEVEL OPTIMIZATION
To resolve the above challenges of regularization based methods, we first present a new perspective that interprets PDECO as bi-level optimization problem (Liu al., 2021), which can facilitate new solver consequentially. Specifically, we solve the following bi-level optimization problem:
min 0
T (w"
s.l.
w arg min,, (w, 0).
In the outer loop, we only minimize with respect to 0 given the optimal value of w*, and optimize PDE losses using PINNs in the inner loop with the fixed 0. The objective & of the inner loop sub-problem
E E / [F un (x)
/ [B(mm, uo (dx
Jan
By transforming the problem in Eq. equation I into a bi-level optimization problem, the optimization of PDEs' state variables and control variables be decoupled, which relieves the headache of of X; in Eq. equation 3.
setting proper
To solve this bi-level optimization problem, design inner loops and outer loops that executed iteratively. As shown in Figure 1, we train PINNs in the inner loop with PDE losses in Eq. equation 5. In the outer loop, we compute the hypergradients based Implicit Function Differentiation inspired
3

---

Published as a conference paper at ICLR 2023
IFT Differentiation
PE
APg Duchi
v
öe
dw
Broyden Iterations
(Low Ranko Vector
min J (w"" ,
Matrix Product
Outer loop
s.t. "² arg min Eur, 0)
Finetume PINNs:
W <
Ow
Inner loop
Figure 1: Illustration of our bi-level optimization framework (BPN) for solving PDE constrained optimization problems. In each iteration, we compute hypergradients of control parameters using IFT Differentiation in the outer loop. We calculate inverse vector-Hessian product based on Broyden's method which uses low rank approximation for acceleration. Then in the inner loop, we fine-tune PINNs using PDE losses only.
by Lorraine et al. (2020). Along this line, we need to calculate highly complex inverse Hessian- Jacobian product. To address this issue, we propose use Broyden's method which provides an efficient approximation particular, we fine-tune the PINNs using the PDE losses each iteration of the inner loop. We further compute the gradients of J with respect to parameters of control variables 9 in the outer loop, which also recognized as hypergradients in bi-level optimization and detailed the following section.
superlinear convergence speed (Rodomanov & Nesterov, 2021). In
3.3 HYPERGRADIENTS COMPUTATION USING BROYDEN ITERATIONS
The upper-level objective J depends on the optimal w" of the lower level optimization, i.e.,
dj OJ
ÖJ Ow"
Ow 00
Thus we need consider the Jacobian of w" with respect to when calculating the hypergradients.
Since w minimizes the lower level problem, we can derive Function Theorem (Lorraine al., 2020) as,
Ow' de
by applying Cauchy Implicit
Proposition (Proof in Appendix B.I). Iffor some (w', 0'), the lower level optimization is solved, i.e. Du (ww ,0) and (wuw )-1 invertible, then there exists a function w² w" (0) surrounding (00)) (0').0') _ = 0 and we have:
8E (w', 0') s.t. Dw
Ow"
E2
OP TwoOT
00 le
dwowT
(w)())))))
By Proposition 1, we could compute the hypergradients analytically as
dJ
03
OJ
02E
OPE TwooT
(8)
d0
00
Ow
@wouu
(w-(()))))
However, computing the inverse Hessian matrix, i.e.,
a²g Twoh
is intractable for parameters
of neural networks. To handle this challenge, we can first compute is also called the inverse vector-Hessian product. Previous works (Lorraine et al., 2020) use the Neumann series to approximate z" with linear convergence speed. However, in practice this approach is usually coarse and imprecise estimation of the hypergradients (Grazzi et al.. 2020). Here employ more efficient and effective approach to compute the inverse oector-Hessiae product which enjoys superlinear convergence speed (Rodomanov & Nesterov, 2021). ). It equals finding the root
which
Dw
Fwow
4

---

Published as a conference paper at ICLR 2023
2" for the following linear equation
ô 50
a5
(x) Ow
Ow 0.
Ow
Note that for each evaluation of Ju (z), we only need to compute two Jacobian-vector products with low computational cost, which does not need to create a giant instance of Hessian matrix. Specifically,
we use low rank Broyden's s method (Broyden, 1965; Rodomanov & Nesterov, 2021) to iteratively
approximate the solution z² In each iteration, first approximate the inverse of 8Pg dwōr
E
RJ
= - B:
Uk
(10)
Owow
k=1
We update u, v and z according to the following rules,
2++1
- (zg((zi)
(11)
A2i+1 B,Agh+1
u4+1 (Az++1) B,Agi+1
(12)
0/11
B,A
(13)
where A241 Fi+1 Ri, Agi41 96+1 gi. and a is the step size (usually set to or using line search). In summary, the inversion of the Hessian could updated by
00 Owow
B B, Az44 B,Agi+1 (Az++1) B,Agi+1 A
(14)
After m iterations, we use 2m as the approximation of z², We store uy 22 in low rank matrices which displayed as two red thin matrices in Figure Since we use a low rank approximation of the inverse of Hessian matrix, we longer need to store the whole matrix and only need to record and ukk where = 1... K and K is a tunable parameter depending on the memory limit. We run Broyden iterations until the maximum iteration is reached or the error is within threshold.
Our method is named Bi-level Physics-inie informed Neural networks with Broyden's 's hypergradients (BPN), of which the pseudo code is outlined in Algorithm 1. Given PDECO problem, we initialize PINNs with random parameters Wp and a guess on control parameters fg. First, we train PINNs under initial control 6 for Nw epochs as warm up (influence of hyperparameters will be diseussed in Appendix C). Then, we compute the hypergradients for using Broyden's method and update with gradient descent. After that, fine-tune PINNs under for N, epochs. These steps iteratively executed until convergence.
Algorithm Bi-level
informed Neural networks with Broyden's 's hypergradients (BPN).
Input: PINNs wwe with parameters Wg. initial guess 6. loss functions E and J, warmup and finetune epochs No. Ny. learning rates E1, c2 for PINNs and respectively
Output:optimal control parameters 0
1: Train uhe under control Ou for Nw epochs
2: for = 1, Niter do
3: Compute hypergradients V0, using Broyden's 's method in Eq (9) and Eq (14).
4: 0++1 = 02 62 VO,-
5: Train u, under control f,41 for Ny epochs.
6: if converged then
7: Break the cycle.
8: end if
9:
end for
4 THEORETICAL ANALYSIS AND DISCUSSION
Error analysis of hypergradients approximation. We now present analyses of the hypergradients approximation using BPN. Due to the complexity of the problem and computational resources, it extremely difficult to calculate the hypergradients exactly (Grazzi al., 2020). Our BPN
5

---

Published as a conference paper at ICLR 2023
using Broyden's 's method based on Implicit Function Differentiation is also an approximation for hypergradients. Then the approximation error or the convergence rate for hypergradients is one of the key factors that determine the performance. Here show that the approximation error could be bounded under mild assumptions using Broyden's 's method. Inspired by the idea of (Grazzi et al., 2020; Ji al., 2021), we prove that the approximation error could be decomposed into two parts, the former error term is induced by the inner loop optimization and the latter term is caused by the linear solver. Moreover, since our hypergradients are approximated by solving the linear equations using Broyden's 's method, it enjoys superlinear convergence rate (Rodomanov & Nesterov, 2021), which is superior compared with other methods such as Neumann Series (Lorraine et al., 2020). Specifically, we state the assumptions to be satisfied as below.
Assumption 1. We denote V Ow
2 ae d1 du d2 de for all (w, 0).
ViE is invertible and | p.
ViE and V VI V2E are Lipschitz continuous with constants p and T respectively.
VIJ. V2J, E and V32E are Lipschitz continuous with constant L.
The inner loop optimization is solved by an iterative algorithm that satisfies 1| w|12 I Pt |||lz l2 and < Pt + 0 as t + 0o where is the number of iterations.
The above assumptions hold for most practical PDECO problems. For the linear equation in Eq. equation 9, we use m step Broyden's method to solve it. Denote = L and d2Jr, d2J are the computed and real hypergradients for 0 respectively, we have the following Theorem.
Theorem (Proof in Appendix B.2). If Assumption holds and the linear equation equation 9 is solved with Broyden's 's method, there exists constant M > O, such that the following inequality holds at iteration I,
d22 12 [1(+K++
PM 122
TM r
pe pellwo12 ll + Mr (1 (1-
+
(15)
mK
The theorem above provides a rigorous theoretical guarantee that the hypergradients are close to real hypergradients when the inner loop optimization, i.e., PINNs' training, could be solved. We observe that the convergence rate relies on the inner loop optimization pa which is usually locally linear (Jiao et al., 2022) for gradient descent and superlinear term introduced by Broyden's 's method (Rodomanov & Nesterov, 2021). The high convergence speed could nearly eliminate this term in a few iterations, so the main error is usually dominated by the first term which also verified in our experiments.
Connections to traditional adjoint methods. Our BPN and adjoint methods (Herzog & Kunisch, 2010) share a similar idea of solving a transposed linear system in Eq. equation 9 to reduce the computational cost. But bi-level optimization is more general framework compared with constrained optimization in adjoint methods. It allows more flexible design choices like replacing FEM solvers with PINNs. Additionally, Eq. (9) is different from the adjoint equation (Herzog & Kunisch, 2010). Our Eq (9) solves a system in the parameter space of neural networks but the traditional adjoint equation corresponds real field defined on admissible space U or meshes.
EXPERIMENTS
In this section, we conduct extensive experiments on practical PDE constrained optimization problems to show the effectiveness of our method.
5.1 EXPERIMENTAL SETUP AND EVALUATION PROTOCOL
Benchmark problems. We choose several classic and challenging PDE constrained optimization problems containing both boundary, domain and time distributed optimization for both linear and non-linear equations. More details of problems are listed in Appendix A.
(1) Poisson's Equations. We solve this classic problem on irregular domain which could be viewed as a prototype of the heat exchanger (Diersch et al., 201 1) that is widely used in multiple domains. We denote this problem as Poisson's 2d CG in Table 1.

---

Published as a conference paper at ICLR 2023
(2) Heat Equations. A dime-distributed control problem for two dimensional heat equation similar to (Wang et al., 2021a) (denoted by Heat 2d).
(3) Burgers Equations. A time-distributed control problem for burgers equations (Burgers 1d) which are used in nonlinear acoustics and gas dynamics.
(4)~(6) Navier-Stokes Equations. They are highly nonlinear equations characterizing flow of fluids. We solve one shape optimization problem similar to (Wang et al., 2021a) (denoted by NS Shape in Table 1) and two boundary control problems (denoted by NS 2inlets and NS backstep in Table 1) (Mowlavi & Nabi, 2021).
Baselines. To demonstrate the effectiveness and superiority of method, compare it with several neural methods for PDECO recently proposed.
(1)~(2) hPINN (Lu et al., 2021): It trains PINNs with PDE losses and objective functions jointly with adaptive weights as regularization term. There are two strategies, i.e. the penalty method and the augmented Lagragian method update weights and use hPINN-P and hPINN-A denote them respectively.
(3) PINN-LS (Mowlavi & Nabi, 2021): It also treats the objective function as regularization term with adaptive weights. And it uses line-search rules to find the maximum tolerated weights for the objective function.
(4) PI-DeepONet (Wang et al., 2021a): As two-stage method, it first trains a DeepONet using physics-infee informed losses and then optimizes the objective using the operator network.
Apart from the methods mentioned above, we also implement several other bi-level optimization algorithms and measure their performance on these PDECO tasks.
(5) Truncated Unrolled InterrentiaioonTRRMI (Shaban et al., 2019): It assumes that the inner loop optimization solved by several gradient updates and uses reverse mode backpropagation to calculate hypergradients.
(6) T1 T2 (Luketina et al., 2016): It computes the hypergradients using Implicit Differentiation Theorem and uses an identity matrix instead of calculating the exact inversion of Hessian matrix.
(7) Neumann Series (Neumann) (Lorraine et al., 2020): It proposed to use the Neumann series to iteratively approximates the inverse of Hessian and uses the vector-Jacobian product and the vector-Hessian product to avoid instantiation of Hessian matrices.
Beyond these baselines, we also choose adjoint method (Herzog & Kunisch, 2010) which is traditional PDECO method as a reference. We run adjoint methods based on high fidelity finite element solvers (Mitusch et al., 2019) to calculate the reference solutions of these tasks. For linear problems, it could be viewed as the ground truth solution. However, it is computationally expensive and sensitive to initial guesses since solves the whole system in each iteration.
Hyperparameters and Evaluation Protocol. We use multi- perceptrons (MLP) with width from 64~1288 and depth from ~ 5 for different problems, and train them using Adam optimizer (Kingma & Ba, 2014) with a learning rate 10-3 Since the accuracy of PINNs of regularization based methods could not be guaranteed, we resort to finite element methods (FEM) to evaluate the performance. The evaluation metric is the objective function for each problem which is specified in Appendix A. We save and interpolate the control variables and solve the system using FEM and calculate the objective function numerically in every validation epoch. Other details and method-specific hyperparameters are reported in Appendix F. We run experiments on single 2080 Ti GPU.
5.2 MAIN RESULTS
Based on the experimental results for all PDECO tasks in Table 1, we have the following observations. First, our BPN achieves state-of-the-arr performance compared with all baselines. It shows that bi-level optimization is an effective and scalable way for large-scale PDECO problems. Second, we observe that our BPN reaches nearly the global optimal result for linear equations (Poisson's 's and Heat equation). For these equations, the problem is convex and reference values provided by adjoint methods could be viewed as the ground truth. For non-linear problems, the reference values provide locally optimal solutions and our BPN sometimes outperforms adjoint methods. Third, the update strategy of loss weights is critical for regularization based methods like hPINN-P, hPINN-A and
7

---

Published as a conference paper at ICLR 2023
Objective (J)
Poisson's & 2d (2d)
Hent (1d)
Burgers (1d)
NS Shape (V)
NS Zinlets (ld)
NS hackstep (Id)
Initial guess
0,400
0.142
0.0974
1,78
0.0830
0.133
hPINN-P
0,373
0.132
0,0941
0.0732
0.0586
hPINN-A A
0,352
0.138
0,0895
0.0867
0.0696
PINN-1.
0.0634
0.0831
Pl-DeepONet
0,392
0,0501
0.0710
1.27
0.0850
0.0670
Ours
0.160
0.0379
0.0662
0.0369
0.0365
Reference
0,0378
0.0392
0.0382
Table I: Main results for performance comparison of different algorithms several PDECO tasks.
Lower score means better performance. We bold best results across baselines except from reference values. "_" means that this method cannot solve the problem. "2d"/"1d"'"V" means the control variable a 2d/1d function or a vector.
Objective (3)
Poisson's ' 2d CG
Heal 2d
Burgers ld
NS Shape
NS 2inlets
NS backstep
Initial guess
0.400
0.142
0.0974
1.77
0.0830
0.138
TRMD
0.224
0.0735
0.0671
0.0985
0.0391
TI 72
0.425
0.0392
0.0865
1.68
0.115
0.0879
Neumann
0.225
0.0782
0.0703
1.31
0.0962
0.0522
Broyden(Ours)
0.160
0.0379
0.0662
1.27
0.0369
0.0365
Table Performance comparison different strategies of computing hypergradients, Lower score means better performance. We bold the best results across all methods.
PINN-LS which limits their performance without theoretical guidelines. A possible reason is that the balance between the objective function and PDE losses is sensitive and unstable for complex PDEs.
5.3 COMPARISON WITH OTHER BI-LEVEL OPTIMIZATION STRATEGIES
We list the performance of all bi-level optimization methods in Table 2. First, we observe that our method using Broyden's 's method for hypergradients computation achieves the best results compared with other methods. This is a natural consequence since our method gives the most accurate approximation for the response gradients. Second, could see that all bi-level optimization methods are effective on solving PDECO problems. Third, the results are better if the hypergradients are more accurate in most cases. For example, Neumann series and TRMD use better approximation for inverse Hessian and they perform better compared with T1 T2.
5.4 EXPERIMENTS ON ITERATION EFFICIENCY
To show that our method is computationally efficient, we conduct efficiency experiments by comparing the iterations required. Since PI-DeepONet a two-stage method and FEMs are based on gradient descent, we only choose three regularization based methods as baselines. Note that for our BPN we count total inner loop iterations for fairness. We plot values of objective functions J in Figure 22. We found that our BPN is much more efficient with nearly linear convergence speed for linear problems. This shows that the hypergradients provided by bi-level optimization are more stable and effective compared with regularization based methods. Besides, we found that PINN-LS is more efficient compared with hPINNs. However, it is a common drawback that regularization based methods are not stable which might take lot of effort to find suitable loss weights.
5.5 FIDELITY OF HYPERGRADIENTS AND ABLATION STUDIES
Fidelity of hypergradients compared with other methods. In this experiment, we aim show how accurate the hypergradients are for these bi-level optimization methods. However, all tasks in our main experiments do not have closed form solution for response gradients. Here we conduct this experiment on toy task of one dimensional Poisson's 's equation with an analytical solution in Appendix A. We first compare our method with several other bi-level optimization methods. We between computed shown in Figure 3. Note that the data
2-y measure the cosine similarity, which is defined as | for any two vectors
hypergradients and analytical hypergradients and the results are collected in the first 75 outer iterations since all methods converge fast on this problem. The left part of the figure shows that Broyden's method gives the most accurate hypergradients with similarity
8

---

Published as a conference paper at ICLR 2023
DP%2000 WINN7 IPNNA
0.65
0.2
0.6
0.18
0.35
0.16
0.45
0.35
0.08
0.25
0.06
0.04
0.15
(min)
(mam)
Figure 2: Results of efficiency experiments on Poisson 2d CG (left) and Heat 2d (right) problem.
Gradient Similarity
Gradient Similarity
1.00
1.0
0.75
0.9
0.50
0.8
0.25
0.7
0.00
-0.25
0.6
0.50
0.5
TRMD
T1-T2
Neumann
Broyden
8
16
32
64
Method
Bruyden Iterations
Figure 3:
Cosine similarity of hypergradients for different methods (left) and number of iterations (right) on toy task of Poisson's 's equation.
close to 1. The other methods also provide positive approximation that helps to optimize the control variable. However, the stability Neumann series not good and some iterations it provides negative estimation of hypergradients.
Fidelity of hypergradients using different number of Broyden iterations. Since Broyden's method an iterative algorithm for solving the optimization problem, there a trade-off between efficiency and performance. Here we compare the fidelity hypergradients using different numbers of iterations for Broyden's 's method. We also use box plots for this experiment and the results are in the right part of Figure 3. We observe that as the number of Broyden iterations increases, the hypergradients become more accurate. And the median cosine similarity of hypergradients is more than 0.9 even we only use iterations which shows the efficiency of Broyden's 's method. We also find that the increment is minor after 16 iterations for this simple problem. This shows that with the number of iterations increases, Broyden's 's method converges fast and the rest the error is dominated by the term caused by inexactness of inner loop optimization in Theorem and Eq. equation 35.
We conduct more ablation studics including the impact of hyperparameters in Broyden's 's method and PINNs optimization in Appendix C, the comparison between our BPN and continuous adjoint method with adjoint PDE solved with PINNs in Appendix G. comparison of running time in Appendix G, and provide more visualization results for these PDECO problems in Appendix E.
CONCLUSIONS
In this paper, we proposed a novel bi-level optimization framework named Bi-level Physics-informed Neural networks (BPN) for solving PDE constrained optimization problems. We used PINNs for solving inner loop optimization and Broyden's 's method for computing hypergradients. Experiments on multiple PDECO tasks, including complex geometry or non-line coneiinear Haiver-Stokes equations, verified the effectiveness our method. As for potential negative social impact, the interpretability of PINNs not comparable with traditional numerical solvers, which is left for future work.

---

Published as a conference paper at ICLR 2023
7 REPRODUCIBILITY STATEMENT
We ensure the reproducibility our paper from three aspects. (1) Experiment: The implementation of experiment is described in Sec. 5.1. Ablation study for our experiments is in See. 5.5. Further details are Appendix A and Appendix C. (2) Code: Our code is included in supplementary materials, (3) Theory and Method: A complete proof of the theoretical results described is provided in Appendix B.
8 ETHICS STATEMENT
PDE constrained optimization has a wide range of real-world applications in science and engineering including physics. fluids dynamics, heat engineering and aerospace industry, etc. Our BPN is general framework for PDECO and thus might accelerate the development of these fields. The potential negative impact is that methods based neural networks like PINNs lack theoretical guarantee and interpretability. Accident investigation becomes more difficult if these unexplainable deployed in risk-sensitive areas. A possible solution mitigate this impact is develop more explainable and robust methods with better theoretical guidance or corner case protection when they are applied sessisens areas.
ACKNOWLEDGEMENT
This work supported by the National Key Research and Development Program of China (2020AA01006000, 2020AAA0106302, 2021YFB3201000), NSFC Projects (Nos. 62061 136001, 62076147, U19B2034, U1811461, U19A2081, 61972224), BÑRist (BNR2022RC01066) Tsinghua Institute for Guo Qiang. and the High Performance Computing Center, Tsinghua University. J.Z was also supported by the XPlorer Prize.
REFERENCES
Shaojie Bai, Vladlen Koltun, and J Zico Kolter. Multiscale deep equilibrium models. Advances in Neural Information Processing Systems, 33:5238-5250, 2020.
Fan Bao, Guoqiang Wu, Chongxuan Li, Jun Zhu, and Bo Zhang. Stability and generalization of bilevel programming in hyperparameter optimization. Advances in Neural Information Processing Systems, 34, 2021.
Alex Beatson, Jordan Ash, Geoffrey Roeder, Tianju Xue, and Ryan P Adams. Learning composable energy surrogates for pde order reduction. Advances in neural information processing systems, 33: 338-348, 2020.
Charles G Broyden. A class of methods for solving nonlinear simultaneous equations. Mathematics of computation, 1909)):57-5593 1965.
Siddhartha P Chakrabarty and Floyd B Hanson. Optimal control of drug delivery to brain tumors for a distributed parameters model. In Proceedings of the 2005, American Control Conference, 2005.,
pp. 973-978. IEFE, 2005.
Qun Chen, Moran Wang, Ning Pan, and Zeng- Yuan Guo. Optimization principles for convective heat transfer. Energy, 199-11066 2009.
Ross M Clarke, Elre Oldewage, and José Miguel Hermnndez-Lobate Scalable one-pas, optimisation of high-dimensional ueight-update hyperparameters by implicit differentiation. arXiv preprint arXiv:2110.10460 2021.
Juan Carlos De los Reyes and PioolaBBbbianea Schönlieb. Image denoising: learning the noise model via nonsmooth pde-constrained optimization. Inverse Problems Imaging, 7(4):1183, 2013.
H-JG Diersch, D Bauer, W Heidemann, Wolfram Rühaak, and Peter Schätzl. Finite element modeling of borehole heat exchanger systems: Part 2. numerical simulation. Computers & Geosciences, 37 (8):1136-114) 2011.
10

---

Published as a conference paper at ICLR 2023
Riccardo Grazzi, Luca Franceschi, Massimiliano Pontil, and Saverio Salzo. On the iteration complexity of hypergradient computation. In International Conference on Machine Learning, pp- 3748-3758 PMLR, 2020.
Roland Herzog and Karl Kunisch. Algorithms for pet- optimiiztionn GAMM-Mittillnnenn 33(2):163-176, 2010.
Raymond M Hicks and Preston Henne. Wing design by numerical optimization. Journal of Aircraft. 15(7):407-412, 1978.
Kaiyi Ji. Junjie Yang, and Yingbin Liang. Bilevel optimization: Convergence analysis and enhanced design. In International Conference on Machine Learning, pp. 4882 4892. PMLR, 2021.
Yuling Jiao, Yanming Lai, Dingwei Li, Xiliang Lu, Fengru Wang, Jerry Zhijian Yang, et al. A rate of convergence of physics informed neural networks for the linear second order elliptic pdes. Communications Computational Physics, (((4):1222-1255 2022.
Carl T Kelley. Iterative methods for linear and nonlinear equations. SIAM, 1995.
Diederik P Kingma and Jimmy Ba. Adam: method for stochastic optimization. arXiv preprint arXiv:1412.6980, 2014.
Samuel Lanthaler, Siddhartha Mishra, and George E Karniadakis. Error estimates for deeponets: A deep learning framework in infinite dimensions. Transactions of Mathematics and Its Applications, 6(1):taac001, 2022.
Hanxiao Liu, Karen Simonyan, and Yiming Yang. Darts: Differentiable architecture search, arXiv 2018.
preprint 1806
Risheng Liu, Jiaxin Gao, Jin Zhang, Deyu Meng, and Zhouchen Lin. Investigating bi-level optimization for learning and vision from a unified perspective: A survey and beyond. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2021.
Jonathan Lorraine, Paul Vicol, and David Duvenaud. Optimizing millions of hyperparameters by implicit differentiation. In International Conference Artificial Intelligence and Statistics, pp. 1540-1552. PMLR, 2020.
Lu Lu, Pengzhan Jin, and George Em Karniadakis. Deeponet: Learning nonlinear operators for identifying differential equations based on the universal approximation theorem of operators. arXiv preprint arxiv: 1910.0910000.. 10000339,, 2019.
Lu Lu, Raphael Pestourie, Wenjie Yao, Zhicheng Wang, Francesc Verdugo, and Steven Johnson. Physics-ifformed neural networks with hard constraints for inverse design. SIAM Journal Scientific Computing, 43(6): 101-11551111 B1122, 2021.
Jelena Luketina, Mathias Berglund, Klaus Greff, and Tapani Raiko. Scalable gradient- based tuning of continuous regularization hyperparameters, In International conference machine learning. pp. 2952-2960. PMLR, 2016.
Dougal Maclaurin, David Duvenaud, and Ryan Adams. Gradient-baad based hyperparameter optimization through reversible learning. In International conference machine learning, pp. 2113-2122.
PMLR, 2015.
Sebastian K Mitusch, Simon W Funke, and Jørgen S Dokken. dolfin-adjoint 2018. 1: automated adjoints for fenics and firedrake. Journal of Open Source Soflware, 4(38):1292, 2019.
Saviz Mowlavi and Saleh Nabi. Optimal control of pdes using physics-informed neural networks. arXiv preprint arXiv:2111.0988) 2021.
Yatin Nandwani, Abhishek Pathak, and Parag Singla. A primal dual formulation for deep learning with constraints. Advances in Neural Information Processing Systems, 32, 2019.
James Ng and Stevan Dubljevic. Optimal boundary control of diffusion-conn convection-reaction pde model with time-dependent spatial domain: Czochralski crystal growth process. Chemical engineering science, 1111)111-111,
2012. 2012.
11

---

Published as a conference paper at ICLR 2023
Fabian Pedregosa.
Hyperparameter optimization with approximate gradient. In International conference on machine learning, Pp. 737-746. PMLR, 2016.
Nestor V Queipo, Raphael T Haftka, Wei Shyy, Tushar Goel, Rajkumar Vaidyanathan, and P Kevin Tucker. Surrogate-based analysis and optimization. Progress in aerospace sciences, 41(1):1-28, 2005.
Maziar Raissi, Paris Perdikaris, and George 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.
Aravind Rajeswaran, Chelsea Finn, Sham M Kakade, and Sergey Levine. Meta-learning with implicit gradients. Advances in neural information processing systems, 32, 2019.
Anton Rodomanov and Yurii Nesterov. Greedy quasi-newton methods with explicit superlinear convergence. SIAM Journal on Optimization, 31(1):755-811,00 2021.
Amirreza Shaban, Ching-An Cheng, Nathan Hatch, and Byron Boots. Truncated back-propagation for bilevel optimization. In The 22nd International Conference Artificial Intelligence and Statistics, pp. 1723-1732. PMLR, 2019.
Xingyuan Sun, Tianju Xue, Szymon Rusinkiewicz, and Ryan P Adams. Amortized synthesis of constrained configurations using a differentiable surrogate. Advances in Neural Information Processing Systems, 34, 2021.
Sifan Wang, Mohamed Aziz Bhouri, and Paris Perdikaris. Fast ple-constrained optimization via sulf-supervised operator leamiing arXiv preprint arXiv:2110.13290 2021a.
Sifan Wang, Yujun Teng, and Paris Perdikaris. Understanding and mitigating gradient flow pathologies in physics-informed neural networks. SIAM Journal Scientific Computing, 4365):43055-3381, 2021b.
Sifan Wang, Hanwen Wang, and Paris Perdikaris. Learning the solution operator of parametric partial differential equations with physics-informed deeponets. Science advances, 7(40):eabi8605, 2021c.
Tianju Xue, Alex Beatson, Sigrid Adriaenssens, and Ryan Adams, Amortized finite element analysis for fast pde-constrained optimization. In International Conference on Machine Learning, pp. 10638-10641 PMLR, 2020.
Olek C Zienkiewiez, Robert I. Taylor, and Jian Z Zhu. The finite element method: its basis and fundamentals. Elsevier, 2005.
12

---

Published as a conference paper at ICLR 2023
A DETAILS OF PDECO TASKS
In this section, we give a mathematical description of these PDECO tasks in detail.
A.1 2D POISSON'S 'S EQUATION WITH COMPLEX GEOMETRY (POTSSON'S 2D CG)
In this problem, optimize wwo-dimensional Poisson's equation defined complex geometry which is a prototype of heat exchanger (Diersch et al., 2011). We solve the equation on a rectangle area ! = [-4, 4 minus four circles 22 { y) (z: zz)² + (y y₂)² & rr} where æ Ma are 22,,, 0.8. We define a field variable u(a, y) E R on domain Q 52" U 29. The control source [2; y) IR distributed circle X { v) x² + y² 1.62), Our goal is to solve the following optimization problem,
[ / Iu 1/
min J
(16)
s.t. Au -flfx E x), Ω
u = 1, E as
0, an
We could visualize the geometry of this problem in Figure The source function is distributed in circle at the origin with radius of 1.6 (not displayed here).
A.2 TIME DISTRIBUTED CONTROL OF 2D HEAT EQUATION (HEAT 2D)
In this problem, our goal is to solve an optimal control task of a system governed by the heat equation. The temperature field uæx y, f) IR is defined on rectangle domain Ω [0, 1]² with time E [0. 2] and the control signal is a function f(() E depends on time but does not depend on spatial coordinates and Our goal is to make u close to target function ú(z, y, (),
1 /
min J = - /
2 [a
Iu i/2dx
(17)
Ou vAu = (z, E x [0.2]
8.
Ot
u(I, = 0, 02 x [0.2]
u(x, y,
(, EΩ
The coefficient is 0.001 and the target function is chosen as â 32x(1 - x)y(1 - sin nt. We choose f(t) = 0.1 as the initial guess for the problem. The solution at timestep t = 2.0 is shown in Figure 6.
A.3 TIME DISTRIBUTED CONTROL OF ID BURGERS EQUATION (BURGERS ID)
Burgers equation is a nonlinear PDE widely used various areas like applied mathematics, fluid mechanics, gas dynamics, traffic flow and nonlinear acousties. Here as an example which is a dissipative system. The field variable u(x, t) is defined on ! = 1, 1] and t E [0, 1]. The control variable f() depends only time. We aim solve the following optimization problem.
viscous Burgers equations
1 1u(x,1)- - 2/2
min = f
(18)
Ou
Du
vAu
fxrel
Øt
dx
u(I,t) t)
0, E 00
u(x, 0) sin(TTr)e 2n² -2n
The diffusion coefficient 0.01 and the target function is chosen
à (x- (x+
(19)
13

---

Published as a conference paper at ICLR 2023
Figure 4: Visualization of geometry shape of Poisson's 2d CG.
Figure Init guess solution of Poisson's ' 2d CG. (Simulated using high fidelity FEM)
61
608
Figure 6: Initial guess solution for heat2d problem at final time / = 2.0. The result is simulated by high fidelity FEM.
14

---

Published as a conference paper at ICLR 2023
15
24
e3
ea
24
Figure 7: Visualization of the initial guess solution for Burgers 1d problem. The horizontal axis is spatial coordinates x and the vertical axis time t.
which is symmetric wave. We use ftt) = 0 the initial guess for this problem. The visualization of the initial solution is shown the Figure 7
A.4 INLET FLOW CONTROL FOR 2D STEADY NAIVER-STOKES EQUATIONS (NS 2INLETS)
Naiver-Stokes equations are one of the most important equations in fluid mechanics, aerodynamics and applied mathematics which are notoriously difficult to solve due to the high non-linearity. In this problem, we solve NS equations in a pipeline and aim to find the best inlet flow distribution fyy to make the outlet flow as uniform as possible. The flow velocity field is = (u, v) and the pressure field is p and they are defined on a rectangle domain 00, 1.5) x [0, 1.0]. We have two inlets and two outlets and several walls for this domain,
ria ((0,)) : 0 0.5)
(20)
ria [(x,0) : 0.5 <
Dot
nout = ((x, <
r an MMaTTurruu urp "I
The whole problem follows,
min J
Iu 2/
(21)
f
8.1(u V)u
v Vp
Re
V.u 0
u (f(y), 0), (I, 3) E r
u (0,v2(z)), (()()))),,,))) E T Ur
(, E
= (x E r
The target function is parabolic function a(x, y) = 4y((1 - y) and the velocity field on the second inlet and outlet is Uz(x) 18(x - 0.5)(1 z). The Reynold number is set to 100 in this problem. We initialize the solution f(y) the same with the target function, iee. f(y) = 4y(1 - y). In this case, the outlet velocity and the whole velocity field is shown in the following pictures
15

---

Published as a conference paper at ICLR 2023
02
67-
0.5
0.4-
01
02
01
1.2
3.2
Figure 8: Visualization of the initial guess solution NS 2inlet. The top two images are velocity fieldd of the X and Y directions. The bottom one is the pressure field.
16

---

Published as a conference paper at ICLR 2023
2.00
1.75
1.50
0.50
0.5
0.00
Figure 9: Geometry shape and collocation points for NS 2d backstep problem.
A.5 DRAG MINIMIZATION OVER AN OBSTACLE OF NS EQUATIONS (NS SHAPE)
This problem is shape optimization task that is to find the best shape of the obstacle that minimizes the drag forces from the flow. The inlet the left side of the area and the outlet is the right side of the area,
r {(0,)) {8)
(22)
rou = ((8,)) : & y & 8
T" = ( r
The flow field is defined on a domain Ω = [0, 8 and Ω is the obstacle. The shape of the obstacle is ellipse parameterized by parameter 2.5). The goal is minimize the following objective.
1 JaY
Yu 2
(v
2
(vv
min J = f
+
+
+
dx
(23)
Ox
Oy
By
8.l.(u ( V)u
V2u Vp
Re
V.u O
(inn( Ty ).0) (z, 8) plu
0, (Ix y) r
= 0, (I, 3) E rou
A.6 INLET FLOW CONTROL FOR BACKSTEP FLOW (NS BACKSTEP)
The backstep flow is a classic example that might exhibits turbulent flow. In this problem, we also aim find the best inlet flow distribution f f2) make the outlet flow symmetric and uniform. The geometry of backstep could be viewed as a union of two rectangles Ω = [0, 1] [0, 0.5] (0, I]. The inlet is the left side of the area and the outlet is the right side of the area,
in
{(0.v) 0 0 & 0.5)
(24)
pout
((2.)) : & 1) ( r
r
17

---

Published as a conference paper at ICLR 2023
We going optimize the velocity fields of the outlet,
min J =
1u -i
(25)
Pqm
8. V)u
v Vp
Re
V u =
(f(y), 0), (x,,)) E pin
- 0, (z, E l
p (z, v) E Nout
The target velocity field is ((r v) 3y(1 y) and initialize the inlet velocity ('v) to be 8y(0.5 y) The Reynold number is also set to 100.
The visualization of the geometry of this problem is shown in Figure 9. A solution of the initial guess is shown in Figure 10.
12
025
03
625
03
0.0
025
Figure 10: Visualization of the initial guess solution for NS 2d backstep. The top two images velocity fields the X and Y directions. The bottom one is the pressure field.
A.7 DISTRIBUTED FORCE CONTROL OF 2D MAVIER-STOKEE EQUATIONS
In previous examples of Naiver-Stokes equations, the control variable is either 1d function (inlet flow) vector (shape parameters). In this problem, consider more challenging setting control a system governed by NS equation by optimizing distributed force (or source) terms over a domain. The dimension of the control variable (function) in this problem is 2. The geometry and boundary
18

---

Published as a conference paper at ICLR 2023
condition is the same with NS2inlet problem. The flow velocity field is (u, v) and the pressure field is and they are defined on a rectangle domain = [0, 1.5] x [0, 1.0]. We have two inlets and two outlets and several walls for this domain,
I {(0,y) :0 0.5]
(26)
2 [(z,0):
I {(1.5,9) y 0.5)
2 ((x,0.5) 0.5 <5 <
r ( Ur ur! u
Ω = (c, : v(x 0.5) + (y 0.5)3&00.2 < <0.2
The whole problem as follows,
/ Iu -lydy
min J
(27)
f
Jqqu
s.t.(u V)u
Vuu Up + F(z.y)lfze
Re
V
0
u = (f))).), (x,y) E E
(0, (2)))))) (x, E Din UT2
u =
O, (z,y) E T
P 0, (2,y) E T
A.8 POISSION'S ' ID Toy PROBLEM FOR EVALUATING HYPERGRADIENTS SIMILARITY.
This is simple problem with an analytic solution which is a boundary control for one-dimensional Poission's 's equation. We use this problem for evaluating the fidelity of hypergradients,
min J 0..0:
(y x² 22)dr
(28)
d'y 2
s.t. dz
v(0) do
y11 = 01
The exact solution of this Poission's 's equation is,
y 2 (0, Ou 1)x + Ou-
(29)
And the objective function J is quadratic form about On and 01,
J = (0) + + 0,0 - 20, 0o + 1).
(30)
The gradient for this function is
28 +1 -1 8g + 20, -
V
(31)
Thus we could compare the approximated hypergradients with the real gradients to study the effectiveness of different bi-level optimization strategies. The optimal solution of this problem
0" - ( (
and the corresponding objective function J = 0.
19

---

Published as a conference paper at ICLR 2023
B PROOF OF THEOREMS
In this section, we prove theorems the main text.
B.1 PROOF OF THEOREM
Theorem 1. 1 If for some (w', ,0), the lower level optimization is solved, i.e. (wo (w",,r and regularity conditions satisfied, then there exists a function u" w² (0) surrounding (w' 0') s.f. DE
(10))))) 0 and we have:
Ow"
00 OwowT
(32)
00
TwoAT (w'')0))0")
Proof. By Cauchy's Implicit Function Theorem, if f(w', ,0) = 0 is satisfied and its Jacobian Zu l or) is invertible then there exists an open set U including O' and continuous differentiable function w" satisfying that w" (6') = w' and ((w"(0), 0, 40 E U. Moreover, the partial derivatives could be computed by
Ow"
(01
05 Ow
0f
0B (0')
(w (00).0))
(33)
00
Here we take f(w', 0')
8e om
(w',0), and we get the following result
Ow" 02E
02 Two0T
(w- (0').0)))
(34)
00
dwdwT
Thus have proven the theorem.
B.2 PROOF OF THEOREM 1
Theorem 1. If assumptions hold and the linear equation is solved with Broyden method, there exists constant M >> O, such that the following inequality holds,
PM
TM
||d2J: - d2Jl12 S (L(1+K+
>pellwlla + Mx(1
(35)
1
Proof. First, since we assume the inner loop optimization is contraction map,
||w, w = pr|lw/| ||w|||
(36)
Then for all f. the weights are in compact set W = {w' ||w0|| < 2||w//]- Note that J is continous and differentiable, then there exists constant M and J is Lipschitz continous with M, i.e. for all W, E W we have
|VIJ & M.
(37)
We define the following notations to simplify the proof,
ViJi = VIJ 0)
(38)
V,E, VIE(Wi, 0)
V3E. = V2E(W,, 0)
2t VIJi (V7E))-
VIJ (Vî8)-1
and 2t.m is the approximation after m iterations of Broyden method. The emror of hypergradients are follows.
|[ - dJll =
|V2J + 2,m VI2&, - V2J - 12
(39)
|VJ, V2JÍ || V VIgE; vV 12
(40)
VVJr - V2JII + || V1 - 2, VIZE + 2, V VI2E -
(41)
IV2J. V2JII VI2E 2122 || VI2E, V
(42)
20

---

Published as a conference paper at ICLR 2023
Then if we could bound these terms respectively, the total error could be controlled. For the first and second term, it is caused by the inexactness of the inner loop optimization, since V2J is Lipschitz with constant L
V2J V2JII
L L| wil
Lpelluol-
(43)
And VIE, V2 is Lipschitz continuous with constant L.
|| 128 - 12
< L|| - =|1
LIIV (VVE)((VV-)))-))))))))()) VIJ (V78)-11
L (II (V/&)- VIJ. VIJI- (Vi8) (V/8) | (V78)-1 VIJ
VII.(VV8)-1) (1)8)-111
< IIVI [1(V))) -1
+ VIJII 1((78)
(44)
And these two terms could be bounded by,
V |[7E 138 |
1 wil
2 Pc ||1l
(45)
IV VIJII
L|lur wll Lpe||u|/-
(46)
Thus we have,
(DML
+
L2 + 1
V12E - 2122 &
Pelluil-
(47)
112
The last term is mixed term caused by both the approximation error when solving the linear system and the inexactness of the inner loop optimization,
kl. V2 - V2V12111 = 1 V V12E, - 2, 12E + 2,V V 12E, 2, 12
22mm VIZE, - 24 V1EE
2 V V - 24 12
V V12Ee1 || zell + || 122 - V12E||
1| =zlll + V (V( 1/ 12
LI - =1ll + M
E2 1
/
MT
L 2l
21|l +
|
wl
Mr
Zl - 24 ||
+
Pt
|||
/
since Broyden method enjoys a superlinear convergence rate (Rodomanov & Nesterov, 2021), we have the following convergence bound,
L 24.m
zell
x
L
||
mK
2 ||V (V7&)-11
1
L
mK
1
L
|V 11 [[(
mA
1
L
M
mR
1
Mr
(48)
mk
21

---

Published as a conference paper at ICLR 2023
Iters
8
16
32
64
Rank
0.0392
0.0390
0.0397
0.0382
16
0.0394
0.0379
0.0384
32
0.0379
0.0386
64
0.0379
Table 3: Ablation Study on Heat2d problem. We report the performance of BPN with different numbers of maximum ranks and iterations for Broyden method.
Thus we sum these error terms together and get,
|dJ. - dJIl
V2
2|ze V12 - V122
Vi2&
V
12
12)
T
MT
Lpe||w||+ +
+
pellwoll + Mr
pt||w||
42
MK
PM
TM
+
Pellu[l + MR
(49)
n1²
mR
Thus we have proven the main theorem for error analysis of the hypergradients computation using Broyden method.
C SUPPLEMENTARY ABLATION EXPERIMENTS AND DETAILS.
In this section, we provide the results of some ablation studies about our BPN and baselines. We choose Heat 2d and NS backstep as examples for these experiments.
C.I PERFORMANCE WITH THE DIFFERENT NUMBER OF BROYDEN ITERATIONS AND MAXIMUM MEMORY STEPS.
In this subsection, investigate the influence of different iterations and memory limit of Broyden method (Broyden, 1965). We conduct experiments on Heat 2d problem and NS backstep problem. Note that the maximum memory limit should be less than the maximum number of iterations.
The results of Heat 2d equation are shown in Table 3. We see that roughly the performance improves (the objective function decreases) with the increase maximum memory limit and maximum iterations. In this problem, the rank for approximating the inverse Hessian more than 16 is able to handle this problem. There is no significant gain we use more than iterations.
The results of NS backstep problem are shown in Table 4. The trend is similar to Heat 2d problem and roughly performance improves with the increase of Broyden iterations. However, we see that this problem is much more difficult and the performance deteriorates if we only use 8 memory steps even though we use many iterations.
C.2 PERFORMANCE WITH DIFFERENT NUMBERS OF FINETUNING BPOCHS.
In this subsection, we investigate the influence of different finetune epochs for the inner loop of PINNs, i.e. N, in Algorithm We choose Heat 2d problem for this experiment and the results are listed Table 5. We also plot the results Figure 11. We see that this hyperparameter has a considerable influence the convergence speed from about 0.2k 0.6k. However, its impact the final performance is minor within a broad range. We also observe that a small Ny might deteriorate final performance. Thus choosing moderate N, is efficient and effective.
22

---

Published as a conference paper at ICLR 2023
Iters
8
16
32
64
Rank
8
0.136
0.066
0.056
0.071
16
0.048
0.037
0.037
32
0.037
0.039
64
0.036
Table 4: Ablation Study on NS backstep problem. We report the performance of BPN with different numbers of maximum ranks and iterations for Broyden method.
Outer loop iterations
Finetune epochs
100
200
300
400
500
600
700
800
900
1000
32
0.0939
0.0856
0.0721
0.0563
0.0416
0.0384
0.0389
0.0400
0.0398
0.0410
64
0.0951
0.0845
0.0703
0.0574
0.0444
0.0400
0.0389
0.0384
0.0384
0.0379
128
0.0870
0.0697
0.0510
0.0430
0.0403
0.0389
0.0396
0.0392
0.0390
0.0386
256
0.0964
0.0924
0.0862
0.0715
0,0525
0.0428
0.0396
0.0380
0.0380
0.0379
512
0.0935
0.0836
0.0662
0,0477
0.0418
0.0401
0.0392
0.0387
0.0380
0,0380
Table 5: Results of ablation study the influence of finetuning epochs of PINNs.
C.3 PERFORMANCE WITH DIFFERENT NUMBERS OF WARMUP BPOCHS.
In this subsection, investigate the influence of different finetune epochs for the inner loop of PINNs, i.c. N in Algorithm 1. We also choose Heat 2d problem for this experiment. The results with different Na are shown in Table 6. The observations for this experiment ablation study on Ny. We see that convergence speed might be different for different optimization schedules of PINNs in the inner loop, but the final performance does not change lot. In summary,
similar to the
BPN is robust to the choice of these two parameters on Heat 2d problem.
D BPN FOR SHAPE OPTIMIZATION.
In this subsection, we show that BPN could also used for a special class of PDECO tasks which called shape optimization. The main challenge for this is the PDE losses & and optimization target I does not directly depend on the shape parameters, but they depend on the collocations points when sample from the region in practice. Thus calculating gradients for shape parameters is challenge for shape optimization tasks. Another viewpoint is that the shape parameters contribute to the integral domain but not the objective function. This means that the gradients with respect to 0 will not appear in the computational graph. Specifically, if J is also an integral over the domain, we have
3=/,0 Saco) (x)dd.
J J=
(50)
Outer loop iterations
Warmup epochs
100
200
300
400
500
600
700
800
900
1000
0.0935
0.0836
0.0662
0.0477
0,0418
0.0401
0.0392
0.0387
0.0384
0.038T
2000
0.0939
0.0856
0.0721
0.0563
0.0416
0.0384
0.0389
0.0387
0.0386
0.0380
6000
0.0951
0.0845
0.0703
0.0574
0.0444
0.0400
0.0389
0.0384
0,0384
0.0379
10000
0.0870
0.0697
0.0510
0.0430
0,0403
0.0389
0.0386
0.0382
0.0380
0.0385
20000
0.0964
0.0924
0.0862
0.0715
0.0525
0.0428
0.0396
0.0380
0.0381
0.0380
Table 6: Results of ablation study on the influence of warmup epochs of PINNs.
23

---

Published as a conference paper at ICLR 2023
32
50
64
2000
10"
128
6000
Objective
256
Objective
512
20000
10
Iterations (k)
Iterations (k)
Figure 11: Ablation study on the influence of finetuning steps (left) and warmup epochs (right) on heat 2d problem.
And the PDE losses could be written as (here we assume the control parameters only influence the geometric shape (2)))).
= / (z) +
(B)
(51)
(E())
Jan(o)
We need compute the gradients of 07 28g not contained in the computational graph. Inspired by the idea from traditional shape optimization, we have the following rules,
30 ThauT which rely on the gradients with respect 0 that is
aj
0J 0e
J(x) Ox 00
+
rdss
(52)
d0
(o
where n is normal direction for boundary 02(0). For the PDE loss, we have
3p a0
- On (IF(UW)))+++
Ox
(B(r)) (B)w))
-nds, On
(53)
Janee)
where ε is the culvature of the boundary shape for two dimensional problems. However, in order to compute the hypergradients, we need to compute the following rector-Hessiee product,
OJ
Ow
E2
Ow
30
SwOOT
i
30
(54)
Dw
OAT
need to compute the directional derivatives for a vector function 80 Direct computation of this term requires times of backpropagation which is extremely inefficient. Here small trick that only requires a constant computational cost by transforming this into a vector-Jacobian product,
JJ Ou
0E
00 )e
V
Ow-
a0
Ov
Ow
We could that ① (v
0E could be computed using single vector-Jacobian product and Bu Bun (u 39 ))) could then computed by two vector-Jacobian products. This trick is easily implemented and efficient in practice.
E VISUALIZATION OF OPTIMIZATION RESULTS
In this subsection, we show the results of the PDECO problems. We also compare our optimized result with the reference results solved by high fidelity adjoint solvers based on FEMs. To make sure the evaluation is fair and accurate, we discretize the control variables into meshes. Then load them into the FEM solver and interpolation the control variable using an interpolation curve or surface. We solve these systems using fidelity FEM solvers.
24

---

Published as a conference paper at ICLR 2023
E. POISSION'S 2D EQUATION WITH COMPLEX GEOMETRY.
We show the results of the control variable of Poisson's 2d CG problem solved using BPN and adjoint method (reference) in Figure 12.
24
60
4.2
Visualization optimization results of Poisson's 's 2d CG problem. The left picture is the control variable using BPN and the right picture is the control variable using adjoint method.
Figure 12:
And the solution of the PDEs under the optimal control variable shown in Figure 13.
25
35
45
Figure 13: Visualization optimization results of Poisson's 's 2d CG problem. The left picture is the solution field using BPN and the right picture is the solution field using adjoint method.
We plot the distance between final solution y and the target function in the Figure 14.
observe that the optimization result using BPN looks different with using adjoint solvers. Although the objective function is nearly the same in this problem, the control variables are different. The results of BPN are not symmetric and smoother compared with adjoint solvers. A possible reason is that randomness in training PINNs leads to breaking of symmetry. And the inductive biases of representing control variables using neural networks tend to find smoother solutions.
E.2 TEMPORAL CONTROL FOR HEAT EQUATIONS.
We show the results of the control variable of Heat 2d problem solved using BPN and adjoint method (reference) in Figure 15. The solution field at time 1.0 is shown Figure 16. And the distance to the target function is shown in Figure 17.
25

---

Published as a conference paper at ICLR 2023
40
17
455
Figure 14: Visualization of error field (u ú)² using BPN and adjoint method.
52
42
230
2.4
Figure 15: Visualization of control variables.
65
03
03
Figure 16: Visualization of solution fields at time f 1.0.
26

---

Published as a conference paper at ICLR 2023
4:5
07
04
633
0.2
Figure 17: Visualization of error (u û) time 1.0.
E.3 NON-LINEAR NAVIER-STOKER EQUATIONS.
Here we only choose the NS backstep problem as example for visualization of PDECO on Navier- Stokes equations. This problem more complex due to the non-linearity of Navier-Stokes equations. Both the results of our BPN and adjoint solvers might not be the ground truth solution. We visualize the velocity field u, v of the solution in Figure 18 and Figure 19. The pressure field is shown in 20. We that both methods allocate higher velocity the lower part of the inlet. But the details different since there are many degrees of freedom for this problem. It also possible there are many local minima for this problem.
14-
13
12-
05
02
06
04
44
02
Figure 18: Visualization of solution field u.
OTHER DETAILS AND HYPERPARAMETER, FOR BASELINES.
In this section, we specify some details and hyperparameters for baselines.
F.I HYPERPARAMETERS REGULARIZATION BASED METHIODS.
For these methods, we need to adjust the weights of the Lagrangian multipliers after several finetune epochs of PINNs. We choose the finetune steps from {1000, 2000, 4000, 6000}. For hPINN, the weights are set to 0.001 for initialization and they amplified by the ratio of 2 after each finetuning stage. The maximum weights are le3. For the line search strategy of PINN-LS, the maximum loss threshold is selected from (10-4, 10-3 10-2) for different problems.
27

---

Published as a conference paper at ICLR 2023
37
223
4.55
42
46
Figure 19: Visualization of solution field v.
14
14
24
12
12
010
0.6
66
3.55
004
022
602
-.00
6.04
Figure 20: Visualization of solution field p.
28

---

Published as a conference paper at ICLR 2023
F.2 HYPERPARAMETER. FOR PI-DEEPONET.
We use modified MLPs as network architectures for Pl-DeepONets proposed by (Wang et al.. 202lb;a). The number of layers is selected from 3~ 8 and the layer width is selected from {32, 64, 128, 256). We use a Gaussian random field for generating the training functions (Wang et al., 2021a). We generate another 1000 input functions also from the GRF as the test functions. We pretrain the PI-DeepONet until convergence. In the second stage, we replace the PDE constraints with the PI-DeepONet and optimize the control variables using the Adam optimizer with a learning rate of 10
F.3 DETAILS OF OTHER BI-LEVEL OPTIMIZATION METHODS.
TRMD. It uses truncated reverse-mode automatic differentiation for computing the hypergradients. If we use the SGD optimizer with L steps as the inner loop update algorithm,
00
W W/ O2 Ow;
(55)
We need compute the following Jacobian matrix,
Owi+1 Sw,
E.
0
(56)
00
00
Ow:000
Instead of create an instance of Jacobian matrix, we calculate the following Jacobian-vectoe vector product,
0J @wi+1
0J ow; DJ 03
(57)
dwL
00
dwr
a0 OwL
Tw;00T
We could iteratively calculate it with a truncated history of m steps. The steps chosen from (1,5, 10, 20) in experiments.
Neumann Series. It also approximates the hypergradients based on Implicit Function Differentiation. But the difference is that it uses the following rule to approximate the inverse Hessian matrix,
E2 Owowo
02e
k
(58)
TwouT T
k=U
+2 awowT T
(59)
And the inverse Hessian-Jacobian product is approximated by.
0J Ow
828 Owoww
L
aJ
(1-
828 Owow
(60)
Dw
A=
The convergence speed is decided by two hyperparameters and K. We choose from {10-5,
10-1 10-3, 10-2) and K from {1, 2, 4, 8, 16, 32). In theory, a should satisfy 10-5, 10
|l2. We also experimentally find that large a and large might lead divergence depending the problems. For too small it might degenerate to the method of T1 T2.
#2g Furtho
G SUPPLEMENTARY EXPERIMENTS
In this section, we provide supplementary experimental results to show the effectiveness of our method.
29

---

Published as a conference paper at ICLR 2023
DP%20v0 WINN7
0.65
0.6
0.1
0.35
0.1
0.45
0.35
0.08
0.25
0.06
0.04
(min)
(mai))
Figure 21: Results of efficiency experiments on Poisson 2d CG (left) and Heat 2d (right) problem. The horizontal axis represents the GPU clock time.
BPN((Our) hPINN-P hPINN-A
PINN-LS
0.16
0.14
0.45
0.35
0.3
0.25
0.2
0.15
Iteratious (k)
Tterations (k)
Figure 22: Results of efficiency experiments on Poisson 2d CG (left) and Heat 2d (right) problem. The horizontal axis represents total number of iterations.
G.1 EFFICIENCY EXPERIMENTS
In this subsection, we conduct efficiency experiments to show that our BPN is more efficient compared with baselines. We measure objective functions varying with both time and number of iterations for baselines on Poisson's 's CG problem and Heat 2d problem. The results are shown in the following two figures.
G.2 RUN TIME ANALYSIS AND ERROR BARS
Here report experimental results measuring error bars and wall clock runtime for these methods on POisson's 2d CG problem and Heat 2d problem.
We run experiments with random seeds and compute the mean and std of these results.
We measure all wall clock time for comparison. For neural networks based methods, we use single 2080 Ti GPU. The reference FEM are runned on laptop with threads of CPU (11-th Gen Intel Core i 1300 H @3.10 GHz ). The mesh size of reference FEM is roughly 500x500 for Poisson's 2d CG and 400x400x20 for Heat 2d. We omit P1-DeepONet for time comparison because it uses pre-trained operator network.
30

---

Published as a conference paper at ICLR 2023
Poisson's 8 2d CG
Heat 2d
Objectivemeenn/ss/tinee
mean
std
time(min)
mean
std
time
hPINN-P
0.368
0.025
0.149
0.015
12
hPINN-A
0.356
0.017
14
0.147
0.023
17
hPINN-LS
0.245
0.021
9
0.107
0.017
15
PI-DeepONet
0.376
0.035
0.0512
0.0012
T1-T2
0,412
0.031
8
0.0407
0.008
14
Neumann
0.217
0.019
15
0.0750
0.023
18
TRMD
0.239
0.041
28
0.0711
0.027
42
Ours
0.163
0.005
18
0.0382
0.004
19
Ref FEM
0.159
12
0.0378
15
Table 7: Experiments for run time analysis and enror bars.
Poisson's 's 2d CG
Heat 2d
Iterations 50
Objective (3)mean/std)
Total Time (min)(maan/std)
Ohjective (J)
Total Time (min)
0.232(0.017)
9
0. 0.2((00011)
13
100
0.225(0.013)
11
0. (06(0.016))
15
500
0.201(0.015)
15
0.0892(0.017)
17
1000
0.189(0.016)
20
0.0592(0.013)
2000
0.180(0.008)
28
0.0590(0.007 (000))
31
Broyden
0.163(0.005)
18
0.0382(0.004)
19
Table 8: Experimental results for comparsion with computing hypergradients with conjugate gradient method.
G.3 COMPARISON WITH CONJUGATE GRADIENT METHOD.
We compare the time cost and the performance on Poisson's 2d CG and Heat 2d problems using the conjugate gradient for computing hypergradients. The results shown in the following table. We use 50 iterations for Broyden's method.
We have two observations from these results. First, the performance using Broyden's method is better than using conjugated gradients, We found the main reason is that the hypergradients are more accurately approximated than using the conjugate gradient. Second, the performance of conjugated gradients improves with the increase in the number of iterations of CG. However, it is not efficient enough to use conjugated descent without preconditioning. But since we cannot store the whole hessian matrix, we are not able to use effective preconditioning techniques like incomplete Cholesky decomposition or Gauss-Seidel preconditioning. In summary, using a quasi-Netwon method like Broyden's 's method is more efficient for solving such a high dimensional and dense linear equation.
G.4 COMPARISON WITH ADJOINT SIMULATION WITH PINNS.
Adjoint PDEs for all cases in our main experiments, exist and they can be analytically derived. We also provide supplementary experiments comparing our algorithm with adjoint method simulated by using PINNs. We denote the method computes gradients from adjoint simulation of PDEs using PINNs as PINN-adjoint.
We found that PINN addoint is also effective for linear problems. However, by checking residual loss of PDEs and adjoint PDEs, we found that adjoint PDEs are more difficult to solve with PINNs, especially for complex cases like nonlinear Naiver-Stoker equations. In summary, that BPN still competitive and outperforms PINN-adjoint average.
Poisson's 's 2d CG
Heat
NS backstep
Objective
Time
Objective
Time
Objective
Time
PINN-adjoint
0.163
24min
0.0378
16min
0.0720
>30min
BPN (Ours)
0.160
18min
0.0379
19min
0.0365
55min
Table 9: Experimental results for comparsion with solving adjoint PDE with PINNs.
31

---

Published as a conference paper at ICLR 2023
(ethod(NS2inlet force control)
Objective
Initial guess
8.37e-2
hPINN-P
4.27e-2
hPINN-A
5.04e-2
PINN-LS
9.15e-3
TRMD
6.19e-3
T1-T2
5.22e-2
Neumann
7.20e-3
Ours
3.24e-3
Reference
3.17e-3
Table 10: Experimental results for NS2inlet force problem.
G.5 SUPPLEMENTARY EXPERIMENTS ON NS2INLET FORCE PROBLEM.
We have conducted a supplementary experiment on controlling a system governed by NS equation by optimizing a force (source) distributed over a 2d domain. The results are shown in the following Table 10.
G.6 DETAILS OF EVALUATION PROTOCOL.
Here specify more details of the evaluation protocol. For the main experiments, measure the performance using the final prediction when the algorithm stops rather than the best solution it founds. The motivation for this metric is that for practical problems, estimating the ground truth objective is expensive. For efficiency experiments, we measure their performance in each validation epoch to compare the efficiency for these methods.
32


Reviews:
[{"review_id": "Xge_eCokOKT", "reviewer": "Reviewer_yrSS", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "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": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "This paper focuses on an emerging and interesting research topic, Physics-informed neural networks (PINNs) for PDE constrained optimization. Different from existing works that use the regularization-based paradigm which is hard to set a proper weights to balance the optimization targets and regularization terms, this paper, for the first time, transforms the constrained optimization problem into a bi-level optimization problem. In addition, by leveraging Broyden’s method to approximate the IHVPs, the authors can more precisely estimate the hypergradient for the outer-loop hypergradient. Extensive experiments on several constrained PDE optimization problems demonstrate the effectiveness of the proposed reformulation.", "summary_of_the_review": "This paper is generally well-written and the authors focus on an interesting research problem. Although the novelty seems to be marginal, I'd like to accept this paper to encourage more attempts in the area of deep learning for science.", "strength_and_weaknesses": "Strength:\n1. The idea of the proposed reformulation is straightforward and extensive experiments verify the effectiveness of this succinct idea.\n2. The paper is well-written, and the main idea is easy to follow.\n3. The authors conduct a series of ablation studies, especially on the hypergradient approximation, to show the superiority of the recently-proposed Broyden’s method.\n\nWeaknesses:\n1. The idea is hard to say exciting, and the contributions can only be summarized as the simple reformulation and the adoption of the recently-proposed IHVPs method.\n2. The theoretical analysis mainly inherits from the cited paper [1] with less new findings.\n\n[1] Anton Rodomanov and Yurii Nesterov. Greedy quasi-newton methods with explicit superlinear convergence. SIAM Journal on Optimization, 31(1):785–811, 2021.", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}, {"review_id": "Nf4Zc5_vgc", "reviewer": "Reviewer_1sqe", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 5, "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": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "The paper proposed bi-level PINNs for solving PDE-constrained optimization problems. Bi-level PINNs have an inner loop optimization and outer loop optimization. In the inner loop, PINNs are used to solve PDE constraints, and in the outer loop, the control variables are optimized using Broyden’s hypergradients. Numerical experiments are performed to demonstrate the effectiveness of bi-level PINNs, and better results are achieved.", "summary_of_the_review": "New method is proposed in the paper, and the results look good. However, the tested problems are still simple, and more challenging problems should be tested.", "strength_and_weaknesses": "Strengths\n- The idea of bi-level optimization is not new, and the main contribution here is applying Broyden’s\nmethod to efficiently approximate the hypergradients in the outer loop. Also, an error bound is derived for the approximation error of the hypergradients, although this is relatively straightforward from existing results.\n- Several numerical examples are tested, and the method achieved good accuracy and speed.\n\nWeaknesses\n- Poisson and heat equations are simple, and NS is more challenging. However, in all the examples, the control variables are either scalers or 1D function, which are still simple. Problems with control variables of 2D functions should be tested. (Also, it is worth mentioning in Table 1 what type of the control variables is, scalers or 1D functions.)\n- The paper only chose deep learning methods as the baseline. Traditional methods should also be used for comparison.\n- In Introduction, “this is the first attempt that solves general PDECO problems using a bi-level optimization framework that enjoys scalability and theoretical guarantee” is not correct. There are traditional bi-level optimization methods for general PDECO problems. Probably this is first deep learning based bi-level optimization framework, but I am not sure.\n- In Introduction, “We conduct extensive experiments and achieve state-of-the-art results on several challenging PDECO problems with complex geometry or non-linear Naiver-Stokes equations.” “state-of-the-art results” is not correct. As shown in the paper, traditional adjoint method is used as the Reference method and is better than the proposed method in many cases.\n", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "Not applicable"}, "extra_text": {}}, {"review_id": "Caa5EhWOiex", "reviewer": "Reviewer_RMiZ", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "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": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "This paper formulates the PDE-constrained optimization problem as a bi-level optimization problem. The gradient of the bilevel problem (using implicit function theorem) involves an inversion of the Hessian. This paper use Broyden's hyper gradient to approximate the inverse of the Hessian.", "summary_of_the_review": "Overall this is a novel paper definitely.  But this paper lacks literature review of bi-level optimization and PDE constrained optimization at this point.", "strength_and_weaknesses": "Strength:\n- The bilevel formulation is interesting and novel.\n- The experiment is solid compared to other learning baselines.\n\nWeakness:\n- No discussion about algorithms used in bi-level optimization. The reviewer also has no idea of popularly used algorithm for bi-level optimization. I think it would be great to see popularly used bi-level algorithms as a baseline. PINN doesn't use large network. Inverse a Hessian will not be so expansive in my mind. \n- No baseline using traditional methods provided. I still doesn't know the benefit of using a NN for PDE in low dimension.\n- number of parameter and time of computation is not provided in this version.", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "4: The contributions are significant, and do not exist in prior works.", "empirical_novelty_and_significance": "4: The contributions are significant, and do not exist in prior works."}, "extra_text": {}}, {"review_id": "6Hm_9-YcayN", "reviewer": "Reviewer_isLd", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "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": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "This paper proposes a bi-level optimization method to solve constrained PDE. The key is to compute a hyper-gradient with a high-efficiency and accuracy. The paper has some unclear notations, making it hard for me to understand the main results. ", "summary_of_the_review": "REVISED: I now understand that the evaluation metric is used to show how well how control variable is found by the proposed method. The main results (table 1) show that it work better state-of-the art methods. Therefore I think the paper could be considered to be published. However, I still think that the evaluation on how well the pde is solved for a given control variable, should also be evaluated (maybe in a future work). For this reason, I raise my score to 6 . ", "strength_and_weaknesses": "The proposed method uses an interesting way to compute the hyper-gradients which is not so standard literature (i.e. to solve eq (8) by eq (9)), but this does not seem to be new. The strength of this paper is to show that the proposed method works well on a wide range of constrained PDE. ", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}]

Author rebuttal:

Decision: Accept: poster