An operator preconditioning perspective on training in physics-informed machine learning

1. (continued) Finally, we have added a new section SM Section C.3, where we consider another general simple preconditioning strategy which consists of using an approximate inverse of the matrix $\mathbb{A}$ (Eqn. (2.6)) as a preconditioner for the gradient descent step, i.e., pre-multiplying $(\mathbb{A} + \varepsilon I)^{-1}$ to the gradient descent update. This provides an additional general preconditioning strategy and we believe further strengthens the empirical aspects of our paper.
2. (W2+W3) Some parts of the paper are not clearly written and the last two sections, in my opinion, should be consolidated into one section as related work. We thank the reviewer for this excellent suggestion. In the process of incorporating all the 6 reviewers' multiple suggestions, we will consider a restructuring of these sections with your suggestion in mind in a CRV, if accepted.
3. We thank the reviewer for pointing out additional relevant references for hard BCs and will incorporate them in a CRV, if accepted.

We sincerely hope that we have addressed all the reviewer's concerns, particularly about the scope of empirical vs. theoretical contributions in our paper, to your satisfaction. We kindly request the reviewer to upgrade your assessment of our paper accordingly.

We thank the reviewer for appreciating the theoretical contributions of our paper. Below, we address the questions raised by the reviewer and thank the reviewer in advance for their patience in reading our reply.

1. (W1) Regarding the reviewer's concern that the The experimental section of the paper is quite weak, we would like to, at length, frame the contributions of our paper. To this end, we take the liberty to repeat the relevant portion of our reply to all the reviewers above here: At the outset, it is amply clear in the literature that a major issue with physics-informed machine learning models such as PINNs is the difficulty of training them. This has been widely recognised and several ad-hoc strategies have been proposed to improve training. However, it is also clear that a suitable theoretical framework that precisely formulates the problem of optimizing physics-informed loss functions was missing. We believe that a lack of precision in defining what the problem actually is also impedes finding suitable measures for ameliorating the problem itself.

   With this background in mind, our main aim in the paper was precisely to propose a suitable theoretical framework in which the issues affecting training of physics-informed machine learning models can be formulated and analyzed. To this end, we have rigorously connected the rate of convergence of gradient descent algorithms for physics-informed loss functions to the condition number of an operator, which involves the Hermitian-square of the underlying differential operator. This contribution is novel to the best of our knowledge and has been acknowledged as such by most of the reviewers. A corollary of this observation is the fact that preconditioning (now defined rigorously as reducing the condition number of the underlying operator) is necessary to accelerate training of physics-informed machine learning models. Again, this provides a mathematically precise definition of preconditioning for these models and stands in stark contrast to the existing literature where the words preconditioning and condition numbers are used very loosely and imprecisely. These two aforementioned aspects constitute our main contribution and has been emphasised as such repeatedly in our paper.

   Given the precise notion of preconditioning that we propose, it was incumbent upon us to also analyse several existing strategies for improving training in physics-informed machine learning through the perspective of our framework of operator preconditioning. To this end, we have in Section 3 and in the SM, analyzed via a combination of theoretical analysis and empirical evaluation, several such strategies (Choice of balance parameter $\lambda$, Hard Boundary conditions, Second-order optimizers, Domain decomposition and Casual learning, we will also include Sobolev Training) through the lens of operator preconditioning, providing a novel interpretation of these strategies and identifying conditions under which they might or might not work. We had gone further and analyzed an additional strategy that essentially boils down to rescaling Fourier features, whether with respect to linear basis functions or nonlinear neural networks. These are secondary contributions and we feel that it is not fair to judge the article based on secondary, rather than primary contributions.

   Moreover, the fact that only linear PDEs was considered in the empirical evaluation was based on the observation that if physics-informed machine learning models do not work for simple linear PDEs, there is little hope of them working for complicated nonlinear PDEs. Furthermore, we follow widely cited papers such as Krishnapriyan et al, 2021, which flagged the problems of training PINNs by focussing considerably on simple linear PDEs. Finally, the page limit also automatically limits the scope of an article. Given these considerations, we believe that linear PDEs already provide a reasonably representative arena for analyzing the issue of training PINNs and related models and lessons learnt there can be applied to nonlinear PDEs in future works. In particular, at some level, nonlinear PDEs have to be linearized and the lessons learnt from linear PDEs will be germane. Finally in this context, we would like to mention that there is a tradition of analyzing training in machine learning through linear models (see ArXiv:1810.02281 and references therein) as analyzing gradient descent even with linear models is highly non-trivial and representative of the difficulties faced in training non-linear models.

We thank the reviewer for their detailed feedback on our paper. Below, we address the questions raised by the reviewer and thank the reviewer in advance for their patience in reading our reply.

1. (W1): The reviewer has pointed to Missing related work. We would like to thank the reviewer for identifying very interesting related work, especially on Sobolev training and competitive PINNs that we will certainly include in a CRV, if accepted. Although we have taken care to mention previous work on preconditioning in the field of numerical analysis, we apologize having missed the stated work in ML and will include them in a CRV, if accepted.

