## Section: Research Program

### Approximation

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

#### Best 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}$, to recover a holomorphic function from its values on a subset $K$ of the boundary of $D$. For the discussion it is convenient to normalize $D$, which can be done by conformal mapping. So, in the simply connected case, we fix $D$ to be the unit disk with boundary unit circle $T$. We denote by ${H}^{p}$ the Hardy space of exponent $p$, which is the closure of polynomials in ${L}^{p}\left(T\right)$-norm 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 $g$ in $D$ matching some measured values $f$ approximately 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
thereof. The value of $p$ reflects the type of
stability which is sought and how much one wants to smooth out the data.
The choice of ${L}^{p}$ classes is suited to handle pointwise
measurements.

To fix terminology, we refer to $\left(P\right)$ as
a *bounded extremal problem*.
As shown in [39], [42],
[48],
the solution to this convex
infinite-dimensional optimization problem can be obtained
when $p\ne 1$ upon iterating with respect to a Lagrange parameter
the solution to spectral equations for
appropriate Hankel and Toeplitz operators.
These spectral equations involve the solution to the special case
$K=T$ of $\left(P\right)$, which is a standard extremal problem
[63]:

(${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$ is more or less open.

Various modifications of $\left(P\right)$ can be tailored to meet specific needs. For instance when dealing with lossless transfer functions (see Section 4.4), one may want to express the constraint on $T\setminus K$ in a pointwise manner: $|g-\psi |\le M$ a.e. on $T\setminus K$, see [43]. In this form, the problem comes close to (but still is different from) ${H}^{\infty}$ frequency optimization used in control [66], [74]. One can also impose bounds on the real or imaginary part of $g-\psi $ on $T\setminus K$, which is useful when considering Dirichlet-Neumann problems, see [68].

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 to recover a harmonic function on
the inner boundary from Dirichlet-Neumann data on the
outer boundary (see Sections 3.2.1, 4.3, 5.1.5). 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 solution.
In this case, the Lagrange parameter indicates
how to deform the inner contour in order to improve
data fitting.
Similar topics are discussed in Section 3.2.1 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 no longer constant.
Then, the Hardy spaces in Problem $\left(P\right)$
are those of a so-called conjugate Beltrami equation:
$\overline{\partial}f=\nu \overline{\partial f}$ [67],
which are studied for
$1<p<\infty $ in
[4], [12], [34] and [58].
Expansions
of solutions needed to constructively handle such issues in the specific
case of linear fractional conductivities (occurring for instance
in plasma shaping)
have been expounded in [60].

Another instance of problem $\left(P\right)$ in with $p=2$ and additional pointwise interpolation constraints inside a simply connected domain (disk) $D$ was studied and solved in [11], Part I, and [15]. Such pointwise interpolation constraints could be of practical interest for inverse Cauchy type problems in cases where interior information is also available or to model uncertainty on boundary data.

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 some ${\mathbb{R}}^{n}$-valued vector field $V$ 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 substitute for holomorphic Hardy spaces is provided by the Stein-Weiss Hardy spaces of harmonic gradients [77]. Conformal maps are no longer available when $n>2$, so that $\Omega $ can no longer be normalized. More general geometries than spheres and half-spaces 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. Given $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$,
Problem $\left({P}_{1}\right)$ was solved in [2] as well as
its analog on a shell, when the tangent component of $V$ is a gradient
(when $O$ is Lipschitz the general case follows easily from this).
The solution extends
the work in [39] to the 3-D case,
using a generalization of Toeplitz operators. The case of the shell was motivated by applications to the processing of EEG data.
An important ingredient is a refinement of the Hodge
decomposition, that we call the *Hardy-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}^{p}\left(B\right)$, a vector field in
${H}^{p}({\mathbb{R}}^{n}\setminus \overline{B})$,
and a tangential divergence free vector field on $S$; the space of such
divergence-free fields
is denoted by $D\left(S\right)$.
If $p=1$ or $p=\infty $,
${L}^{p}$ must be replaced by the real Hardy space
or the space of functions with bounded mean oscillation.
More generally this decomposition, which is valid on any sufficiently
smooth surface (see
Section 5.1), seems to play a fundamental role
in inverse potential problems. In fact, it was first introduced
formally on the plane to describe
silent magnetizations supported in ${\mathbb{R}}^{2}$
(*i.e.* those generating no field
in the upper half space) [36].

Just like solving problem $\left(P\right)$ appeals to the solution of problem $\left({P}_{0}\right)$, our ability to solve problem $\left({P}_{1}\right)$ will depend on the possibility to tackle the special case where $O=S$:

$\left({P}_{2}\right)$ Let $1\le p\le \infty $ and $V\in {L}^{p}\left(S\right)$ be a ${\mathbb{R}}^{n}$-valued vector field. Find a harmonic gradient $G\in {H}^{p}\left(B\right)$ such that ${\parallel G-V\parallel}_{{L}^{p}\left(S\right)}$ is minimum.

Problem $\left({P}_{2}\right)$ is simple when $p=2$ by virtue of the Hardy Hodge decomposition together with orthogonality of ${H}^{2}\left(B\right)$ and ${H}^{2}({\mathbb{R}}^{n}\setminus \overline{B})$, which is the reason why we were able to solve $\left({P}_{1}\right)$ in this case. Other values of $p$ cannot be treated as easily and are still under investigation, especially the case $p=\infty $ which is of particular interest and presents itself as a 3-D analog to the Nehari problem [73].

Companion to problem $\left({P}_{2}\right)$ is problem $\left({P}_{3}\right)$ below.

$\left({P}_{3}\right)$ Let $1\le p\le \infty $ and $V\in {L}^{p}\left(S\right)$ be a ${\mathbb{R}}^{n}$-valued vector field. Find $G\in {H}^{p}\left(B\right)$ and $D\in D\left(S\right)$ such that ${\parallel G+D-V\parallel}_{{L}^{p}\left(S\right)}$ is minimum.

Note that $\left({P}_{2}\right)$ and $\left({P}_{3}\right)$ are identical in 2-D, since no non-constant tangential divergence-free vector field exists on $T$. It is no longer so in higher dimension, where both $\left({P}_{2}\right)$ and $\left({P}_{3}\right)$ arise in connection with inverse potential problems in divergence form, like source recovery in electro/magneto encephalography and paleomagnetism, see Sections 3.2.1 and 4.3.

#### Best meromorphic and rational approximation

The techniques set forth in this section are used to solve step 2 in Section 3.2 and they are instrumental to approach inverse boundary value problems for the Poisson equation $\Delta u=\mu $, where $\mu $ is some (unknown) measure.

##### Scalar meromorphic and rational approximation

We put ${R}_{N}$ for the set of rational functions with at most $N$ poles in $D$. By definition, meromorphic functions in ${L}^{p}\left(T\right)$ are (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 $f$ continuous is it known how to solve $\left({P}_{N}\right)$ in semi-closed form. The unique solution is given by AAK theory (named after Adamjan, Arov and Krein), which connects the spectral decomposition of Hankel operators with best approximation [73].

The case where $p=2$ is of special importance for it
reduces to rational approximation. Indeed,
if we write the Hardy decomposition $f={f}^{+}+{f}^{-}$ where ${f}^{+}\in {H}^{2}$ and
${f}^{-}\in {H}^{2}(\u2102\setminus \overline{D})$,
then ${g}_{N}={f}^{+}+{r}_{N}$ where ${r}_{N}$ is a best approximant to ${f}^{-}$ from ${R}_{N}$
in ${L}^{2}\left(T\right)$.
Moreover, ${r}_{N}$ has no pole outside $D$,
hence it is a *stable* rational
approximant to ${f}^{-}$. However, in contrast to the case
where $p=\infty $, this best approximant may *not* be unique.

The former Miaou project (predecessor of Apics) designed a dedicated
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. This gradient algorithm proceeds
recursively with respect to $N$ on a compactification
of the parameter space [32].
Although it has proved to be
effective in all applications carried out so far
(see Sections 4.3, 4.4),
it is still unknown 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 3.4.4).

In order to establish global convergence results, Apics has undertaken a
deeper study of the number and nature of critical points
(local minima, saddle points...), in which
tools from differential topology and
operator theory team up with classical interpolation theory
[45], [47].
Based on this work,
uniqueness or asymptotic uniqueness of the approximant
was proved for certain classes of functions like
transfer functions of relaxation
systems (*i.e.*
Markov functions) [49] and more
generally Cauchy integrals over hyperbolic geodesic arcs [52].
These are the only results of this kind. Research by Apics on this topic
remained dormant for a while by reasons of opportunity,
but revisiting the work [30] in higher dimension is
a worthy and timely endeavor today. Meanwhile,
an analog to AAK theory
was carried out for $2\le p<\infty $ in [48].
Although not as effective
computationally, it was recently used
to derive lower bounds [3].
When $1\le p<2$, problem $\left({P}_{N}\right)$ is still quite open.

A common
feature to the above-mentioned problems
is that critical point equations
yield non-Hermitian orthogonality relations for the denominator
of the approximant. This stresses connections with interpolation,
which is a standard way to build approximants,
and in many respects best or near-best rational approximation
may be regarded as a clever manner to pick interpolation points.
This was exploited in [53], [50],
and is used in an essential manner to assess the
behavior of poles of best approximants to functions with branched
singularities,
which is of particular interest for inverse source problems
(*cf.* Sections 3.4.2
and 5.1).

In higher dimensions, the analog of Problem $\left({P}_{N}\right)$ is best approximation of a vector field by gradients of discrete potentials generated by $N$ point masses. This basic issue is by no means fully understood, and it is an exciting field of research. It is connected with certain generalizations of Toeplitz or Hankel operators, and with constructive approaches to so-called weak factorizations for real Hardy functions [59].

Besides,
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 (*i.e.* a Schur function). In particular, Schur interpolation
lately received renewed attention
from the team, in connection with matching problems.
There, interpolation data are subject to
a well-known compatibility condition (positive definiteness of the so-called
Pick matrix), and the main difficulty is to put interpolation
points on the boundary of $D$ while controlling both the degree and the
extremal points (peak points for the modulus) of the interpolant.
Results obtained by Apics in this direction generalize
a variant of contractive interpolation
with degree constraint as studied in [64],
see Section 5.2.
We mention that contractive interpolation with nodes approaching the boundary
has been a subsidiary research topic by the team in the past,
which plays an interesting role in the
spectral representation of certain non-stationary
stochastic processes [35], [38]. The subject is
intimately connected to orthogonal polynomials on the unit circle,
and this line of investigation has been pursued
towards an asymptotic study of orthogonal polynomials on planar domains,
which is today an active area in approximation theory with application to
quantum particle systems, spectra of random matrices, and Hele-Shaw flows, see
Section 5.6.

##### Matrix-valued rational approximation

Matrix-valued approximation is necessary to handle systems with several
inputs and outputs but it generates additional difficulties
as compared to scalar-valued approximation,
both theoretically and algorithmically. In the matrix case,
the McMillan degree (*i.e.* the degree of a minimal realization in
the System-Theoretic sense) generalizes the usual notion of degree
for rational functions. For instance when poles are simple, the McMillan degree is the sum of the ranks of the residues.

The basic problem that we consider now goes as follows:
*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 derived in
[32] and mentioned in
Section 3.3.2.1
generalizes to
the matrix-valued situation [62]. The
first difficulty here is to parametrize
inner matrices (*i.e.* matrix-valued functions
analytic in the unit disk and unitary on the unit circle) of
given McMillan degree degree $n$.
Indeed, inner matrices play the role of denominators
in fractional representations of transfer matrices (using the so-called
Douglas-Shapiro-Shields factorization).
The set of inner matrices of given degree is
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
(local parametrizations) and to handle changes of charts in the course of the algorithm. Such parametrization can be obtained using
interpolation theory and Schur-type algorithms, the parameters of which
are vectors or matrices
( [28], [65], [71]). Some of
these parametrizations are also interesting to compute
realizations and achieve filter synthesis
( [65], [71]). The
rational approximation software “RARL2” developed
by the team is described in Section 3.4.4.

Difficulties relative to multiple local minima of course 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 degree $n$ or small perturbations thereof (the consistency problem) was solved in [46]. Matrix-valued Markov functions are the only known example beyond this one [31].

Let us stress that RARL2 seems the only algorithm handling rational approximation in the matrix case that demonstrably converges to a local minimum while meeting stability constraints on the approximant. It is still a working pin of many developments by Apics on frequency optimization and design.

#### Behavior of poles of meromorphic approximants

Participant : Laurent Baratchart.

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. Normalizing the contour to be the unit circle $T$, we are back to Problem $\left({P}_{N}\right)$ in Section 3.3.2.1; invariance of the latter under conformal mapping was established in [5]. Research so far has focused on functions whose singular set inside the contour is polar, meaning that the function can be continued analytically (possibly in a multiple-valued manner) except over a set of logarithmic capacity zero.

Generally speaking in approximation theory, assessing the
behavior of poles of rational approximants is essential
to obtain error rates as the degree goes large, and to tackle
constructive issues like
uniqueness. However, as explained in Section 3.2.1,
the original twist by Apics is to consider this issue also as a means
to extract information on
singularities of the solution to a
Dirichlet-Neumann problem.
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.3). It can be used
as a computationally cheap
initial condition for more precise but much heavier
numerical optimizations which often do not even converge
unless properly initialized.
As regards crack detection or source recovery, this approach
boils down to
analyzing the behavior of best meromorphic
approximants of given pole cardinality to a function with branch points, which is the prototype of
a polar singular set.
For piecewise analytic cracks, or in the case of sources, we were able to
prove ([5], [6], [37]),
that the poles of the
approximants accumulate, when the degree goes large,
to some extremal cut of minimum weighted
logarithmic capacity connecting
the singular points of the crack, or the sources
[40].
Moreover, the asymptotic density
of the poles turns out to be the Green equilibrium distribution
on this cut in $D$, therefore it charges the
singular points if one is able to approximate in
sufficiently high degree (this is where the method could fail, because
high-order approximation requires rather precise data).

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

It is remarkable that inverse source problems inside a sphere or an ellipsoid in 3-D can be approached with such 2-D techniques, as applied to planar sections, see Section 5.1. The technique is implemented in the software FindSources3D, see Section 3.4.2.