## Section: Scientific Foundations

### Identification and approximation

Identification typically consists in approximating experimental data by the prediction of a model belonging to some model class. It consists therefore of two steps, namely the choice of a suitable model class and the determination of a model in the class that fits best with the data. The ability to solve this approximation problem, often non-trivial and ill-posed, impinges on the effectiveness of a method.

Particular attention is payed within the team to the class of stable linear time-invariant systems, in particular resonant ones, and in isotropically diffusive systems, with techniques that dwell on functional and harmonic analysis. In fact one often restricts to a smaller class —e.g. rational models of suitable degree (resonant systems, see section  4.3 ) or other structural constraints— and this leads us to split the identification problem in two consecutive steps:

1. Seek a stable but infinite (numerically: high) dimensional model to fit the data. Mathematically speaking, this step consists in reconstructing a function analytic in the right half-plane or in the unit disk (the transfer function), from its values on an interval of the imaginary axis or of the unit circle (the band-width). We will embed this classical ill-posed issue (i.e. the inverse Cauchy problem for the Laplace equation) into a family of well-posed extremal problems, that may be viewed as a regularization scheme of Tikhonov-type. These problems are infinite-dimensional but convex (see section 3.1.1 ).

2. Approximate the above model by a lower order one reflecting further known properties of the physical system. This step aims at reducing the complexity while bringing physical significance to the design parameters. It typically consists of a rational or meromorphic approximation procedure with prescribed number of poles in certain classes of analytic functions. Rational approximation in the complex domain is a classical but difficult non-convex problem, for which few effective methods exist. In relation to system theory, two specific difficulties superimpose on the classical situation, namely one must control the region where the poles of the approximants lie in order to ensure the stability of the model, and one has to handle matrix-valued functions when the system has several inputs and outputs, in which case the number of poles must be replaced by the McMillan degree (see section 3.1.2 ).

When identifying elliptic (Laplace, Beltrami) partial differential equations from boundary data, point 1. above can be recast as an inverse boundary-value problem with (overdetermined Dirichlet-Neumann) data on part of the boundary of a plane domain (recover a function, analytic in a domain, from incomplete boundary data). As such, it arises naturally in higher dimensions when analytic functions get replaced by gradients of harmonic functions (see section  4.2 ). Motivated by free boundary problems in plasma control and questions of source recovery arising in magneto/electro-encephalography, we aim at generalizing this approach to the real Beltrami equation in dimension 2 (section 6.3.3 ) and to the Laplace equation in dimension 3 (section 6.3.1 ).

Step 2. above, i.e., meromorphic approximation with prescribed number of poles—is used to approach other inverse problems beyond harmonic identification. In fact, the way the singularities of the approximant (i.e. its poles) relate to the singularities of the approximated function is an all-pervasive theme in approximation theory: for appropriate classes of functions, the location of the poles of the approximant can be used as an estimator of the singularities of the approximated function (see section 6.3.2 ).

We provide further details on the two steps mentioned above in the sub-paragraphs to come.

#### Analytic approximation of incomplete boundary data

Keywords : extremal problems, inverse problems, Hardy spaces, harmonic functions, Beltrami equations.

Participants : Laurent Baratchart, Slah Chaabi, Yannick Fischer, Juliette Leblond, Jean-Paul Marmorat, Jonathan Partington, Stéphane Rigat [ Univ. Aix-Marseille I ] , Emmanuel Russ [ Univ. Aix-Marseille III ] , Fabien Seyfert.

Given a planar domain D , the problem is to recover an analytic function from its values on a subset of the boundary of D . It is convenient to normalize D and apply in each particular case a conformal transformation to meet a “normalized” domain. In the simply connected case, which is that of the half-plane, we fix D to be the unit disk, so that its boundary is the unit circle T . We denote by Hp the Hardy space of exponent p which is the closure of polynomials in the Lp -norm on the circle if 1p< and the space of bounded holomorphic functions in D if p = . Functions in Hp have well-defined boundary values in Lp(T) , which make it possible to speak of (traces of) analytic functions on the boundary.

A standard extremal problem on the disk is  [68] :

(P0 )  Let 1p and fLp(T) ; find a function gHp such that g-f is of minimal norm in Lp(T) .

When seeking an analytic function in D which approximately matches some measured values f on a sub-arc K of T , the following generalization of (P0 ) naturally arises:

(P )  Let 1p , K a sub-arc of T , fLp(K) , and M>0 ; find a function gHp such that and g-f is of minimal norm in Lp(K) under this constraint.

Here is a reference behavior capsulizing the expected behavior of the model off K , while M is the admissible error with respect to this expectation. The value of p reflects the type of stability which is sought and how much one wants to smoothen the data.