2. (W2). The reviewer has addressed the issue of Sobolev Training and its relation to preconditioning at great length. We apologize for having missed this work on Sobolev training in discussing related work. Following the excellent suggestions of the reviewer, we have investigated Sobolev training in some detail and would like to explain our point of view on it below,

   • First, in our reading of references [1-4] mentioned by you, the Sobolev training refers to using $H^r$ norms with $r \geq 1$ for the PDE residual, rather than the $L^2$ norm, as the physics-informed loss function. These papers have also shown that using this loss can improve training. However, we were unable to see how exactly these papers consider preconditioning in a precise sense as it is unclear to us what exactly was the underlying condition number. Thus, there is a clear difference in our approach where we precisely define the condition number of the underlying operator that composes the Hermitian square of the differential operator with the Tangent Kernel operator.

   • It turns out that we can also analyze Sobolev training within the operator preconditioning framework that we have formulated in our paper. To this end, we revisit the example considered in Theorem 3.1 of our paper concerning the Laplacian in one space dimension as calculations can be explicitly performed in this case. We simply repeat the steps of the proof to account for the $H^r$-Sobolev norm. The matrix $\mathbb{A}$ is now given as $\mathbb{A}^* = \mathbb{D}^* + \lambda vv^\top$ with $\mathbb{D}^* := \sum_{\ell=0}^r \langle \mathcal{D} \partial_x^\ell \phi, \mathcal{D} \partial_x^\ell \phi\rangle_{L^2}$. One can then calculate that $\mathbb{D}^*$ is diagonal with $\mathbb{D}^*_{kk} = \sum_{\ell=0}^r k^{4+2\ell}$. An analogous calculation to the one in Theorem 3.1 then shows that the lower bound is $\kappa(\mathbb{A}^*) \geq K^4 \cdot \tfrac{1}{r+1}\sum_{\ell=0}^r K^{2\ell}$, which is clearly a worse lower bound than without Sobolev training, i.e., with $r=0$. For instance, for $H^1$ we find that $\kappa(\mathbb{A}^*) \geq K^6/2$ and for $H^2$ we find that $\kappa(\mathbb{A}^*) \geq K^8/3$, which is considerably worse.

   • We also ran numerical experiments, for both Fourier features as well as MLPs using Sobolev training. The comparative results with/without Sobolev training are provided in the Table below:

     Model	Epochs	Total loss	Relative $L^1$ error	Condition number
     ff_no_precond_no_sobolev	2000	0.01537	0.02386	31520.74
     ff_no_precond_sobolev	2000	3750.37	0.15199	1536899
     ff_precond_sobolev	20	$\approx 3.37 \times 10^{-32}$	$\approx 5.29 \times 10^{-16}$	1.16
     mlp_no_precond_no_sobolev	2000	0.05171	0.02653	-
     mlp_no_precond_no_sobolev	2000	5360540	1.29184	-

     Here, Sobolev training uses the $H^1$-norm. As we see from the table, for both Fourier features and MLPs, Sobolev training significantly worsens the training landscape as well as the resulting test errors and the results are consistent with the theory elaborated above. Hence, even for this simple problem, Sobolev training does not seem to pay off nor can it be viewed as a preconditioning technique, when preconditioning is defined as in our framework.

3. (continued) Given the precise notion of preconditioning that we propose, it was incumbent upon us to also analyse several existing strategies for improving training in physics-informed machine learning through the perspective of our framework of operator preconditioning. To this end, we have in Section 3 and in the SM, analyzed via a combination of theoretical analysis and empirical evaluation, several such strategies (Choice of balance parameter $\lambda$, Hard Boundary conditions, Second-order optimizers, Domain decomposition and Casual learning, we will also include Sobolev Training) through the lens of operator preconditioning, providing a novel interpretation of these strategies and identifying conditions under which they might or might not work. We had gone further and analyzed an additional strategy that essentially boils down to rescaling Fourier features, whether with respect to linear basis functions or nonlinear neural networks. These are secondary contributions and we feel that it is not fair to judge the article based on secondary, rather than primary contributions.

   In response to the reviewer's points about lack of general preconditioning strategies, we highlight Sec 3 where many existing strategies are examined through the lens of operator preconditioning. Moreover, in *SM Section C.3, we have additionally explored another general simple preconditioning strategy which consists of using an approximate inverse of the matrix $\mathbb{A}$ (Eqn. (2.6)) as a preconditioner for the gradient descent step, i.e., pre-multiplying $(\mathbb{A} + \varepsilon I)^{-1}$ to the gradient descent update. This strategy is also applicable to any PDE and any linear or nonlinear model.

4. (W3) We thank the reviewer for their careful inspection of our paper and have edited all typos pointed out by you. In addition, we will expand the bibliography following their suggestions in the CRV, if accepted.

5. (Q1) We have discussed the main difference between the preconditioning method proposed in this paper and sobolev training and Is there any numerical example that your method outperforms (higher accuracy, lower computational cost, etc) sobolev training? in some detail above.

