Neural network augmented inverse problems for PDEs
Jens Berg Kaj Nyström
stat.ML math.NA

















outdir=./







 
equationsection

comment
redKommentar


theoremTheorem[section]
proposition[theorem]Proposition
corollary[theorem]Corollary
lemma[theorem]Lemma

remark
remarkRemark[section]

definition
definitionDefinition

definition
*definition*Definition





 
Department of Mathematics, Uppsala University

S-751 06 Uppsala, Sweden





Neural network augmented inverse problems for PDEs
    Kaj Nyströmkaj.nystrom@math.uu.se
    
==================================================




In this paper we show how to augment classical methods for inverse problems with artificial neural networks. The neural network acts as a prior for the coefficient to be estimated from noisy data. Neural networks are global, smooth function approximators and as such they do not require explicit regularization of the error functional to recover smooth solutions and coefficients. We give detailed examples using the Poisson equation in 1, 2, and 3 space dimensions and show that the neural network augmentation is robust with respect to noisy and incomplete data, mesh, and geometry.


 



§ INTRODUCTION

In this paper we study the classical coefficient approximation problem for partial differential equations (PDEs). The inverse problem consists of determining the coefficient(s) of a PDE given more or less noisy measurements of its solution. A typical example is the heat distribution in a material with unknown thermal conductivity. Given measurements of the temperature at certain locations, we are to estimate the thermal conductivity of the material by solving the inverse problem for the stationary heat equation.

The problem is of both practical and theoretical interest. From a practical point of view, the governing equation for some physical process is often known, but the material, electrical, or other properties are not. From a theoretical point of view, the inverse problem is a challenging often ill-posed problem in the sense of Hadamard. Inverse problems have been studied for a long time, starting with <cit.> in the 40's, <cit.> and <cit.> in the 60's, to being popularized by <cit.> in the 70's by their work on regularization techniques. Today there is a vast amount of literature available on inverse, and PDE constrained optimization problems in general, and we refer the reader to the many survey papers and textbooks, and the references therein, for an overview. See for example <cit.>.

As the PDE cannot be solved analytically in general, a numerical method is required. A common choice is the finite element method (FEM) which have been frequently used with success for a large class of inverse problems. See for example <cit.>. Recently, there has been a lot of activity in a branch of numerical analysis known as probabilistic numerics, with successful applications in inverse problems and PDE constrained optimization in general. The probabilistic numerical methods are similar, yet different, to the approach taken in this paper. In the probabilistic setting, one places a statistical prior, typically a Gaussian prior, on the unknown coefficient and uses Bayesian inference to compute statistics of the unknown coefficient. As more data is used in the inference, the variance of the coefficient is reduced and one ends up with a confidence interval for the true coefficient. See for example <cit.>.

In this work, we use the finite element method as supplied by the software package  <cit.> to solve the underlying PDE. To compute the gradient we need for the PDE constrained optimization, we use the software  <cit.> which works together with  to automatically derive and solve the adjoint PDE. Finally, we augment  and  with an artificial neural network to represent the unknown coefficient in the PDE. The neural networks we consider in this paper are simple feedforward neural networks with sigmoid activation functions in the hidden layers, and linear activations in the output layer. Such a neural network defines a smooth mapping ℝ^N →ℝ which can approximate any continuous function to arbitrary accuracy <cit.>. The similarity with the probabilistic method lies in that a model for the coefficient is selected a priori, and as more data is analyzed, the model improves. The probabilistic method uses a statistical non-parametric model, and here we use neural networks which is a deterministic parametric prior.

The model problem we are working with in this paper is the scalar Poisson equation in a subset Ω⊂ℝ^N, N = 1, 2, 3, given by

    -∇· (q(x) ∇ u(x)) = f(x),