To fix terminology we generically refer to (P ) as a bounded extremal problem . The solution to this convex infinite-dimensional optimization problem can be obtained upon iteratively solving spectral equations for appropriate Hankel and Toeplitz operators, that involve a Lagrange parameter, and whose right hand-side is given by the solution to (P0 ) for some weighted concatenation of f and . Constructive aspects are described in  [43] , [45] , [87] , for p = 2 , p = , and 1<p< , while the situation p = 1 is essentially open.

Various modifications of (P) have been studied in order to meet specific needs. For instance when dealing with loss-less transfer functions (see section 4.3 ), one may want to express the constraint on in a pointwise manner: |g-|M a.e. on , see [29] , [47] for p = 2 and = 0 .

The above-mentioned problems can be stated on an annular geometry rather than a disk. For p = 2 the solution proceeds much along the same lines [23] . When K is the outer boundary, (P ) regularizes a classical inverse problem occurring in nondestructive control, namely to recover a harmonic function on the inner boundary from overdetermined Dirichlet-Neumann data on the outer boundary (see sections  4.2 and 6.3 ). Interestingly perhaps, it becomes a tool to approach Bernoulli type problems for the Laplacian, where overdetermined observations are made on the outer boundary and we seek the inner boundary knowing it is a level curve of the flux (see section  6.3.3 ). Here, the Lagrange parameter indicates which deformation should be applied on the inner contour in order to improve the fit to the data.

Continuing effort is currently payed by the team to carry over bounded extremal problems and their solution to more general settings.

Such generalizations are twofold: on the one hand Apics considers 2-D diffusion equations with variable conductivity, on the other hand it investigates the ordinary Laplacian in . The targeted applications are the determination of free boundaries in plasma control and source detection in electro/magneto-encephalography (EEG/MEG, see section 6.3.2 ).

An isotropic diffusion equation in dimension 2 can be recast as a so-called real Beltrami equation  [73] . This way analytic functions get replaced by “generalized” ones in problems (P0 ) and (P ). Hardy spaces of solutions, which are more general than Sobolev ones and allow one to handle Lp boundary conditions, have been introduced when 1<p<   [46] . The expansions of solutions needed to constructively handle such problems have been preliminary studied in  [64] , [65] . The goal is to solve the analog of (P ) in this context to approach Bernoulli-type problems (see section 6.3.1 ).

At present, bounded extremal problems for the n -D Laplacian are considered on half-spaces or balls. Following [88] , Hardy spaces are defined as gradients of harmonic functions satisfying Lp growth conditions on inner hyperplanes or spheres. From the constructive viewpoint, when p = 2 , spherical harmonics offer a reasonable substitute to Fourier expansions [13] . Only very recently were we able to define operators of Hankel type whose singular values connect to the solution of (P0 ) in BMO norms. The Lp problem also makes contact with some nonlinear PDE's, namely to the p -Laplacian. The goal is here to solve the analog of (P ) on spherical shells to approach inverse diffusion problems across a conductor layer.

#### Meromorphic and rational approximation

Keywords : meromorphic approximation, rational approximation, critical point theory, orthogonal polynomials.

Participants : Laurent Baratchart, José Grimm, Martine Olivi, Edward Saff, Herbert Stahl [ TFH Berlin ] .

Let as before D designate the unit disk, T the unit circle. We further put RN for the set of rational functions with at most N poles in D , which allows us to define the meromorphic functions in Lp(T) as the traces of functions in Hp + RN .

A natural generalization of problem (P0 ) is

(PN )  Let 1p , N0 an integer, and fLp(T) ; find a function gNHp + RN such that gN-f is of minimal norm in Lp(T) .

Problem (PN ) aims, on the one hand, at solving inverse potential problems from overdetermined Dirichlet-Neumann data, namely to recover approximate solutions of the inhomogeneous Laplace equation u = , with some (unknown) distribution, which will be discretized by the process as a linear combination of N Dirac masses. On the other hand, it is used to perform the second step of the identification scheme described in section 3.1 , namely rational approximation with a prescribed number of poles to a function analytic in the right half-plane, when one maps the latter conformally to the complement of D and solve (PN ) for the transformed function on T .

Only for p = and continuous f is it known how to solve (PN ) in closed form. The unique solution is given by the AAK theory, that allows one to express gN in terms of the singular vectors of the Hankel operator with symbol f . The continuity of gN as a function of f only holds for stronger norms than uniform, [85] .

The case p = 2 is of special importance. In particular when , the Hardy space of exponent 2 of the complement of D in the complex plane (by definition, h(z) belongs to if, and only if h(1/z) belongs to Hp ), then (PN ) reduces to rational approximation. Moreover, it turns out that the associated solution gNRN has no pole outside D , hence it is a stable rational approximant to f . However, in contrast with the situation when p = , this approximant may not be unique.