We sincerely hope that we have addressed all the reviewer's concerns, particularly about the relation between our paper and the literature on Sobolev training, to your satisfaction. We kindly request the reviewer to upgrade your assessment of our paper accordingly.

2. (continued)

   • Next, the papers like [5, 8-10] also (essentially) suggest to use an $H^r$-norm during training but now with $r<0$, in contrast with [1-4]. We agree with the reviewer that this could be a valid method of preconditioning. Again, we analyze it from the preconditioning perspective suggested by us. Repeating the analysis for the Poisson equation above by applying $(I-\Delta)^{-1}$ to the operator $\mathcal{D}$, the matrix $\mathbb{A}$ is now given as $\mathbb{A}^* = \mathbb{D}^* + \lambda vv^\top$ with $\mathbb{D}^* := \langle (I-\Delta)^{-1}\mathcal{D} \phi, (I-\Delta)^{-1}\mathcal{D} \phi\rangle_{L^2}$. One can then calculate that $\mathbb{D}^*$ is diagonal with $\mathbb{D}^*_{kk} = k^4/(1+k^2)^2$. An analogous calculation to the one in Theorem 3.1 then shows that the lower bound is $\kappa(\mathbb{A}^*) \geq 4K^4/(1+K^2)^2$. This is a major improvement over the Fourier model without any preconditioning, but it can not achieve condition numbers close to 1 (like in Theorem 3.1) as for large $K$ the lower bound will be approximately 4.

     Nevertheless, we would like to point that in practice, computing Sobolev norms with negative orders (or computing $(I-\Delta)^{-1}(\mathcal{D} u_\theta-f))$ can be very hard, even of comparable difficulty as solving the PDE. Hence, it is unclear to us if Sobolev training with negative Sobolev norms is going to be useful in solving PDEs.

   Given that above discussion, we would like to emphasize that while Sobolev training might be a viable strategy in some contexts, we cannot support the reviewer's contention that Sobolev training is equivalent to preconditioning, let alone a preconditioning strategy that works in general. To function effectively as a preconditioner, it is crucial that the preconditioning approach adapts when dealing with different PDEs, as detailed in the paper. For instance, the operator effective in preconditioning the Poisson equation is unlikely to be suitable for the advection equation. This aspect is key and as Sobolev regularization is fixed it may not behave well for any given PDE.

   Hence, we suggest that we will certainly mention Sobolev training and include the suggested references as another example of a preconditioning strategy in a CRV, if accepted. If the reviewer so wishes, we are happy to include an elaborate discussion on Sobolev training with the above results in the CRV. We thank the reviewer again for pointing out this avenue which we believe further strengthens the contributions of our paper.

3. Regarding the reviewer's question about our contributions, we would like to explicitly point out the context and contributions of the paper. To this end, we take the liberty to repeat the relevant portion of our reply to all the reviewers above here: At the outset, it is amply clear in the literature that a major issue with physics-informed machine learning models such as PINNs is the difficulty of training them. This has been widely recognised and several ad-hoc strategies have been proposed to improve training. However, it is also clear that a suitable theoretical framework that precisely formulates the problem of optimizing physics-informed loss functions was missing. We believe that a lack of precision in defining what the problem actually is also impedes finding suitable measures for ameliorating the problem itself.

   With this background in mind, our main aim in the paper was precisely to propose a suitable theoretical framework in which the issues affecting training of physics-informed machine learning models can be formulated and analyzed. To this end, we have rigorously connected the rate of convergence of gradient descent algorithms for physics-informed loss functions to the condition number of an operator, which involves the Hermitian-square of the underlying differential operator. This contribution is novel to the best of our knowledge and has been acknowledged as such by most of the reviewers. A corollary of this observation is the fact that preconditioning (now defined rigorously as reducing the condition number of the underlying operator) is necessary to accelerate training of physics-informed machine learning models. Again, this provides a mathematically precise definition of preconditioning for these models and stands in stark contrast to the existing literature where the words preconditioning and condition numbers are used very loosely and imprecisely. These two aforementioned aspects constitute our main contribution and has been emphasised as such repeatedly in our paper.

We thank the reviewer for their detailed feedback on our paper. Below, we address the questions raised by the reviewer and thank the reviewer in advance for their patience in reading our reply.

1. (W1): The reviewer has pointed to Missing related work. We would like to thank the reviewer for identifying very interesting related work, especially on Sobolev training and competitive PINNs that we will certainly include in a CRV, if accepted. Although we have taken care to mention previous work on preconditioning in the field of numerical analysis, we apologize having missed the stated work in ML and will include them in a CRV, if accepted.