subject to Dirichlet boundary conditions. The inverse problem is to compute the coefficient q, given a series of more or less noisy measurements û_i of u_i, for i = 0, …, M. It is clear that the inverse problem is ill-posed whenever ∇ u = 0 on an open set, as then there is no information available about the coefficient q. Moreover, well-posedness of the inverse problem depends on the norm under consideration. For example, consider the Poisson equation in one space dimension, x ∈ (0, π), q = 1/2, u = x^2, and a sequence of coefficients and observations q̂_N = (2 + cos(Nx))^-1, û_N = u + δ_N with

    δ_N = x/Nsin(Nx) + 1/N^2cos(Nx).

In the L^∞-norm given by

    ||u||_∞ = max_x ∈Ω |u(x)|

we have

    ||û_N - u||_∞≤c/N→ 0,        ||q̂_N - q||_∞ = 1/2↛0,

and we can clearly see that q is not a continuous function of the data in the sense of Hadamard <cit.>. One can, however, show that this problem is well-posed in the Sobolev norm H^1. However, as PDE constrained optimization in the H^1-norm is in general not practically applicable, since one does not have measurements of ∇û in general, we will consistently use the standard L^2-norm.



§ THE GENERAL INVERSE PROBLEM

We consider stationary partial differential equations of the general form

    F(u, q)    = 0,        x ∈Ω, 
    
    Bu    = g,        x ∈∂Ω,

where Ω⊂ℝ^N is a bounded domain and ∂Ω its boundary. Moreover, u = u(x) is the solution variable, B is a boundary operator, g is the boundary data, and q = q(x) is an a priori unknown parameter.

The deviation between u and û is defined by the error functional

    J(u, q) = 1/2||u - û||^2 = 1/2∫|u - û|^2dx,

and the goal is to compute q^* such that

    q^* = min_q J(u, q)

subject to u satisfying (<ref>). To efficiently solve the minimization problem (<ref>), we require the gradient of the error functional (<ref>) with respect to q. To compute the gradients, we will use the adjoint approach.



 §.§ The adjoint equations

Note that the parameter q defines the mapping

    q ↦ u(q),

which can be computed by solving (<ref>) by analytical or numerical means. This turns (<ref>), and (<ref>) into pure functions of q,

    F(u, q)    ↦ F(u(q), q), 
    
    J(u, q)    ↦ J(u(q), q),

and the PDE constrained problem (<ref>) becomes an unconstrained problem as the constraint is built-in to the solution of (<ref>). The total derivative of the error functional can be computed as

    J(u(q), q)/ q = ∂ J(u(q), q)/∂ u(q)∂ u(q)/∂ q + ∂ J(u(q), q)/∂ q,

where ∂ J(u(q), q)/∂ u(q) and ∂ J(u(q), q)/∂ q are in general easily computed. To compute the troublesome solution Jacobian ∂ u(q)/∂ q, we formally differentiate the PDE itself to obtain

    F(u(q), q)/ q = ∂ F(u(q), q)/∂ u(q)∂ u(q)/∂ q + ∂ F(u(q), q)/∂ q,

or equivalently

    -∂ F(u(q), q)/∂ u(q)∂ u(q)/∂ q = ∂ F(u(q), q)/∂ q.

The relation (<ref>) is known as the tangent linear system related to the functional (<ref>). In the rest of the derivations, we will write F(u(q), q) = F, J(u(q), q) = J and u(q) = u to ease the notation. The tangent linear operator ∂ F/∂ u is the linearization of the differential operator around the solution u which acts on the solution Jacobian ∂ u/∂ q. In the case of a linear PDE, the differential operator in (<ref>) and the tangent linear operator coincides. Assuming that the tangent linear system is invertible, we can use (<ref>) and (<ref>) to solve

    ∂ u/∂ q = -( ∂ F/∂ u)^-1∂ F/∂ q.

Substituting (<ref>) into (<ref>) and taking the Hermitian transpose gives

    J^*/ q = -∂ F^*/∂ q( ∂ F^*/∂ u)^-1∂ J^*/∂ u + ∂ J^*/∂ q.

