Y Favennec
email: yann.favennec@univ-nantes.fr

IRDL P Le Masson

Function estimation in inverse heat transfer problems

Keywords: optimization, convexity, zero-order method, deterministic method, stochastic method, gradienttype method, conjugate gradient, BFGS, Gauss-Newton, Levenberg-Marquardt, gradients, direct di erentiation, adjoint-state

This lecture presents some commonly-used numerical algorithms devoted to optimization, that is maximizing or, more often minimizing a given function of several variables. The goal is function estimation. At rst, some general mathematical tools are presented. Some gradient-free optimization algorithms are presented and then some gradient-type methods are pointed out with pros and cons for each method. The gradient of the function to be minimized is presented according to three distinct methods: nite di erence, forward di erentiation and the use of the additional adjoint-state problem. The last part presents some practical studies where some tricks are given, along with some numerical results.

Introduction

This lecture is devoted to the solution of inverse problems in heat transfer, speci cally when function are to be recovered. Usually, such problems are non-linear and may fall into the category of large-scale inverse problems, so that speci c optimization tools are to be developed. The lecture rst presents some basic examples of IHCP (Inverse Heat Conduction Problems) and points out the distinction between estimation of parameters on the one hand, and functions in the other hand. Indeed, as a simple example, we have the distinction between estimating i) λ as a parameter, ii) λ(x) as a function of the space x (x = (x 1 , x 2 ) t for instance) and, iii) λ(T ) as a function of the state T . The lecture then presents the most usual optimization tools for the solution of di erent kinds of inverse problems. It rst gives notions on the functional to be minimized, and convexity. It gives de nitions of constraints (equality and inequality) added to the functional to be minimized, the added constraints being related to either the state or the parameter/functional. Then, before tackling the detailed iterative optimization algorithms, the most usual stopping criteria are presented. Zero-order, rst-order and quasi-second order optimization methods are brie y presented with pros and cons for each of them. Concerning zero order methods, both deterministic and stochatic methods are very brie y presented with some speci c examples (Simplex, PSO, and GA). Within the frame of rst-order methods, one presents the steepest-descent method with and without linesearch, then the conjugate gradient method for quadratic and arbitrary functions. Some quasi-Newton algorithms are then presented: the BFGS, the Gauss-Newton and the Levenberg-Marquardt methods. A comparison is given in terms of gradient needed for all previously presented methods along with the convergence rate, if possible. The next part presents the computation of the functional gradient: through the nite di erence method, through the direct di erentiation of the PDEs (partial di erential equations), and through the use of the adjoint-state problem. Several ways to access the adjoint-state problems are given. A comparison of gradient computations is given through examples to emphasize the di erences. Note that this lecture has been prepared with some well-known books such as [START_REF] Allaire | Analyse numérique et optimisation[END_REF][START_REF] Minoux | Mathematical Programming, Theory and Applications[END_REF][START_REF] Nocedal | Numerical Optimization[END_REF][START_REF] Gill | Practical Optimization[END_REF]. These books being considered as "standard" popular books, some parts of this lecture are taken from these references. Note also that this lecture is being continuously improved, starting from its very rst version in 2005 [START_REF] Favennec | Optimization, sensitivity and adjoint states[END_REF].

2 Estimation in heat transfer -Optimization

Parameter and function estimation

The modeling of a physical system is based on several requirements. In addition to the physical modeling equations that include some physical parameters (e.g. conductivity coe cients), the initial state and the sources are also to be known if the physical problem is to be solved. If all these data is known, then the so-called direct problem -or forward problem -can be solved. Now, if some of the previously expressed quantities are missing, the physical problem cannot be solved any longer, but some inversion procedure may evaluate the missing quantity, tting the model output with some real ones (i.e. obtained through experiments). The evaluation of such missing quantities needs an inverse problem -or a backward problem -to be solved.

Depending on the nature of the missing quantity, the estimation is performed on parameters or on functions.

These last years, a debate took place within the heat transfer community about the di erence and the meaning of, on the one hand, parameter identi cation and, on the other hand, function estimation. According

Parameter and function estimation

Function estimation in IHTP to the author, both are very di erent, though some similarities exist between both of them.

Let us work on following examples of physical properties estimation to back up our methodology.

i) If a single material conductivity λ is to be identi ed, then the problem clearly belongs to the category of parameter estimation. In such a case, the number of unknowns (the parameters) is very low: only one for a uniform isotropic medium, and only six at maximum for a uniform orthotropic medium. Due to the low dimensionality of the inverse problem, any optimizer can be used (either gradient-free or gradient-type). Moreover, such problems are likely to be well-posed, and the use of regularization tools may not be necessary. These parameter estimation problems are not di cult both from mathematical and computational points of view. The same comments can be drawn if di erent non-varying thermal conductivities are to be estimated in di erent locations (for example dealing with the case of a multi-layer medium).

ii) If several physical properties are to be estimated, for example a thermal conductivity λ [W m -1 K -1 ]), a heat capacity C p [J K -1 ]), and a convective heat transfer coe cient h [W m -2 K -1 ]), then we consider a collection of elements that can be put together into a vector, such that classical optimizers can solve this parameter estimation problem. However, taking a norm of such a vector would not make any sense in a physical point of view. This is one of the reasons why some priors are used (in this speci c case λ 0 , C 0 p and h 0 ), and the estimation is performed on adimensionalized parameters (in this speci c case λ = λ /λ 0 , Cp = Cp /C 0 p and h = h /h 0 ). Doing so, norms (on the collection of adimensionalized parameters) are understandable by both mathematicians and physicists. Note that another reason why it is preferable to adimensionalize parameters is that it usually slightly attenuates the ill-posed character of the inverse problem, and thus the process of adimensionalization can be seen, somehow, as the very rst regularization tool.

iii) If a physical property now depends continuously on the state, (e.g. temperature-dependent conductivity λ(T )), then one may think that the problem of conductivity estimation falls into the category of function estimation. However, a parameterization of this function is anyway necessary to use numerical algorithms, and the type of parameterization can make the di erence between parameter and function estimation. If -for example -a polynomial expansion is used, say λ(T ) ≈ N i α i T i , then the collection of the N coe cients α i is to be estimated, and, therefore, such a problem eventually falls into the category of a parameter estimation problem (the parameters are the polynomial coe cients). Moreover, because the number of unknowns is likely to be low (say less than ten), the choice of the optimizer does not matter much. (Note however that this choice of polynomial expansion is unlikely to be a good candidate for the parameterization; the one presented in the following item iv) is likely to be much better.) iv) If a space-dependent physical property is to be estimated, for example a thermal conductivity λ (x), then the estimation is performed on a function. As in the previous case, a parameterization of this function is anyway necessary to use any numerical algorithm. Building a basis {ξ i } N i=1 and using it to project the function, i.e. with λ (x) = λ i ξ i (T ), the estimation in the end is performed on discrete parameters {λ i } N i=1 , all of these having the same unit, say [W m -1 K -1 ]. At this stage, one may think we face again a parameter estimation problem. However, most often, the function has to satisfy some regularity properties. For example, the conductivity is nite and varies continuously in space, so

λ(x) ∈ H 1 (D) = {λ ∈ L 2 (D), λ ∈ L 2 (D)}.
Because such a regularity property is be satis ed, this problem falls into the category of a function estimation problem, and speci c regularization tools are to be used to enforce the function to satisfy these constraints of regularity. Added to that, the dimension of the discrete unknown, N , is very likely to be big. (As an example, a property de ned in a cube discretized with only 100 voxels per side gives 10 6 unknowns.) Therefore, some speci c optimization algorithms have been designed to cope with such high dimensions.

The function to be minimized

Function estimation in IHTP

It could be seen from previous examples how function estimation is di erent from parameter estimation. Main di erences between both of them come from, on the one hand, the regularity of the functions to be estimated, and, on the other hand, the high dimensionality of the optimization problem due to the parameterization. The regularity issue demands speci c regularization tools and a special care on the parameterization, and the high dimensionality demands powerful optimization algorithms.

The function to be minimized

In an inversion process, one usually minimizes some errors between some experimental data (say u d ) and related model data (say u). The cost function (also called somewhere discrepancy function or objective function) is often expressed as the square of a norm of the di erence between u and u d . The most often, one uses the L 2 (•) norm if some "quasi-"continuous u and especially u d are available (i.e. u - 2 dx but, when data u d is given only on speci c locations (in space and/or time) , then the squared euclidean norm is to be used:

u d 2 L 2 (S) = S (u -u d )
u -u d 2 2 := i (u(x i ) -u d (x i )) 2 = S δ j i (u -u d ) 2 dx where δ j i = δ(x i -x j ).
Often, some function of the state and of the measure are used, for instance state derivation, integration, weighted summation, etc. Moreover, some selection process is, most of the time considered. So, in order to write down a general form for the cost function to be minimized, we use :

J (u) = u -u d 2 X (1)
without specifying any choice for the norm on X at this early stage. Though the cost function is explicitly given in terms of the state u, the cost function is actually to be minimized with respect to what it is searched, i.e. the parameters ψ. Hence we write the equality (by de nition):

j(ψ) := J (u, ψ) (2)
where the function j is often called the reduced cost function, as opposed to J which is the calculated cost function. One actually computes the cost function in terms of the state (by eq. ( 1) for instance), but the cost function is to be minimized with respect to another quantity, say ψ.

Elements of minimization

The function denoted j in eq. ( 2) is de ned on K with values in R. K is a set of admissible elements of the problem. In some cases, K de nes some constraints on the parameters or functions. The minimization problem is written as: inf φ∈K⊂V j(φ).

According to [START_REF] Allaire | Analyse numérique et optimisation[END_REF], if the notation "inf" is used for a minimization problem, it means that one does not know, a priori, if the minimum is obtained, i.e. if there exists φ ∈ K such that Let us now recall basic de nitions needed for mathematical optimization [START_REF] Allaire | Analyse numérique et optimisation[END_REF]:

j (φ) = inf ψ∈K⊂V j(ψ).
De nition 1. ψ is a local minimum of j on K if and only if ψ ∈ K and ∃δ > 0, ∀φ ∈ K, φ -ψ < δ → j(φ) ≥ j(ψ).

Optimality conditions

Function estimation in IHTP

De nition 2. ψ is a global minimum of j on K if and only if ψ ∈ K and j(φ) ≥ j(ψ) ∀φ ∈ K.

De nition 3.

A minimizing series of j in K is a series (ψ n ) n∈N such that

ψ n ∈ K ∀n and lim n→+∞ j (ψ n ) = min φ∈K j(φ). De nition 4. a set K ∈ V is convex if, for all ψ, φ ∈ K and ∀θ ∈ [0, 1], the element (θψ + (1 -θ)φ) is in K (see gure 1).

De nition 5.

A function j is said to be convex when de ned on a non-null convex set K ∈ V with values in R if and only if

j (θψ + (1 -θ) φ) ≤ θj (ψ) + (1 -θ) j (φ) ∀ψ, φ ∈ K, ∀θ ∈ [0, 1] .
Moreover, j is said to be strictly convex if the inequality is strict when ψ = φ and θ ∈ ]0, 1[ (see g. 2).

Ending, if j if a convex function on K, the local minimum of j on K is the global minimum on K.

Optimality conditions

For convex functions, there is no di erence between local minima and global minimum. In the following we are more interested in minimizing a function without specifying whether the minimum is local or global. It

Stopping criteria

Function estimation in IHTP will be seen in next sections that some gradient-free optimization algorithms may nd the global minimum even if the cost function contains local minima. Let us derive here the minimization necessary and su cient conditions. These conditions use the rst-order derivatives (order-1 condition), and second-order derivatives (order-2 condition) on the cost function j. Using gradient-type algorithms, the rst-order condition is to be reached, while the second-order condition requires assuming a local convexity hypothesis, and then make a distinction between minima, maxima and optima.

Let us assume that j(ψ) is continuous and has continuous rst partial derivatives (∂j/∂ψ i )(ψ) and second partial derivatives (∂ 2 j/∂ψ i ∂ψ j )(ψ). Then a necessary condition for ψ to be a minimum of j (at least locally) is that: i) ψ is a stationary point, i.e. ∇j( ψ) = 0, ii) the Hessian ∇ 2 j( ψ) = ∂ 2 j/∂ψ i ∂ψ j ( ψ) is a positive semi-de nite matrix, i.e. ∀y ∈ R n , ∇ 2 j( ψ)y, y ≥ 0 where (., .) is a scalar product in R n (we have dim(ψ) = n).