The former Miaou project (predecessor of Apics) has designed an adapted steepest-descent algorithm for the case p = 2 whose convergence to a local minimum is guaranteed; it seems today the only procedure meeting this property. Roughly speaking, it is a gradient algorithm that proceeds recursively with respect to the order N of the approximant, in a compact region of the parameter space  [40] . Although it has proved rather effective in all applications carried out so far (see sections 4.2 , 4.3 ), it is not known whether the absolute minimum can always be obtained by choosing initial conditions corresponding to critical points of lower degree (as done by the Endymion software section  5.5 and RARL2 software, section 5.2 ).

In order to establish convergence results of the algorithm to the global minimum, Apics has undergone a long-haul study of the number and nature of critical points, in which tools from differential topology and operator theory team up with classical approximation theory. The main discovery is that the nature of the critical points (e.g. local minima , saddles...) depends on the decrease of the interpolation error to f as N increases [48] . Based on this, sufficient conditions have been developed for a local minimum to be unique. This technique requires strong error estimates that are often difficult to obtain, and most of the time only hold for N large. Examples where uniqueness or asymptotic uniqueness has been proved this way include transfer functions of relaxation systems (i.e., Markov functions) [49] , the exponential function, and meromorphic functions [8] . The case where f is the Cauchy integral on a hyperbolic geodesic arc of a Dini-continuous function which does not vanish “too much” has been recently answered in the positive, see section  6.7 . An analog to AAK theory has been carried out for 2p< [9] . Although not computationally as powerful, it has better continuity properties and stresses a continuous link between rational approximation in H2 and meromorphic approximation in the uniform norm, allowing one to use, in either context, techniques available from the other(When 1p<2 , problem (PN ) is still fairly open.).

A common feature to all these problems is that critical point equations express non-Hermitian orthogonality relations for the denominator of the approximant. This is used in an essential manner to assess the behavior of the poles of the approximants to functions with branched singularities which is of particular interest for inverse source problems (cf.  sections 6.3.2 , 6.7 ).

In higher dimensions, the analog of problem (PN ) is the approximation of a vector field with gradients of potentials generated by N point masses instead of meromorphic functions. The issue is by no means understood at present, and is a major endeavor of future research problems.

Certain constrained rational approximation problems, of special interest in identification and design of passive systems, arise when putting additional requirements on the approximant, for instance that it should be smaller than 1 in modulus. Such questions have become over years an increasingly significant part of the team's activity (see sections  4.3 , 6.6 , 6.7 , and 6.10 ). When translated over to the circle, a prototypical formulation consists in approximating the modulus of a given function by the modulus of a rational function of degree n . When p = 2 this problem can be reduced to a series of standard rational approximation problems, but usually one needs to solve it for p = . The case where |f| is a piecewise constant function with values 0 and 1 can also be approached via classical Zolotarev problems [86] , that can be solved more or less explicitly when the pass-band consists of a single arc. A constructive solution in the case where |f| is a piecewise constant function with values 0 and 1 on several arcs (multiband filters) is one recent achievement of the team. Though the modulus of the response is the first concern in filter design, the variation of the phase must nevertheless remain under control to avoid unacceptable distortion of the signal. This is an important issue, currently under investigation within the team under contract with the CNES, see section  6.10 .

From the point of view of design, rational approximants are indeed useful only if they can be translated into physical parameter values for the device to be built. This is where system theory enters the scene, as the correspondence between the frequency response (i.e., the transfer-function) and the linear differential equations that generate this response (i.e., the state-space representation), which is the object of the so-called realization process. Since filters have to be considered as dual modes cavities, the realization issue must indeed be tackled in a 2×2 matrix-valued context that adds to the complexity. A fair share of the team's research in this direction is concerned with finding realizations meeting certain constraints (imposed by the technology in use) for a transfer-function that was obtained with the above-described techniques (see section  6.8 ).

#### Behavior of poles of meromorphic approximants and inverse problems for the Laplacian

Keywords : singularity detection, free boundary inverse problems, meromorphic and rational approximation, orthogonal polynomials, discretization of potentials.

Participants : Laurent Baratchart, Edward Saff, Herbert Stahl [ TFH Berlin ] , Maxim Yattselev.

We refer here to the behavior of the poles of best meromorphic approximants, in the Lp -sense on a closed curve, to functions f defined as Cauchy integrals of complex measures whose support lies inside the curve. If one normalizes the contour to be the unit circle T , we are back to the framework of section 3.1.2 and to problem (PN ); the invariance of the problem under conformal mapping was established in [6] . The research so far has focused on functions whose singular set inside the contour is zero or one-dimensional.