We define the adjoint variable λ by

    λ = -( ∂ F^*/∂ u)^-1∂ J^*/∂ u,

or equivalently

    -∂ F^*/∂ uλ = ∂ J^*/∂ u.

The above equation (<ref>) is the adjoint, or dual, equation associated with the forward problem (<ref>) and the error functional (<ref>). By solving the adjoint equation (<ref>) and substituting into (<ref>), we can compute the total derivative of the error functional as

    J/ q = λ^* ∂ F/∂ q + ∂ J/∂ q.

With the total derivative of the error functional given by (<ref>), any gradient based optimization method can be used to solve the minimization problem (<ref>).



§ THE NUMERICAL METHOD

In general, the PDE (<ref>) cannot be solved analytically and hence neither the mapping (<ref>) nor the gradient (<ref>) can be computed. In practice, a numerical method is used to solve both the PDE (<ref>) and the adjoint equation (<ref>). A common method is to use FEM. To effectively use FEM, one selects a finite element space V for the solution, and another space Q for the coefficient. Typically, Q is the space of piecewise constants to reduce the number of degrees of freedom. As the space V consist of local basis functions, one needs regularization of the error functional to reconstruct a smooth coefficient. Typically, one adds Tikhonov regularization to the error functional of the form

    J̃ = 1/2∫ |u - û|^2 dx + α/2∫ |q|^2 dx,

where α > 0 is a regularization parameter. If, for some reason, an estimate of the coefficient and noise level is known a priori, one can further improve the results by adding generalized Tikhonov regularization of the form

    J̃ = 1/2∫ |u - û|^2 dx + α^δ/2∫ |q - q_*|^2 dx,

where q_* is the a priori estimated coefficient and the regularization parameter α^δ depends on the noise level. The actual value of regularization parameters α and α^δ depends on the problem at hand, and finding an optimal value is a non-trivial task. The optimal generalized Tikhonov regularization parameter α^δ can be computed by using, for example, the Morozov discrepancy principle or heuristic L-curves. These methods are, however, of limited practical use in a PDE constrained optimization context as they become extremely expensive. We will, however, discuss optimal regularization of the 1D example in appendix <ref>.

In this section, we do not use any regularization in the error functional (<ref>). We will use only the most basic settings and avoid any fine tuning of hyperparameters, or any other parameters, to get as much of a black box solution strategy as possible. This means that the comparisons we make with pure FEM are not entirely fair as finding the optimal function spaces and regularizations are beyond the scope of this paper. The main purpose of this paper is to present neural networks as possible augmentations to classical methods for inverse problems. Adding explicit regularization will improve the neural network augmentation as well as can be seen in appendix <ref>.



 §.§ The finite element method and automatic differentiation

In this work, we have used the FEM software  to solve the PDE (<ref>), together with  to automatically derive and solve the adjoint equation (<ref>). The automatic differentiation in  works by overloading the solver functions in  and recording a graph of solver operations. Each node in the graph has an associated gradient operation and the total derivative can be computed by reversed mode differentiation, i.e. backpropagation. The reversed mode differentiation has the benefit that the gradient computation is exact with respect to the FEM discretization, as long as it is smooth enough. The only source of error in the computation of (<ref>) is thus in the numerical solution of the forward problem[At least as long as the adjoint equation can be solved by a direct method. For large enough systems, an iterative method must be used which introduces additional errors. These are, however, usually small compared to the discretization errors of the forward problem.].



 §.§ Neural network representation of the coefficient

Rather than representing the unknown coefficient in a finite element space, we can represent the coefficient by a feedforward artificial neural network. That is, we let q = q(x; W, b) be parameterized by the weights and biases of a neural network, here denoted by W and b. This approach yield some immediate benefits. A neural networks is a global, smooth, function approximator which can approximate any continuous function to arbitrary accuracy by having a sufficient amount of parameters <cit.>. Since it is a global function, it can be cheaply evaluated anywhere in the domain without first searching for the correct mesh element and performing interpolation, which is an expensive operation.