A point ψ which satis es condition item i) is called a stationary point. It is important to point out that stationarity is not a su cient condition for local optimality. For instance the point of in exion for cubic functions would satisfy the condition i), while there is no optimum. Hence the Hessian is not positive de nite but merely positive semi-de nite.

The su cient condition for ψ to be a minimum of j (at least locally) is that i) ψ is a stationary point, i.e. ∇j( ψ) = 0, ii) the Hessian ∇ 2 j( ψ) = ∂ 2 j/∂ψ i ∂ψ j ( ψ) is a positive de nite matrix, i.e. ∀y ∈ R n , y = 0, ∇ 2 j(ū)y, y > 0.

We remark that the condition item ii) amounts to assuming that j is strictly convex in the neighbourhood of ψ.

Stopping criteria

Since the convergence of the iterative algorithms is, in general, not nite, a stopping criterion must be applied. Here below are given some commonly used criteria. We denote ψ p the vector parameter ψ at the optimization iteration p. ∇j(ψ p ) ≤ ε;

(3)

j(ψ p ) -j(ψ p-1 ) ≤ ε; (4)
ψ p -ψ p-1 ≤ ε; (5)
j(ψ p ) ≤ v (6)
For each of the above criteria, it may be judicious to demand that the test is satis ed over several successive iterations. The three rst above-presented criteria are convergence criteria applied on the cost function gradient, on the cost function evolution, or on the parameter evolutions. These criteria are very commonly-used when dealing with optimization and optimal control problems.

The last criterion is, in one sense, more speci c to inverse problems: when the cost function reaches a critical value that depends on the variance of measurement errors, then the optimization algorithm should stop [START_REF] Av Goncharskii | A generalized discrepancy principle[END_REF][START_REF] Vladimir | Methods for solving incorrectly posed problems[END_REF][START_REF] Alifanov | Inverse Heat Conduction, Ill-posed problems[END_REF]. It can be shown that the consequence of lowering the cost function below v a ects the result in dramatically highliting its inherent noise. This criterion is the "maximum discrepancy principle".

Often, the maximum discrepancy principle eq. ( 6) is used together with the other criteria and also with a maximum number of iterations.

Classi cation of optimization methods

Function estimation in IHTP

Classi cation of optimization methods

The solution of the optimization problem may be performed in numbers of ways. Among numerous methods found in the litterature, the classi cation of methods given below (see g. 3) is based on our experience. First, one can distinguish gradient-free methods from methods relying on gradients. Among gradient-free methods, there are those deterministic and those stochastic (the latter introducing random in the search of the optimum). Among gradient-based methods, one can distinguish between rst and second-order methods, and those in between. 3 Zero-order n-dimensional optimization algorithms

n-D Optimization Methods

Gradient

Zero-order methods, also called "derivative-free optimization" (DFO) or "gradient-free methods" are based on a global vision of the cost function value j on the search space. The main interest of using such methods is when the cost function gradient is not available, or when the cost gradient is not easy to compute, or when the cost function presents local minima. There is an increasing number of computation tools to solve optimization problems with no gradient [START_REF] Onwubolu | New optimization techniques in engineering[END_REF]. In the sequel, we restrict our-self in very brie y presenting, among the enormous number of existing methods, one deterministic algoritm which is the so-called simplex method, and one probabilistic method which is the particle swarm optimization method.

Simplex

We present here the Nelder-Mead simplex method (1965). This method is popular and simple to code. Moreover, there exists a large number of freeware that can be used to minimize a function using such an algorithm. Let a simplex S 0 be a set of n + 1 "points" linearly independent (n = dim ψ) with S 0 = ψ I , I = 1, . . . , n + 1 . One iteration of the simplex optimization algorithm consists in generating a new simplex closer to the minimum eliminating the point with the higher cost function value. The basic operations of n = 2 are given in g. 4: let ψ the isobarycenter of ψ I , I = 1, . . . , n, (without ψ h = arg I=1,...,n max j(ψ I )), let the ordering so that j(ψ 1 ) ≤ j(ψ 2 ) ≤ . . . ≤ j(ψ n+1 )

and let ψ = arg I=1,...,n min j(ψ I ). At each iteration, the simplex improvement is performed in three steps:

1. [Re ection] One builds ψ R symmetry of ψ h with respect to the segment [ ψ, ψ ]. According to the value of the cost j(u R ) with respect to j(ψ ), the parametric space is then extended (step 2), or contracted (step 3);

PSO

Function estimation in IHTP 2.

[Extension] if j(ψ R ) < j(ψ ), one searches a new point in the same direction. The point ψ E is such that

ψ E = γψ R + (1 -γ) ψ with γ > 1. If j(ψ E ) < j(ψ R ), ψ h is replaced by ψ R , otherwise ψ h is replaced by ψ E ; 3. [Contraction] If j(ψ R ) > j(ψ ), the point ψ C such that ψ C = γψ h + (1 -γ) ψ, γ ∈]0, 1[ is created. If j(ψ C ) < j(ψ R
), ψ h is replaced by ψ C otherwise the simplex is contracted (inside contraction) in all directions replacing ∀I = L ψ I by (ψ I + ψ )/2.

• ψ h • ψ l • × ψ ψ R • ψ h • ψ l • × ψ ψ R ψ E • ψ h • ψ l • × ψ ψ R ψ C • ψ h • ψ l • × ψ

PSO

The particle swarm optimization is a stochastic algorithm described by Kennedy and Eberhart in 1995. One considers an initial set of individuals (particles) located randomly. Each particle moves within the space K interacting with other particles on their best locations. From this information, the particle shall change its position ψ i and its velocity δψ i . The general formulation for this behavior is given by:

δψ i = χδψ i + λ 1 rand 1 φ g -ψ i + λ 2 rand 2 φ i -ψ i ψ i = ψ i + δψ i (7)
where ψ i is the position of the particle i, δψ i is its velocity, φ g is the best position obtained in its neighborhood, and φ i is its best position (see g. 5). χ, λ 1 and λ 2 are some coe cients weighting the three directions of the particule [START_REF] Onwubolu | New optimization techniques in engineering[END_REF]:

• how much the particle trusts itself now;

• how much it trusts its experience;

• how much it trusts its neighbours.

Next, rand 1 and rand 2 are random variables following a uniform distribution in [0, 1]. There are several con guration parameters for the method, see [START_REF] Clerc | L'optimisation par essaims particulaires[END_REF]:

• swarm size, usually between 20 and 30;

• initialization of both the position of the particles and their velocity ∼ U[0, 1];

• neighborhood topology such that a particule communicates with only some other particles;

• inertial factor χ which de nes the exploration capacity of the particules;

• con dence coe cients λ 1 and λ 2 which are constriction coe cients;

• stopping criterion which is usually the maximum of iterations, or the critical value of the cost function j(ψ).

Usually, a circular neighborhood topology is used, along with χ = 0.72 and λ 1 = λ 2 = 1.46. A large number of free software are available, see for instance [START_REF] Clerc | L'optimisation par essaims particulaires[END_REF].

• ψ i φ i φ g δψ i

One-dimensional unconstrained opimization -line search algorithm

In order to nd the optimum of a function j of n variables, we shall describe in section 5 a number of iterative methods which require, at each step, the solution of an optimization problem in one single variable, of the type:

Find ᾱ = arg min α>0 g(α) = j (ψ p + αd p ) , (8)
where ψ p = (ψ p 1 . . . ψ p n ) t is the obtained point at iteration p and where d p = (d p 1 . . . d p n ) t is the direction of descent (see section 5). As a matter of fact we have the problem of nding the optimum of the function j, starting from the guess ψ 0 in the direction of descent d 0 . Since this problem must be solved a great number of times, it is important to design e cient algorithms that deal with it.In any case, one has to keep in mind that the main objective is not to solve eq. ( 8) but to nd the minimum of j(ψ). Thus one has to design e cient tools for the one-dimensional algorithm that nds the minimum of g(α), or approach it, in a not so expensive way. Note that we always assume that g (0) = (∇j(ψ p ), d p ) < 0, which means that d p is indeed a descent direction.

The dichotomy method

This method halves, at each step, the length of the interval which contains the minimum, by computing the function g in two new points. By carrying out n computations of the function g, the length of the initial interval [a 0 , b 0 ] is reduced in a proportion of 2 (n-3)/2 . The general procedure is the following: starting from the interval [a 0 , b 0 ], and taking the midpoint c 0 = a 0 + b 0 /2, and the two points d 0 = a 0 + c 0 /2, and e 0 = c 0 + b 0 /2, one obtains ve equidistant points of length δ 0 = (b 0a 0 )/4; computing the cost function values at these points, two of the four sub-intervals may be eliminated, while the two adjacent sub-intervals remain; the same procedure is repeated within the selected interval [a 1 , b 1 ], and so on. Since the step length is divided by 2 at each iteration, the dichotomy method converges linearly to the minimum [START_REF] Minoux | Mathematical Programming, Theory and Applications[END_REF].

The Newton-Raphson method

Let us assume that the function g(α) is twice continuously di erentiable. The search for a minimum of g(α) is carried out by looking for a stationary point, i.e. ᾱ satisfying the possibly nonlinear relationship g (ᾱ) = 0. If α q is the point obtained at stage q, then the function g (α) is approximated by its tangent, and the next point α q+1 is chosen to be at the intersection of this tangent with the zero-ordinate axis. The relationship to pass from one step to the next comes from g α q+1 = g (α q ) + g (α q ) × α q+1α q = 0 which gives:

α q+1 = α q - g (α q ) g (α q ) . ( 9
)

The secant method

Function estimation in IHTP It is of interest that this method has the property of nite convergence when applied to quadratic functions. This is an interesting feature because any function which is su ciently regular (at least twice continuously di erentiable) behaves as a quadratic function near the optimum [START_REF] Minoux | Mathematical Programming, Theory and Applications[END_REF]. On the other hand, the main drawback of this method is that it requires the computation of the rst and of the second derivative of g at each stage. That is the reason why the secant method (next section) is also widely used, especially when there is no way for computing the second order derivative, or when the exact second derivative is complicated to compute or too time consuming.

g (α 0 ) α 0 α 1 α 2

The secant method

The second-order derivative g (α) is approximated by nite di erences so that the Newton-Raphson's equation initially given by eq. ( 9) becomes eq. ( 10):

α q+1 α q -g (α q ) α q -α q-1 g (α q ) -g (α q-1 ) . ( 10
)
This method is the so-called secant method. Applied to the search of g (α) = 0, this method consists in searching the intersection between the zero-ordinate axis and the straight line passing by the points α q-1 , g (α q-1 ) and [α q , g (α q )].

The quadratic interpolation

By comparison of those of sections 4.2 and 4.3, this method has the advantage of not requiring the computation of rst or second order derivatives of the function . Let three points

α 1 ≤ α 2 ≤ α 3 such that g(α 1 ) ≥ g(α 2 ) ≤ g(α 3
) and let us approximate the function g on the related interval by a quadratic function g with the same values as those of g at the points α 1 , α 2 and α 3 . The minimum of g is obtained at the new point α 4 satisfying:

α 4 = 1 2 r 23 g(α 1 ) + r 31 g(α 2 ) + r 12 g(α 3 ) s 23 g(α 1 ) + s 31 g(α 2 ) + s 12 g(α 3 ) , (11)
where r ij = α 2 iα 2 j and s ij = α iα j . This procedure may be repeated again with the three new selected points. Under some regularity hypothesis, the convergence rate of this method is super-linear [START_REF] Minoux | Mathematical Programming, Theory and Applications[END_REF].

Another approach consists in di erentiating the cost function towards the direction of descent with a Taylor expansion, andin neglecting second order derivatives:

g (α) = d dα u ψ k + αd k -u d 2 X = d dα u ψ k -αu ψ k ; d k -u d 2 X = 0 (12)
, u -u d X + α u , u X = 0 (13) α = - (u , u -u d ) X (u , u ) X (14)
This latter method -which is widely used in the heat transfer community -can give easily an accurate step size α when the cost function j is close to quadratic, i.e. when the state u varies almost linearly with ψ.

Other methods -Inexact line-search

A great number of other one-dimensional optimization methods may be found in the literature. These methods may be more or less complicated and some of them may be much more optimal than the abovepresented methods. In practice the Fibonacci method,the golden section search method and the cubic interpolation method are also very widely used in practice (the reader may refer to [START_REF] Gill | Practical Optimization[END_REF][START_REF] Minoux | Mathematical Programming, Theory and Applications[END_REF] for more details). All these methods can be quite CPU-time consuming, and in fact, the convergence of some of the methods presented afterwards in Section 5 (typically the BFGS method) can be reached without getting a point very close to satisfying g (α) = 0. Well-accepted conditions used to build inexact line-search algorithms are based on the two rules: a) α must not be too large in order, for instance, to avoid oscillations, b) α must not be chosen too small in order to prevent from premature convergence.