Generally speaking, the behavior of poles is particularly important in meromorphic approximation to obtain error rates as the degree goes large and also to tackle constructive issues like uniqueness. However, the original motivation of Apics is to consider this issue in connection with the approximation of the solution to a Dirichlet-Neumann problem, so as to extract information on the singularities. The general theme is thus how do the singularities of the approximant reflect those of the approximated function? The approach to inverse problem for the 2-D Laplacian that we outline here is attractive when the singularities are zero- or one-dimensional (see section  4.2 ). It can be used as a computationally cheap preliminary step to obtain the initial guess of a more precise but heavier numerical optimization.

For sufficiently smooth cracks, or pointwise sources recovery, the approach in question is in fact equivalent to the meromorphic approximation of a function with branch points, and we were able to prove [4] , [6] that the poles of the approximants accumulate in a neighborhood of the geodesic hyperbolic arc that links the endpoints of the crack, or the sources [44] . Moreover the asymptotic density of the poles turns out to be the equilibrium distribution on the geodesic arc of the Green potential and it charges the end points, that are thus well localized if one is able to compute sufficiently many zeros (this is where the method could fail). The case of more general cracks, as well as situations with three or more sources, requires the analysis of the situation where the number of branch points is finite but arbitrary, see section  6.7 ). This are outstanding open questions for applications to inverse problems (see section  6.3 ), as also the problem of a general singularity, that may be two dimensional.

Results of this type open new perspectives in non-destructive control, in that they link issues of current interest in approximation theory (the behavior of zeroes of non-Hermitian orthogonal polynomials) to some classical inverse problems for which a dual approach is thereby proposed: to approximate the boundary conditions by true solutions of the equations, rather than the equation itself (by discretization).

Let us point out that the problem of approximating, by a rational or meromorphic function, in the Lp sense on the boundary of a domain, the Cauchy transform of a real measure, localized inside the domain, can be viewed as an optimal discretization problem for a logarithmic potential according to a criterion involving a Sobolev norm. This formulation can be generalized to higher dimensions, even if the computational power of complex analysis is then no longer available, and this makes for a long-term research project with a wide range of applications. It is interesting to mention that the case of sources in dimension three in a spherical or ellipsoidal geometry, can be attacked with the above 2-D techniques as applied to planar sections (see section  6.3 ).

#### Matrix-valued rational approximation

Keywords : rational approximation, inner matrix, reproducing kernel space, realization theory.

Participants : Laurent Baratchart, Martine Olivi, José Grimm, Jean-Paul Marmorat, Bernard Hanzon, Ralf Peeters [ Univ. Maastricht ] .

Matrix-valued approximation is necessary for handling systems with several inputs and outputs, and it generates substantial additional difficulties with respect to scalar approximation, theoretically as well as algorithmically. In the matrix case, the McMillan degree (i.e., the degree of a minimal realization in the System-Theoretic sense) generalizes the degree.

The problem we want to consider reads: Let and n an integer; find a rational matrix of size m×l without poles in the unit disk and of McMillan degree at most n which is nearest possible to in (H2)m×l . Here the L2 norm of a matrix is the square root of the sum of the squares of the norms of its entries.

The approximation algorithm designed in the scalar case generalizes to the matrix-valued situation  [67] . The first difficulty consists here in the parametrization of transfer matrices of given McMillan degree n , and the inner matrices (i.e., matrix-valued functions that are analytic in the unit disk and unitary on the circle) of degree n enter the picture in an essential manner: they play the role of the denominator in a fractional representation of transfer matrices (using the so-called Douglas-Shapiro-Shields factorization).

The set of inner matrices of given degree has the structure of a smooth manifold that allows one to use differential tools as in the scalar case. In practice, one has to produce an atlas of charts (parametrization valid in a neighborhood of a point), and we must handle changes of charts in the course of the algorithm. Such parametrization can be obtained from interpolation theory and Schur type algorithms, the parameters being interpolation vectors or matrices [33] , [10] , [11] . Some of these parametrizations have a particular interest for computation of realizations [10] , [11] , involved in the estimation of physical quantities for the synthesis of resonant filters. Two rational approximation codes (see sections 5.2 and 5.5 ) have been developed in the team.

Problems relative to multiple local minima naturally arise in the matrix-valued case as well, but deriving criteria that guarantee uniqueness is even more difficult than in the scalar case. The already investigated case of rational functions of the sought degree (the consistency problem) was solved using rather heavy machinery  [7] , and that of matrix-valued Markov functions, that are the first example beyond rational function has made progress only recently [39] .

Let us stress that the algorithms mentioned above are first to handle rational approximation in the matrix case in a way that converges to local minima, while meeting stability constraints on the approximant.

Logo Inria