The coefficient q is computed by feed forwarding

    q(x)    = σ(z^L) 
    
    z^L    = W^Lσ_L-1(z^L-1) + b^L 
    
    z^L-1   = W^L-1σ_L-2(z^L-2) + b^L-1
       [=]⋮
    
    z^2    = W^2σ_1(z^1) + b^2 
    
    z^1    = W^1x + b^1,

where W^l_ij = w^l_ij are the weight matrices connecting layers l and l-1 and b^l are the vectors of biases for each layer, and σ_l is the activation function of layer l. To determine weights and biases, we seek W^*, b^* such that

    W^*, b^* = min_W, b J.

As there is one extra layer of parameters, we need to apply the chain rule once more to compute the total derivative of the error functional. Let p denote any of the weight or biases parameters, then

    J/ p = ∂ J/∂ u∂ u/∂ q∂ q/∂ p + ∂ J/∂ q∂ q/∂ p = ( ∂ J/∂ u∂ u/∂ q + ∂ J/∂ q) ∂ q/∂ p.

The first factor on the right hand side of (<ref>) is computed as before by using  and . The last factor is the gradient of the network output with respect to the network parameters, which can be computed exactly by using backpropagation or automatic differentiation.

Note that the neural network augmentation is not restricted to FEM. Any numerical method which efficiently can solve the forward problem (<ref>), and the adjoint problem (<ref>), can be used in combination. The neural network is simply a prior for the coefficient with good function approximation properties, and efficient means for computing its gradient. The application of the chain rule in (<ref>) factors the discretization of the forward and backward problems from the computation of the gradient of the network. As a side effect, the discretization of the PDE can be refined without increasing the number of parameters in the neural network. In the simplest FEM case, the coefficient is represented as a constant on each mesh element. The number of optimization parameters is then equal to the number of mesh elements, which quickly become infeasible if the mesh is refined.

In the numerical examples which follow, we have consistently used a neural network design with as few parameters as possible. A low number of parameters is beneficial both for the numerical optimization as well as for the the implicit regularization. The designs are still, however, based on heuristics and experience. We found during this work that surprisingly small networks perform better than larger networks as long as the number of parameters is sufficient to approximate the coefficient. Once a sufficient design has been found, the results start to deteriorate with increasing network complexity due to overfitting and less efficient implicit regularization.



§ MODERATELY ILL-POSED EXAMPLES
 
In the numerical examples, we have used an exact solution and coefficient to compute a forcing function in (<ref>). The exact solution is used to generate the data in the error functional. In the case of added noise, we add a normally distributed value δ r, r ∈𝒩(0, 1), for some noise level δ, to each of the interior data points. The data on the boundary is always exact.

The results differ somewhat depending on the noise and random initialization of the network's weights. In the figures and tables, we use 's random number generators with , for reproducibility, to add noise and initialize the weights. We use the same initial guess for the coefficient and settings for the optimizer in both the FEM and network cases.

The equation we use in the examples is the Poisson equation. The equation requires that the coefficient to be estimated is positive. To ensure that the coefficient is initially positive, we initialize the network's weights from a uniform distribution 𝒰[0, 1) and all the biases are initialized to zero. Note that the endpoint is excluded by default in 's random number generators. It is of no practical importance.

The examples which follow are called moderately ill-posed since the coefficients to be estimated are smooth and the measurement data is available in the whole domain. The purpose of the examples is to show the applicability and ease of use of the method due to the implicit regularization by the neural network prior.



 §.§ One-dimensional Poisson equation
 
The one-dimensional Poisson equation with Dirichlet boundary conditions is given by

    -(qu_x)_x    = f,        x ∈ (0, 1), 
    
    u(0)    = g_0,        u(1) = g_1,