Among the large number of inexact line-search algorithms, one is based on the Goldstein rules (see Figure 7) which rst ensures condition a) by satisfying ( 15) choosing m 1 ∈ [0, 1], and second ensures condition b) satisfying [START_REF] Favennec | Space-dependent sobolev gradients as a regularization for inverse radiative transfer problems[END_REF]

choosing m 2 ∈ [m 1 , 1]. g (α) ≤ g (0) + m 1 αg (0) (15)
g (α) ≥ g (0) + m 2 αg (0) (16)
Other rules can be stated in similar ways. For instance, the Armijo's method is a variant of the Golstein method. Related algorithms are very simple and can be found in any book on optimization. Algorithm 1: Typical method based on the Goldstein rules input :

0 g (0) m 2 g (0) m 1 g (0) a b c d
α min = 0, α max = ∞, ψ = ψ k , ∇j, d = d k output : ᾱ 1 Give some initial value to α; Compute g (0) = (∇j, d); 2 Compute g(α) = j(ψ + αd); if g(α) ≤ g(0) + m 1 αg (0) then go to 3) else set α max = α and go to 5 end 3 Compare g (α) and g (0) + m 2 αg (0); if g(α) ≥ g(0) + m 2 αg (0) then END else
go to 4 end 4 Set α min = α; 5 Look for new value in ]α min , α max [ and return to 2

Gradient-type n-dimensional optimization algorithms

Since in all cases, the stationarity of j is a necessary optimality condition, almost all unconstrained optimization methods consist in searching the stationary point ψ where ∇j( ψ) = 0. The usual methods are iterative and proceed this way: one generates a sequence of points ψ 0 , ψ 1 ,. . . ψ p which converges to a local optimum of j. At each stage p, ψ p+1 is de ned by ψ p+1 = ψ p + α p d p where d p is a displacement direction which may be either the opposite of the gradient of j at ψ p (i.e. d p = -∇j(ψ p )), or computed from the gradient, or chosen in any another way, provided that it is a descent direction, i.e. satisfying (∇j(ψ p ), d p ) < 0.

5.1 1st order gradient methods 5.1.1 The gradient with prede ned steps method (1st order method) At each iteration step p, the gradient ∇j(ψ p ) gives the direction of the largest increase of j. The procedure consists in computing the gradient, and in nding the new point according to the prede ned strictly positive step size α p as:

ψ p+1 = ψ p -α p ∇j(ψ p ) ∇j(ψ p ) . ( 17
)
It may be shown that this iterative scheme converges to ψ provided that α p → 0 (p → ∞) and ∞ p=0 α p = +∞. One can choose for instance α p = 1/p. The main drawback of this method is its very low convergence rate.

The steepest descent method (1st order method)

In this frequently used method, α p is chosen at each iteration p so as to minimize the function g(α) = j (ψ p -α∇j(ψ p )) on the set of α ≥ 0. The algorithm is thus the following. One chooses a starting point ψ 0 and set p = 0. At each iteration p, one computes the gradient and set d p = -∇j(ψ p ). One then solves the one-dimensional problem (see section 4) and set ψ p+1 = ψ p + α p d p . This procedure is repeated until a stopping test is satis ed (see section 2.5). The main disadvantage of the steepest descent method is the fact that the convergence can still be very slow. As a matter of fact, since α p minimizes g(α) = j (ψ p + αd p ) then g (α p ) = (d p , ∇j (ψ p + αd p )) = d p , ∇j ψ p+1 . Hence d p , d p+1 = 0. This means that two successive displacements are strictly orthogonal. As a direct consequence, the number of steps to minimize elongated valley-type functions for instance may be very high (see g. 8 and then g. 10d on page 22).

ψ 0 d 0 ψ 1 d 1 ψ 2 d 2
Figure 8: When the steepest descent method is used, the two consecutive directions are orthogonal.

The conjugate gradient method for quadratic functions (1st order method)

In this section we shall rstly assume that the cost function is quadratic. The case of arbitrary functions shall be dealt with in section 5.1.4. Let the quadratic functional be of the form:

j(ψ) = 1 2 (A ψ, ψ) , (18)
and let us recall the de nition for two conjugate vectors. Let A be a given symmetric matrix (operator). Two vectors x 1 and x 2 are A -conjugate if (A x 1 , x 2 ) = 0. The general method to optimize j is the following. Let us start with a given ψ 0 and choose d 0 = -∇j(ψ 0 ). One may remark that for quadratic functions, the one-dimensional minimization procedure may be analytically solved. Recalling that the minimization of g(α) along the direction d 0 should lead to the fact that this current direction (d 0 ) would be orthogonal to the next gradient ∇j(ψ 1 ), one has:

d 0 , ∇j ψ 1 = 0. (19)
Using the relationship ∇j(ψ) = A ψ given by the di erentiation of (18) and the reactualization formulation ψ 1 = ψ 0 + α 0 d 0 , (19) becomes:

d 0 , ∇j ψ 1 = d 0 , A ψ 1 = d 0 , A ψ 0 + α 0 d 0 = d 0 , A ψ 0 + α 0 d 0 , A d 0 . (20)
Equaling (20) to zero gives the step size α 0 :

α 0 = - d 0 , A ψ 0 (d 0 , A d 0 ) . (21)
Next, at stage p, we are at the point ψ p and we compute the gradient ∇j(ψ p ). The direction d p is obtained by combining linearly the gradient ∇j(ψ p ) and the previous direction d p-1 , where the coe cient β p is chosen in such a way that d p is A -conjugate to the previous direction. Hence:

d p , A d p-1 = -∇j (ψ p ) + β p d p-1 , A d p-1 = -∇j (ψ p ) , A d p-1 + β p d p-1 , A d p-1 . ( 22
)
5.1 1st order gradient methods

Function estimation in IHTP

Next, choosing β p such that the previous equation equals zero yields to:

β p = ∇j(ψ p ), A d p-1 (d p-1 , A d p-1 ) . ( 23
)
The algorithm based on the above relationships is given in algorithm 2. Also, it is proved, see [START_REF] Minoux | Mathematical Programming, Theory and Applications[END_REF], that the conjugate gradient method applied to quadratic functions converges in at most n iterations, where n = dim ψ.

Algorithm 2: The conjugate gradient algorithm applied to quadratic functions 1. Let p = 0, ψ 0 be the starting point, compute the gradient and the descent direction,

d 0 = -∇j(ψ 0 ), compute the step size α 0 = - d 0 , A ψ 0 (d 0 , A d 0 ) ; 2.
At step p, we are at the point ψ p . We de ne ψ p+1 = ψ p + α p d p with:

• the step size

α p = - (d p , ∇j(ψ p )) (d p , A d p ) • the direction d p = -∇j(ψ p ) + β p d p-1
• where the coe cient needed for conjugate directions:

β p = ∇j(ψ p ), A d p-1 (d p-1 , A d p-1 ) ; 3.
Stopping rule (see Section 2.5). If satis ed: End, otherwise set p ← p + 1 and return to step [START_REF] Minoux | Mathematical Programming, Theory and Applications[END_REF].

ψ 0 d 0 ψ 1 d 1 ψ 2
Figure 9: When the conjugate gradient method is used, the two consecutive directions are conjugate instead of orthogonal. Applied to a quadratic function, the method converges in at most n iterations (in this gure two iterations are needed since dim ψ = 2).

5.1.4

The conjugate gradient method for arbitrary (non quadratic) functions (1st order)

Before presenting the application of conjugate gradient methods on arbitrary functions, let us give some properties inherent to quadratic functions. Di erentiating eq. ( 18), and taking into account of the re-actualization relationship, one has:

∇j(ψ p ) -∇j(ψ p-1 ) = A ψ p -ψ p-1 = A ψ p-1 + α p-1 d p-1 -ψ p-1 = α p-1 A d p-1 , (24)
which also gives the following relationship:

1 α p-1 ∇j(ψ p ), ∇j(ψ p ) -∇j(ψ p-1 = ∇j(ψ p ), A d p-1 . ( 25
)
On the other hand, substituting (25) into (23) gives

β p = ∇j(ψ p ), A d p-1 (d p-1 , A d p-1 ) = ∇j(ψ p ), ∇j(ψ p ) -∇j(ψ p-1 ) (d p-1 , ∇j(ψ p ) -∇j(ψ p -1)) . ( 26
)
Next, expanding the descent direction d p-1 , (26) becomes:

β p = ∇j(ψ p ), ∇j(ψ p ) -∇j(ψ p-1 ) (-∇j(ψ p-1 ) + β p-1 d p-2 , ∇j(ψ p ) -∇j(ψ p-1 )) ; (27)
β p = ∇j(ψ p ), ∇j(ψ p ) -∇j(ψ p-1 ) (-∇j(ψ p-1 ) -β p-1 ∇j(ψ p-2 ) + Λ, ∇j(ψ p ) -∇j(ψ p-1 )) , ( 28
)
where Λ is the series given from the re-actualizations. All the gradients being orthogonal to each other, (28) becomes:

β p = ∇j(ψ p ), ∇j(ψ p ) -∇j(ψ p-1 ) (∇j(ψ p-1 ), ∇j(ψ p-1 )) , (29)
and also:

β p = (∇j(ψ p ), ∇j(ψ p )) (∇j(ψ p-1 ), ∇j(ψ p-1 )) . ( 30
)
It is pointed out that in the neighborhood of the optimum, non-quadratic functions may always be approximated by quadratic functions. The Fletcher and Reeves' method consists in applying (30) to access β p while the Polak and Ribiere's method consists in applying (29) to access β p . Taking into account of above remarks, the conjugate gradient algorithm applied to arbitrary functions is given in Algorithm 3. It is important to note that the global convergence of the presented methods is only ensured if a periodic restart is carried out . The restart d n = -∇j(u n ) is usually carried out every n iterations, at least.

Algorithm 3:

The conjugate gradient algorithm applied to arbitrary functions 1. Let p = 0, ψ 0 be the starting point, d 0 = -∇j(ψ 0 ); 2. At step p, we are at the point ψ p ; we de ne ψ p+1 = ψ p + α p d p with:

• the step size α p = arg min α∈R + g(α) = j (ψ p + αd p ) with:

• the direction d p = -∇j(ψ p ) + β p d p-1 where • the conjugate condition β p satis es either (29) or (30) depending on the chosen method; 3. Stopping rule (see subsection 2.5). If satis es: End, otherwise, set p ← p + 1 and return to step (2).

The Newton's method (2nd order)

Function estimation in IHTP

The Newton's method (2nd order)

Let us assume that the cost function j(ψ) is now twice continuously di erentiable and that second derivatives exist. The idea is to approach the next cost function gradient by its quadratic approximation through a Taylor development:

∇j(ψ p+1 ) = ∇j(ψ p ) + ∇ 2 j(ψ p ) δψ p + O (δψ p ) 2 , (31)
and equaling the obtained approximated gradient to zero to get the new parameter ψ p+1 = δψ p + ψ p :

ψ p+1 = ψ p -∇ 2 j (ψ p ) -1 ∇j(ψ p ). (32)
Note that while using second-order optimization algorithms, the direction of descent as well as the step size are obtained from (32) in one go. Another interesting point is the fact that the algorithm converges to ψ in a single step when applied to strictly quadratic functions. However, for arbitrary functions, O (δψ p ) 2 may be far from zero in eq. ( 31); yielding to some errors in the displacement δψ p , and thus in the new point ψ p+1 . As a consequence, if the starting point ψ 0 is too far away from the solution ψ, then the Newton method may not converge. On the other hand, since the approximation of j(ψ) by a quadratic function is almost always valid in the neighborhood of ψ, then the algorithm should converge to ψ if the starting point ψ 0 is chosen closely enough to the solution. Moreover, it is very common to control the step size this way. One rst calculates the direction d p = -∇ 2 j(ψ p ) -1 ∇j(ψ p ) and control the step size through an iterative one-dimensional minimization problem of the kind min g(α) = j (ψ p + αd p ) before the actualization ψ p+1 = ψ p + αd p . One limitation of the Newton's method is when the Hessian ∇ 2 j(u p ) is not positive de nite. In these cases, the direction given by d p = -∇ 2 j(ψ p ) -1 ∇j(u p ) may not be a descent direction, and the global convergence of the algorithm may not be guaranteed any more. Moreover, and above all, the Hessian is usually very di cult to compute and highly time consuming. To overcome these di culties, one should prefer using one of the numerous quasi-Newton methods detailed afterwards.

Quasi-Newton methods

Quasi-Newton methods consist in generalizing the Newton's recurrence formulation (32). Since the limitation of the Newton's method is the restriction of the Hessian ∇ 2 j(u p ) to be positive de nite, the natural extension consists in replacing the inverse of the Hessian by an approximation to a positive de nite matrix denoted H p . Obviously, this matrix is modi ed at each step p. There is much exibility in the choice for computing the matrix H p . In general, the condition given by (33) is imposed:

H ∇j(ψ p ) -∇j(ψ p-1 ) = ψ p -ψ p-1 . ( 33
)
Various corrections of the type

H p+1 = H p + Λ p (34)
may be found in the literature [START_REF] Minoux | Mathematical Programming, Theory and Applications[END_REF]. Depending on whether ∆ p is of rank 1 or 2, we shall speak of a correction of rank 1 or 2.

Rank 1 correction

The point is to choose a symmetric matrix H 0 and to perform the corrections so that they preserve the symmetry of the matrices H p . The rank 1 correction matrix consists in choosing ∆ p = α p v p v p where v p is a vector and α p is a scalar such that, from a symetric matrix H 0 , the correction preserves the symetry of matrices H p . Denoting

δ p = ψ p+1 -ψ p (35)
γ p = ∇j(ψ p+1 ) -∇j(ψ p ) (36)
one chooses α p and v p such that H p+1 γ p = δ p , thus:

H p + α p (v p v pt ) γ p = δ p , (37)
and

γ p H p γ p + α p γ p v p v p γ p = γ p δ p , (38)
thus α p v p γ p 2 = γ p (δ p -H p γ p ) . (39)
Using the identity

α p v p v p = α p v p v p γ p α p v p v p γ p α p v p γ p 2 , ( 40
)
and using (37) and (38) to get

α p v p v p γ p = δ p -H p γ p , (41)
α p v p γ p 2 = γ p (δ p -H p γ p ) , (42)
one obtains the correction (of rank 1) of the inverse Hessian:

H p+1 -H p = α p v p v p = (δ p -H p γ p ) (δ p -H p γ p ) γ p (δ p -H p γ p ) . (43)

The rank 2 Davidon-Fletcher-Powell (DFP) algorithm

The Davidon-Fletcher-Powell algorithm (in short DFP) consists in modifying the inverse Hessian with the correction formulation of rank 2:

H p+1 = H p + δ p (δ p ) (δ p ) γ p - H p γ p (γ p ) H p (γ p ) Hγ p (44)
where we have de ned above δ p = ψ p+1ψ p and γ p = ∇j(ψ p+1 ) -∇j(ψ p ), and where the new point ψ p+1 is obtained from ψ p through the displacement

d p = -H p ∇j(ψ p ). (45)
The global DFP method is presented in algorithm 4.

The rank 2 Broyden -Fletcher -Goldfarb -Shanno (BFGS) algorithm

The Broyden -Fletcher -Goldfarb -Shanno algorithm (in short BFGS) developped in 1969-1970 uses a rank 2 correction matrix for the inverse Hessian that is derived from eq. ( 44). It can be shown [START_REF] Minoux | Mathematical Programming, Theory and Applications[END_REF] that the vectors δ p and γ p can permute in eq. ( 44) and in the relationship H p+1 γ p = δ p . The correction eq. ( 44) can thus also approximate the Hessian itself, and the correction for the inverse Hessian H p+1 can thus be given from H p through the correction formulation:

H p+1 = H p + 1 + γ pt H p γ p δ p t γ p δ p (δ p ) t (δ p ) t γ p - δ p γ pt H p + H p γ p δ pt δ p t γ p . ( 47
)
When applied to a non purely quadratic function, one has, as for the conjugate gradient method and the DFP method, to carry out a periodic restart in order to ensure the convergence [START_REF] Gill | Practical Optimization[END_REF][START_REF] Culioli | Introduction à l'optimisation[END_REF]. It is known that the BFGS algorithm is superior than the DFP algorithm is the sense that it is much less sensitive on the line-search inaccuracy, allowing the use of economical inexact line-search algorithms [START_REF] Minoux | Mathematical Programming, Theory and Applications[END_REF].

Algorithm 4: The Davidon -Fletcher -Powell (DFP) algorithm 1. Let p = 0, ψ 0 be the starting point. Choose any positive de nite matrix H 0 (often the identity matrix);

2. at step p, compute the displacement direction d p = -H p ∇j(ψ p ), and nd ψ p+1 at the minimum of j(ψ p + αd p ) with α ≥ 0;

3. set δ p = ψ p+1ψ p and compute γ p = ∇j(ψ p+1 ) -∇j(ψ p ) to actualize:

H p+1 = H p + δ p (δ p ) t (δ p ) t γ p - H p γ p (γ p ) t H p (γ p ) t Hγ p ; (46)
4. Stopping rule (see section 3.4). If satis es: End, otherwise, set p ← p + 1 and return to step item 2.

Algorithm 5: The BFGS algorithm 1. Let p = 0, ψ 0 be the starting point. Choose any positive de nite matrix H 0 (often the identity matrix);

2. at step p, compute the displacement direction d p = -H p ∇j(ψ p ), and nd ψ p+1 at the minimum of j(ψ p + αd p ) with α ≥ 0;

3. set δ p = ψ p+1ψ p and compute γ p = ∇j(ψ p+1 ) -∇j(ψ p ) to actualize:

H p+1 = H p + 1 + γ pt H p γ p δ p t γ p δ p (δ p ) t (δ p ) t γ p - δ p γ pt H p + H p γ p δ pt δ p t γ p (48)
4. Stopping rule (see section 3.4). If satis es: End, otherwise, set p ← p + 1 and return to step item 2.

Gauss-Newton

When the cost function is explicitly a square norm of the error between the prediction and the state, that is of the form j(ψ)

:= J (u) = u -u d 2 X , (49)
then the Gauss-Newton method or some derivatives or it (e.g. Levenberg-Marquardt) may be interesting to deal with, especially if the number of parameters is small. Before going deeper into the cost function gradient computation (see section 6), de ning u (ψ; δψ) as the derivative of the state at the point ψ in the direction δψ as:

u (ψ; δψ) := lim →0 u(ψ + δψ) -u(ψ) , (50)
then the directional derivative of the cost function writes out as:

j (ψ; δψ) = u -u d , u (ψ; δψ) X , (51)
where j (ψ; δψ) = (∇j(ψ), δψ). In the analogue way, the second derivative of j(ψ) at the point ψ in the directions δψ and δφ is given by: j (ψ; δψ, δφ) = uu d , u (ψ; δψ, δφ) X + u (ψ; δψ), u (ψ; δφ) X .

(52)

Elements of comparison between some presented methods

Function estimation in IHTP

Neglecting the second-order term (this is actually the Gauss-Newton approach), we have:

j (ψ; δψ, δφ) ≈ u (ψ; δψ), u (ψ; δφ) X . (53)
In order to build up the gradient vector of cost function and the approximated Hessian matrix, one has to choose the directions for the whole canonical basis of ψ. Doing so, one can use the so-called sensititivity matrix S which gathers the derivatives of u in all directions δψ i , i = 1, . . . , dim ψ, and the product (u (ψ; δψ i ), u (ψ; δψ j )) involved in (53) is the product of the so-called sensitivity matrix with its transposed. The Newton relationship is thus approximated as:

S t Sδψ k = -∇j(ψ k ).
(54)

The matrix system S t S is obviously symmetric and positive de nite with a dominant diagonal yielding thus to interesting features (Cholesky factorization, etc.). Though the Gauss-Newton system eq. ( 54) presents inherent interesting features (it almost gives in one step the descent direction and the step size), the matrix S t S is likely to be ill-conditionned. One way to decrease signi cantly the ill-condition feature is to "damp" the system, using:

S t S + I δψ k = -∇j(ψ k ), (55)
or better:

S t S + diag(S t S) δψ k = -∇j(ψ k ). (56)
Note that → 0 yields the Gauss-Newton algorithm while bigger gives an approximation of the steepest descent gradient algorithm. In practice, the parameter may be adjusted at each iteration.

Elements of comparison between some presented methods

Some of the presented methods are below tested on the well-known Rosenbrock function:

f (x, y) = (x -α) 2 + β(x 2 -y) 2 . ( 57
)
For the considered case, the chosen parameters are α = 1 and β = 100, so that the optimum is at (1, 1). Figure 10a on page 22 presents the function. This function presents a long elongated valley where the function gradient is very low. Next, the PSO algorithm is the one from [START_REF] Clerc | L'optimisation par essaims particulaires[END_REF].

The deterministic simplex method from the GSL library starting from the point x 0 = -1, y 0 = 1 needs 64 evaluations of the cost function. The stopping criterion is based on the simplex characteristic size equal to 10 -2 . The PSO algorithm taken from [START_REF] Clerc | L'optimisation par essaims particulaires[END_REF] with 20 particles with 3 informed particles, φ = 4.14,

χ = 2 φ-2+ √ φ 2 -4φ
, λ 1 = λ 2 = 0.5χφ. The stopping criterion is based on the cost function equal to 10 -5 . With these parameters, around 6,000 evaluations of the cost function is needed for the minimization. For the Steepest descent, the conjugate gradient and the BFGS algorithms, the stopping criterion is based either on the gradient norm equal to 10 -3 , or on a maximum number of iterations equal to 10,000. For the steepest descent method, the maximum of iteration criterion is achieved. For the conjugate gradient, and the BFGS method, 49 and 11 iterations are needed, respectively.

To sum up about this numerical optimization test case in which we were searching the minimum of the Rosenbrock function, we can give the following comments and conclusions:

• The PSO method, which is a stochastic zero-order method -as genetic algorithms are also -does converge to the minimum, but at a huge expense. In fact, usually, such stochastic methods are even able to nd the minimum of non-convex functions, which is their most important advantage, but anyway at the price of being very expensive.

• When the function is likely to be convex (which is not the case of the Rosenbrock function), one should prefer less expensive deterministic optimization algorithms. Among those, the simplex method is also a gradient-free (zero-order optimizer) so it nds the minimum of the convex function at a more moderate expense, because it is deterministic. But we are here -in this simple example -handling a function of only two parameters, which is very few, and tens of function evaluations are necessary. With more parameters, say hundreds or thousands (this is at least what one usually has in function estimation), zero-order gradient-free are still too expensive and thus cannot be used in practice; gradient-based optimizers should be prefered.

• When the function is likely to be convex, and when the cost function depends on some states -solution of partial di erential equations -, then the model itself is likely to be di erentiated. In such cases, gradient-based optimizers are to be chosen. Among those, with respect to the most basic steepest descent algorithm, the numerical e ort of implementing the conjugate gradient, or better the BFGS, is highly recommended.

• The example presented here, on the only two-dimensional Rosenbrock function, has demonstrate this result. Such conclusions are of course much solid when it comes to function estimation where higher dimensions are encountered.

Cost function gradient

We recall here that the function to be minimized is the cost function J (u), expressed in terms of the state u, but minimized with respect to the parameters ψ. We thus have the equality (by de nition) between the cost function and its reduced version: j(ψ) := J (u). The state u is related to the parameters ψ through an operator (which may be linear, or not) that combines the partial di erential equations along with the boundary conditions, initial conditions, etc. This operator is denoted as S for the state problem. To be concise, one writes down S(u, ψ) = 0,

where we have the mapping ψ → u(ψ). Often, the space (and time) is discretized so that the state operator S is approximated in some matrix formulation. In this case, we have R(u, ψ) = 0, with dim R = dim u. Note that u involved in (58) is continuous while u involved in R(u, ψ) = 0 is likely to be already discretized (using nite di erence, nite elements, etc.). We now need the de nition of the directional derivative of j(ψ) in the direction δψ (see de nition 6). Other kinds of derivatives can also be used, such as the Gâteaux or Fréchet derivatives, see [START_REF] Allaire | Analyse numérique et optimisation[END_REF] for technical de nitions.

De nition 6 (Directional derivative). Let a point ψ ∈ K and a direction φ ∈ K. One de nes (t) := ψ + tφ and the function J (t) := j( (t)). The directional derivative of j at the point ψ in the direction φ is:

j (ψ; φ) := J (0) = lim t→0 t>0 j(ψ + tφ) -j(ψ) t . (59)
It has been seen before (see eq. ( 51)) that we have the equality

j (ψ; δψ) = u -u d , u (ψ; δψ) X , (60)
and, because of linearity of both u (ψ; δψ) and j (ψ; δψ) in δψ:

j (ψ; δψ) = (∇j(ψ), δψ) Z . ( 61
)
where Z is most of the time chosen equal to Y but it can be chosen di erently for regularization purposes.

Finite di erence

Function estimation in IHTP

Finite di erence

The nite di erence approach consists in approaching the cost function gradient through a substraction of the cost function with a perturbed cost function for the whole canonical base of ψ, that is δψ = δψ 1 , δψ 2 , . . . , δψ dim ψ . For the i th component, we have:

(∇j(ψ)) i = (∇j(ψ), δψ i ) Z ≈ j(ψ + δψ i ) -j(ψ) . (62)
Usually, in order to perform the same relative perturbation on all components ψ i , one rather uses i ← ε|ψ i |, where the positive scalar ε is xed. The very simple related algorithm is described in algorithm 6.

Algorithm 6:

The nite di erence algorithm to compute the gradient of the cost function Set the length ε > 0; At iteration p, compute the state u(ψ p ), compute j(ψ p ); foreach i = 1, . . . , dim ψ do Compute the cost j(ψ p + ε|ψ i |δψ i );

Set the gradient (∇j(ψ)) i ← j(ψ p + ε|ψ i |δψ i )j(ψ p ) ε|ψ i | end Integrate the gradient within the optimization methods that do not rely on the sensitivities (conjugate gradient or BFGS for instance among the presented methods)

In practice, the tuning parameter ε has to be chosen within a region where variables depend roughly linearly on ε. Indeed for too small values, the round-o errors dominate while for too high values one gets a nonlinear behavior. Even though the nite di erence method is easy to implement, it has the disadvantage of being highly CPU time consuming. Indeed, the method needs as many integrations of the model S(u, ψ) = 0 as the number of parameters, dim ψ. The gradient computed this way can be integrated to the previously presented optimization methods that do not rely on u , such as the conjugate gradient methods, or better the BFGS.

When performing the nite di erenciation with respect to ψ i , one also accesses the approximated perturbed state u (ψ; δψ i ). This way, one can use again the conjugate gradient methods or the BFGS method for instance, but also the Gauss-Newton-type methods based on matrix inversion and which do rely on the sensitivities u (ψ; δψ i ), i = 1, . . . , dim ψ. Doing so, the related optimization is given in algorithm 7.

Forward di erentiation

The forward di erentiation approach consists in computing u (ψ; δψ i ) di erentiating the state equations S(u, ψ) = 0, to get:.

S u (u, ψ)u + S ψ (u, ψ)δψ = 0. ( 63
)
As in the previous section, the gradient computation needs one integration of eq. ( 63) per parameter ψ i ; so one needs dim ψ integrations in total to access the full gradient ∇j(ψ). However, in this case, eq. ( 63) is linear, while eq. ( 58) was not linear.

As for the nite di erence approach, one may use the sensitivities u and integrate them into the Gauss-newton-type methods, or simply use the cost function gradient, and then use the methods that do not rely on the sensitivities.

When compared to the nite di erence approach, the forward di erence method leads to exact cost function gradient components. Moreover, though S is likely to be a nonlinear operator, the system given by eq. ( 63) is linear, thus yielding to much less CPU time. Another singularity is that the discrete version of S u (u, ψ), i.e. R u , is the tangent matrix that is to be used anyway for solving the "forward" problem

Forward di erentiation

Function estimation in IHTP Algorithm 7: The nite di erence algorithm to compute the gradient of the cost function and the sensitivities Set the step ε > 0; At iteration p, compute the state u(ψ p ), compute j(ψ p ); foreach i = 1, . . . , dim ψ do Compute the perturbed state u(ψ p + ε|ψ i |δψ i ) and the cost j(ψ p + ε|ψ i |δψ i );

Set the state sensitivity u (ψ;

δψ i ) ← u(ψ p + ε|ψ i |δψ i ) -u(ψ p ) ε|ψ i | ;
Set the gradient (∇j, δψ i ) with either (uu d , u (ψ; δψ i )) or as in previous algorithm with

j(ψ p + ε|ψ i |δψ i ) -j(ψ p ) ε|ψ i | .
end Integrate the gradient within the optimization methods that do not rely on the sensitivities (conjugate gradient or BFGS among the presented methods) or within optimization methods that do rely on the sensitivities (Gauss-Newton or Levenberg-Marquardt).

S(u, ψ) = 0. The computation of this linear tangent matrix is most often the task that takes the longer time in solving S(u, ψ) = 0. The optimized procedure is thus the one given in algorithm 8.

Algorithm 8:

The forward di erentiation algorithm to compute the cost gradient and the sensitivities At iteration p, solve iteratively S(u, ψ p ) = 0, compute j(ψ p ) and save the discrete version of the linear tangent operator S u (u, ψ p ); foreach i = 1, . . . , dim ψ do Solve S u (u, ψ)u + S ψ (u, ψ p )δψ i = 0; Set (∇j, δψ i ) Z = (uu d , u (ψ; δψ i ) X ; end Integrate the gradient within the optimization methods that do not rely on the sensitivities (conjugate gradient or BFGS among the presented methods) or within optimization methods that do rely on the sensitivities (Gauss-Newton or Levenberg-Marquardt). Note: the linear tangent matrix which is to be assembled for the solution of the "forward" model can be re-used for all canonical components δψ i .

Remark. Equation (63) is often called the sensitivity equation.

Example. Let us consider the unsteady heat conduction equation, with known heat capacity C, conductivity λ, volume source term f , initial condition u 0 , and Dirichlet condition on a part of the boundary, e.g. u 0 on ∂D D . The unknown ψ is the ux φ on the rest of the boundary, i.e. on ∂D N = ∂D\∂D D . The unperturbed and perturbed models are:

S(u, ψ) ≡        C ∂u ∂t -∇ • λ∇u = f in D u(x, t = 0) = u 0 in D u(x, t) = u 0 on ∂D D -∇u(x, t) • n = φ on ∂D N ; S(u + , ψ+ δψ) ≡        C ∂u + ∂t -∇ • λ∇u + = f in D u + (x, t = 0) = u 0 in D u + (x, t) = u 0 on ∂D D -∇u + (x, t) • n = φ + δψ on ∂D N
(64) Substracting the equations involved in these two models, dividing by , and searching the limit when → 0 gives:

lim →0 S(u + , ψ + δψ) -S(u, ψ) ≡            lim →0 C ∂ u + -u ∂t -∇ • λ∇ u + -u = 0 in D lim →0 u + -u (x, t = 0) = 0 in D lim →0 u + -u (x, t) = 0 on ∂D D lim →0 -∇ u + -u (x, t) • n = δψ on ∂D N (65)
that gives:

S u (u, ψ)u + S ψ δψ ≡        C ∂u ∂t -∇ • λ∇u = 0 in D u (x, t = 0) = 0 in D u (x, t) = 0 on ∂D D -∇u (x, t) • n = δψ on ∂D N (66)

Adjoint state

In this section we present the use of an additional problem -the so-called adjoint-state problem -that gives also the exact cost function gradient, but in a computational cheap way. We present one method based on the identi cation procedure (section 6.3.1), and another one that uses the Lagrange function (section 6.3.2). For the latter method, the model equation is treated as an equality constraint for the optimization. Both methods can deal with either the continuous equations or the discrete ones. One has to keep in mind that when the continuous method is used, all the obtained equations have later on to be discretized. Both strategies are equivalent in usual, but if the cost is computed through the integration of some discretized equations, then we consider that the discretized equations have to be di erented (it is the so-called "discretize-thendi erentiate" method). The other way is to deal with the continuous equations, then discretize the state model, etc. (it is the so-called "di erentiate-then-discretize" method). Some examples of adjoint derivation will be given in the last sections.

Identi cation method

In this rst part, we derive the adjoint-state method using the identi cation method. From the de nition of the functional gradient, one writes the gradient:

(∇j, δψ) Z = j (ψ; δψ) = u -u d , u (ψ; δψ) X . (67)
One then introduces a new variable (the adjoint-state variable u * ) such that the gradient equation given by eq. (67) also satis es the "easier-to-compute":

j (ψ; δψ) = S ψ (u, ψ)δψ, u * U . (68)
On the other hand, since we have the relationship S(u, ψ) = 0, then

S u (u, ψ)u + S ψ (u, ψ)δψ = 0 (69)
and thus, we have:

j (ψ; δψ) = -S u (u, ψ)u , u * U . (70)
Identifying eq. (67) and eq. ( 70), we obtain the adjoint-state problem that must satisfy the equality:

-S u (u, ψ)u , u * U = u -u d , u (ψ; δψ) X . (71)
Next, if the adjoint problem given by eq. ( 71) is satis ed (it means that we accessed the adjoint state u * ), then the cost function gradient is very simply given by eq. ( 68). We then use the inner product property The Lagrange function introduced in this section is a function of three variables, namely the state u , the parameter to be to be identi ed ψ, and the adjoint state variable u * . This means that both variables u and ψ are somehow considered to be independent, even though there exists, at least implicitly, the relationship S(u, ψ) = 0 that maps ψ to u. Moreover, since u is the solution of the forward model, then the Lagrange function L is always equal to the cost function J (u), and the constraints -which represent the partial di erential equations of the forward problem -are always satis ed. We now show that a necessary condition for the set ψ to be solution of the optimization problem eq. ( 1) is that there exists a set (u, ψ) such that (u, ψ, u * ) is a saddle point (stationary point) of L . Indeed, let us show that the necessary condition j (ψ; δψ) = 0, ∀δψ, is equivalent to:

∃ (u, u * , ψ) | L u (•) δw = 0; L u * (•) δw = 0; L ψ (•) δw = 0, (75)
for all directions δw taken in appropriate spaces (u , δu * and δψ). First, since the state is satis ed, then:

L u * = S(u, ψ) = 0.
Moreover, since we have also S u (u, ψ)u + S ψ (u, ψ)δψ = 0, we get:

L ψ (•) δψ = S ψ (u, ψ)δψ, u * U = -S u (u, ψ)u , u * U . (76)
In another hand, the di erentiation of the Lagrange function with respect to the state gives:

L u (•) u = u -u d , u X + S u (u, ψ)u , u * U . (77)
So far, the choice for the adjoint variables u * has not been xed yet. However, choosing the adjoint variable such that L u (•) u = 0 ∀u considerably simpli es the relationship between the di erentiated lagrangian with respect to ψ and the cost function gradient. One actually chooses u * such that it satis es the adjoint-state equation:

S u (u, ψ)u , u * U + u -u d , u (ψ; δψ) X = 0. (78)
This way we obtain the cost function gradient:

L ψ (•) δψ = u -u d , u (ψ; δψ) X = j (ψ; δψ) = (∇j, δψ) Y (79)
The adjoint-state equation is thus:

S * (u, ψ)u * + (u -u d ) = 0, (80)
and the gradient is given by:

∇j = S ψ (u, ψ), u * Y . (81)
Summarizing, the minimum of the cost function is to be found at the stationary point of the Lagrange function eq. ( 74). When the adjoint-state equation eq. ( 80) is satis ed, then the components of the cost function gradient are simply given through the inner product eq. (81).

Examples

In the examples presented below, we do not specify what the parameters are. We just give the form of the adjoint-state problem related to the "forward" state problem form.

The global optimization algorithm

Function estimation in IHTP

Case of parabolic PDE The discretization of the space and time dependent di usive model yields to the so-called parabolic problem. It is somehow the union between both just above presented cases:

S (u, ψ) = C u -∆u -B = 0 for t ∈ I, x ∈ D u = u 0 for t = 0, (87)
with associated boundary conditions. Injecting the di erentiated operators involved in eq. ( 87) into the adjoint eq. ( 78) gives:

C u , ψ U -∆u , ψ U + u -u d , u X = 0 with (a, b) U = T D ab dx dt.
Transposing all operators through integration by parts (once in time and twice space) gives:

-u , C u * U + u , C u * D T 0 - I [...] ∂D dt + u -u d , u (ψ; δψ) X = 0
Eventually, the adjoint problem becomes:

-C u * -∆u * + J (u) = 0 for t ∈ I u * = 0 for t = t f . (88)
along with associated spatial boundary conditions.

Remark. A more detailed example given later on in section 9.2 provides the full calculation of the boundary condition for a similar case.

The global optimization algorithm

The general algorithm is given in algorithm 10. The global procedure described in this algorithm is run until (at least) one of the stopping criteria presented in section 2.5 is reached.

Continuous gradient and discretized continuous gradient

In previous examples as well as in the derivation of both forward di erentiated and adjoint-state models, all derivations were performed on continuous equations. To be solved, such equations will have later on to be discretized, for example with nite elements or any other method. This ordinary process yields the so-called continuous gradient. For example, refering to g. 12, the continuous state equation for the continuous variable u is rst di erentiated, yielding a continuous di erentiated state u , solution of a continuous di erentiated partial di erential equation. Then, after discretization (which is an approximation process), one has the discretized di erentiated state (u ) h , such that the discretized continuous gradient j | DCG can be computed.

The other way round consists in rst discretizing, then di erentiating. All partial di erential equations are discretized, so that the forward state is u h , and the corresponding cost function, based on this approximated state is j h . This approximated model can then be di erentiated so that the derivative we get is (u h ) , which gives the discrete gradient j | DG .

Both gradients j | DCG and j | DG are di erent because the approximations are not performed on the same operators. The cost function is always j h because a numerical solver is used to compute u h . The minimum of j h corresponds to j | DG = 0 but at this minimum, it is likely that we have j | DCG = 0. This means that the discrete gradient is compatible with the cost function which is calculated while the discretized continuous gradient is not. However, one has to keep in mind that errors coming from discretization are likely to be negligible when compared to measurement errors of the inverse problem. As such, the computation of discretized continuous gradients is -according to the author -the better strategy because all derivations are performed on partial di erential equations, and di erentiated models are also partial di erential equation very similar to model equation, and such similarity is the easiest way to go: similar equation, re-use the forward solver for the di erentiated model or for the adjoint-state, etc.

u 1. discretize u h → j h 2. di erentiate (u h ) → j | DG 1. di erentiate u 2. discretize (u ) h → j | DCG

Elements of comparison

We give in this section some elements of comparison between the previously presented optimization algorithms and between the di erent gradient computation strategies.

Convergence speed

The optimization algorithms presented in section 5 yield to a series {ψ k } k≥1 that converges to ψ. Hereafter are some some convergence rate de nitions [START_REF] Nocedal | Numerical Optimization[END_REF][START_REF] Minoux | Mathematical Programming, Theory and Applications[END_REF].

De nition 7. The convergence rate of the series {ψ k } k≥1 is said to be linear if

ψ k+1 -ψ k ψ k -ψ ≤ τ, τ ∈ (0, 1). ( 89
)
This means that the distance to the solution ψ decreases at each iteration by at least the constant factor τ .

De nition 8. The convergence rate of the series {ψ k } k≥1 is said to be superlinear in n steps if

lim k→∞ ψ k+n -ψ k ψ k -ψ = 0. ( 90
)

Gradient computation cost

Function estimation in IHTP De nition 9. The convergence rate of the series {ψ k } k≥1 is said to be quadratic if

ψ k+1 -ψ k ψ k -ψ 2 ≤ τ, τ > 0. (91)
Quasi-Newton methods usually converge super-linearly and the Newton method converges quadratically. The steepest descent method converge linearly. Moreover, for ill-posed problems, this method may converge linearly with a constant τ close to 1. Next, the conjugate-gradient method converges superlinearly in n steps to the optimum [START_REF] Minoux | Mathematical Programming, Theory and Applications[END_REF].

Thus the quasi-Newton methods convergence-rate is much higher than the conjugate gradient methods convergence-rate which need approximately n times more steps (n times more line-search) at the same convergence behavior. However, for the quasi-Newton method, the memory place is proportionnal to n 2 .

Gradient computation cost

Let S(u, ψ) = 0 the state problem that maps ψ → u, R being possibly nonlinear (for highlighting di erences between the distinct strategies), and dim ψ the number of parameters to be evaluated. We compare the number of times the model S, the di erentiated model and/or the adjoint-state model are computed to access the full gradient of the cost function.

Finite di erence method:

(dim ψ + 1) nonlinear computation of S(u, ψ) = 0.

Forward di erentiation method:

1 nonlinear computation of S(u, ψ) = 0, dim ψ linear computation of S u (u, ψ)u + S ψ (u, ψ)δψ = 0.

Adjoint state method:

1 nonlinear computation of S(u, ψ) = 0, 1 linear computation of S * (u, ψ)u * + uu d = 0.

Thus, the nite di erence method is very time consuming, though it is easy to use. Next, comparing the two latter methods, the operator involved in the adjoint-state method is almost the same as the one involved in the forward di erentiation method, though the adjoint-state method yields to higher algorithmic complexity (backward time integration, memory, etc.). When dim ψ is high (even if dim ψ is bigger than say 100), the use of the direct di erentiation method becomes cumbersome and computationally expensive; the adjoint-state method is, in fact, the only acceptable method.

Gradient computation needs

We recall in the following table the way (the required needed steps) one computes the cost function gradient.

Steepest, conjugate-gradients, Newton

Gauss-Newton, BFGS, DFP, . . . Levenberg-Marquardt, . . .

u ← S(u, ψ) = 0 u ← S(u, ψ) = 0 u ← S(u, ψ) = 0 j ← u j ← u j ← u ∇j ←    Forward di . or Adjoint state ∇j ←    Forward di . or Adjoint state ∇j ← S t S ← S ← u (Forward di .) ∇ 2 j (complicated)
When the inverse problem is ill-posed (which is likely to be the case in real cases, especially when the control space dimension is big), regularization is needed and sometimes compulsary. Regarding function estimation, for instance space-dependent physical properties or sources, speci c regularization strategies di erent from the ones used in parametric estimation are required. Regularization may be viewed as adding a priori information, but other means can also be used, including (see [START_REF] Dubot | Regularization opportunities for the di use optical tomography problem[END_REF] for elements of comparison on applications of optical tomography):

• choose of speci c X -norm for the cost function expression according to the prior knowledge of the unknown (use for instance the L 1 (D)-norm instead of the ordinary L 2 (D)-norm).

• add prior information through Tikhonov penalization, If some Tikhonov-type regularization terms are added to the cost function, the cost function j (ψ) := J (u) + J + (ψ) is the one to be minimized, with:

J + := Dψ 2 Y ( 92
)
where D is often a di erential operator acting on the function ψ and • Y is another norm to be de ned according to the chosen control space.

• choose an appropriate Z-norm for extracting the cost function gradient. The use of speci c inner products when extracting the cost function gradient is a recent regularization tool. In order to present this regularization strategy, let us work on the example where a space-dependent physical property in D is to be estimated. In such a case it is usual to use the ordinary L 2 (D)-inner product, i.e., one uses j (ψ; φ) = (∇j, φ) L 2 (D) ; this gives the ordinary L 2 (D) cost function gradient, denoted here as ∇ L 2 j(ψ). Besides, the Sobolev inner product can give much better (smoother) results when the noise has propagated to the adjoint-state variable and then to the cost function gradient. Even better, the weighted version has recently proven to give excellent results. This one de ned as:

(u, v) Z = (u, v) H 1( ) (D) := D uv + 2 ∇u • ∇v dx (93)
is used in the cost function gradient extraction relationship j (ψ; φ) = (∇j, φ) H 1( ) (D) in order to compute the weighted Sobolev cost function gradient ∇ H 1( ) j(ψ).

• choose an appropriate functional space for the control space parameterization. In practice, the control space must be approximated in order to be nite. Often, the nite element method is used so that one searches ψ that belongs to a nite dimensional subspace, say V. Let us consider a triangulation M of the computational domain D, and let us note n p the number of vertices in M. It has been shown, through numerical means on a speci c OT problem that, among the large number of tested possibilities, the piecewise linear continuous functions (dim ψ = n p ) are the most appropriate for the estimation of space-dependent functions.

• choose an appropriate dimension dim ψ of the control space parameterization. Usually the nite element space used to solve the forward model (58) has to be ne enough to ensure that numerical errors stay small enough. Most often, the triangulation chosen for the control space is the one chosen for the state. It has been shown again, through numerical means, that both the convergence and the quality of the reconstructions are much improved when dim ψ is lowered, up to a certain limit, at least for quasi-Newton algorithms.

• Multi-scales approaches is also a fabulous opportunity to regularize solutions and in the same time accelerate the convergence and avoid converging to local minima. Coupled with wavelets on one side, and the BFGS in the other side, this method relies on a reformulation of the original inverse problem into a sequence of sub-inverse problems of di erent scales using wavelet transform, from the largest scale to the smallest one. Successful applications of this method include the estimation of space-dependent absorption and scattering coe cients in optical tomography [START_REF] Dubot | A wavelet multi-scale method for the inverse problem of di use optical tomography[END_REF].

Examples

Parametric conductivities in a transient heat conduction problem

This rst simple example deals with the estimation of uniform conductivity coe cients in di erent subdomains. Heat transfer is considered. Initial temperature is assumed to be known and equal to T 0 . T 0 is also the Dirichlet temperature for positive time on the whole boundary ∂D. The domain has the shape of a head with two eyes, one nose and one mouth. The geometry being known, as well as the initial and boundary conditions, the heat capacity and the time-dependent heat source, the inverse heat conduction problem consists in estimating, through infra-red like temperature measurements on D × I, the set of conductivities λ i , i = 1, • • • , 4 (1, 2, 3, 4 corresponding to the left eye, the right eye, the noze and the mouth, respectively). Noisy (1 % white noise) synthetic data was generated with conductivities equal to 20, 30, 40 and 50, respectively. Guessed conductivities were all equal to 10.

If optimization methods based on sensitivities are chosen, one will have to compute, sucessively:

ρC ∂T ∂t -∇ • λT = ∇ • λ ∇T (94)
with null initial and boundary conditions. In the sensitivity model, the direction λ equals 0 or 1 depending on the location for the four considered sensitivity problems. Corresponding sensitivities are presented in g. [START_REF] Larrouturou | Méthodes mathématiques pour les sciences de l'ingénieur: optimisation et analyse numérique[END_REF].

With so few parameters to identify (4 in total in this example), it is not really necessary to use the adjoint-state method to compute the cost gradient. We however give in g. 13 the evolution of the adjointstate variable which is solved backwardly from nal time to initial time, while integrating along time the errors integrated within the cost function (this rst application was actually chosen for this purposes : small number of unknowns, and possible visualization).

From the knowledge of these temperature sensitivities, one can compute the sensitivity matrix S such that s i×k,j gathers for instance the sensitivity of temperature on the ( nite element) node i at time k with respect to λ j . In the same manner, the error vector e i×k gathers the error (di erence between the predictions and the measurements) at the (FE) node i, and at time t k . Consequently, the cost gradient is computed straightforwardly through ∇j = S e, and the Gauss-Newton algorithm can be used without any regularization because this parametric problem is not ill-posed. Very few Gauss-Newton iterations are needed to converge as can be seen in g. 14.

Space-dependent convection coe cient in a transient heat conduction problem

In this section, we consider an application of a nonlinear transient heat transfer inverse problem arising in thermal treatment for instance. D being an open bounded set of R 2 and I =]0, t f ], the modeling equation in D × I is where the temperature-dependent physical properties are considered. We also consider the following initial and boundary conditions with ∂D 1 ⊕ ∂D 3 ⊕ ∂D 3 forming a partition of ∂D:

C ∂T ∂t -∇ • (λ∇T ) = f for (x, t) ∈ D × I (95)
λ i , i = 1, • • • , 4 j(λ i ) Iterations λ 1 λ 2 λ 3 ] λ 4 j(λ)
T = T 0 at t = 0 -λ∇T • n = h (T -T ∞ ) for x ∈ ∂D 1 ∇T • n = 0 for x ∈ ∂D 2 -λ∇T • n = εσ T 4 -T 4 ∞ for x ∈ ∂D 3 (96)
The estimation of the heat transfer function h(x, t), x ∈ ∂D 1 , t ∈]0, t f ] is performed through the minimization of the cost function:

j(h) := J (T ) = t f 0 N j=1 (T (x j , t) -T d (x j , t)) 2 dt ( 97
)
where T (x j , t) and T d (x j , t) represent respectively the predicted and measured temperatures at N various locations x := (r j , z j ) in the material. For such application, the minimization can be carried out by using conjugate gradients or better quasi-Newton methods. In any case, the optimization is based on the gradient computation.

The cost function gradient is obtained for all values x ∈ ∂D 1 , t ∈]0, t f ] by the following relationship:

∇j(h) = T * × (T -T ∞ ) (98)
where T * is the solution of the adjoint problem:

-C ∂T * ∂t -∇ • (λ∇T * ) = j (T -T d ) × δ (x -x j ) (99)
with the condition T * = 0 at nal time t f and the conditions on the boundaries:

-λ∇T * • n = hT * for x ∈ ∂D 1 ∇T * • n = 0 for x ∈ ∂D 2 -λ∇T * • n = 4εσT 3 T * for x ∈ ∂D 3

• n + h (T -T ∞ ) , γ) U 1 + (∇T • n, ξ) U 2 + λ∇T • n + εσ T 4 -T 4 ∞ , U 3 (102)
The di erentiated Lagrange function with respect to h in the direction δh is:

(L h (•), δh) = (T -T d , T ) X + C ∂(T ) ∂t -∇ • λ∇T , T * U + (λ∇T ) • n + hT + δh(T -T ∞ ), γ) U 1 + (∇T • n, ξ) U 2 + λ∇T • n + 4εσT 3 T , U 3 (103)
We then use the following integrations by parts to express some particular terms: Choosing γ = T * on ∂D 1 , ξ = λT * on ∂D 2 and = T * on ∂D 3 , the adjoint problem can eventually be written as:

           -C ∂T * ∂t -λ∆T * = -j (T -T d ) × δ(x -x j ) x ∈ D, t ∈ I T * = 0 x ∈ D, t = t f -λ∇T * • n = hT * x ∈ ∂D 1 , t ∈ I ∇T * • n = 0 x ∈ ∂D 2 , t ∈ I -λ∇T * • n = 4 σT 3 T * x ∈ ∂D 3 , t ∈ I (106)