2. (W2). The reviewer has addressed the issue of Sobolev Training and its relation to preconditioning at great length. We apologize for having missed this work on Sobolev training in discussing related work. Following the excellent suggestions of the reviewer, we have investigated Sobolev training in some detail and would like to explain our point of view on it below,

   • First, in our reading of references [1-4] mentioned by you, the Sobolev training refers to using $H^r$ norms with $r \geq 1$ for the PDE residual, rather than the $L^2$ norm, as the physics-informed loss function. These papers have also shown that using this loss can improve training. However, we were unable to see how exactly these papers consider preconditioning in a precise sense as it is unclear to us what exactly was the underlying condition number. Thus, there is a clear difference in our approach where we precisely define the condition number of the underlying operator that composes the Hermitian square of the differential operator with the Tangent Kernel operator.

   • It turns out that we can also analyze Sobolev training within the operator preconditioning framework that we have formulated in our paper. To this end, we revisit the example considered in Theorem 3.1 of our paper concerning the Laplacian in one space dimension as calculations can be explicitly performed in this case. We simply repeat the steps of the proof to account for the $H^r$-Sobolev norm. The matrix $\mathbb{A}$ is now given as $\mathbb{A}^* = \mathbb{D}^* + \lambda vv^\top$ with $\mathbb{D}^* := \sum_{\ell=0}^r \langle \mathcal{D} \partial_x^\ell \phi, \mathcal{D} \partial_x^\ell \phi\rangle_{L^2}$. One can then calculate that $\mathbb{D}^*$ is diagonal with $\mathbb{D}^*_{kk} = \sum_{\ell=0}^r k^{4+2\ell}$. An analogous calculation to the one in Theorem 3.1 then shows that the lower bound is $\kappa(\mathbb{A}^*) \geq K^4 \cdot \tfrac{1}{r+1}\sum_{\ell=0}^r K^{2\ell}$, which is clearly a worse lower bound than without Sobolev training, i.e., with $r=0$. For instance, for $H^1$ we find that $\kappa(\mathbb{A}^*) \geq K^6/2$ and for $H^2$ we find that $\kappa(\mathbb{A}^*) \geq K^8/3$, which is considerably worse.

   • We also ran numerical experiments, for both Fourier features as well as MLPs using Sobolev training. The comparative results with/without Sobolev training are provided in the Table below:

     Model	Epochs	Total loss	Relative $L^1$ error	Condition number
     ff_no_precond_no_sobolev	2000	0.01537	0.02386	31520.74
     ff_no_precond_sobolev	2000	3750.37	0.15199	1536899
     ff_precond_sobolev	20	$\approx 3.37 \times 10^{-32}$	$\approx 5.29 \times 10^{-16}$	1.16
     mlp_no_precond_no_sobolev	2000	0.05171	0.02653	-
     mlp_no_precond_no_sobolev	2000	5360540	1.29184	-

     Here, Sobolev training uses the $H^1$-norm. As we see from the table, for both Fourier features and MLPs, Sobolev training significantly worsens the training landscape as well as the resulting test errors and the results are consistent with the theory elaborated above. Hence, even for this simple problem, Sobolev training does not seem to pay off nor can it be viewed as a preconditioning technique, when preconditioning is defined as in our framework.

2. (continued) Given the precise notion of preconditioning that we propose, it was incumbent upon us to also analyse several existing strategies for improving training in physics-informed machine learning through the perspective of our framework of operator preconditioning. To this end, we have in Section 3 and in the SM, analyzed via a combination of theoretical analysis and empirical evaluation, several such strategies (Choice of balance parameter $\lambda$, Hard Boundary conditions, Second-order optimizers, Domain decomposition and Casual learning) through the lens of operator preconditioning, providing a novel interpretation of these strategies and identifying conditions under which they might or might not work. We had gone further and analyzed an additional strategy that essentially boils down to rescaling Fourier features, whether with respect to linear basis functions or nonlinear neural networks. These are secondary contributions and we feel that it is not fair to judge the article based on secondary, rather than primary contributions.

   Regarding the reviewer's specific point about Spectral methods for the Poisson Equation, we would like to point that preconditioning for spectral methods is indeed well-established. However, there are essential differences between our approach and using a spectral method. Because, we minimize the loss in an $L^2$-sense, we are confronted with having to analyze the spectrum of the Hermitian-square of the underlying differential operator, where as in a conventional spectral method, only the differential operator itself has to be pre-conditioned. Following the reviewer's excellent suggestion, we will explicitly mention the differences with a spectral method in a CRV, if accepted.

   In response to the reviewer's points about lack of general preconditioning strategies, we highlight Sec 3 where many existing strategies are examined through the lens of operator preconditioning. Moreover, in *SM Section C.3, we have additionally explored another general simple preconditioning strategy which consists of using an approximate inverse of the matrix $\mathbb{A}$ (Eqn. (2.6)) as a preconditioner for the gradient descent step, i.e., pre-multiplying $(\mathbb{A} + \varepsilon I)^{-1}$ to the gradient descent update. This strategy is also applicable to any PDE and any linear or nonlinear model.