where f, g_0,1 are known, and q = q(x) : ℝ→ℝ_+ is the coefficient to be estimated. We discretize the domain Ω into M=101 continuous piecewise linear elements where the solutions are represented. We solve the minimization problem using both standard FEM, with the coefficient represented in the solution space, as well as represented by a neural network. Here, we use the BFGS <cit.> optimizer from  <cit.> as the number of optimization parameters is rather small. We iterate until the norm of the gradient of the error functional is less than 10^-6.

The neural network is a tiny feedforward network with one input neuron, one linear output neuron, and one hidden sigmoid layer with 3 neurons. We choose the solution u to be

    u(x) = sin^2(2π x)

and we will compute the inverse problem for a few different exact coefficients q̂ and noise levels. The neural network has 10 parameters which will be optimized, compared to the FEM problem which has 101 (same as the number of elements). We consider the exact coefficients q̂(x) = 1, q̂(x) = 1+x, q̂(x) = 1+x^2, and q̂(x) = 1 + 0.5 sin(2 π x) with noise level δ=0 and δ=5*10^-2. The results can be seen in Figures <ref>–<ref> and Table <ref>. It is clear that the network representation outperforms FEM in this case. The network representation is insensitive to noise and always produces smooth solutions and coefficients. FEM does not converge to smooth solutions and coefficients due to the lack of regularization in the error functional. We can also see that the number of iterations increase in the network case as the coefficient becomes more oscillatory. A deeper network with higher capacity might help to mitigate the increase as was seen in <cit.>.

For one-dimensional problems the number of mesh elements is usually rather limited. This means that the network is at a high risk of overfitting. In particular with noisy data, a network with too high capacity will fit the noise which causes a severe increase in the number of optimization iterations, and of course also lack of generalization. This is the reason we have used a tiny network with only three neurons in the hidden layer. A larger network would require weight regularization, dropout and/or  any other method which reduces overfitting. See for example <cit.>.










1.2




 §.§ Two-dimensional Poisson

The higher-dimensional Poisson equation with Dirichlet boundary conditions is given by

    -∇· (q ∇ u)    = f,        x ∈Ω, 
    
    u    = g,        x ∈∂Ω,

where Ω⊂ℝ^N is the computational domain, ∂Ω its boundary. The functions f, g are known, and q = q(x, y): ℝ^N →ℝ_+ is the a priori unknown coefficient to be estimated. Here, we choose the solution to be

    u = sin(π x) sin(π y),

for the exact coefficients q̂(x, y) = 1, q̂(x, y) = 1 + x + y, q̂(x, y) = 1 + x^2 + y^2, and q̂(x, y) = 1 + 0.5 sin(2 π x) sin(2 π y). The network is a single hidden layer feedforward network with 10 neurons in the hidden layer. The optimization is performed using BFGS and we iterate until the norm of the gradient of the error functional is less than 10^-7.



  §.§.§ Unit square with uniform mesh

Here, we let Ω = [0 ,1] × [0, 1] be the unit square discretized by 101 × 101 piecewise linear elements elements.  The coefficient represented by the network has 41 parameters, while it has 10201 when represented in the finite element space. Due to the very different nature of the optimization problems, a direct comparison is no longer meaningful. The BFGS method scales quadratically in complexity with the number of optimization parameters, and is only useful for small-scale optimization problems, where it is very efficient. The FEM minimization problems can instead be solved by, for example, the memory limited BFGS method <cit.>. By doing so, we see the same pattern as in the 1D case. FEM does not converge to smooth solutions or coefficients with added noise due to the lack of regularization.

We plot the difference between the exact and computed coefficient in Figures <ref>–<ref> and summarize the results in Table <ref>. We can see the same pattern as in the 1D case. The neural network representation is insensitive to noise and we always recover a smooth coefficient even without regularization. Note that the errors are larger at the corners of the domain. This is probably due to the fact that the optimization is performed in the L^2 norm which allows larger errors in isolated points. The errors can probably be reduced by performing the optimization under the H^1 norm instead. See for example <cit.> and references therein.