and the cost gradient is written as:

∇j = -(T -T ∞ ) T * . ( 107
)
From the integration of the adjoint-state, the cost function gradient is computed. From the knowledge of the cost function gradient, the direction of descent is computed, for instance with the conjugate gradient method, or with any other faster method if a ne parameterization for h is required. It is also to be pointed out that the temperature state being varying almost linearly with the convection property, the line-search equation can be for instance given by the solution of (14).

Adjoint RTE

Function estimation in IHTP

Adjoint RTE

This last example aims at developping adjoint-state equation of the radiative transfer equation (RTE). The main objective behind this developments is the solution of optical tomography problems, in which the problem is the reconstruction of radiative properties (κ(x) and σ(x)) within the medium, input intensity being prescribed on the boundary, and measurements being also performed on boundaries. Some of the di culties for solving such problems include: i) the dimension of the discrete control space is likely to be high, in 2-D and especially in 3-D. This means that e cient optimizers such as the ones based on the gradient-type BFGS are the only ones to be used. Others such as the conjugate gradients for instance may be too slow and gradient-free are not appropriate at all;

ii) the RTE is integro-di erential, so appropriate inner products must be used through all the derivations. Therefore, mathematical developments for the derivation of the adjoint-state as well as for the cost function gradient must be undertaken very carefully. In the same spirit, numerical algorithm and implementation must be undertaken very carefully;

