Section: Research Program

Approximation of boundary data

Participants : Laurent Baratchart, Sylvain Chevillard, Juliette Leblond, Martine Olivi, Dmitry Ponomarev, Elodie Pozzi, Fabien Seyfert.

The following people are collaborating with us on these topics: Bernard Hanzon (Univ. Cork, Ireland), Jean-Paul Marmorat (Centre de mathématiques appliquées (CMA), École des Mines de Paris), Jonathan Partington (Univ. Leeds, UK), Ralf Peeters (Univ. Maastricht, NL), Edward Saff (Vanderbilt University, Nashville, USA), Herbert Stahl (TFH Berlin), Maxim Yattselev (Purdue Univ. at Indianapolis, USA).

Best constrained analytic approximation

In dimension 2, the prototypical problem to be solved in step 1 of Section  3.1 may be described as: given a domain $D\subset {ℝ}^2$, we want to recover a holomorphic function from its values on a subset of the boundary of $D$. Using conformal mapping, it is convenient for the discussion to normalize $D$. So, in the simply connected case, we fix $D$ to be the unit disk with boundary the unit circle $T$. We denote by $H^p$ the Hardy space of exponent $p$ which is the closure of polynomials in the $L^p$-norm on the circle if $1\le p<\infty$ and the space of bounded holomorphic functions in $D$ if $p=\infty$. Functions in $H^p$ have well-defined boundary values in $L^p\left(T\right)$, which makes it possible to speak of (traces of) analytic functions on the boundary.

To find an analytic function in $D$ approximately matching measured values $f$ on a sub-arc $K$ of $T$, we formulate a constrained best approximation problem as follows.

$\left(P\right)$  Let $1\le p\le \infty$, $K$ a sub-arc of $T$, $f\in L^p\left(K\right)$, $\psi \in L^p(T\setminus K)$ and $M>0$; find a function $g\in H^p$ such that ${\parallel g-\psi \parallel }_{L^p(T\setminus K)}\le M$ and $g-f$ is of minimal norm in $L^p\left(K\right)$ under this constraint.

Here $\psi$ is a reference behavior capturing a priori assumptions on the behavior of the model off $K$, while $M$ is some admissible deviation from them. The value of $p$ reflects the type of stability which is sought and how much one wants to smoothen the data. The choice of $L^p$ classes is well-adapted to handling point-wise measurements.

To fix terminology we refer to $\left(P\right)$ as a bounded extremal problem. As shown in [43] , [45] , [51] , for $1<p\le \infty$, the solution to this convex infinite-dimensional optimization problem can be obtained upon iterating with respect to a Lagrange parameter the solution to spectral equations for some appropriate Hankel and Toeplitz operators. These equations in turn involve the solution to the standard extremal problem below [64] :

($P_0$)  Let $1\le p\le \infty$ and $\varphi \in L^p\left(T\right)$; find a function $g\in H^p$ such that $g-\varphi$ is of minimal norm in $L^p\left(T\right)$.

The case $p=1$ of $\left(P_0\right)$ is essentially open.

Various modifications of $\left(P\right)$ have been studied in order to meet specific needs. For instance when dealing with loss-less transfer functions (see Section  4.5 ), one may want to express the constraint on $T\setminus K$ in a point-wise manner: $|g-\psi |\le M$ a.e. on $T\setminus K$, see [47] . In this form, it comes close to (but still is different from) $H^{\infty }$ frequency optimization methods for control [66] , [73] .

The analog of Problem $\left(P\right)$ on an annulus, $K$ being now the outer boundary, can be seen as a means to regularize a classical inverse problem occurring in nondestructive control, namely recovering a harmonic function on the inner boundary from Dirichlet-Neumann data on the outer boundary (see Sections  3.2.1 , 4.2 , 6.1.1 , 6.2 ). It may serve as a tool to approach Bernoulli type problems where we are given data on the outer boundary and we seek the inner boundary, knowing it is a level curve of the flux. Then, the Lagrange parameter indicates which deformation should be applied on the inner contour in order to improve data fitting.

This is discussed in Sections  3.2.1 and  6.2 for more general equations than the Laplacian, namely isotropic conductivity equations of the form $\mathrm{div}\left(\sigma \nabla u\right)=0$ where $\sigma$ is non constant. In this case, the Hardy spaces in Problem $\left(P\right)$ are those of a so-called conjugate or real Beltrami equation $\overline{\partial }f=\nu \overline{\partial f}$ [67] , which were studied for $1<p<\infty$ in [13] , [4] . Expansions of solutions needed to constructively handle such issues have been carried out in  [61] .

Though originally considered in dimension 2, Problem $\left(P\right)$ carries over naturally to higher dimensions where analytic functions get replaced by gradients of harmonic functions. Namely, given some open set $\Omega \subset {ℝ}^n$ and a ${ℝ}^n$-valued vector $V$ field on an open subset $O$ of the boundary of $\Omega$, we seek a harmonic function in $\Omega$ whose gradient is close to $V$ on $O$.