1.2




  §.§.§ Mesh independent convergence

It was shown in <cit.> that the number of optimization iterations grows polynomially with the ratio between the volumes of the smallest and largest mesh elements, h_max/h_min. The reason is that most gradient based optimization libraries do not take the underlying function space inner product into account when computing the step size and direction. As the neural network augmentation is mesh independent, we expect that the number of iterations remains constant when locally refining the mesh, as long as the mesh is fine enough to accurately compute the solution of the forward problem.

To test the mesh independent convergence, we consider a simple test case. We solve the inverse problem for the Poisson equation (<ref>), on the unit square, with coefficient q̂ = 1 and the exact solution given by <ref> without any added noise. We locally refine the mesh in concentric circles around the center of the domain with radius r(n) = 0.5/n, for n = 1, …, 7, as can be seen in Figure <ref>. The results are summarized in Table <ref>. We can clearly see that the number of iterations required until convergence as the mesh is refined becomes constant, and a stable minimum is found.


1.2




 §.§ Three-dimensional Poisson

As a final example, we consider the three-dimensional Poisson equation (<ref>) on a complex geometry from <cit.>. We take the exact coefficient q̂ = 1 and the analytical solution

    u = sin(π x) sin(π y) sin(π z)

with no noise and compute the inverse problem. We use a network with 10 neurons in the single hidden layer which gives a total of 51 optimization parameters. Convergence was achieved after 41 BFGS iterations with gradient norm tolerance 10^-7, which is consistent with the result in Table <ref>. This indicates that the neural network augmentation method is also independent of the geometry. The domain and difference between the exact and computed network coefficient can be seen in Figure <ref>.




The extra BFGS iteration required in the 3D case is because the forward problem could not be solved by a direct method due to memory requirements. We used the GMRES iterative method with an incomplete LU preconditioner to solve the forward problem. The solution to the forward problem is thus slightly less accurate due to the tolerance settings of the iterative method, which in turn means that the computed error functional gradients are also less accurate and more BFGS iterations are required.




§ SEVERELY ILL-POSED EXAMPLES
 
In the previous section we discussed moderately ill-posed problems where we have measurements of the smooth coefficients in the whole domain. Here we will show some examples of severely ill-posed problems where we have discontinuous coefficients or incomplete measurements.



 §.§ 1D discontinuous coefficient
 
Here we repeat the examples in section <ref> for the Poisson problem (<ref>) with

    q = 
    0.5,     0 ≤ x < 0.5 
    
    1.5,     0.5 ≤ x ≤ 1
    
    ,       
    f = 10,        g_1 = g_2 = 0.

In this case, the analytical solution is not known and we compute the reference solution and data using FEM. We used the same small network with a single hidden layer with 3 neurons and the sigmoid activation function. The results can be seen in Figure <ref> and Table <ref>. The network captures the discontinuity perfectly for noiseless data, and the smoothing effect of the implicit regularization is clearly seen when noise is added.


1.2




 §.§ 1D incomplete data
 
We repeat the examples in section <ref> for the constant coefficient case, q=1, under the assumption that we only have access to measurement data at a distance d from the boundaries. In this case, the error functional (<ref>) we want to minimize is reduced to

    J(u, q) = 1/2∫_0^d |u - û|^2 dx + 1/2∫_1-d^1 |u - û|^2 dx.

We let, as before, the neural network representing the unknown coefficient have a single hidden layer with 3 neurons and the sigmoid activation function. The result can seen in Figures <ref>–<ref> for d=0.4, 0.2, and 0.1, respectively. We can clearly see that the neural network is able to reconstruct the coefficient and solution in the whole domain despite the data being known only close to the boundaries.









 §.§ Remark on multiple discontinuities and multiscale coefficients