iii) the state being non linear with respect to the physical properties and overall the nonlinear inverse problem being ill-posed, several regularization strategies must be used and combined together.

Let the radiative transfer equation (RTE) being written as, ∀ (x, s) ∈ Dπ:

(s • ∇ + κ + σ) I(x, s) = σ 4π I(x, s)Φ(s, s ) dω(s ), (108)
where s is the considered direction of propagation, Φ(s, s ) is the phase function representing the probability that a photon arriving from the direction s is scattered to the direction s, and κ and σ are the absorption and di usion space-dependent functions, respectively. On a part of the boundary, there is a prescribed Dirichlet condition: I(x, s) = Ī for x ∈ ∂D s and s • n < 0.

Also, let a cost function measuring the mis t between predictions and measurements somewhere on the boundary, i.e., ∂Ω d ⊂ ∂Ω, the mis t being expressed (it is actually a norm) in terms of the radiance, e = I(x, s d ) -I d (x, s d ) for x ∈ ∂D d and s d • n > 0.

(110)

In order to make the derivation of the adjoint-state, the tools described in previous sections are used. Additionally, the state variable I being de ned in (x, s) ∈ D × 4π, the inner product U de ned in eq. ( 68) and in following equations is:

(u, v) U = 4π D uv dxds. (111)
After integration by parts and some -technical -manipulations in the inner products, one nds the adjoint RTE to be: This continuous version of the components of the cost function gradient is then projected onto the basis used to parameterize the control-space. A very detailed derivation of the adjoint RTE can be found in [START_REF] Favennec | Space-dependent sobolev gradients as a regularization for inverse radiative transfer problems[END_REF], for instance.