3. (Q3) Regarding the reviewer's question on Lemma 2.2, is there a way of controlling the error $\epsilon_k$ with respect to $m$ for large $m$? Our analysis for particular Lemma builds on the results from Wang et. al, which do not control the error in terms of $m$. However, a careful examination of the proof reveals that this precise control is likely to be possible, yet very tedious. We will attempt to perform this analysis and include it in a CRV, if accepted and thank the reviewer for pointing it to us.

4. (Q6) In response to the reviewer's question on Fig. 3 (left): I don't really see an improvement of the approach as there are still a large number of 0 eigenvalues with the MLP + preconditioned Fourier features approach, suggesting that the system is ill-conditioned, we would like to point out that he spectrum of the preconditioned model is much more localized, which is often correlated to better convergence as seen in Figure 3 (right).

5. We have addressed all the other minor concerns of the reviewer in the modified version of the paper.

We sincerely hope that we have addressed all the reviewer's concerns to your satisfaction. We kindly request the reviewer to upgrade your assessment of our paper accordingly.

We start by thanking the reviewer for their appreciation of the strengths of our theoretical and empirical strengths of our paper. Below, we address the questions raised by the reviewer and thank the reviewer in advance for their patience in reading our detailed reply.

1. (W1) Regarding the reviewer's point that they neglect all higher order terms in an error term $\epsilon_k$, which they assume to be small, we want to highlight that we do not just assume that this error term is small, but that rather we show that it is identically zero for linear models and that we prove rigorously in Lemma 2.2 that this error term is small in the NTK regime for neural networks. Furthermore, the reviewer asks Is that equivalent to performing a local analysis close to a local minima: we stress that this is not the setting of our analysis, but rather our analysis requires the smallness of the error term in Lemma 2.1. This error term is exactly zero for linear models and can be made small (Lemma 2.2) in the NTK regime for neural networks. There are no further assumptions. Hence, our analysis also covers the pre-asymptotic behavior of gradient descent as we explicitly consider the discrete time setting with finitely many parameters. This is precisely the setting where condition numbers make sense.

2. (W2) Regarding the reviewer's question on However, to perform preconditioning, one should instead derive an upper bound on the condition number of the Gram matrix, we fully agree with the reviewer that this general upper bound does not prove that a well-conditioned matrix implies a well-conditioned operator $\mathcal{A}\circ TT^*$. Fortunately, the equality holds for many practical cases. For linear basis functions the Gram matrix is trivially invertible, and also for randomly initialized neural networks (e.g. using normal distribution) it will be invertible as the $\phi_i$ will be independent functions with probability 1, and hence equality will also hold. Finally, in practice it is the conditioning of the matrix $\mathbb{A}$ that determines the convergence of gradient descent, and not directly that of the operator $\mathcal{A}\circ TT^*$, so the condition number of the matrix $\mathbb{A}$ is much more important in preconditioning in practice.

3. (W3) Regarding the reviewer's comment that The practical section 3 is the main weakness of the paper, we would like to explicitly point out the context and contributions of the paper. To this end, we take the liberty to repeat the relevant portion of our reply to all the reviewers above here: At the outset, it is amply clear in the literature that a major issue with physics-informed machine learning models such as PINNs is the difficulty of training them. This has been widely recognised and several ad-hoc strategies have been proposed to improve training. However, it is also clear that a suitable theoretical framework that precisely formulates the problem of optimizing physics-informed loss functions was missing. We believe that a lack of precision in defining what the problem actually is also impedes finding suitable measures for ameliorating the problem itself.

   With this background in mind, our main aim in the paper was precisely to propose a suitable theoretical framework in which the issues affecting training of physics-informed machine learning models can be formulated and analyzed. To this end, we have rigorously connected the rate of convergence of gradient descent algorithms for physics-informed loss functions to the condition number of an operator, which involves the Hermitian-square of the underlying differential operator. This contribution is novel to the best of our knowledge and has been acknowledged as such by most of the reviewers. A corollary of this observation is the fact that preconditioning (now defined rigorously as reducing the condition number of the underlying operator) is necessary to accelerate training of physics-informed machine learning models. Again, this provides a mathematically precise definition of preconditioning for these models and stands in stark contrast to the existing literature where the words preconditioning and condition numbers are used very loosely and imprecisely. These two aforementioned aspects constitute our main contribution and has been emphasised as such repeatedly in our paper.

At the outset, we would like to sincerely thank the reviewer for reading our response and providing constructive feedback on it. We are grateful to hear that most of the reviewer's concerns are addressed and we thank them for the clarification of the final question. Regarding Q6, we respond below,

