An operator preconditioning perspective on training in physics-informed machine learning
Authors: Tim De Ryck
Florent Bonnet
Siddhartha Mishra
Emmanuel de Bezenac
Venue: ICLR 2024 poster (2024)
Keywords: physics-informed machine learning
operator preconditioning
deep learning
neural network training
Primary area: applications to physical sciences (physics, chemistry, biology, etc.)

Abstract: In this paper, we investigate the behavior of gradient descent algorithms in physics-informed machine learning methods like PINNs, which minimize residuals connected to partial differential equations (PDEs). Our key result is that the difficulty in training these models is closely related to the conditioning of a specific differential operator. This operator, in turn, is associated to the Hermitian square of the differential operator of the underlying PDE. If this operator is ill-conditioned, it results in slow or infeasible training. Therefore, preconditioning this operator is crucial. We employ both rigorous mathematical analysis and empirical evaluations to investigate various strategies, explaining how they better condition this critical operator, and consequently improve training.

Published as a conference paper at ICLR 2024
AN OPERATOR PRECONDITIONING PERSPECTIVE ON TRAINING IN PHYSICS-INFORMED MACHINE LEARNING
Tim De Ryck
Florent Bonnet
Seminar for Applied Mathematics,
Institute of Intelligent Systems and Robotics, Extrality,
ETH Zürich, Switzerland
Sorbonne Université, France
Siddhartha Mishra
Emmanuel de Bézenac'
Seminar for Applied Mathematics, ETH AI Center,
Seminar for Applied Mathematics, ETH Zürich. Switzerland
ETH Zürich, Switzerland
ABSTRACT
In this paper, we investigate the behavior of gradient descent algorithms in physics-infrrmed machine learning methods like PINNs, which minimize resid- uals connected to partial differential equations (PDEs). Our key result is that the difficulty in training these models closely related to the conditioning of a specific differential operator. This operator, in turn, is associated to the Hermi- tian square of the differential operator of the underlying PDE. this operator is ill-conditioned, it results in slow or infeasible training. Therefore, precondition- ing this operator is crucial. We employ both rigorous mathematical analysis and empirical evaluations investigate various strategies, explaining how they better condition this critical operator, and consequently improve training.
INTRODUCTION
Partial Differential Equations (PDEs) are ubiquitous as mathematical models of interesting phe- nomena in science and engineering (Evans 2010). Traditionally, numerical methods such finite difference, finite element etc (Quarteroni & Valli] 1994) are used simulate PDEs. However, given the prohibitive cost of these methods for variety of PDE problems such those with multiple scales, in high dimensions or involving multiple calls to the PDE solver like in UQ control and in- verse problems, machine learning based alternatives are finding increasing traction as efficient PDE simulators, see Karniadakis et al. (2021) and references therein.
Within the plethora of approaches that leverage machine learning techniques to solve PDEs, models which directly incorporate the underlying PDE into the loss function are widely popular. promi- example of this framework, often referred as physics-informed machine learning, are physics- informed neural networks or PINNs (Dissanayake & Phan- Thiin 1994 Lagaris et al.) 2200aa (b] Raissi al.) 2019). which minimize the PDE residual within the ansatz space of neural networks. Related approaches in which the weak or variational form of the PDE residual is minimized in- clude Deep Ritz (E & Yu 2018), neural Galerkin (Bruna et al. 2022 variational PINNs (Kharazmi et al. 2019) and weak PINNs (De Ryck et al. 2022). Similarly, PDE residual minimization methods other ansatz spaces such as Gaussian processes (Raissi & Karniadakis 2018). Fourier features (Tancik et al. 2020), random features (Ye et al. 2023) etc have also been considered.
Despite the considerable success of PINNs and their menioneed variants in solving numerous types of PDE forward and inverse problems (see Karniadakis et al. (202T); Cuomo et al. (2022) and references therein for extensive reviews), significant problems have been identified with physics- informed machine learning. Arguably, the foremost problem lies with the training these frameworks
Thhee authors contributed equally this work.

---

Published as a conference paper at ICLR 2024
with (variants of) gradient descent methods (Krishnapriyan al. 2021 Moseley et al.) 2021| Wang
et al. 2021a 2022b). It has been increasingly observed that PINNs Is and their variants are slow, even infeasible, to train even certain model problems (Krishnapriyan et al. (02T), with the training process either not converging or converging to unacceptably large loss values.
What is the reason behind the issues observed with training physics-informed machine learning algorithms? Empirical studies such as Krishnapriyan et al. (2021) attribute failure modes to the non-convex loss landscape, which is much more complex when compared to the loss landscape of supervised learning. Others like Moseley et al. (202)) Dolean et al.] (2023) have implicated the well-known spectral bias (Rahaman et al.) (2019) of neural networks as being a cause for poor training whereas Wang et al. (2021a by used infinite-wide
NTK theory to propose that the subtle balance between the PDE residual and supervised components of the loss function could explain and possibly ameliorate training issues. Nevertheless, it is fair to say that there is a paucity of principled analysis of the training process for gradient descent algorithms in the context of physics-informed machine learning. This provides the context for the current work where aim to rigorously analyze gradient descent based training in physics-informed machine learning, identify a potential cause of slow training and provide possible strategies to alleviate it. To this end, our main contributions are,
We derive precise conditions under which gradient descent for a informed loss function can be approximated by simplified gradient descent algorithm, which amounts to the gradient descent update for a linearized form of the training dynamics.
Consequently, we prove that the speed of convergence of the gradient descent is related to the condition number of an operator, which in turn is composed of the Hermitian square (D*D) of the differential operator (D) of the underlying PDE and kernel integral opera- tor, associated to the tangent kernel for the underlying model.
This analysis automatically suggests that preconditioning the resulting operator is neces- sary alleviate training issues for physics-informed machine learning.
By combination of rigorous analysis and empirical evaluation, we examine how differ- ent preconditioning strategies can overcome training bottlenecks and also investigate how existing techniques, proposed the literature for improving training, can be viewed from this new operator preconditioning perspective.
2 ANALYZING TRAINING FOR PHYSICS-IFFORMED MACHINE LEARNING IN TERMS OF OPERATOR CONDITIONING.
Setting. Our underlying PDE is the following abstract equation,
Du(x) = f(x), x x Ω u(x) = g(z), x x an.
(2.1)
Here, Q C IRd is an open bounded subset of either space or space-time, depending on whether the PDE depends on time or not. The PDE (2.I is specified in terms of the differential operator D and the boundary conditions given by g. Specific forms of the differential operator 'D are presented later on, whereas for simplicity, we fix Dirichlet-type boundary conditions in (2.1). while other types of boundary conditions can be similarly treated. Finally, we consider the solution u + IR as a scalar for simplicity although all the considerations below also for apply to the case of a vector
Physics-informed machine learning relies on an ansatz space of parametric functions, u(- Q IR for all 0 E Σ C R". This ansatz space could consist of linear (affine) combinations of basis functions Σx=
1 Oedx. with possible basis functions as trigonometric functions or finite-clement type piecewise polynomial functions it could consist of nonlinear parametric functions such as neural networks (Goodfellow et al. 2016) or Gaussian processes (Rasmussen) 2003).
The aim is to find parameters 0 E Σ such that the resulting parametric function u(= ()) = u, ap- proximates the solution u of the PDE (2.1). In contrast to supervised learning, where the parameters 0 would be chosen fit (possibly noisy) data u(z;) with s; D. the key ingredient in physics-
2

---

Published as a conference paper at ICLR 2024
informed machine learning is to consider the loss function
/ [((x) g(x)) da(x),
I(O)= = 2 (uu(x) - f(x)]² dx +
(2.2)
2
02
R(0)
B(O)
with PDE residual R, supervised loss B at the boundary and a parameter À > 0 that relatively weighs the two components of the loss function. In practice, the integrals in the loss function (2. need to replaced by suitable quadratures, but as long as the number of quadrature (sampling) points sufficiently large, the corresponding generalization errors (Mishra & Molinaro 2020 De Ryck & can be made arbitrarily small.
Mishra 2021
Characterization of Gradient Descent for Physics-informed Machine Learning. Physics- informed machine learning boils down to minimizing the physics-infommed loss 2.2 i... to find,
Ø argmin L(0). BEE
(2.3)
Once such an (approximate) minimizer At is obtained, one appeals to theoretical results such as those in Mishra & Molinaro (2020): De Ryck & Mishra (2021): De Ryck al. (202I) to show that u(- ; approximates the solution u of the PDE 2.1 to high accuracy. Moreover, explicit error estimates terms of the training error L(Bt) can also be obtained (Mishra Molinaro 2020)]
De Ryck & Mishra (021). As is customary in machine learning (Goodfellow et 201,)) the non- convex optimization problem (2.3) solved with (variants of) gradient descent algorithm which takes the following generic form,
B = 8 - 1) (gL()))
(2.4)
with descent steps > 0, learning rate v) >> 0, loss (2.2 and the initialization On chosen randomly.
Our aim here is to analyze whether this gradient descent algorithm (2.4) converges as k -> oo to a minimizer of (2.3). Moreover, we want investigate the rate of convergence to ascertain the computational cost of training. As the loss L (2.4) is non-convex, it is hard to rigorously analyze the training process in complete generality. One needs to make certain assumptions on (2.4) to make the problem tractable. To this end, fix step l in (2.4) and start with the following Taylor expansion,
T u( (2; 0k)) u1x; 0o) Vou 0o) T(p, 0o) (0x 0o) (x: 0o)- (2.5)
Here, Ha (x) == (es Ox + (1 Tx)) is the Hessian of u(-, O) evaluated at intermediate values, with 0 S Tk < 1. Now introducing the notation (x) = 0g (x;; On), and assuming that Díb, E (Q2), we define the matrix A Raxx and the vector B E IR" as,
Aij = (Dé,, (o,)) + Mos L.)
(2.6)
IB = (f - DDug, D (2(S) + A(u uso (.) L2(8))-
Substituting the above formulas in the GD algorithm 2.4). we can rewrite it identically as,
0x+1 = 0z VoLk) = (1 - n4)0, + (A06 B) nEkt
(2.7)
where ER is an term that collects all terms that depend on the Hessians H₃ and DH A full definition and further calculations can be found in SM A.T From this characterization of gradient descent (.. clearly see that (2.4) is related to simplified version of gradient descent given by,
Ú (I nA)0 + vTAO + B), 0, = Ou.
(2.8)
modulo the error term E defined in (2.
In the following Lemma (proved in SM A.2 we show that this simplified GD dynamics (.8 approximates the full GD dynamics 2.4 desired accuracy as long as the error term 5x is small.
Lemma 2.1. Let 8 > 0 be such that max; ||ex||2 S 5. If A is invertible and n c/ max, (A)
for some 0 < c < 1 then it holds for any k E N that,
10 0x112 S 8/min (A)))
(2.9)
3

---

Published as a conference paper at ICLR 2024
The key assumption in Lemma 2.1 is the smallness of the error term EK 2.7 for all k. This is trivially satisfied for linear models Ug( = En ByOr as E = 0 for all in this case. From the definition of Eg (SM (A.I), we see that more general sufficient condition for ensuring this smallness is to ensure that the Hessians of ug and Dug (resp. Hą and DH₃ in (2.5) are small during training. This amounts requiring approximate linearity of the parametric function u(- ;)) near the initial value 9o of the parameter 0. For any differentiable parametrized function fa. its linearity is equivalent to the constancy of the associated tangent kernel @[fe] y) Vofo(x)) (V) (Liu al. 2020). Hence, it follows that if the tangent kernel associated to tig and Dug is (approximately) constant along the optimization path, then the error term EK will be small.
For neural networks this entails that the neural tangent kernels (NTK) θ[ue] and [[Dug] stay ap- proximately constant along the optimization path. The following informal lemma, based on Wang et al. (2022b), confirms that this is indeed the case for wide enough neural networks. A rigorous version of the result and its proof can be found in SM[A.3]
hidden layer of width and linear differ- ential operator DD it holds that limm-+oo @[uo,] = limm-soc @luo,] and limm- @Duu,, = lim (Duu,,) for all k. Consequently, the error term Ek 2.7 is small for wide neural net-
Lemma 2.2. For a neural network us with
works, limm-+so max; ||=x||2 =0.
Given the much simpler struc- ture of (2.8). when compared to (2.4). we can study the corresponding gradient descent dynamics explicitly and obtain the following convergence theorem (proved in SMA.4).
Convergence of Simplified Gradient Descent Iterations 2.8
Theorem 2.3. Lel A in 2.8 be invertible with condition number K< (AA),
x(A)= Amax (A)/ Amiu (A) = max /A)\ min [A)(
(2.10)
5
and let 0 < <<< Ser = c/ Amax (A) and 6" =0 + A BB. It holds for any k E N that,
1 0* 1|2 s (1-c/ c//((A)))) & 110 1-
(2.11)
An immediate consequence of the quantitative convergence rate 2.11 is as follows: to obtain an error of size E, i.c., |10g- - 112 S E, we can readily calculate the number of GD steps (s) as, ((e) In (n(=////// (e/110- 0r /22 In(1 - c/x(A)) = O (m)
(2.12)
Hence, for a fixed value c, large values of the condition number x((A) will severely impede conver- gence of the simplified gradient descent 2.8) by requiring much larger number of steps.
Operator Conditioning. So far, we have established that, under suitable assumptions, the rate of convergence of the gradient descent algorithm for ahysics-informed machine learning boils down to the conditioning of the matrix A (2.6). However, at first sight, this matrix is not very intuitive and we want to relate it to the differential operator DD from the underlying PDE (2.I). To this end, first introduce the so-called Hermitian square AA given by A = D'D, in the sense of operators, where D the adjoint operator for the differential operator D. Note that this definition implicitly assumes that the adjoint DD' exists and the Hermitian square operator A is defined on an appropriate function space. As an example, consider as differential operator the Laplacian, i.e., Du -u,, defined for instance u E H' (22) then the corresponding Hermitian square is Au Auu, identified as the bi-Laplacian that is well defined on u E (?))..
Next for notational simplicity, we set X = 0 in (2.2) and omit boundary terms in the following. Let H. be the span of the functions Ox := 00:uu(;;8 60). Define the maps T R" + H. v + Σx= Ukok and T* : ((2) + RR:, f - {(01. f) f) We define the following scalar product on T, (12),
(f.g'= ==
= f," g)...
(2.13)
Note that the maps T, T provide a correspondence between the continuous space (I.2) and discrete space (H) spanned by the functions Φ2. This continoous-discrtte correspondence allows us to relate the conditioning of the matrix A (2.6) to the conditioning of the Hermitian square operator = D through the following theorem (proved in A.5..
Theorem 2.4. It holds for the operator A o TT (22 + (2) that k<(s) 2 x(A T Moreover, if the Gram matrix (db, 6)x H is invertible then equality holds, ie. x(A) K(A TT*).
4

---

Published as a conference paper at ICLR 2024
Thus, we show that the conditioning of the matrix A that determines the speed of convergence of the simplified gradient descent algorithm (2.8) for pepssics-infomed machine learning is intimately tied with the conditioning of the operator TTT This operator, in turn, composes the Hermitian square of the underlying differential operator of the PDE (2.1. with the so-called Kernel Integral operator TT", associated with the (neural) tangent kernel θ [ug) Theorem 2.4 implies in particular that if the operator A o TT is ill-conditioned, then the matrix A is ill-conditioned and the gradient descent algorithm (2.8) for physics-informed machine learning will converge very slowly.
Remark 2.5. One can readily generalize Theorem 2.4 to the setting with boundary conditions (i.e.. with X > 0 the loss (2.2g) In this case one can prove for the operator A A= = DD an Id, and its corresponding matrix A (as in (2.6)) that k(A) 2 x(A TT*) in the general case and x(A) = x(A T if the relevant Gram matrix is invertible. The proof is given in SM
an
Remark 2.6. It is instructive to compare physics-informeet machine learning with standard super- vised learning through the prism of the analysis presented here. It is straightforward to see that for supervised learning, i.e., when the physics-informed loss in (.2) is replaced with the supervised loss 2llu- =0 by simply setting 'D = Id, the corresponding operator in Theorem 2.4 is simply the kernel integral operator TT associated with the tangent kernel as A = Td. Thus, the complexity in training physics-informed informed machine learning models is entirely due to the spectral properties of the Hermitian square A of the underlying differential operator D.
3 PRECONDITIONING AND IMPROVING TRAINING IN PHYSICS-INFORMI MACHINE LEARNING.
Having established in the previous section that, under suitable assumptions, the speed of training physics-informed machine learning models depends on the condition number of the operator A TT" or, equivalently the matrix A (2.6). we now investigate whether this operator ill-conditioned and if so, how can we better condition it by reducing the condition number. The fact that AA T (equiv. A) is very poorly conditioned for most PDEs of practical interest will be demonstrated both theoretically and empirically below. This makes preconditioning. i.e., strategies to improve (reduce) the conditioning of the underlying operator (matrix), a key component in improving training for physics physics-informet informed machine learning models.
Intuitively, reducing the condition number of the underlying operator T can amount to finding new maps T, 7 for which the kemel integral operator 7T ze A- i.e. choosing the architecture and initialization of the parametrized model ug such that the associated Kernel Integral operator TT is an (approximate) Green's 's function for the Hermitian square of the differential operator D. For an operator AA with well-defined eigenvectors 1/s and eigenvalues wk. the ideal case TT = A-1 is realized when T ox. = wk thk-
Explicit preconditioning by linearly transforming parameters. The above ideal case can be achieved by transforming é (in (2.6) linearly with a (positive definite) matrix IP such that (IP Φ)x Vop thk which corresponds to the change of variables Puo == upo- Assuming the invertibility of
(é, din, Theorem 2.4 then shows that k(A o TT*) k(A) for new matrix A, which can be computed as,
À = puro, DV DVouo La()) = (DP Vo Vougo DPP V Vouoo) = ATAP (3.1)
This implies a general approach for preconditioning, namely linearly transforming the parameters of the model, i.e. considering Pug := upe instead of ug. which corresponds to replacing the matrix A by its preconditioned variant Ā = P
AP. The new simplified GD update rule is then O.+1 = 02 - (0, 0b) + B. Hence, finding TT R A-1, which is the aim of preconditioning, reduces to constructing a matrix P such that 1 R x(A) < n(A). We emphasize that TT need not serve as the exact inverse of A: even an approximate inverse can lead to significant performance improvements, this is the foundational principle of preconditioning.
Explicit preconditioning by linearly transforming the gradients. Given that any positive def- inite matrix can be written as PPT linearly transforming the parameters is equivalent to precondi-
5

---

Published as a conference paper at ICLR 2024
tioning the gradient of the loss by multiplying with positive definite matrix, in the sense:
a++1 = P0 P0, - 7PP (gL(POk)) = 0g - nPP voL(Ok).
(3.2)
which corresponds to performing gradient descent using the transformed parameters 02 := IPO Hence, parameter transformations are all that are needed in this context.
Analysis of the impact of preconditioning for the Poisson equation. As an example, we start with linear parametrized models of the form uz = Σx 0x0x (z), where $1,..., On are any smooth functions. A corresponding preconditioned model, as explained above, would have the form ät (a) = Ex (PO) (I), where P E R is the preconditioner. We motivate the choice of this preconditioner with simple, yet widely used example.
Faurier Features Number
Epochs for Different Fourier Features
10
NR Fourier Foats.
10"
Condition Number
10
25
10-1
10a
55
10-5
10
75
10.
10
Precond
102
Precond.
100
NP of Fourter Features
Epochs
Figure 1: Poisson equation with Fourier features. Left: Optimal condition number vs. Number of Fourier features. Right: Training for the unpreconditioned and preconditioned Fourier features.
Our differential operator the one-dimensional Laplacian D de defined on the domain (- for simplicity with periodic zero boundary conditions. Consequently, the corresponding PDE 2 is the Poisson equation. As the machine learning model, we choose (x())) = Σ b= K Oudi (x)). with (2) - d- v* cos(kz) and Φu (x) - sin(kx) for 1 S k K. This model corresponds to the widely used learnable Fourier Features in the machine learning literature (Tancik) jet al.) 2020) or spectral methods in numerical analysis (Hesthaven et al.) (007). We can readily verify that the resulting matrix A. (2.6) is given by A = D + Aeu where ID is a diagonal matrix with DD 1 and == ((x). Preconditioning solely based on D'D would correspond to finding a matrix IP such that PDP = ld. However, given that Doo = O, this is not possible. We therefore set IP R = 1/k² for ls + 0 and PPo - R. The preconditioned matrix is therefore
d=
V2n
y
7A, PDP APo(Po)
(3.3)
The conditioning of the unpreconditionedd and preconditioned matrices considered above are sum- marized in the theorem (proved in B. below,
Theorem 3.1. The following statements hold for all K N:
1. The condition number of the unpreconditioned matrix above satisfies k(A(A)) 2 K4.
2. There exists a constant C(A,7) 7) >> that is independent of K such that Ks (A,, 7) S C.
3. It holds that x² ((AA(1//1) (22 7)) + 0(1/7) and hence lim xs (2 (2m /?2 ()) = 1.
We observe from Theorem 3.1 that (i) the matrix A, which governs gradient descent dynamics for approximating the Poisson equation with learnable Fourier features, is very poorly conditioned and (ii) we can (optimally) precondition it by rescaling the Fourier features based on the eigenvalues of the underlying differential operator (or its Hermitian square).
These conclusions are also observed empirically. In Figure I (left), we plot the condition number of the matrix A, minimized over À (see SM Figure 8 and SM C for details), as a function of maximum K+, predicted by the Theorem 3.1 Consequently as shown in Figure (right), where we plot the loss function in terms of increas- ing training epochs, the underlying Fourier features model is very hard to train with large losses
frequency K and verify that this condition number increases

---

Published as a conference paper at ICLR 2024
(particularly for higher values of K), showing very slow decay of the loss function the number of frequencies is increased. On the other hand, in Figure I (left), we also show that the condition number (minimized over X) of the preconditioned matrix (3.3) remains constant with increasing fre- quency and is very close to the optimal value of 1, verifying Theorem [ As a result, we observe from Figure (right) that the loss in the preconditioned case decays exponentially fast the num- ber of epochs are increased. This decay is independent of the maximum frequency of the model. The results demonstrate that the preconditioned version of the Fourier features model can learn the solution of the Poisson equation efficiently. in contrast to the failure of the unpreconditioned model do so. Entirely analogous results are obtained with the Helmholtz equation (see SMC)
B vs Condition Number
Loss Epochs Different
104
104
12m
10
Conditon Number
18m
204
10"
30m
3б
10-3
10a
48
10-4
10"
10a
Precond,
10-22
15
20
25
30
Precand.
150
Epochs
Figure 2: Linear advection equation with Fourier features. Left: Optimal condition number vs. B. Right: Training for the unpreconditioned and preconditioned Fourier features.
As a different example, we consider the linear advection equation U + Bu, 0 on the one- dimensional spatial domain x [0, 2n] and with 2 2r-periodic solutions in time with t E [0, 1]. As in Krishnapriyan et al. (202T), our focus in this case to study how physics-informed machine learning models train when the advection speed / > O is increased. To empirically evaluate this example, we again choose learnable time-dependent Fourier features as the model and precondition the resulting matrix A. (2.6) as described in SM B.2.2 see also SM In Figure [ (left), we see that the condition number of A(3) ~ 8² grows quadratically with advection speed. On the other hand, the condition number of the preconditioned model remains constant. Consequently as shown in (riure((riight), the unpreconditioned model trains very slowly (particularly for increasing values of the advection speed (3) with losses remaining high despite being trained for large number of epochs. In complete contrast, the preconditioned model trains very fast, irrespective of the values of the advection speed B. Further details including visualizations of the resulting solutions and comparison with a MLP are presented in SM B.2.2 and Figure [3 In particular, we show that the preconditioned Fourier model readily outperforms the MLP. Other additional experiments can be found in (M
Viewing available strategies for improving training in physies-infrrmdd machine learning models through the lens of operator (pre-Jconditioning. Given the difficulties encountered in training physics-informed machine learning models, several ad-hoc strategies have been proposed in the recent literature to improve training. It turns out that many of these strategies can also be interpreted using the framework of preconditioning that we have proposed. We provide succinct summary below while postponing the details to the SM.
Choice of X. The parameter X in the loss (.2) plays a crucial role as it balances the relative contributions of the physics-informer loss R and the supervised loss at the boundary B. Given our framework, it is natural to suggest that this parameter should be chosen as X == min, (A)))), in order to obtain the smallest condition number of A. and accelerate convergence. In SM B.2 we present X for the 1-D Poisson equation with learnable Fourier features to find that A (K) ~ K2, with K being the maximum frequency. In turns out that finding suitable values of À has been widely proposed as strategy, see for instance Wang et al. (202 2022b) which propose algorithms to iteratively learn À during training. It turns out that applying these strategies leads to different scalings of À with respect to increasing K for the Fourier features model (see SM B.2 for details), distinguishing our approach for selecting X.
7

---

Published as a conference paper at ICLR 2024
Hard boundary conditions. From the very advent of PINNs (Lagaris et al., 2000a b). sev- eral authors have advocated modifying machine learning models such that the boundary conditions in PDE (2.1) can be imposed exactly and the boundary loss in (.2 is zero. Such hard imposition of boundary conditions (BCs) has been empirically shown to aid training, e.g. Moseley et al.] (2021); Dolean et al. (2023) and references therein. In SM [B.3] we present example where the linear advection equation is solved with learnable Fourier features and show that imposing hard BCs reduces the condition number of A, when compared to soft BCs. Thus, hard BCs can improve training by better conditioning the gradient descent dynamics, at least in some cases.
Second-order optimizers. There are many empirical studies which demonstrate that first-order optimizers such as (stochastic) gradient descent or ADAM are not suitable for physics-informed machine learning and one needs to use second-order (quasi-NNewton type optimizers such as L-BGFS in order to make training of physics-informed machine learning models feasible. In models and show that as the Hessian of the loss is identical to the matrix A (2.6) in this case, (quasi-JNewton methods automatically compute an (approximate) inverse of the Hessian and hence, precondition the matrix A. relating the use of (quasi-)Newtoo type optimizers to preconditioning operators in this context.
SM B.4] we examine this issue for linear physics-informe
Domain decomposition. Domain decomposition (DD) is widely used technique in numer- ical analysis to precondition linear systems that arise out of classical methods such as finite elements (Dolcan et al., 2015). Recently, there have been attempts to use DD-inspired methods within physics-informed machine learning, see Moseley et al. (202T): Dolean et al.] (2023) and references therein, although no explicit link with preconditioning the models was established. In SM B.2.2. we re-examine the case of linear advection equation with learnable Fourier features to demonstrate that increasing the number of Fourier features in time by decomposing the time domain simply amounts to changing the effective advection speed / and reducing the condition number, leading to a better--onditioned model. Moreover, in this case, this algorithm also correlates the causal learning based training of PINNs (Wang et al., 2022a). which also can be viewed as improving the condition number.
4 DISCUSSION.
Summary. Physics-informed machine learning models are notoriously hard to train with gradient descent methods. In this paper, we aim for a rigorous explanation of the underlying causes as well as examining possible strategies to mitigate them. To this end, under suitable assumptions that coincide with approximate linearity of models, we prove that gradient descent with physies-informe informed losses is approximated well by novel simplified gradient descent procedure, whose rate of convergence can be completely characterized in terms of the conditioning of an operator, composing the Hermi- tian square of the underlying differential operator with the Kernel integral operator associated with the underlying tangent kernel. Thus, the ill-conditioning of this Hermitian square operator can ex- plain issues with training of physics-informed learning models. Consequently, preconditioning this operator (equivalently the associated matrix) could improve training. By a combination of rigorous examine strategies with view of how one can precondition the associated operators. In particular, we find that rescaling the model parameters, as dictated by the spectral properties of the underlying differential operator, was effective in significantly improving training of physics-informa informed models for the Poisson, Helmholtz and linear advection equations.
analysis and empirical evaluation,
Related Work. While many studies explore the mathematical aspects of PINNs, the majority fo- cus on approximation techniques or generalization properties (De Ryck & Mishra 202T Doumèche et al. 2023). Few works have targeted training error and training dynamics, even though it stands as a significant source of overall error (Krishnapriyan et al., 2021). Some exceptions include Diang et al. (2023). who examine global convergence for linear elliptic PDEs in the NTK regime. However equations are derived in continuous time, thereby sidestepping ill-conditioning (which is intrinsi- cally linked to discrete time) and thus potential training issues. Wang et al. (202 identified that PINNs might converge slowly due to a stiff gradient flow ODE. Our work allows to interpret their proposed novel architecture, which reduces the maximum eigenvalue of the Hessian, as way pre- condition TT*, as the Hessian of the loss equals A (SM (S.M..4). thereby improving the convergence rate (Theorems 2.3 and 2.4 Wang et al. (20225) derive a continuous-tima evolution equation exclu-
8

---

Published as a conference paper at ICLR 2024
sively for the residual during training, leaving out direct exposition of the Hermitian square term, and contrasting our discrete evolution equation in parameter space, as opposed to function space. Wang et al. (2021a; 2022b) also propose algorithms to adjust the À multiplier between boundary and residual loss terms, which we assess within the context of operator preconditioning in SM SM Works aiming improve convergence of PINNs based domain decomposition strategies include Jagtap & Karniadakis (2020); Jagtap et al. (2020); Wang et al. 2022a); Koopnnčááooá et al. 20233) some which can be reinterpreted as methods to precondition A by changing A or T7"
Limitations and Future Work. In this work, our examples for elucidating the challenges in train- ing innfores-nformee machine learning models focussed on linear PDEs. Nevertheless, the analysis already revealed the key role played by equation-dependent preconditioning. Extending our results to nonlinear PDEs is a direction for future work. Moreover, while highlighting the necessity of preconditioning, the current work does not claim provide universal preconditioning strategy, particularly for nonlinear models such as neural networks.
We strongly believe that the complications arising from ill-conditioning merit further scrutiny from the scientific computing community, such as those specializing in domain and operator precondition-
ing. There is much work in this domain (Mardal & Winther 2011 Hiptmair 2006) and references therein, providing a fertile ground for innovative approaches, including the potential application of non-linear preconditioning techniques commonly used in these fields. However, extending our work to these settings exceeds the scope of this paper and remains a direction for future inquiry.
Another aspect worth discussing pertains to our linearized training dynamics (NTK regime), in which feature learning is absent (Chizat et al. (019)). For low-dimensional problems typical in many scientific settings (1-3D), the lack feature learning may not be a significant handicap, as one discretize the underlying domains. Extensive evidence in this paper has shown that the linear bases often outperform nonlinear models. However, neural networks might still outperform linear models high- sigh-dimensional problems (Mishra & Molinaro 2021), highlighting the significance of deviations from the lazy training regime.
Finally, we would like point that our analysis be readily extended to cover physics-informet operator learning models such as those considered in [L et al.] (2023): Goswami et al. (2022) by adopting the theoretical framework of representative neural operators (Bartolucci et al.
2023
Raonić et al. 2023).
REFERENCES
Francesca Bartolucei, Emmanuel de Bézenac, Bogdan Raonić, Roberto Molinaro, Siddhartha Mishra, and Rima Alaifari. Are neural operators really neural operators? frame theory meets operator learning, 2023.
J. Bruna, B. Peherstorfer, and E. Vanden-EEjnden. Neural galerkin scheme with active learning for high-dimensional evolution equations. arXiv preprint arXiv:2203.01354 2022.
Lénaïe Chizat, Edouard Oyallon, and Francis R. Bach. On lazy training in differen- tiable programming. In Hanna M. Wallach, Hugo Larochelle, Alina Beygelzimer, Flo- rence d'Alché- Buc, Emily B. Fox, and Roman Garnett (eds.), Advances in Neural In- formation Processing Systems 32: Annual Conference on Neural Information Process- ing Systems 2019, NeurIPS 2019, December 8-14, 2019, Vancouver, BC, Canada, pp. 2933,2943, 2019, URL
neu.isss....e................... CC cc/pppprr/019/haa.. html
Salvatore Cuomo, Vincenzo Schiano Di Cola, Fabio Giampaolo, Gianluigi Rozza, Maziar Raissi, and Francesco Piccialli. Scientific Machine Learning Through Physics-Informed Neural Net- works: Where we are and What's Next. Jourmal of Scientific Computing, 92(3): jul 2022. doi: 10. 10.10077/1109115002--019999 1tt.:/1111110/11111111001 020-0011--22--11--- -1122-0013922
ISSN 15737691.
URL httpp::/ springer.
Com
T De Ryck, S. Mishra, and R. Molinaro. wpinns: Weak physics informed neural networks for ap- proximating entropy solutions of hyperbolic conservation laws. arXiv preprint arXiv:2207.08483, 2022.

---

Published as a conference paper at ICLR 2024
Tim De Ryck and Siddhartha Mishra. Error analysis for physies informed neural networks (PINNs) approximating Kolmogorov PDEs. arXiv preprint arXiv:2106. 2021.
Tim De Ryck, Ameya D. Jagtap, and Siddhartha Mishra. Error analysis for PINNs approximating the Navier-Stokes equations. In preparation, 2021.
MWMG Dissamayake and N Thien. Neural-network-basee approximations for solving partial differential equations. Communications in Numerical Methods in Engineering, 1994.
Dolean, A. Heinlein, Mishra, and B. Moseley. Multilevel domain decomposition- based archi- tectures for physics-informan informed neural networks. arXiv preprint WrXii:23000004886 2023.
Victorita Dolean, Pierre Jolivet, and Frédéric Nataf. An introduction to domain decomposition meth- ods. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2015. ISBN 978- 1-611974-05-8. URL http: /dx. doi. rithms, theory, and parallel implementation.
ch1 Algo-
Nathan Doumèche, Gérard Biau, and Claire Boyer. Convergence and error analysis of pinns, 2023.
Weinan E and Bing Yu. The Deep Ritz Method: A Deep Learning-Bai -Bandd Numerical Algorithm for Solving Variational Problems. Communications in Mathematics and Statistics, 6(1):1-12, March 2018. ISSN 2194-671X. doi: 10.1007/140334-018--077--
Lawrence C Evans. Purtial differential equations, volume 19. American Mathematical Soc., 2010.
Behrooz Ghorbani, Shankar Krishnan, and Ying Xiao. An investigation into neural net opti- mization via hessian eigenvalue density. In Kamalika Chaudhuri and Ruslan Salakhutdinov (eds.), Proceedings of the 36th International Conference on Machine Learning, ICML 2019, 9- 15 June 2019, Long Beach, California, USA, volume 97 of Proceedings of Machine Learning Research, pp. 2232-22413 PMLR, 2019. URL http: hattp://pren http::///ccceeeeeeeeeeoee mir. wtt...press/v/. ghorbani 9b.9b... html
Gene H Golub. Some modified matrix eigenvalue problems. SIAM review, 15(2):1183335 1973.
Ian Goodfellow. Yoshua Bengio, and Aaron Courville. Deep learning. MIT press. 2016.
Somdatta Goswami, Aniruddha Bora, Yue Yu, and George Em Kamiadakis. Physics-informed deep neural operator networks, 2022.
Jan S Hesthaven, Sigal Gottlieb, and David Gottlieb. Spectral methods for time-dependent problems, volume 21. Cambridge University Press, 2007.
R. Hiptmair. Operator preconditioning. Computers & Mathematics with Applications, 52(5)6699- 706, 2006. ISSN 1221. doi: http::/lbo.org/10.0001....aaaa.. 10.008. URL https:
www. www.sssceeeeeireettt....ee com/ http:///iiieerrreeee///////22221 Topics in Applied and Industrial Mathematics.
Hot
Ameya Jagtap and George Em Karniadakis. Extended physics-informed neural networks (XPINNs): A generalized space-time domain decomposition based deep learning framework for nonlinear partial differential equations. Communications in Computational Physics, 28(5):2002-
2020. 2020.
Ameya D Jagtap, Ehsan Kharazmi, and George Em Karniadakis. Conservative physics-informed neural networks on discrete domains for conservation laws: Applications to forward and inverse problems. Computer Methods in Applied Mechanics and Engineering, 365:1 13028, 2020.
Deqing Jiang, Justin Sirignano, and Samuel N Cohen. Global convergence of deep galerkin and pinns methods for solving partial differential equations. arXiv preprint arxiv:2305.06000 2023.
George Em Karniadakis, Ioannis G. Kevrekidis, Lu Lu, Paris Perdikaris, Sifan Wang, and Liu Yang. Physics informed machine leaming. Nature Reviews Physics, pp. 1-19, may 2021. doi: 10. 1031/ 542254-021-00314: URL www. nature. ctm/ hevphhyys
Kharazmi, Z Zhang, and G. Em Karniadakis. Variational physics informed neural networks for solving partial differential equations. arXiv preprint arXiv: 2019.
10

---

Published as a conference paper at ICLR 2024
Alena Kopaničáková, Hardik Kothari, George Em Karniadakis, and Rolf Krause. Enhancing training of physics-innooe informed neural networks using domain cecompiiccoomooittionn based preconditioning strate- gies. arXiv preprint arXiv:2306. 2023.
Aditi S. Krishnapriyan, Amir Gholami, Shandian Zhe, Robert M. Kirby, and Michael W. Mahoney. Characterizing possible failure modes in physics-informed neural networks. In Mare Aurelio Ranzato, Alina Beygelzimer, Yann N. Dauphin, Percy Liang, and Jennifer Wortman Vaughan (eds.), Advances in Neural Information Processing Systems 34: Annual Conference ral Information Processing Systems 2021, NewrlPS 2021, December 6-14, 2021, virtual, pp. 26548 26560, 2021. URL https::/proceed:: haa.se1tbbsstttt..... .ttrractttmmmm
Neu-
I. E Lagaris, A Likas, and D. Fotiadis. Artificial neural networks for solving ordinary and partial differential equations. IEEE Transactions on Newral Networks, 95):987-10000 2000a.
Isaac E Lagaris, Aristidis Likas, and Papageorgiou G. D. Neural-network methods for bound- ary value problems with irregular boundaries. IEEE Transactions on Neural Networks, 11:1041- 1049, 2000b.
Zongyi Li, Hongkai Zheng, Nikola Kovachki, David Jin, Haoxuan Chen, Burigede Liu, Kamyar Azizzadenesheli, and Anima Anandkumar. Physics-informed neural operator for learning partial differential equations, 2023.
Chaoyue Liu, Libin Zhu, and Misha Belkin. On the linearity of large non-linear models: when and why the tangent kernel constant. Advances Neural Information Processing Systems, 33: 15954-15964, 2020.
Kent-Andre Mardal and Ragnar Winther. Preconditioning discretizations of systems of partial dif- ferential equations. Numerical Linear Algebra with Applications, 18(1):1-40, 2011. doi: https: //bbi.org/10.1002/nla... URL 1002/nla.17 716
wwwewwiwly............ com
Siddhartha Mishra and Roberto Molinaro. Estimates on the generalization error of physics informed
neural networks (PINNs) for approximating PDEs. arXiv preprint 210.20066.6144, 2020.
Siddhartha Mishra and Roberto Molinaro. Physies informed neural networks for simulating radiative transfer. Journal of Quantitative Spectroscopy and Radiative Transfer, 270:107705, 2021.
Ben Moseley, Andrew Markham, and Tarje Nissen-Meyer. Finite Basis Physics-Informed Neu- ral Networks (FBPINNs): scalable domain decomposition approach for solving differential equations. arXiv, jul 2021. URL http::/aawwa...c/////////20/20/1
//arxiv. wttp://r//aa////uuuee//s 2ttpob/221777778.8
A. Quarteroni and A. Valli. Numerical approximation of Partial differential equations, volume 23. Springer, 1994.
Nasim Rahaman, Aristide Baratin, Devansh Arpit, Felix Draxler, Min Lin, Fred A. Hamprecht, Yoshua Bengio, and Aaron Courville. On the spectral bias of neural networks. In 36th In- fernational Conference Machine Learning, ICML 2019, volume 2019-June, pp. 9230-9239. International Machine Learning Society (IMLS), jun 2019. ISBN 9781510886988. URL
fttp:///rxivve org/ waw/gg/bbb. abs/ 1806. 08734
M. Raissi, P. Perdikaris, and G. E. Karniadakis. Physics-informed neural networks: A deep learn- ing framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational Physics, 378:686-707, 2019.
Maziar Raissi and George Em Karniadakis. Hidden physics models: Machine learning of nonlinear partial differential equations. Journal of Computational Physics, 357:125-141, 2018.
Bogdan Raonić, Roberto Molinaro, Tim De Ryck, Tobias Rohner. Francesca Bartolucci. Rima Alai- fari, Siddhartha Mishra, and Emmanuel de Bézenac. Convolutional neural operators for robust and accurate learning of pdes, 2023.
11

---

Published as a conference paper at ICLR 2024
Carl Edward Rasmussen. Gaussian processes in machine learning. In Summer School Machine Learning. pp. 63-71. Springer, 2003.
Matthew Tancik, Pratul Srinivasan, Ben Mildenhall, Sara Fridovich-Keil, Nithin Raghavan, Utkarsh Singhal, Ravi Ramamoorthi, Jonathan Barron, and Ren Ng. Fourier features let networks learn high frequency functions in low dimensional domains. Advances in Neural Information Process- ing Systems, 33:7537-7547, 2020.
N. Trefethen and M. Embree. Spectra and Pseudo-spectra. Princeton University Press, 2005.
Sifan Wang, Yujun Teng, and Paris Perdikaris. Understanding and mitigating gradient flow patholo- gies in physics-informed neural networks. SIAM Journal on Scientific Computing, A3355 A3081, 2021a.
Sifan Wang, Hanwen Wang, and Paris Perdikaris. Learning the solution operator of parametric partial differential equations with phyysics-infommed deeponets. arXiv preprint arXiv:2103.1 10974, 2021b.
Sifan Wang. Shyam Sankaran, and Paris Perdikaris. Respecting causality is all you need for training physics-informed neural networks. arXiv preprint MxXiv:2003.07004, 2022a.
Sifan Wang, Xinling Yu, and Paris Perdikaris. When and why pinns fail to train: A neural tangent kernel perspective. Journal of Computational Physics, 449: 10768, 2022b.
Yinlin Ye, Yajing Li, Hongtao Fan, Xinyi Liu, and Hongbing Zhang. Slenn-elm: shifted legen- dre neural network method for fractional delay differential equations based on extreme learning machine, NHM, 18(1):4944555 512,4511222 2023.
12

---

Published as a conference paper at ICLR 2024
Supplementary material for:
An operator preconditioning perspective on training in informed machine learning
A DETAILS FOR SECTION 2
A.1 DERIVATION OF EQUATION
We provide detailed calculations for how to obtain (2.7 First of all, plugging in (2.5) into (2.2) and calculating yields,
V VgR(Or) = (Du(x; 0b) f())) ((a) Oo)z + / DV DVou(x; (o) Bb) T (0x- 8o )zz
1
(Du(i;; 0x) (())) 00)da + ./ (0s (0s- 0o) DHz (0x- -00) pu(a; 0o)
2
(u(x; 80) -(x)))())uaa))) Vou(x; 0o)dr / Vou(z; Vou 0 (0s 0o)dz
VoB(Ok)
80
(m(:;; 0k) - u(x)) (0x- - 0o 0o)da + 2/ (0x - 00) HA (Ox - 6o) vou 0o)
2 aa
an
Now introducing the notation (x) 0g u(;; ; 0o). we define the vector Ek E R as,
Ek -= (Duo, - -f 0o)) (2(()) + N(uo. u, H0 0o)) 1.(
+ ((0x DD) Bu), DVe.e,) 12(0) + A((0g 0o) Ha 00), Voue))ra(00))
(A.1)
Combining this definition with those of A and gives rise to equation 2.7
A.2 PROOF OF LEMMA 2.1
Proof. If maxe ||ex|l2 < then the error one makes compared to the simplified GD update rule is bounded by
10: 0:112
Σ nA)* TEk
Sn
11
hAll
hA
(A.2)
1 1
max A)(()m)1 -
E - c/n())(()))
(A.3)
k=
k=
x(A)o
(A.4)
miny (4)
A.3 PROOF OF LEMMA
We first provide rigorous version of the statement of Lemma 2.2 which includes all technical de- tails. To do so, we follow the setting of (Wang et al. 2022b Section 4). The result is generalization of (Wang al. 20226 Theorem 4.4) that is suggested Wang et al.) (2022b) itself.
Lemma A.1. Let D be linear differential operator, a smooth function and m, K E N. Let Wo Rixm W22 Rmxx1 the weights and to E Rmm and b² E R the biases that constitute the parameters = (W', W1 bb,b²) of the neural network un R - IR defined by,
1 _
ua(x) u(x;; 0)
W/ (W -
+ +b'
(A.5)
Vm
Let (8x)k be the series of parameters obtains from performing gradient descent the loss L(0), where all entries of the parameters are initialized as iid N 1) random variables. Furthermore, assume that,
13

---

Published as a conference paper at ICLR 2024
1. The parameters stay uniformly bounded during training, iee there is constant C independent of such that maxisksK 110x llo I C.
2. There exists a constant > 0 such that Σ L(O
3. If is an n-th order differential operator, then there exists a constant C > such that
(000) (x)) S for all and 0 < 2
Under these conditions it holds for all k that,
lim O[uoy] = lim O Ouuaol
and lim @(Dun,] = lim O[Duopl. (A.6)
Hence, the error term
A./ is small for wide neural networks, limm- maxk ||a 0.
Proof. In (Wang et al. 2022b Theorem 4.4) they proved that limm-+00 U[ue,] limm θ [ue,]] and limm-soo [[u, lim [Du,,,] the case that D 0 and state that can be generalized
general linear differential operators. We provide the steps to how this more general result can be proven, by summarizing all necessary changes in the proof of (Wang et al. 20226 Theorem 4.4).
Changes to Wang et al. 2022b Lemma D.I). Instead unnn 0) one needs to provide a formula for Du(a; (x; 0). This can be done by noting that,
Σ W" (W))),
(A.7)
where we slightly abuse notation by letting W and bj. be the k-th components of the respec- tive vectors, i.e. they do not correspond the k-th gradient descent update 0g. The bounds (2; then follow in a similar way from the boundedness of both 0 and a and its deriva- tives.
Changes to Wang el al. 2022b Lemma D.2). Here we replace the gradient flow by the gradient descent formula 0k+1 Os the only difference that we consider a problem that is discrete (training) time (as we are using gradient descent).
Vo The rest of the calculations in the same spirit, with
Changes to Wang et al. 2022b Lemma D.4). The calculations are completely analogous, only that one needs to replace the Formulas for Vouga (x; (x;0) 0) by formulas for VoDu(a; 0). Again in our simplified case that D = OF we find that,
@Du(x;;) 0)
(W2)"-1W3 W
no (n) (WW.WWWW b0) + fut ( + b2)
(A.8)
aw
vm
(Du(x;0) 0)
(n) (n) (WW. (Wox+ (2)
(A.9)
OW1
Vm
(uu(x; 0)
w Wilo (n+)(++)((()))) + (
(A.10)
06g Vm
(u;(;; 0)
ab 0.
(A.11)
The rest of the proof is a matter of calculations that are analogous to those in Wang et al. 2022b) One can conclude that in the limit vn - so the tangent kernel of both up and Dup constant dur- ing training. As (approximate) NTK constancy is equivalent to (approximate) linearity Liu et al. 2020), find that the Hessians Hş and DHg, and consequently the error term Eg must be (ap- proximately) zero. More concretely, the error term E (.I)) small for wide neural networks, lim max; || = 0.
14

---

Published as a conference paper at ICLR 2024
A.4 PROOF OF THEOREM 2.3
Proof. We calculate that,
02 = (1 - 1-74)00011 + +7(8800 + 1B)
(A.12)
4-1
= (1 - An)*
-1 Anj (AA0 + 13) 1=
(A.13)
= [I An) 0o + (Am)) -I (I -nn) )y + B))
(A.14)
= (I An) (n)000 00 (1 -1-(1- -
(A.15)
and hence
02 0" = An) (0p 8*)
(A.16)
We can then take the norm, use that 1j = c/ Amax (A) and calculate,
10x 0² 2l = 11 - nAll 1|2 I max || - (A)) $100 - 1-
(A.17)
cAmin (A) Amax (A)
1100 -0 (I-x(A)
m((a) ) 1100 00 1-
(A.18)
A.5 PROOF OF THEOREM 2.4
Proof. Using that A = TTAT AT we find
Amax (A) max T Av 2 ||=||=1
sup
(T (T
(A.19)
SELP:IT--IILL
sup
(T" fe1 ||T* f||=1
sup (o fh
(A.20)
FE1?:
[T*
sup fEL2: [T-11=1
(f.A.TT' fix = Amax (A TT*)
(A.21)
The last equality holds by the minmax theorem since T is a self-adjoint operator. Indeed, this operator is self-adjoint since & is symmetric and given that for f.g H. we find that,
(f TT 9) -T* f) ((T g).
(A.22)
One can do similar calculation for Amia (with the reverse inequality). which concludes the proof.
Now let M- (Ø) o) be the Gram matrix, i.c. M,j = (Ø,, (y, and we assume is invertible. For any parameter vector we can define fo TM As result we have the chain of inequalities,
7 L2(S).. One can verify that T fo v.
sup IE
(T f) A(T*))
f) < max TTAN TAv |+||=1
(A.23)
max (T" So) So) ATT f.) ||||=
(A.24)
<
sup
(T (TT))T(TT' f),
(A.25)
SEL: 1| ʃ||=1
and hence suP [ [|= (T
ma|x T Av. Doing analogous calculation
the one above, we find that s(A o TT*) =AA)
A.6 PROOF OF REMARK 2.5
Let H. be the span of the de == 0o ((:;; Define the maps T : R= - H. + Ex_ &R and
T* (Ω)) - R + ((0x, S) + Mob, (') k I...... We define the following scalar
product on (Ω)::
(f.9) := (f.TT" 9) L2(0) + Xf.TT" 9) (I g).
(A.26)
15

---

Published as a conference paper at ICLR 2024
We define the operator A
A= 1DDDDD al ld,
(A.27)
and the matrix A as,
Aut (0k, (D)DD)) + N(px: 01) L2(0))-
(A.28)
Let M = (Ø, 1) h be the Gram matrix, i.e. M,j (di, k5) If it invertible then for any parameter
vector e we can define fo TM- E H C (9t). One can also verify that T fc. v.
We give the proof in the case where M is invertible. The general case similar to the proof of Theorem 2.4
Proof. We calculate that for any vector w E R" it holds that, (dx, (Tw) L (0x, D"DTT)) A(0x Tw) 12(80)
(A.29)
(A.30)
Σ (0r, D'D) X(xx, Ø1) (0)))) wi
(A.31)
(Aw))
(A.32)
where we used that the measure of aΩ zero with respect to the Lebesgue measure on f and that f is zero on an. Using this and the definition of the scalar product on H we find that,
X max (A)
max Av
(A.33)
||]| =1
sup (I" fo)TA(I AATT fo)
(A.34)
|||==
sup (T" fo) (T ATT fa)
(A.35)
|M|
sup (fu, ferH
(A.36)
|+|=1
Amaxx (A TT*)
(A.37)
This concludes the proof, as a similar calculation can be done for Amin (A).
B DETAILS FOR SECTION 3
B. PROOF OF THEOREM 3.1
We first the state the following auxiliary result.
Lemma B.1 (Section 5, Golub (173)). Let À 2 0, let u E R", let D E Ruxa be a diagonal matrix with eigenvalues d, such that d, S d++1 and let C D Aua The eigenvalues wi...., Wn of (in ascending order) are the roots of
p(w) II (d. - w)
Σ
2
II (d, w) 0.
(B.1)
1=1
j=1jjii
and it holds that d, S Wi S d:+1 and dn ww & dn
+ Au T u.
We now prove Theorem values of X, y and K.
3.1 In Figure we also plot the condition number of A(), 7) for various
Proof. We first prove statement (1). Denote by w; resp. d, the i-th eigenvalue (in ascending order)
of A resp. D. It follows from Lemma B.1 that w1 S d2 = 1 and w2K+1 +1 2 d2x+1 = K Hence we find that for any À > 0 it holds that,
W2K-1
d2K+1 K
x(A(A))
(B.2)
Wi
d2
16

---

Published as a conference paper at ICLR 2024
We continue with statement (2). It is easy to verify that PDP has eigenvalues 1 (multiplicity 2K)
and 0 (multiplicity 1) and that the vector = (((T) has entries w_k (-1)* /TT 1/v2m,
ea for K. Define v Σ 1/k and note that 1 1 म / 1.08.
Now, by Lemma B.1 the eigenvalues of are given by the roots of,
1 w)2k mww w)2k-1
pak+1(w) -w(1 w)2k +
(B.3)
2
1
(1 - w)2K-1
(-2Tw(( (1 - w) + >> -w) -2XwT)
(B.4)
2m
1
w)2k (2Tw² - (2 X7² + 2\n) + X7?).
(B.5)
2
We can already see that 1 is an eigenvalue with multiplicity at least 2K - The other two eigen- values are given by.
w+ (1)
2 + X7² 22n = v (2x + Xy2 + 2Xn)² -ATT?)
(B.6)
4T
Now since
< x² /90 there exists constant C independent of K such that,
max w[+ min 1, w- SC
VKEN: m(A(A, c))
(B.7)
Now suppose that sets À 2m//2 Then we find that
w1( n/7 + 2n/7² 1/2/77
(B.8)
and hence,
+ 2 x(A(A,))) 7)) = 1
0(1/)
for 7 + 00.
(B.9)
As result, <(A(2m/77:7)) 11
y-10.0
1.35
2=20.0
2-00
40.0
conditio runEa
10a
y-60.0
2.25
1-000
1.0
y=00.0
y-100.0
10
10--1
Figure 3: Evolution of condition number of preconditioned matrix for different y and A (left) and for different for the optimal À (right).
B.2 CHOOSING
As mentioned in the main text, it seems natural to suggest that the parameter À in loss 2.2 should
be chosen as
X = min (((()))
(B.
in order to obtain the smallest condition number of & and accelerate convergence, provided that be calculcated in numerically stable way.
17

---

Published as a conference paper at ICLR 2024
B.2.1 POISSON EQUATION
As first example, we revisit the setting of Theorem 3.1 where considered solving the Poisson equation in one dimension by learning the coefficients of basis of Fourier features (see formulas in main text).
In Figure 8 was already shown how the condition number of (A) changes for various values of and K. particular, there is very clear minimum in condition number for every K. In Figure we can clearly see that the optimal X scales as K
10
optimal
10r
1ptmna
107
10F
condition number optimal condition number with A g22 number
10
10
Figure 4: Left: optimal X in terms of K. Right: Evolution of condition number of A(A) in terms of model size K for multiple choices of X.
We now compare our choice of as equation B. 10 with related work. We first compare with Wang al. (2021a). In this work, the authors propose learning rate annealing algorithm to iteratively update throughout training. At each iteration, they define
maxe (V R(O)) (2K + 1)-1 Σx (V)B()))
X
(B.11)
In our example it holds that (VoR(0)): K00 for + l and (VoR(0))e e (Op 1). We find that max; (V, ((0)) K at initialization. Next we calculate that (given that initialization holds that 6m (0, iid),
K
K
E (VoB()))
Σ E
1mm(11kktm
mm
k+m
E
0m
(K
1)3//2
(B.12)
k=
m==
k=
m==
This brings us to the rough estimate that (2K + 1)-1 EE (VgB(0)) VK. So in total we find that K3.5.
Next, we compare with the proposed algorithm of Wang et al.
2022b).
They propose to set
X (n)) Tr(Kuu
Jg VoDu, Vo Jan
2x EK K
2k=1h b4
K
(B.13)
Voug Vouo
K+1
So, both works propose a choice of À that increases a lot faster in than the optimal choice in terms of the condition number (X* K2). In Figure (right) we compare the obtained the condition numbers for the various choices of X* and we observe that the relative difference is however small for small K but increases for larger K. This phenomenon be explained by the fact that x(A(A)) has relatively wide minimum for larger K.
B.2.2 LINEAR ADVECTION EQUATION
As a second example, we revisit the linear advection equation that was studied in Section B The PDE is given by 44 + Bu, 0 on (-æ, )²² and we periodic boundary conditions. As model
18

---

Published as a conference paper at ICLR 2024
us uz t) Ekeza IIEIl.IINIl (x t), where for instance (x(2) sin(k1 + ket) when &r and ex (x) = cos(kko + k2t) when 22 2 0.
One can calculate that Atm (k2 3.k1) Aōk where the boundary term now comes form initial condition for u2 0), and where we rescaled A and À to get rid of any multiplicative unnormalized Fourier features. Re-indexing from matrix å with indices in 2² to a matrix with indices in such that A + Km,+m = Agm
factors stemming from the fact that
then gives the formula A A'()) ld & Ca 23C & C C @ ld AId @ (1 1..11 1. where C is diagonal matrix with Cr L.
In Figure we plot the condition number of A' for various values of À and B. We see that the condition number increases with increasing and that for each 3 there is a clear optimal X In Figure 6 we verify that X indeed scales as B², which was already observed in Figure
Condition number versus and for K
Condition number and
0.01
10P
1.0
0014800 runber
10.0
1003
100
10000
100002.0
109
100
10
15
20
Figure 5: Evolution of condition number of A for soft boundary conditions (full lines) and hard boundary conditions (dotted lines) in terms of 33 and for K K=3
This analysis also raises the questions what happens if T increased. This depends whether the model size K is changed along with K or not.
T' is increased but the number of used Fourier basis functions is unchanged. The first part of A', i.e. the one coming from D' is independent of X and scales linear in T. As result, rescaling T corresponds to rescaling optimal condition number will not depend on T'
Thus, the optimal X" will change but the
T is increased and the number of Fourier basis functions increased accordingly in time. This essentially corresponds changing and rescaling the measure. Moreover, this is intrinsically connected to what happens for domain decomposition in the time dimension,
35000
mincond(a.A .)
quadralic growth
30000
25000
conditiin number
20000
15000
20000
Figure 6: For each 3 we have minimized the condition number terms of and displayed it. Growth is approximately quadratic for large B.
19

---

Published as a conference paper at ICLR 2024
which on its can be form of causal learning strategy (Wang et al. 2022a Since the condition number scales as 3 one can argue that splitting the time domain in N parts will reduce the condition number of each individual submodel by factor of N
B.3 HARD BOUNDARY CONDITIONS
In the previous section, have how varying the weighting parameter of the boundary term can influence the condition number of A. In view of Theorem 2.4 this is due to the fact that one changes the operator A. In this section, argue that different choices to implement hard boundary conditions can also change the condition number of A, through changing the operator TT* We also give some examples where the condition number for hard boundary conditions is strictly smaller than for soft boundary conditions i.e., ((haard BC) < min, K AAgtt (A)).
Demonstrative examples. We first consider a toy model to demonstrate these phenomena in straightforward way. We consider again the Poisson equation -u sin with zero Dirichlet boundary conditions. We choose as model ug (z) 8 cos(x) + 6u + 81 sin(x) and compare the condition number between various options.
Soft boundary conditions. Using the approach discussed in Section B.2 we find that the condition number for the optimal X is given by + 2/2 ss 5.83.
Hard boundary conditions variant A first common method implement hard boundary conditions is to multiply the model with function n(z)) so that the product exactly satisfies the boundary conditions, regardless of ug. our setting, we could consider 17(x)) jud (x) with m)(=)) 0. For sin the total model is given by
((x)ue(z) 0-1 cos(22t) (xx)
0-1
01 + 0o (insin())) + (2r),
(B.14)
and gives rise condition number of 4. Different choices of will inevitably lead to different condition numbers.
Hard boundary conditions variant 2. Another option would be to subtract u=(x) from the model so that the boundary conditions are exactly satisfied. This corresponds to the model,
(x( ug u (cos(((2)) + 1) + 01 sin(x)
(B.15)
Note that this implies that one can discard 0g as parameter, leaving only two trainable parameters. The corresponding condition number is 1.
Linear advection equation. With these instructive examples in place, we once again revisit the linear advection equation (as in SM B.2.2). If uar, 0) = sin(ar) then we could define our model as,
(x(x)
Σ Ox(ex(x) (x, 0)) sin(ax),
(B.16)
and correspondingly.
Due(x) =
Σ
0x((ks + 3k1))x(x)) 3k1ex(z,0) + 3a cas()x),
(B.17)
&
and hence,
Akm (12 +
+bk+Bk1)²0 (2,Bk)))0mn 20km 3km Bkg (k2 + Bkg )k, 88, + 8m,o] + (8k1)² 8k.m, (8kk))
(B.18)
If we define the matrix E 1 1.e by Eas 6ep then we can write the re-indexed matrix A' as,
A' = Id & C2 + 23C @ C+ c @ Id 3(3C @ C + 30 & Id)((Id @ (E+ (E+E + (1-1T).
(B.19)
We find that for this specific case the condition number of A' with hard boundary conditions is always smaller than that with soft boundary conditions, for any or X (see Figure
20

---

Published as a conference paper at ICLR 2024
B.4 SECOND-ORDER OPTIMIZERS
Second-order optimizers distinguish themselves from first-order optimizers, such as gradient de- scent, by using information of the Hessian in their update formula. In the case of Newton's method, the Hessian is used to precondition the gradient. We explore the connection of preconditioning based on our matrix & (2.6) with Newton's 's method.
Given any loss [0), Newton's 's method's gradient flow is given by d0(l)
(B.20)
dt where H is the Hessian. In our case, we consider the physics-informed loss L(O) = R(0) + (B
(2.)) and consider the model (z(x) = E, (x) which corresponds linear model or a neural network in the NTK regime (e.g. Lemma 2.2
Using that 0,,0o, uo 0, we calculate,
0g,00 0o, L(0) = (Co,, uo(x)) Ldd, uo(x)dr + f do uo(o)(2) do, us (x))d
22
0n
(D)(()) ( X A 2((x) - (00)),0(1)))) (x)dr
(B.21)
Aaj-
We conclude that is identical the Hessian of the loss. Hence, using Newton's method auto- matically leads to perfect preconditioning. However, the high potential cost of computing the exact Hessian inhibits the use of these class of optimizers in practice.
C EXPERIMENTAL DETAILS AND ADDITIONAL RESULTS
of Fourier Features vs Condition Number
Loss vs Epochs for Different of Fourler Features
NO Fourier Feats,
10f
10a
10a
Conditoo Number
25
10a
45
10
10a
75
10
10.*
No Precond Precond.
100
MP of Fourler Features
Epochs
Figure 7: Helmholtz equation with Fourier features. Left: Optimal condition number vs. No. Fourier
features. Right: Training curves for the unpreconditioned and preconditioned Fourier features.
Throughout all the experiments, from torch. func. The matrix A is computed using the formula
torch version 2.0.1 and incorporate functional routines
Aij (Dé;, Dōj) (2)) + X(ss, (2)
Here, d& 3g 0o) is initially obtained through autograd. Subsequently, the differential opera- tor 'D is applied, also leveraging autodiff. For the scalar product, we employ Monte Carlo approx- imation and can utilize the same input coordinates used for training. scenarios where the matrix becomes too harge-often the case for large networks-numerical approximations of the D operator can be employed to reduce both computational and memory load.
For training, in order avoid random errors we fix the grid to equispaced grid of size large enough to resolve the neural network/lincar function.
Once the matrix A is computed, for the linear models we use routine compute its condition number and find the optimal À before training using the black-box gradient free 1d optimisation method from scipy opt inttimiz ze. golden. The learning rate is then chosen as 1/Amar. For models with MLPs this no longer possible because of the zero eigenvalues that are plaguing the matrix; resort grid search sweeping wide ranges of learning rates and À values.
21

---

Published as a conference paper at ICLR 2024
Condition Number vs À for Different Ne of Fourier Features With Preconditioning
Without Preconditioning
Condition Number
Condition Number
20
10*
202
107
10a
10
10--
Figure 8: Sensitivity of the condition number on parameter X, the multiplier between boundary and residual loss. the non-preconditioned and preconditioned case (left and right, respectively). has a non non- negligiiee impact on the condition number. However, in the preconditioned case, this impact always the same for different Fourier features, and thus plots for different Fourier features overlap. This implies that the optimal À will not have to be adjusted when changing this hyperparameter, as opposed to the non-preconditioned case.
C.1 POISSON AND HELMHOLTZ, ID
We solve the following boundary value problems on the domain [-т, T.. For the experiments with Fourier features the network is already periodic and thus there is no need enforce these conditions.
Poisson equation, 1d.
The solution is given by u(x) sin (kx).
u' u²(x) 2 sin (kxx)
u(-r)= O, u(T) = O.
(C.1)
Helmholtz equation, 1d.
The solution given by u(x) = cos(war).
u²(x) + ++(u(z)) 0,
u(0) 1, v0) = =0, u(-=)= = u(T), u'(-T)= u'()).
(C.2)
For the Fourier features approximating Helmholtz equation, we use exactly the same formulation as for the Poisson equation in the main text and a similar preconditioning matrix as in (3.3] but with diagonal entries scaling as - The results for Helmholtz equation, presented in Figure 7, are entirely analogous to the observations for the Poisson equation in Main text Figure | All the observations done on those cases in the main paper and the appendices looking at the mean squared error, see Figure and In Table we report the computational time for one epoch the preconditioned and unpreconditioned cases where, for fairness, we computed the preconditioning at each step instead of computing it once for all at the beginning of the training.
also valid when
Architectural Details. For MLP we used a 3 hidden neural network with a hidden dimension of 64 and hyperbolic tangent tanh activation function. All models are optimized using SGD.
C.2 LINEAR ADVECTION EQUATION, 2D
We solve the following boundary value problems on the domain Ω = [0, 2Tl] x 10, 1]:
Du Øt + 3 dr
= 0
u(x, 0) = sin(z) [x),
x E [0, 2m te [0
u(0,t) u(2 t),
(C.3)
22

---

Published as a conference paper at ICLR 2024
solution 8 6
Fourier features, parameter preconcitioiing
-0.5
-0.5
0.4
0.6
0.5
1.0
-1.0
0.2
0.4
0.6
0.8
10
-1.0
Fourier features, no preconditioning
MUP, no preconditioring
1.0
1.0
5
0.5
2-
8.0
0.2
0.4
0.6
0.8
1.0
-100
8.0
0.2
0.4
0.6
0.8
1.0
-1.0
Time
Time
(a)
6x
Exact solution
precondiiioning
0.0
0.0
0.5
0.5
1
80
0.0
Fourier features, no precenditioning
MLP: no preconditioriig
1.0
-0.5
-0.5
0.4
0.6
0.5
1.0
-1.0
0.2
0.4
0.6
0.8
1.0
-1.0
(b) B
30m
Exact solution for 6- 60n
Fourier features, preconcitioning.
S
0.5
0.5
2-
8
0.5
8.0
Time
Time
Fourler preconditioning
MUP: preconditioring
0.0
0.0
0.5
0.5
1-
8.
8.
Time
(c) B = 60
Figure Solutions for the linear advection experiment
23

---

Published as a conference paper at ICLR 2024
Residual loss vs Epochs
Mean squared error vs Epochs
NR Fourier Features
10"
10a
15
10a
10a
10""
45
10-1
75
10-31
95
10-10
10
No Precond. Precond.
10
100
150
200
50
1000
150
200
Epochs
Epochs
Figure 10: Training curves and mean squared error for the Fourier basis in the Poisson problem for different numbers of Fourier features using preconditioned SGD and Adam.
Residual loss vs Epochs
Mean squared error vs Epochs
No Fourier Festures
10a
15
10a
25
10-1
10-1
10
45
10
55
75
10 10
05
10
95
10
10
10-46
Precond.
Epochs
Epochs
Figure 11: Training curves and mean squared error for the Fourier basis in the Helmholtz problem for different numbers of Fourier features using preconditioned SGD and Adam.
where B is the speed of advection.
Each model up. with 0 the set of its parameters, trained with Physics-Informed Loss defined as:
R(0) J [ |0:u0 30.00 1 drdt
us (x, 0) - -in(x)]] dr + / (us(0, t) ua(2TT f) ()) dt
B(0) = |0,2x]
L(0) = R(0)= + XB(O)
(C.4)
The solution of such differential equation is (x, + sin(x 36).
Architectural Details. We test models, linear model via two-dimensional Fourier series with bounded spectrum and a MLP with 5 layers, 50 neurons per hidden layer, and tanh activation function as used Krishnapriyan et (202T). More explicitly, the ansatz for the linear model is:
5
30
up(x,t)
Σ
Σ
ORan cos(k:r + 22mm)) gein sin + 2xmd)
(C.5)
V
m==
with 0 (oro Ogin }
k.m.
In the preconditioned case, the preconditioning matrix is set a diagonal matrix with diagonal terms m for (k. m) o (0, 0). For the parameter goos we set the preconditioning factor to I.
Training details. We discretize the domain Ω on grid of size 256 100 for the learning process to be consistent with the configuration proposed in Krishnapriyan et al. (2021) and grid of size 256 2048 to compute the matrix A defined in 2.6 with finite differences.
For the linear model, we use classical batch gradient descent on the full grid with and without preconditioning for 200 epochs. The learning rate and À are automatically defined to minimize the
24

---

Published as a conference paper at ICLR 2024
Table 1: Computational time per step in milliseconds. The preconditioning time refers the param- eters preconditioning of the Fourier basis via the approximate analytical conditioning described in the Appendices.
Problem
Preconditioning
SGD
Poisson ID
14.1 = 2.1
8.9 + 2.2
Helmholtz ID
15.9 3.2
11.2 + 2.8
Advection 2D
17.2 = 1.5
10.9 + 0.7
condition number of A via a Golden section method. For the MLP, we use Adam on the full grid without preconditioning for 10000 epochs. The learning rate is set to 0.0001 and À is set via grid search as we encountered out-of-memory errors when computing A. Training curves for the linear model with and without conditioning for different 3 and a plot of the condition number with respect to 3 is given in Figure 2 Training curves for the MLP for different 3 are given in Figure 13
The behavior of condition numbers and loss functions for the linear model and its precon- ditioned version were already shown in Figure 2 and discussed in the main text. The loss functions with MLPs (trained with ADAM), for different Iss are shown in Figure 13 and we observe that the MLP (trained with ADAM) converged very slowly to an unacceptably large loss function of ampli- tude 10-2, which can be contrasted with the very fast decay of the loss function to 10-10 10 for the preconditioned Fourier model (Figure[(right)) These issues in training severely impact the overall quality of the results with the three models considered (unpreconditioned Fourier, preconditioned Fourier and MLP). In Figure 9 we plot the exact solution (in space-time) and compared it with the solutions obtained with the three models, at three different advection speeds / = 69, 30x, 60v. We see from this figure that the unpreconditioned Fourier model fails completely, even at the slow- est considered speed 3 = 6τ. MLP (trained with ADAM) is better at the slowest speed than the unpreconditioned Fourier model but fails completely at the higher speeds, as already observed in Krishnapriyan et al. (2021). On the other hand, consistent with theory and the low values of the loss function observed in Figure[] (right), the preconditioned Fourier model was able to approximate the solution (in space-time) with high-accuracy, even for the fastest considered advection speed of 3 = 60x. This example clearly demonstrates the superiority of the preconditioned linear models over their unconditioned versions as well as over nonlinear neural network models. All of this ob- servations are also valid when looking at the mean squared error (computed at the same collocation points as the training loss), see Figure T2 In Table T the computational time is reported for each method where, for fairness, we computed the preconditioning at each step instead of computing it once for all at the beginning of the training.
Results.
Residual vs Epochs
squared emor Epochs
8
10
10a
12n
107
10m
107
24
LoGs
10-1
30m
10-1
36/
10
10
40m
10
Sán
10-41
10-1:
10-14
10-1a
150
Precond.
150
Epochs
Epochs
Figure 12: Training curves and mean squared error for the Fourier basis in the advection problem for different 3 using preconditioned SGD and Adam.
25

---

Published as a conference paper at ICLR 2024
Advection, MLP, Adam
Loss vs Epochs for Different Values of
10a
10
6n
12t
10a
24h
10-7
10
420
10
10
60m
10-30
10-12
2000
4000
6000
8000
10000
Epachs
Figure 13: Training curves for the MLP in the advection problem for different / using Adam.
C.3 PRECONDITIONING VIA DIRECT INVERSION OF THE MATRIX
In this section, we wish validate the theory and assess the assumptions made in a practical, non- linear scenario. To this extent, we propose to directly precondition the gradient using the inverse of the matrix defined in (2.6). recalled below:
Aijj = (Dos, D (9()) + X(ghi 05) (2())))
(C.6)
where (x) = 0g, u(x; 0o). Given that the simplified gradient assumption holds, i.c. the smallness of the error term in (.7). conveeeeneeeeeee is entirely governed by the conditioning of the matrix A. This matrix can be readily computed using autodifferentiation, first differentiating w.r.t. A, then w.r.t. the coordinates for the differential operator D. The scalar products then computed using Monte Carlo approximation on the training points.
We then precondition the gradients, using the regularized inverse A-1:= (A +el) the parameter updates given by:
0x41= = 02 - -
(C.7)
Recalling equations (3.2 and 3.1 the conditioning of the problem now governed by the matrix
Ā = A AA z= I, whose condition number is thus approximately optimal. Experiments conducted using a one-hidden layer, 32 neurons, tanh MLP, using a regularization parameter & E [0.0001, 0.004), and double precision.
The preconditioned gradient directly numerically computed via the resolution of a least-square problem of the form:
argmin | VoL(Ok)
(C.8)
We must underline that the regularizing term E/ and the double precision are crucial for numerical stability as the matrix A is highly ill-conditioned for such MI.. Also, we computed the matrix A. for each new set of parameters, solving the least-square problem each epoch instead of solving only at the initialization.
For Poisson, Helmholtz, and Advection, we compare the convergence of the preconditioned MLP to non-precondiitooee MLP, trained with SGD and Adam respectively. In figures 14 report training curves (left) as well as mean squared errors (right) for different PDE parameters. We observe that preconditioning with this strategy not only consistently yields much lower errors, but also much faster training times. Solutions for the advection equation can be visually inspected in Figures 18 and T9 consistent with findings in Krishnapriyan et al. (2021), the ML. trained
15 16
26

---

Published as a conference paper at ICLR 2024
with preconditioning does not converge the true solution, even after indefinite number of timesteps. In contrast, solutions with the preconditioned MLP are visually indistinguishable from the true solution. However, this preconditioning strategy comes at price, computational costs are higher (as opposed to the preconditioning strategy described in the main paper), as reported in Table 2
Residual loss Eoochs
Mean squared error vs Epochs
10a
10a
10
10a
10
10-3
10r
10"1
10
10.4
10"
10
No Precond.
10"9
200
400
800
800
1000
Procond.
200
400
600
800
1000
Epochs
Epochs
Figure 14: Training curves and mean squared error for the shallow MLP in the advection problem for different 3 using preconditioned SGD and Adam.
Residual Epochs
Mean squared Épochs
10a
10-3
10-1
toss
Los
10
10
10
10-1
10
10
10-3
No Precand
10-33
200
400
600
800
1000
Precond.
200
400
600
800
1000
Epachs
Figure 15: Training curves and mean squared error for the shallow MLP in the Poisson problem for different & using preconditioned SGD and Adam.
Residual vs Epochs
squared Epochs
10
10a
10a
107
107
Loss
10-3
10-3
10
10-4
10t4
10
Prezond
200
Precond:
200
Epochs
Epochs
Figure 16: Training curves and mean squared error for the shallow MLP in the Helmholtz problem for different w using preconditioned SGD and Adam.
C.4 PRECONDITIONING NONLINEAR PHYSICS-INFORMET MACHINE LEARNING MODELS
We investigate conditioning of nonlinear models by considering Poisson equation on (- n and learning its solution with neural networks of the form up(x) = Φ with æ E (-T, π) and T being a three hidden layer MLP. We compute the resulting matrix A (2.6) and plot its normalized eigenvalues in Figure 20 (left) to observe that most of the eigenvalues are clustered
zero, in accordance with the fact that Hessians (which are correlated with A) for neural networks have a large number of (near) eigenvalues (Ghorbani al. 2019). Consequently, it is difficult to
27

---

Published as a conference paper at ICLR 2024
Distribution eigenvalluss of
Residual less vs Epochs
207
Precond, maxic) 1.00
Precond.
10
301
2011
10-1
2011
aymaxxim)
Spechs
Figure 17: Left: Normalized Spectrum of the shallow MLP and preconditioned shallow MLP for the Poisson Equation. Right: Training loss for the shallow MLP and preconditioned shallow MLP for the Poisson Equation.
MI.P. gradient proconditianing
no proconsitioning
S
0.5
6.5
0.5
Space
0.5
0.5
0.5
1
0.
80
-1.0
8.0
Time
Time
Time
Figure 18: Comparison of the solutions for the shallow MLP in the advection problem using pre- conditioned SGD and Adam with / 10.
Exact solution for 8 30
1.0
MLP. gradient preconttioinng
1.0
MLP. no preconeitioning
1.0
Space
0.0
5.0
0.0
0.2
0.4
0.5
3.3
-1.0
0.2
0.4
0.0
0.5
-1.0
0.4
0.5
Figure 19: Comparison the solutions for the shallow MLP in the advection problem using pre- conditioned SGD and Adam with B 30.
Table 2: Computational time per step in milliseconds. The preconditioning time refers to the gradient preconditioning of the shallow MLP via inversing the A matrix.
Problem
Preconditioning
Adam
Poisson ID
16.2 + 2.0
8.2 + 1.8
Helmholtz ID
20.5 3.6
9.7 + 2.0
Advection 2D
200.0 + 25.8
5.9 + 1.1
analyze the condition number per se. However, we also observe from this figure that there are only couple of large non-zero eigenvalues of A indicating very uneven spread of the spectrum. It well- known in classical numerical analysis (Trefethen Embree (005)) that such spread-out spectra are very poorly conditioned and this will impede training with gradient descent. This is corroborated in Figure 20 (right) where the onysics-informed MLP trains very slowly. Moreover, turns out that preconditioning also localizes the spectrum (Trefethen & Embree 2005). This is attested in Figure 20 (left) where we see the localized spectrum of the preconditioned Fourier features model considered previously, which is also correlated with its fast training (Figure (right))
28

---

Published as a conference paper at ICLR 2024
Distribution of eigenvalues of
101
fhecondtiooning MLP Preconditoncd Founter Festures
10-1
1.0
0.0
0.2
0.4
0.6
0.5
Figure 20: Poisson equation with MLPs. Left: Histogram of normalized spectrum (eigenvalues multiplied with learning rate). Right: Loss number of epochs.
However, preconditioning & for nonlinear models such as neural networks can, in general, be hard. Here, we consider a simple strategy by coupling the MLP with Fourier features Φa defined above, ie., setting
D (E aubk). Intuitively, one must carefully choose og to control d di (x)
as it will include terms such as (Eiyt (u((k))² oR ) (Ex axox). Hence, such rescaling can better condition the Hermitian square the differential operator. To test this for Poisson's 's equation, we choose d% 1/k² (for A: f 0) in this FF-MLP model and present the eigenvalues in Figure 20 (left) to observe that although there are still quite few (near)-zero eigenvalues, the number of non-zero cigenvalues is significantly increased leading much more even spread in the spectrum, when compared to the unpreconditioned MLP case. This possibly accounts for the fact that the resulting loss function decays much more rapidly, shown in Figure (right,,, when compared to unpreconditioned MLP.
29


Reviews:
[{"review_id": "cYbJBZe5i6", "reviewer": "Reviewer_Wdfd", "summary": "The paper derives that under certain conditions, the gradient descent of the PINNs can be approximated using a simplified gradient descent algorithm, which further provides a linearized version of the training dynamics. The paper demonstrates that the speed of convergence of this simplified training dynamics depends on the condition number of an operator described using the underlying differential operator and the kernel integral operator. The paper provides a preconditioning approach that can almost achieve the ideal condition number of 1. The paper demonstrates the efficiency of preconditioning the operator approach for simple problems like linear Poisson and Advection Equation.", "questions": "**Questions:**\n1. For a typical Taylor series expansion, the Hessian would be computed at the point $\\theta_0$. Can the authors provide some justification on why the Hessian was computed at an interpolated $\\theta$ value between 0 and k?\n2. Can the authors provide some discussions on the connections and differences of their approach and NTK theory, especially considering that the latter also offers a convergence rate analysis for PINNs?\n3. For the Advection Equation with Fourier Features without preconditioning (Figure 2), it seems that the loss does not converge at all (the value of the loss is $10^3$). Can the authors provide some justification behind this? Is it because the model was chosen to be a simple linear one with Fourier Features?\n\n\n**Minor Comments:**\n1. Typo: Equation 2.2 the $\\lambda$ is missing.\n2. The Taylor expansion in Equation 2.5 should contain higher order terms or the paper should mention that the 3rd order and higher terms are ignored. Similarly, the proof in the appendix A.1 should mention that higher order multiplicative terms by substituting 2.5 in 2.2 are ignored. \n3. The notation: $(\\theta_k − \\theta_0)^T H_k(\\theta_k − \\theta_0)$ can be confusing as it seems like the Hessian is computed at $H_k(\\theta_k - \\theta_0)$. The authors can consider using the following: $(\\theta_k − \\theta_0)^T H_k(x) (\\theta_k − \\theta_0)$", "limitations": "", "rating": 8, "confidence": 3, "soundness": 4, "presentation": 3, "contribution": 4, "strengths": "1. The paper derives a simplified gradient descent algorithm that linearizes the training dynamics, allowing easier analysis of the convergence speed.\n2. The paper provides a strong theoretical foundation, and shows a novel result that the convergence of PINNs depends on the condition number of a certain operator.\n3. Some of the results shown in the paper demonstrates super fast convergence (less than 50 epochs), which is quite impressive.", "weaknesses": "1. The paper assumes that the solution can be expressed as a linear combination of basis functions (such as fourier features, which was shown in the experiments). Such an approximation can be quite limiting, and the paper does not show how the preconditioning can be extended for non-linear models. \n2. The evaluation of the proposed preconditioning approach is poor. Firstly, the paper only covers linear PDEs with linear PINN models. Second, only the training loss curves were shown for the two PDEs. It is well-known that the PINNs can learn trivial solutions and achieving very low training losses does not guarantee convergence. See paper [1]. \n\n[1] Daw, Arka, Jie Bu, Sifan Wang, Paris Perdikaris, and Anuj Karpatne. \"Mitigating Propagation Failures in Physics-informed Neural Networks using Retain-Resample-Release (R3) Sampling.\" ICML 2023", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "t1JuyFCyYD", "reviewer": "Reviewer_gVcB", "summary": "This paper proposes a preconditioning technique for solving optimisation problems associated with solving ODE, which can be expressed through the differential operator used and some kernel integral operator which is variable and can be chosen to minimize the impact of conditioning number of the optimisation problem. This approach appears to be practically efficient for model of Fourier features.", "questions": "I assume that there are models for which kernel integral operator cannot be chosen manually. What you propose to do if expressing precodnitioned operator explicitly is not possible? What do you think about applying additional optimisation procedure for training this kernel operator?", "limitations": "", "rating": 8, "confidence": 4, "soundness": 4, "presentation": 3, "contribution": 4, "strengths": "The proposed preconditioning is elegantly inferred from the approximate form of iteration so that the approach seems to be natural. Theoretical framework allows obtaining improved convergence rates with no dependence on conditioning number of original problem, which is the issue even in finite-dimensional optimisation. Practical efficiency of the precodnitioning operation in its characterising property, decreasing conditioning number, was demonstrated.", "weaknesses": "No significant weakness that I can notice as a specialist in finite-dimensional optimisation. What was described is quite earnestly and helpful for future readers.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "Vl5fuQT2vo", "reviewer": "Reviewer_2Mn2", "summary": "The paper investigates behaviors of gradient descent algorithms when the loss term involves minimizations of the residual connected to partial differential equations (PDEs). The paper simplifies the gradient descent rules via utilizing the Taylor expansion and derives the relationship between the rate of convergence and the differential operator. In particular, the paper provides theoretical analysis showing that the rate of convergence depends on the conditioning of an operator denoted by ``Hermitian square’’, which is derived from the differential operator and its adjoint. The paper proposes a method for preconditioning the differential operator and tests the algorithms on benchmark problems.", "questions": "- could the authors provide more information/insight on how to build preconditioner for more general cases? for complex PDEs where the PINN does not employ Fourier feature transform? \n\n- In the implementation, the gradient descent rule has been reimplemented following Eq. (3.2)?", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "- The paper is well-written and theoretically sound (although there are some typos, most notably, there is no $\\lambda$ in Eq. 2.2). \n\n- The paper tackles an important issue in minimizing PDE residual losses, which is not popularized in many scientific machine learning applications.", "weaknesses": "- It is less convincing if the proposed method can be applicable to a wide variety of neural networks, which is also pointed out in the manuscript, “preconditioning A for nonlinear models such as neural networks can, in general, be hard” on page 7 and in the limitation paragraph. In the manuscript, preconditioning strategies are showcased only with the simple cases, approximating solutions of Poisson equations via linear and nonlinear models with Fourier feature mapping. It does not seem that the paper provides some guidelines or insights for preconditioning strategies for general cases.\n\n- Computational aspects are weak. In particular, comparisons with other works improving optimizers of PINNs in more general classes of PDEs, which also include the assessment on computational wall times, would be more appreciated and informative for readers/users to decide which method to use. Even if it’s not computational comparisons, it would be better to have some descriptions on differences from other baselines [2,3] (as shown in page 8, where the connections to [1] is presented in the Choice of lambda paragraph.)\n\n[1] Wang, et al, 2022, JCP (the reference in the paper)\n\n[2] Basir and Senocak, Physics and equality constrained artificial neural networks: Application to forward and inverse problems with multi-fidelity data fusion, 2022, JCP\n\n[3] Kim et al, DPM: A novel training method for physics-informed neural networks in extrapolation, 2021, AAAI.\n\n- Although the problem that is investigated in the paper is an important problem and the authors provide good contributions in the theoretical side, it is less convincing that the paper makes significant contributions in ML or DL perspectives.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "HGsj6d7WTs", "reviewer": "Reviewer_95TT", "summary": "This paper considers the problem of training neural networks to solve partial differential equations (PDEs) using a physics-informed neural networks (PINNs) approach. In the introduction, the authors mention empirical studies suggesting that training might be the bottleneck for achieving robust convergence in PINNs, which motivates their approach of studying the behavior of gradient descent algorithms in a PINNs context. Then, in section 2, the authors study a simplification of a gradient step and show that the rate of convergence depends on the condition number of a Gram matrix associated with the PDE. This motivates the authors approach for preconditioning the operator to achieve better convergence rates in practical applications.", "questions": "- In eq. (2.6), the function space should be defined so that D\\phi_i lies in L^2.\n- I would suggest using \\bm notation for matrices and vectors instead of \\mathbb, which could be confused with the set of complex numbers for \\mathbb{C}.\n- Lemma 2.2, is there a way of controlling the error epsilon_k with respect to m for large m?\n- Typo: Lemma 2.1: \"Let If eta = ...\"\n- Typo: p.7 \"precondiioned Fourier model\"\n- Fig. 3 (left): I don't really see an improvement of the approach as their are still a large number of 0 eigenvalues with the MLP + preconditioned Fourier features approach, suggesting that the system is ill-conditioned.\n- The related work section would be better placed in the introduction.", "limitations": "", "rating": 5, "confidence": 3, "soundness": 2, "presentation": 3, "contribution": 2, "strengths": "- The paper is well-written and focuses on an important problem in PINNs.\n- The authors provide a clear exposition of their approach with both theoretical and empirical results.", "weaknesses": "- On p.3, the authors mention that their aim is to analyze the convergence of gradient descent. However, the analysis in section 2 only considers a simplified gradient step where they neglect all higher order terms in an error term epsilon_k, which they assume to be small. Is that equivalent to performing a local analysis close to a local minima and if so, how does the current analysis relate to the NTK analysis in the context of PINNs (Wang et al., 2022)? It is not clear whether this explains the poor performance of PINNs in practice as the main issue is to understand the pre-asymptotic behavior of gradient descent.\n- Thm. 2.4 shows that the condition number of the Gram matrix is greater than the condition number of a Hermitian square operator A\\circ TT^*. However, to perform preconditioning, one should instead derive an upper bound on the condition number of the Gram matrix. I am aware that Thm 2.4 gives an equality case but it assumes that the Gram matrix is invertible.\n- The practical section 3 is the main weakness of the paper. It seems that the authors approach might only be useful when one already know a really good approximant to the eigenvectors and eigenvalues of the PDE, when it is linear (e.g. the Poisson equation example in Fig. 1). In this case, why would one use a neural network to solve the PDE instead of a spectral method? The authors states that they employ a linear combination of eigenvectors of the Poisson equation as the machine learning model, but this is exactly a spectral method and we know that the performance of an iterative solver (e.g. conjugate gradient, GMRES) for solving the system is related to the condition number of the matrix. The authors give very little details on how one might use their approach in a real context with a complex neural networks to solve a nonlinear PDE.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "EPCHeMs9zT", "reviewer": "Reviewer_E6qc", "summary": "This paper studies the training dynamics of gradient descent (GD) for solving partial differential equations (PDE) with physics-informed neural networks (PINN). Theoretically, the authors analyze a linearized form of the GD update. Let $\\mathcal{D}$ denote the differential operator associated with the underlying PDE. The authors show that the convergence rate of GD depends on the conditional number of the Hermitian square operator $\\mathcal{D}^\\ast\\mathcal{D}$, which implies that preconditioning the operator $\\mathcal{D}^\\ast\\mathcal{D}$ helps accelerate the training process. Empirically, the authors propose preconditioning strategies for several PDEs and provide numerical experiments to exhibit the effectiveness of using preconditioning during training.", "questions": "1. Would it be possible for the authors to discuss the main difference between the preconditioning method proposed in this paper and sobolev training? Is there any numerical example that your method outperforms (higher accuracy, lower computational cost, etc) sobolev training when the models used for parametrizing the PDE solution are the same?\n\n2.  The authors mentioned that there are empirical studies showing that first-order method are less effective compared to second order methods for physics-informed ML. Would it be possible for the authors to include related references here? It seems to me that adaptive methods like Adam can also improve the conditioning as it takes use of diagonal rescaling. Hence, it might be meaningful for the authors to provide some numerical examples to explain why preconditioned GD can outperform Adam/SGD. \n\nA full list of references mentioned in this review is given below:\n\n[1] Lu, Y., Blanchet, J., & Ying, L. (2022). Sobolev acceleration and statistical optimality for learning elliptic equations via gradient descent. Advances in Neural Information Processing Systems, 35, 33233-33247.\n\n[2] Yu, J., Lu, L., Meng, X., & Karniadakis, G. E. (2022). Gradient-enhanced physics-informed neural networks for forward and inverse PDE problems. Computer Methods in Applied Mechanics and Engineering, 393, 114823. \n\n[3] Son, H., Jang, J. W., Han, W. J., & Hwang, H. J. (2020). Sobolev training for the neural network solutions of pdes.\n\n[4] Son, H., Jang, J. W., Han, W. J., & Hwang, H. J. (2021). Sobolev training for physics informed neural networks. arXiv preprint arXiv:2101.08932.\n\n[5] Zeng, Q., Kothari, Y., Bryngelson, S. H., & Schäfer, F. (2022). Competitive physics informed networks. arXiv preprint arXiv:2204.11144.\n\n[6] Czarnecki, W. M., Osindero, S., Jaderberg, M., Swirszcz, G., & Pascanu, R. (2017). Sobolev training for neural networks. Advances in neural information processing systems, 30.\n\n[7] Amari, S. I., Ba, J., Grosse, R., Li, X., Nitanda, A., Suzuki, T., ... & Xu, J. (2020). When does preconditioning help or hurt generalization?. arXiv preprint arXiv:2006.10732.\n\n[8] Zhu, B., Hu, J., Lou, Y., & Yang, Y. (2021). Implicit regularization effects of the Sobolev norms in image processing. arXiv preprint arXiv:2109.06255.\n\n[9] Calder, J., Mansouri, A., & Yezzi, A. (2010). Image sharpening via Sobolev gradient flows. SIAM Journal on Imaging Sciences, 3(4), 981-1014.\n\n[10] Richardson Jr, W. B. (2008). Sobolev gradient preconditioning for image‐processing PDEs. Communications in Numerical Methods in Engineering, 24(6), 493-504.\n\n[11] Soliman, Y., Chern, A., Diamanti, O., Knöppel, F., Pinkall, U., & Schröder, P. (2021). Constrained willmore surfaces. ACM Transactions on Graphics (TOG), 40(4), 1-17.\n\n[12] Yu, C., Brakensiek, C., Schumacher, H., & Crane, K. (2021). Repulsive surfaces. arXiv preprint arXiv:2107.01664.\n\n[13] Raissi, M., Perdikaris, P., & Karniadakis, G. E. (2019). 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.", "limitations": "", "rating": 6, "confidence": 4, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "This paper includes both theoretical analysis and numerical experiments to support the argument that preconditioning is essential for training PINNs to solve PDEs. For numerical studies, the authors provide plots of both the condition number and loss function to exhibit that preconditioning does improve the training efficiency. Moreover, the authors also explain why preconditioning serves as a general framework for improving training efficiency by showing that several training strategies, such as tuning the regularization parameter, enforcing hard boundary conditions, using high-order optimizers and domain decomposition, can be treated as preconditioning the operator in some way.", "weaknesses": "(1) The reviewer's first concern is that this paper doesn't discuss its relation to previous work on preconditioning. Specifically, one line of work that researchers within the ML for PDE community have studied is \"Sobolev Training\", where the main idea is to incorporate the Sobolev gradients of the training objective into the loss function. Examples of theoretical and empirical work in this direction include [1] (theory) and [2, 3, 4] (empirical studies). In fact, sobolev training can be interpreted as a form of operator preconditioning, which we will explain it briefly below: Consider the task of solving the PDE $\\mathcal{D}u = f$ over the domain $\\Omega$ with zero boundary condition via PINN, where $\\mathcal{D}$ is some differential operator. Let $E(u) :=\\frac{1}{2} \\|\\mathcal{D}u - f\\|_{L^2(\\Omega)}^2$ be the standard $L^2$ loss function, whose associated gradient flow is $u_t =  -\\frac{\\delta}{\\delta u}E(u) = -\\mathcal{D}^\\ast(\\mathcal{D}u-f)$. \n\nIn contrast, if we consider solving the same problem under the $H^2$ loss $E_{S}(u) := \\frac{1}{2} \\|(I-\\Delta)(\\mathcal{D}u - f)\\|_{L^2(\\Omega)}^2$  (Sobolev training), we have that new gradient flow is given by $u_t =  -\\frac{\\delta}{\\delta u}E_S(u) = -\\mathcal{D}^\\ast\\mathcal{P}^\\ast\\mathcal{P}(\\mathcal{D}u-f)$, where the differential operator $\\mathcal{P} = I - \\Delta$. This implies that training under a loss function that includes Sobolev gradient is equivalent to preconditioning the Hermitian operator $\\mathcal{D}^\\ast\\mathcal{D}$ via the differential operator $\\mathcal{P}$. Another approach for reducing the condition number is to reformulate the problem of minimizing the $L^2$ loss as a min max optimization problem, which was proposed in recent work [5]\n\nMoreover, the idea of using preconditioning/sobolev training to speed up optimization/training of a prespecified model in machine learning and scientific computing is not new either. A few examples of previous work include using preconditioning in neural network training [6, 7], application of Sobolev preconditioned GD in image processing [8, 9, 10] and graphics [11, 12], etc. Hence, the authors are encouraged to have a subsection discussing how this paper is related to previous work on the usage of preconditioning/sobolev training in ML for PDE and related fields. \n\n(2) Secondly, given that previous work [1, 2, 3, 4, 5] have already studied how sobolev training, which is equivalent to preconditioning, accelerates the optimization process of solving PDEs with PINNs, the reviewer thinks that this paper's novelty seems to be a bit limited. \nIt seems to me that sobolev training is a natural and simple method to implement (The main idea is to use integration by parts to kill any term involving the derivatives of the parametrized PDE solution $u_{\\theta}$). However, for the method proposed in this work, it seems hard to find a good preconditioner $P$ to reduce $\\kappa(A)$, especially when the differential operator $\\mathcal{D}$ and the neural network model $u_{\\theta}$ becomes a bit complicated. Consequently, to complement the numerical experiments exhibited in this paper, it would probably be meaningful to test the proposed method's performance on more complicated nonlinear PDEs, such as the Allen-Cahn equation, the Schrodinger equation and the Navier Stokes equation used in [13]. \n\n(3) Thirdly, the presentation of this paper can possibly be improved from a few aspects: Firstly, the review of ML-based/Physics-Informed PDE solvers in the second paragraph of Section 1 seems to be incomplete. A lot of recent work in this field seems to be missing. The authors are encouraged to take a look at the first part of Subsection 1.1 in [1] for a complete review of related work in this field. Secondly, it's probably better to move the subsection of related work from Section 4 to Section 1 (Introduction), as this might help the readers under the main contributions presented at the end of the introduction in a better way. It would also be helpful to include a subsection of the paper's organization and a subsection discussing the mathematical notations used in this paper within Section 1 (Introduction). Moreover, there seem to be some typos in the proofs presented in this paper. A non-exhaustive list of typos is given below:\n\n1. In (2.5), a factor of $\\frac{1}{2}$ is missing in front of the Hessian term $H_k$. Similarly, the $\\frac{1}{2}$ factor is also missing in the expansion of the two gradients $\\nabla_{\\theta}R(\\theta_k)$, $\\nabla_{\\theta}B(\\theta_k)$ and the error term $\\epsilon_k$ in Appendix A.1. \n\n2. Based on the definition of the vector $C$, the sign in front of $C$ in (2.7), (2.8) should be negative, i.e, $C$ should be replaced by $-C$. Also, the expression of $\\theta$ in Theorem 2.3 should be $\\theta = \\theta_0 - A^{-1}C$. Similarly, $\\eta \\epsilon_k$ should be replaced by $-\\eta \\epsilon_k$ in (2.7)\n\n3. In the proof of Lemma 2.1 presented in Appendix A.2, the term $\\epsilon_k$ should be replaced by $\\eta \\epsilon_k$ in the first line (A.2) of the proof. Therefore, the proof and the main result of Lemma 2.1 also need to be slightly modified.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "du1EznJn71", "reviewer": "Reviewer_Wj3M", "summary": "The paper mainly explores how the convergence of training in Physics-Informed Neural Networks (PINNs) is dependent on the condition number of the operators constituting the Partial Differential Equations (PDEs). From the perspective of novelty, the idea posited in this paper is quite interesting. In terms of conclusions, the paper articulately presents the theoretical result that a lower condition number leads to better convergence of PINNs. Moreover, a preconditioning method for linear PINNs is proposed in the paper, significantly reducing the condition number and accelerating the training process. \n\nOverall, I appreciate the theoretical analysis conducted in this paper, which I believe meets the submission standards for ICLR. However, the experimental section is somewhat weak (as noted in the drawbacks), and I believe that more extensive experiments should be conducted to validate the theoretical analysis presented in the paper. Additionally, some related works concerning hard constraints were overlooked, specifically the following two papers:\n\n**Reference**\n1. A Unified Hard-Constraint Framework for Solving Geometrically Complex PDEs (https://arxiv.org/abs/2210.03526)\n2. PFNN: A Penalty-Free Neural Network Method for Solving a Class of Second-Order Boundary-Value Problems on Complex Geometries (https://arxiv.org/abs/2004.06490)", "questions": "None.", "limitations": "", "rating": 5, "confidence": 5, "soundness": 4, "presentation": 2, "contribution": 3, "strengths": "1. The idea presented in this paper is novel, to the best of my knowledge.\n2. The paper's approach of using condition number analysis to study the convergence of PINNs is innovative, providing fundamental theoretical guarantees for PINN convergence.\n3. The writing of most sections is easy to understand. Overall, the theoretical analysis part of this paper is very engaging to read.", "weaknesses": "1. The experimental section of the paper is quite weak. It solely provides a method for reducing the condition number in linear PDEs using PINNs based on the theoretical analysis presented, which considerably narrows the applicability of the proposed method.\n2. Some parts of the paper are not clearly written, such as the calculation of the preconditioner matrix P, and the discussion regarding complexity appears to be vague.\n3. The logic in the latter experimental and discussion sections is somewhat jumbled. For instance, the last two sections, in my opinion, should be consolidated into one section as related work.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (poster)