Concluding remarks

This lecture was devoted to the presentation of mathematical and numerical algorithms used in the estimation of functions while solving inverse heat transfer problems. Contrarily to parameter estimation problem, the dimension of the "control space" is likely to be big after the process of parameterization, this being an essential issue and reason of why using e cient optimization algorithms. Among those e cient algorithms, the ones based on the cost function gradient as well as on adjoint-states are of rst importance. Functions to be estimated usually contain several regularity requirements and thus, the use of e cient regularization tools are compulsory to cope with the ill-posed character of the inverse problem.

For

  indicating that the minimum is obtained, one should prefer the notations φ = arg min ψ∈K⊂V j(ψ) and j(φ) = min ψ∈K⊂V j(ψ)

Figure 1 :Figure 2 :

 12 Figure 1: Convex and non-convex domaine K

Figure 4 :

 4 Figure 4: Basic operations on a simplex for n = 2. From left to right: re ection, expansion, contraction, and inside contraction.

Figure 5 :

 5 Figure 5: PSO algorithm: a particle displacement.

Figure 6 :

 6 Figure 6: Schematic representation of the Newton-Raphson line search method.

4. 5

 5 Other methods -Inexact line-search Function estimation in IHTP with u ψ, d k = u the derivative of u at the point ψ k and in the direction d k . This equation gives straightforwardly: u