When $\Omega$ is a ball or a half-space, a convenient substitute of holomorphic Hardy spaces is provided by Stein-Weiss Hardy spaces of harmonic gradients [77] . Conformal maps are no longer available in ${ℝ}^n$ for $n>2$ and other geometries have not been much studied so far. On the ball, the analog of Problem $\left(P\right)$ is

$\left(P_1\right)$  Let $1\le p\le \infty$ and $B\subset {ℝ}^n$ the unit ball. Fix $O$ an open subset of the unit sphere $S\subset {ℝ}^n$. Let further $V\in L^p\left(O\right)$ and $W\in L^p(S\setminus O)$ be ${ℝ}^n$-valued vector fields, and $M>0$; find a harmonic gradient $G\in H^p\left(B\right)$ such that ${\parallel G-W\parallel }_{L^p(S\setminus O)}\le M$ and $G-V$ is of minimal norm in $L^p\left(O\right)$ under this constraint.

When $p=2$, spherical harmonics offer a reasonable substitute to Fourier expansions and Problem $\left(P_1\right)$ was solved in [2] , together with its natural analog on a shell. The solution generalizes the Toeplitz operator approach to bounded extremal problems [43] , and constructive aspects of the procedure (harmonic 3-D projection, Kelvin and Riesz transformation, spherical harmonics) were derived. An important ingredient is a refinement of the Hodge decomposition allowing us to express a ${ℝ}^n$-valued vector field in $L^p\left(S\right)$, $1<p<\infty$, as the sum of a vector field in $H\left(B\right)$, a vector field in $H^p({ℝ}^n\setminus \overline{B})$, and a tangential divergence free vector field. If $p=1$ or $p=\infty$, $L^p$ must be replaced respectively by the real Hardy space $H^1$ and the bounded mean oscillation space $BMO$, and $H^{\infty }$ should be modified accordingly. This decomposition was fully discussed in [14] (for the case of the half-space) where it plays a fundamental role.

Problem $\left(P_1\right)$ is under investigation in the case $p=\infty$, where even the case where $O=S$ is pending because a substitute of the Adamjan-Arov-Krein theory [71] is still to be built in dimension greater than 2.

Such problems arise in connection with source recovery in electro/magneto encephalography and paleomagnetism, as discussed in Sections  3.2.1 and  4.2 .

Best meromorphic and rational approximation

The techniques explained in this section are used to solve step 2 in Section  3.2 via conformal mapping and subsequently are instrumental to approach inverse boundary value problems for Poisson equation $\Delta u=\mu$, where $\mu$ is some (unknown) distribution.

Scalar meromorphic and rational approximation

Let as before $D$ designate the unit disk, and $T$ the unit circle. We further put $R_N$ for the set of rational functions with at most $N$ poles in $D$, which allows us to define meromorphic functions in $L^p\left(T\right)$ as traces of functions in $H^p+R_N$.

A natural generalization of Problem $\left(P_0\right)$ is:

($P_N$)  Let $1\le p\le \infty$, $N\ge 0$ an integer, and $f\in L^p\left(T\right)$; find a function $g_N\in H^p+R_N$ such that $g_N-f$ is of minimal norm in $L^p\left(T\right)$.

Only for $p=\infty$ and continuous $f$ it is known how to solve $\left(P_N\right)$ in closed form. The unique solution is given by AAK theory (named after Adamjan, Arov and Krein), that connects the spectral decomposition of Hankel operators with best approximation in Hankel norm  [71] . This theory allows one to express $g_N$ in terms of the singular vectors of the Hankel operator with symbol $f$. The continuity of $g_N$ as a function of $f$ only holds for norms finer than uniform.

The case $p=2$ is of special importance. In particular when $f\in {\overline{H}}^2$, the Hardy space of exponent 2 of the complement of $D$ in the complex plane (by definition, $h\left(z\right)$ belongs to ${\overline{H}}^p$ if, and only if $h(1/z)$ belongs to $H^p$), then $\left(P_N\right)$ reduces to rational approximation. Moreover, it turns out that the associated solution $g_N\in R_N$ has no pole outside $D$, hence it is a stable rational approximant to $f$. However, in contrast with the situation when $p=\infty$, 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; until now it seems to be 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 [38] . Although it has proved effective in all applications carried out so far (see Sections  4.2 , 4.5 ), it is not known whether the absolute minimum can always be obtained by choosing initial conditions corresponding to critical points of lower degree (as is done by the RARL2 software, Section  5.1 ).

In order to establish global convergence results, APICS has undertaken a deeper 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, saddle points...) 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. These conditions are hard to use in practice because they require strong estimates of the approximation error. These are often difficult to obtain for a given function, and are usually only valid for large $N$. Examples where uniqueness or asymptotic uniqueness has been proved this way include transfer functions of relaxation systems (i.e. Markov functions) [52] and more generally Cauchy integrals over hyperbolic geodesic arcs  [54] and certain entire functions  [50] .