First, we would like to make an important point regarding conditioning of the matrix A, the condition number, and its relationship to convergence of gradient descent. The reviewer is absolutely correct that--as a consequence of having near zero eigenvalues in both the MLP and MLP + preconditioned Fourier features case (as shown in Fig. 3) , the minimum and maximum eigenvalue are approximately the same, hence they yield approximately the same condition number. Although a smaller condition number is a sufficient condition for faster convergence, it is not a necessary condition, at least for the pre-asymptotic convergence. Indeed, if we decompose the difference $\widetilde{\theta_0}-\vartheta = \sum_i \alpha_i v_i$ along the eigenspaces of $\mathbb{A}$, we can rewrite equation SM (A.16) as follows: $\widetilde{\theta_k}-\vartheta = (I-\eta \mathbb{A})^k(\widetilde{\theta_0}-\vartheta) = \sum_{i=1}^n (1-\eta\lambda_i)^k \alpha_i v_i,$ where $\lambda_{\min} := \lambda_1 \leq \lambda_2 \leq \ldots \leq \lambda_n =: \lambda_{\max}$ the $n$ eigenvalues of $\mathbb{A}$ and $v_1, \ldots, v_n$ the corresponding eigenvectors. This is a refinement of Theorem 2.3. Hence by recalling that $\eta = c/\lambda_{\max}$ and developing the argument of SM Eqn (A.17) and (A.18), we observe that the more eigenvalues are closer to $\lambda_{\max}$, the more terms in the above sum are closer to 0 and the smaller $\Vert \widetilde{\theta_k}-\vartheta\Vert_2$ will be. In Figure 3 this is the case going from the MLP to the MLP + Preconditioned Fourier features, and even more so when considering the preconditioned Fourier features. We believe this finer grained perspective is very useful when thinking about convergence, and thus we will take care to discuss in the CRV of our paper, if it is accepted.

In addition, regarding the comment I would have expected a similar transformation of the spectrum like the blue eigenvalues to explain the difference in training on the right plot, we point to the fact that we only presented the spectrum resulting from preconditioning the Fourier features within the MLP in this particular figure and did not precondition the MLP as a whole. As a result, the improvement in the spectrum is not as large, and also not expected to be as large, as for the Fourier features without MLP, which is in complete correspondence with the observation that the preconditioned Fourier features (without MLP) achieve much lower errors and converge much faster than the MLP with preconditioned Fourier features (see Figure 1).

However, with your feedback, along with feedback from other reviewers we were able to push preconditioning further in the revised version, and have corroborated our findings in the fully non-linear setting by preconditioning by pre-multiplying $(\mathbb{A} + \varepsilon I)^{-1}$ to the gradient descent update (see Section C.3 of SM). We have now added an additional Figure 18 to the SM where we plot the spread of the eigenvalues and convergence of the loss, in the MLP and preconditioned-MLP case. In this case, the result is much clearer: although there is still quite a few (near) zero eigenvalues (the zero eigenvalues originate from the fact that the PINNs loss with an MLP is highly ill-conditioned), there are much more eigenvalues near $\lambda_{\max}$ for the preconditioned-MLP, and as a consequence of this and following the above argument, training is indeed much faster.

Needless to say, we will add these results and the corresponding discussion to the CRV, if accepted and thank the reviewer again for pointing out this avenue of improvement to us. We hope to have addressed your remaining concern satisfactorily and if so, we would kindly request the reviewer to update their assessment accordingly.

We thank the reviewer for their positive evaluation of our paper. Below, we address the questions raised by the reviewer and thank the reviewer in advance for their patience in reading our reply.

1. Regarding the reviewer's question about What to do if the kernel integral operator cannot be chosen manually: we agree that it is indeed often not possible to chose the kernel integral operator $TT^*$ explicitly. In this case, one can consider the various preconditioning techniques listed and analysed in the last paragraphs of Section 3 as well those in SM, Section C.3, which can always be applied.
2. The reviewer's sugggestion to add another optimization layer to learn the kernel integral operator is excellent. However, given the limited time for rebuttal, we will consider it as an interesting avenue for future research.

We sincerely hope that we have addressed all the reviewer's concerns to your satisfaction.

6. Regarding the reviewer's Q3 about non-convergence of unpreconditioned Fourier features for advection, we start by pointing out that this cannot only be attributed to the linear form of the underlying model as the MLP also barely converges (see SM Figure 11). Rather, the ill-conditioning of the model is responsible and only we precondition the linear model, it converges very fast, as acknowledged by the reviewer. The non-convergence of MLPs in this case was also observed in Krishnapriyan et. al., 2021.
7. We have addressed the other minor comments of the reviewer in the revised version.

We sincerely hope that we have addressed all the reviewer's concerns, particularly about use of linear vs. nonlinear models, to your satisfaction. We kindly request the reviewer to update your assessment of our paper accordingly.

We start by thanking the reviewer for their appreciation of the strength of our theoretical foundation, and the very fast convergence when one uses this theory to properly precondition models, as you are among the few to highlight these contributions which we believe are of great value to the community. Below, we address the questions raised by the reviewer and thank the reviewer in advance for their patience in reading our detailed reply.