Figure 7 :

 7 Figure 7: Set of points satisfying the Goldstein's rules: [a, b] ∪ [c, d].

  PSO algorithm: ≈ 6,000 cost function evaluations. Conjugate gradients descent algorithm: 45 cost function evaluations.

  BFGS descent algorithm: 11 cost function evaluations.

Figure 10 :

 10 Figure 10: Numerical comparison of optimizers on the 2-D Rosenbrock function.

Algorithm 10 :Compute the direction of descent 4 .

 104 The global optimization algorithm 1. Integrate the cost function value through integration of the forward (maybe nonlinear) problem; Store all state variables to reconstruct the tangent matrix (or store the tangent matrix); 2. Integrate the backward linear adjoint-state problem, all matrices being possibly stored or recomputed from stored state variables 3. Compute the cost function gradient; Solve the line-search algorithm

Figure 12 :

 12 Figure 12: Discrete gradient vs discretized continuous gradient

Figure 13 :

 13 Figure 13: Variables needed to solve the parametric inverse transient heat conduction problem. Columns correspond to increasing time from 0 to 40 by steps of 5.The rst raw presents the source term that follows, in this particular case, a lemniscate. The second raw presents the temperature evolution. The four following raws present the evolution of sensitivities with respect to λ

Figure 14 :

 14 Figure 14: Evolution of the parameters and of the cost function with respect to the G-N iterations.

  Remark. The following notations are used:U := L 2 (I; L 2 (D)), U i := L 2 (I; L 2 (∂D i )), ∀i = 1, 2, 3 and U * := L 2 (I; L 2 (∪ i ∂D i )).

9. 2

 2 Space-dependent convection coe cient in a transient heat conduction problem Function estimation in IHTP Proof. The derivative of the state T at the point h and towards δh, T (h; δh) is de ned by: ∂t -∇• λ∇T = 0 x ∈ D, t ∈ I T = 0 x ∈ D, t = 0 λ∇T • n + hT + δh(T -T ∞ ) = 0 x ∈ ∂D 1 , t ∈ I ∇T • n = 0 x ∈ ∂D 2 , t ∈ I λ∇T • n + 4εσT 3 T = 0 x ∈ ∂D 3 , t ∈ I(101)The Lagrange function is formally de ned as:L (T, {T * , γ, ξ, }, h) = J (T ) + C ∂T ∂t -∇ • (λ∇T )f, T * U + (λ∇T

=

  T , -C ∂T * ∂t U + (CT , T * ) D (t = t f ) -(CT , T * ) D (t = 0) (λ∆T , T * ) U = (λ∆T * , T ) U + (λ∇T * • n, T ) U * -(λT * , ∇T • n) U *(104)We bring together similar terms to get:(L h (T, {T * , γ, ξ, }, h), δh) = (T -T d , T ) X + (δh (T -T ∞ ) , γ) U 1 + -C ∂T * ∂t -λ∆T * , T U + (CT * , T ) D (t = t f ) + (λ∇T * • n, T ) U * -(λT * , ∇T • n) U * + (λ∇T ) • n + hT , γ) U 1 + (∇T • n, ξ) U 2 + λ∇T • n + 4εσT 3 T , U 3 (105)

(

  -s • ∇ + κ + σ) I * (x, s) = σ 4π I * (x, s)Φ(s, s ) dω(s )(112)coupled with the Dirichlet boundary condition:I * (x, -s d ) = (s d • n) -1 (I -I d ) (x, -s d ) for x ∈ ∂D d and s d • n > 0(113)The adjoint-state being computed, the cost function gradient, here as an implicit function of x can be computed:∇ κ j(x) =4π I(x, s)I * (x, s) ds ∇ σ j(x) = 4π I(x, s)I * (x, s) -4π I(x, s )Φ(s, s ) ds I * (x, s) ds (114)

(the left eye), λ

(the right eye), λ

(the noze), and λ

(the mouth), respectively. The last raw presents the evolution of the adjoint-state: it is null at nal time and then increases due to cost function integration while going back to initial time.

(A u, v) = (u, A * v) where A * is the transposed conjugate operator of A (adjoint) to modify the adjoint equation given by eq. (71) to:

where S * is the conjugate transposed of the linear tangent operator S u , i.e., we used:

where the term [• • • ] may contain some additional terms coming from some integrations by parts. Figure 11 schematically represents the process of identi cation method. Integrate the gradient within the optimization methods that do not rely on the sensitivities (conjugate gradient or BFGS among the presented methods)

Lagrange formulation

The use of a Lagrange formulation means that the state equations are taken as constraints in the optimization problem. To do so, let us introduce the Lagrange function [START_REF] Céa | Optimisation, théorie et algorithmes[END_REF][START_REF] Larrouturou | Méthodes mathématiques pour les sciences de l'ingénieur: optimisation et analyse numérique[END_REF]:

Function estimation in IHTP Case of ODE Let us start with the case where the state model is simpli ed to a single linear continuous ordinary di erential equations integrated in time I = (0, t f ]. The forward problem thus writes:

where C is an inertial scalar term and B contains the loadings. Injecting the di erentiated time-dependent relationship eq. ( 82) into the adjoint-state relationship eq. ( 78) gives:

C u , u * U + uu d , u X = 0 where the inner must be understood as (a, b) U = I ab dt. One then integrates by part the rst term to get:

Since there is no reason that the initial state depend on the parameters ψ (except if the initial state is searched), then the derivatives u of u at initial time is zero. The adjoint-state problem is eventually:

Remark. There is a minus sign just before the operator C involved in the rst equation. At the same time, the boundary-time condition is given at nal time t f . Therefore, when considering these two points, there is no way to solve the adjoint problem forwardly, i.e., from t = 0 to t f . The trick consists in introducing a new time variable τ = t ft (the dual time). Doing so, the initial condition is given at the initial time τ = 0, and the time-dependent equation eq. ( 83) is solved in the forward way in the dual time variable τ -which corresponds to the backward way in the primal time variable t.

Remark. The loading component (uu d ) involved in eq. ( 83) is non-zero only at times where the cost function j is to be integrated, i.e., in accordance with the de nition of the X -norm.

Remark. Inherently, the adjoint-state problem is linear: even though the forward problem is likely to be nonlinear (it was not the case in the considered exemple), the adjoint-state problem is still linear since the operators do not depend on the adjoint-state variables. An equivalent remark was given for the forward di erentiation method which used the linear tangent operator.

Case of elliptic PDE This second example concerns the case where the state model is simpli ed to a di usive-type continuous partial di erential equation independent of time:

Injecting the di erentiated space-dependent relation eq. ( 84) into the adjoint eq. ( 78) gives:

-λ∆u , u * U + uu d , u X = 0. with (a, b) U = D ab dx. Using twice the Green theorem on the rst integral, one gets:

Owing to be veri ed for all directional derivatives u , the general adjoint-state problem becomes:

Remark. The second term, (. . .) ∂U comes in eq. (85) because of the integration by parts. These terms depend on the boundary conditions associated to eq. (84) that formed the complete forward model. Taking into account of these terms will also complete the de nition of the adjoint-state model, yielding the boundary conditions associated to the adjoint-state equation eq. (86).

Remark. The loading component uu d involved in the space-dependent equation is non-zero only at the selected locations where the cost function j is to be integrated, i.e., in accordance with the de nition of the X -norm.