## 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 {\mathbb{R}}^{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 {\mathbb{R}}^{n}$ and a ${\mathbb{R}}^{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 ${\mathbb{R}}^{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 {\mathbb{R}}^{n}$ the unit ball. Fix $O$ an open subset of the unit sphere $S\subset {\mathbb{R}}^{n}$. Let further $V\in {L}^{p}\left(O\right)$ and $W\in {L}^{p}(S\setminus O)$ be ${\mathbb{R}}^{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 ${\mathbb{R}}^{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}({\mathbb{R}}^{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 $\mathcal{F}\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 $\mathcal{F}$ 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.