The neural network augmentation without explicit regularization works well when there are only a few discontinuities of the same size. As the number of discontinuities grow, or the scales of the coefficient start to differ too much, the low capacity networks are no longer sufficient to approximate the coefficient. To approximate more complicated coefficients, the number of network parameters needs to be increased. However, without explicit regularization, we have been unable to obtain good results in these cases. Multiple discontinuities, highly oscillatory and multiscale coefficients require more advanced networks and/or explicit regularization and will be the topic of future research.



§ SUMMARY AND CONCLUSIONS

We have shown how to augment classical inverse problems with artificial neural networks. We used the finite element method to discretize and solve the forward problem, and automatic differentiation to derive and solve the adjoint problem related to the error functional. The total gradient of the error functional could then be computed exactly by combining the automatic adjoint for the finite element part, and backpropagation/automatic differentiation for the neural network part. By the chain rule, the discretization of the forward and backward problems are factored from the neural network prior in the computation of the gradient of the error functional. With the gradient of the error functional given, any gradient based optimization method can be used.

The neural network acts as a prior for the coefficient to be estimated from noisy data. As neural networks are global, smooth, universal approximators, they do not require explicit regularization of the error functional to recover both smooth solutions and coefficients. The convergence of classical optimization methods, with the neural network augmentation, is independent of both the mesh and geometry in contrast to standard finite elements. As the number of optimization parameters is independent of the mesh, large scale inverse problems can be computed with efficient optimizers, such as BFGS, as long as the number of network parameters is rather small.

The implicit regularization by low capacity network priors allows a smooth reconstruction of discontinuous coefficients as long as there are not too many discontinuities. Since neural networks are globally defined, they allow the reconstruction of the coefficient in the whole domain even when the measurement data is available only close to the boundaries.

For multiple discontinuities and multiscale coefficients, more research is needed on explicit regularization and the choice of neural network design.



§ ACKNOWLEDGMENTS

We would like to thank the anonymous reviewers for their thorough reviews and constructive suggestions on how to improve the quality of the paper.

The authors were partially supported by a grant from the Göran Gustafsson Foundation for Research in Natural Sciences and Medicine.



§ OPTIMAL REGULARIZATION OF THE 1D POISSON EQUATION
 
We return to the problem in section <ref> to discuss optimal regularization for a fair comparison with FEM. We are to minimize the generalized Tikhonov functional (<ref>) which is restated here for convenience,

    J̃ = 1/2∫ |u - û|^2 dx + α^δ/2∫ |q - q_*|^2 dx.

In this case we have a complete description of the noise which is uniformly distributed with a value of δ r added to each measurement, where r ∈𝒩(0, 1) and δ = 0.05. Moreover, since the problem is artificially constructed we know the exact value of q_* in (<ref>). This allows us to compute the optimal value of α^δ using the Morozov discrepancy principle for the various coefficients to be estimated. The Morozov discrepancy principle amounts to finding the unique zero of the monotone function

    f(α^δ) = ∫ |u - û|^2 dx - ||δ||

where ||δ|| denotes the noise level. The computation of u requires that the inverse coefficient problem is solved by minimizing (<ref>). Each evaluation of (<ref>) thus require the solution of the forward problem and the backward problem until convergence in the coefficient q for a given α. If (<ref>) is computed numerically by an iterative method, this procedure is repeated in each iteration making the computation of an optimal regularization extremely expensive. For the 1D Poisson problem, however, we can perform the computations to obtain the optimal regularization parameters for δ = 0.05 as shown in Table <ref>.


The results in section <ref> are repeated here where we have used optimal regularization. It can clearly be seen that FEM outperforms neural networks this time. However, since optimal regularization is rare we propose that neural networks are used to compute the estimated coefficient q_* in (<ref>) which can be used in a generalized Tikohonov regularization.






1.2


abbrv