An analog to AAK theory has been carried out for $2\le p<\infty$ [51] . Although not computationally as powerful, it can be used to derive lower bounds [29] and to analyze the behavior of poles. When $1\le p<2$, Problem $\left(P_N\right)$ 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 makes connection with interpolation theory [55] , [53] and is used in an essential manner to assess the behavior of the poles of the approximants to functions with branchpoint-type singularities, which is of particular interest for inverse source problems (cf. Sections  5.6 and  6.1 ).

In higher dimensions, the analog of Problem $\left(P_N\right)$ is best approximation of a vector field with gradients of potentials generated by $N$ point masses instead of meromorphic functions. This issue is by no means fully understood, and is an exciting line of research. It is connected with spectral properties of certain operators generalizing classical Toeplitz and Hankel ones, and to constructive approaches to so-called weak factorizations of div-curl type for real Hardy functions.

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 attracted significant attention of members of the team (see Section  4.5 ). For instance, convergence properties of multi-point Schur approximants, which are rational interpolants preserving contractivity of a function, were analyzed in [41] . Such approximants are useful in prediction theory of stochastic processes, but since they interpolate inside the domain of holomorphy they are of limited use in frequency design.

In another connection, the generalization to several arcs of classical Zolotarev problems [72] is an achievement by the team which is useful for multi-band synthesis [10] . Still, though the modulus of the response is the first concern in filter design, variation of the phase must nevertheless remain under control to avoid unacceptable distortion of the signal. This specific but important issue has less structure and was approached using constrained optimization; a dedicated code has been developed under contract with the CNES (see Section  5.5 ).

Matrix-valued rational approximation

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 consider is now: let $ℱ\in {\left(H^2\right)}^{m\times l}$ and $n$ an integer; find a rational matrix of size $m\times l$ without poles in the unit disk and of McMillan degree at most $n$ which is nearest possible to $ℱ$ in ${\left(H^2\right)}^{m\times l}$. Here the $L^2$ norm of a matrix is the square root of the sum of the squares of the norms of its entries.

The scalar approximation algorithm [38] , mentioned in Section  3.3.2.1 , generalizes to the matrix-valued situation [63] . The first difficulty here consists 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$. The latter enter the picture in an essential manner as 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 ( [35] , [9] , [11] ). Some of them are particularly interesting to compute realizations and achieve filter synthesis ([9] [11] ). For rational approximation software codes developed by the team, see Section  5.1 .

Difficulties relative to multiple local minima naturally arise in the matrix-valued case as well, and deriving criteria that guarantee uniqueness is even more difficult than in the scalar case. The case of rational functions of sought degree or small perturbations thereof (the consistency problem) was solved in  [49] . The case of matrix-valued Markov functions, the first example beyond rational functions, was treated in [37] .

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.

Behavior of poles of meromorphic approximants

Participant : Laurent Baratchart.

The following persons did collaborate with us on this subject: Herbert Stahl (TFH Berlin), Maxim Yattselev (Purdue Univ. at Indianapolis, USA).

We refer here to the behavior of poles of best meromorphic approximants, in the $L^p$-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.3.2.1 and to Problem $\left(P_N\right)$; invariance of the problem under conformal mapping was established in [5] . 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 to tackle constructive issues like uniqueness. As explained in Section  3.2.1 , we consider this issue in connection with 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? This approach to inverse problem for the 2-D Laplacian turns out to be attractive when singularities are zero- or one-dimensional (see Section  4.2 ). It can be used as a computationally cheap initialization of more precise but heavier numerical optimizations.

As regards crack detection or source recovery, the approach in question boils down to analyzing the behavior of best meromorphic approximants of a function with branch points. For piecewise analytic cracks, or in the case of sources, we were able to prove ([5] , [6] , [40] ), that the poles of the approximants accumulate on some extremal contour of minimum weighted energy linkings the singular points of the crack, or the sources [44] . Moreover, the asymptotic density of the poles turns out to be the Green equilibrium distribution of this contour in $D$, hence puts heavy charge around the singular points (in particular at the endpoints) which are therefore well localized if one is able to approximate in sufficiently high degree (this is where the method could fail).

The case of two-dimensional singularities is still an outstanding open problem.

It is interesting that inverse source problems inside a sphere or an ellipsoid in 3-D can be attacked with the above 2-D techniques, as applied to planar sections (see Section  6.1 ). This is at work in the software FindSources3D, see Section  5.6 .

Miscellaneous

Participant : Sylvain Chevillard.

Sylvain Chevillard, joined team in November 2010. His coming resulted in APICS hosting a research activity in certified computing, centered on the software Sollya of which S. Chevillard is a co-author, see Section  5.7 . On the one hand, Sollya is an Inria software which still requires some tuning to a growing community of users. On the other hand, approximation-theoretic methods at work in Sollya are potentially useful for certified solutions to constrained analytic problems described in Section  3.3.1 . However, developing Sollya is not a long-term objective of APICS.