1. (W1) Regarding the point that the paper assumes that the solution can be expressed as a linear combination of basis functions and the paper does not show how the preconditioning can be extended for non-linear models: we start by pointing out that for our theoretical analysis, we only assume that the tangent kernel of the model is (approximately) constant along the optimization path, which can also hold for nonlinear models as neural networks. Infact, in Lemma 2.1, a precise charecterization of this error has been provided. Moreover, in the last paragraphs of Section 3 we have discussed and analysed many different preconditioning approaches that also hold for nonlinear models. We have also added another simple preconditioning strategy which consists of using an approximate inverse of the matrix $\mathbb{A}$ (Eqn. (2.6)) as a preconditioner for the gradient descent step, i.e., pre-multiplying $(\mathbb{A} + \varepsilon I)^{-1}$ to the gradient descent update (see Section C.3 of SM) which is applicable to any nonlinear model.
2. (W2) Regarding the reviewer's concern that the paper only covers linear PDEs with linear PINN models, we refer the reviewer to our discussion and experiments using neural networks in Section 3 as well as SM section C, and hence we do consider nonlinear models. Moreover, the fact that only linear PDEs was considered in the empirical evaluation was based on the observation that if physics-informed machine learning models do not work for simple linear PDEs, there is little hope of them working for general nonlinear PDEs. Furthermore, we burrowed from widely cited papers such as Krishnapriyan et al, 2021, which flagged the problems of training PINNs and focus considerably on simple linear PDEs. Finally, the page limit also automatically limits the scope of an article. Given these considerations, we believe that linear PDEs already provide a reasonably representative arena for analyzing the issue of training PINNs and related models and lessons learnt there can be applied to nonlinear PDEs in future works. In particular, at some level, nonlinear PDEs have to be linearized and the lessons learnt from linear PDEs will be germane. Finally in this context, we would like to mention that there is a long tradition of analyzing training in machine learning through linear models (see ArXiv:1810.02281) and references therein) as analyzing gradient descent even with linear models is highly non-trivial and representative of the difficulties faced in training non-linear models.
3. Regarding the reviewer's point that only the training loss curves were shown for the two PDEs, we follow the reviewer's excellent suggestion and plot the resulting $L^2$-test errors in the SM Figures 17-19, where the tight coupling between the residual losses and the test errors are demonstrated.
4. Regarding reviewer's (Q1+MC2) about infinite Taylor series and ignoring higher-order terms, we wish to clarify that we did not consider an infinite Taylor series and then ignored the higher-order terms. Instead, we considered a first-order Taylor polynomial with the corresponding exact remainder term (this is where the Hessian comes in and why it is evaluated at an interpolated value due to the mean-value theorem). We again stress that this formula is exact and no terms were ignored herein. We later show that the remainder term is small under suitable conditions (Lemma 2.2).
5. Regarding reviewer's (Q2) on the relation of our analysis with NTK theory, we would like to emphasize that NTK theory was originally formulated for the supervised MSE ($L^2$)-loss and for infinite-width neural networks. In contrast, we consider a physics-informed loss and finite-width neural networks, resulting in significant departures from the conventional NTK framework. Compared to other works that use NTK theory to analyse PINNs e.g. (Wang et al, 2022), we want to note that our analysis is valid in the practical case of discrete time rather than continuous time (as in Wang et. al.) and also we consider a finite number of parameters (corresponding to reality), along with the fact that we rigorously keep track of all error terms. Discrete time and finite parameters are absolutely crucial to our analysis as the condition number only appears in this discrete setting. Moreover, to the best of our knowledge, Wang et. al. do not identify an underlying operator and its conditioning as the key issues in training physics-informed machine learning models.

2. (continued) Moreover, the fact that only linear PDEs was considered in the empirical evaluation was based on the observation that if physics-informed machine learning models do not work for simple linear PDEs, there is little hope of them working for complicated nonlinear PDEs. Furthermore, we follow widely cited papers such as Krishnapriyan et al, 2021, which flagged the problems of training PINNs by focussing considerably on simple linear PDEs. Finally, the page limit also automatically limits the scope of an article. Given these considerations, we believe that linear PDEs already provide a reasonably representative arena for analyzing the issue of training PINNs and related models and lessons learnt there can be applied to nonlinear PDEs in future works. In particular, at some level, nonlinear PDEs have to be linearized and the lessons learnt from linear PDEs will be germane. Finally in this context, we would like to mention that there is a long tradition of analyzing training in machine learning through linear models (see ArXiv:1810.02281 and references therein) as analyzing gradient descent even with linear models is highly non-trivial and representative of the difficulties faced in training non-linear models.

   Finally, we have added a new section SM Section C.3, where we consider an additional and very general simple preconditioning strategy which consists of using an approximate inverse of the matrix $\mathbb{A}$ (Eqn. (2.6)) as a preconditioner for the gradient descent step, i.e., pre-multiplying $(\mathbb{A} + \varepsilon I)^{-1}$ to the gradient descent update.

3. Following the reviewer's excellent suggestion, we have now included **wall-clock times for the methods in SM Table 2. The reviewer's suggestion to discuss Refs. [2] and [3] are well-taken and would be included in a CRV, if accepted. In this context, we would like to comment that the contributions of these papers can be viewed through a preconditioning perspective as done for the contributions of your reference [1] as described in SM section B.2.

4. Regarding the reviewer's concern about it is less convincing that the paper makes significant contributions in ML or DL perspectives, please refer to our detailed framing of the contributions as we have outlined them to you in point 2. above.

5. Regarding the reviewer's question on how the gradient descent rule has been reimplemented following Eq. (3.2), we would like to point out that we have shown that rescaling the model parameters is equivalent to preconditioning the gradient directly. Hence, in our experiments with linear models or MLP with Fourier features, we rescaled our model parameters and therefore did not use (3.2) directly. On the other hand, with the matrix inversion strategy of SM Sec C.3, we use a form of Eqn (3.2) for the gradient update.

We sincerely hope that we have addressed all the reviewer's concerns, particularly about the framing of our contributions, to your satisfaction. We kindly request the reviewer to upgrade your assessment of our paper accordingly.

We thank the reviewer for acknowledging the importance of the problem our paper addresses, and for their other positive comments. Below, we address the questions raised by the reviewer and thank the reviewer in advance for their patience in reading our detailed reply.\\

1. Regarding the reviewer's comment about us explicitly mentioning preconditioning A for nonlinear models such as neural networks can, in general, be hard, we start by clarifying what me meant in this context, namely that, in the neural network setting the operator $TT^*$ appearing in the main theorem 2.4 does not reduce to a simple operator like it does in the case of Fourier features or other linear basis functions. This implies that the condition number of the operator $\mathcal{A} \circ TT^*$ is more difficult to study, but by no means is it impossible as we consider in Section 3 as well as in SM Section C.

2. Regarding the reviewer's concern that It does not seem that the paper provides some guidelines or insights for preconditioning strategies for general cases and the question how to build preconditioner for more general cases, we start by repeating our detailed reply to all the reviewers providing the context of our contributions. At the outset, it is amply clear in the literature that a major issue with physics-informed machine learning models such as PINNs is the difficulty of training them. This has been widely recognised and several ad-hoc strategies have been proposed to improve training. However, it is also clear that a suitable theoretical framework that precisely formulates the problem of optimizing physics-informed loss functions was missing. We believe that a lack of precision in defining what the problem actually is also impedes finding suitable measures for ameliorating the problem itself.

   With this background in mind, our main aim in the paper was precisely to propose a suitable theoretical framework in which the issues affecting training of physics-informed machine learning models can be formulated and analyzed. To this end, we have rigorously connected the rate of convergence of gradient descent algorithms for physics-informed loss functions to the condition number of an operator, which involves the Hermitian-square of the underlying differential operator. This contribution is novel to the best of our knowledge and has been acknowledged as such by most of the reviewers. A corollary of this observation is the fact that preconditioning (now defined rigorously as reducing the condition number of the underlying operator) is necessary to accelerate training of physics-informed machine learning models. Again, this provides a mathematically precise definition of preconditioning for these models and stands in stark contrast to the existing literature where the words preconditioning and condition numbers are used very loosely and imprecisely. These two aforementioned aspects constitute our main contribution and has been emphasised as such repeatedly in our paper.

   Given the precise notion of preconditioning that we propose, it was incumbent upon us to also analyse several existing strategies for improving training in physics-informed machine learning through the perspective of our framework of operator preconditioning. To this end, we have in Section 3 and in the SM, analyzed via a combination of theoretical analysis and empirical evaluation, several such strategies (Choice of balance parameter $\lambda$, Hard Boundary conditions, Second-order optimizers, Domain decomposition and Causual learning) through the lens of operator preconditioning, providing a novel interpretation of these strategies and identifying conditions under which they might or might not work. We had gone further and analyzed an additional strategy that essentially boils down to rescaling Fourier features, whether with respect to linear basis functions or nonlinear neural networks. These are secondary contributions and we feel that it is not fair to judge the article based on secondary, rather than primary contributions.

(see above for the first part )

Moreover, the fact that only linear PDEs was considered in the empirical evaluation was based on the observation that if physics-informed machine learning models do not work for simple linear PDEs, there is little hope of them working for complicated nonlinear PDEs. Furthermore, we follow widely cited papers such as Krishnapriyan et al, 2021, which flagged the problems of training PINNs by focussing considerably on simple linear PDEs. Finally, the page limit also automatically limits the scope of an article. Given these considerations, we believe that linear PDEs already provide a reasonably representative arena for analyzing the issue of training PINNs and related models and lessons learnt there can be applied to nonlinear PDEs in future works. In particular, at some level, nonlinear PDEs have to be linearized and the lessons learnt from linear PDEs will be germane. Finally in this context, we would like to mention that there is a long tradition of analyzing training in machine learning through linear models (see ArXiv:1810.02281) and references therein) as analyzing gradient descent even with linear models is highly non-trivial and representative of the difficulties faced in training non-linear models.

With these general comments, we hope to have addressed a common point of criticism of several of the reviewers and look forward to a healthy discussion.

Yours sincerely,

Authors of An operator preconditioning perspective on training in physics-informed machine learning