Section: Research Program

Approximation

Participants : Laurent Baratchart, Sylvain Chevillard, Juliette Leblond, Martine Olivi, 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 {ℝ}^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  [40], [43],  [49], 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  [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)$.

In the case $p=1$, partial results are known but computational issues remain 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  [44]. In this form, the problem comes close to (but still is different from) $H^{\infty }$ frequency optimization used in control  [67], [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.

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, 6.1.3). 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 (i.e., varies in the space). 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}$  [68], which are studied for $1<p<\infty$ in [6],  [31],  [35] and  [59]. 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  [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 some ${ℝ}^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  [78]. 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 {ℝ}^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. 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  [40] 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 ${ℝ}^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({ℝ}^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 6.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 ${ℝ}^2$ (i.e. those generating no field in the upper half space)  [37].

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 ${ℝ}^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({ℝ}^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 ${ℝ}^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(ℂ\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 Miaou project (predecessor of Apics) already designed a dedicated steepest-descent algorithm for the case $p=2$ whose convergence to a local minimum is guaranteed; the algorithm ha evolved over years and still 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  [33]. 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.5).

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  [46], [48]. 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)  [50] and more generally Cauchy integrals over hyperbolic geodesic arcs  [53]. 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  [29] 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  [49]. Although not as effective computationally, it was recently used to derive lower bounds [5]. 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  [54], [51], 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.3 and 6.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  [60].

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  [65]. 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  [36], [39].

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 $ℱ\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 derived in  [33] and mentioned in Section 3.3.2.1 generalizes to the matrix-valued situation  [63]. 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 ( [27], [66], [71]). Some of these parametrizations are also interesting to compute realizations and achieve filter synthesis ( [66], [71]). The rational approximation software “RARL2” developed by the team is described in Section 3.4.5.

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  [47]. Matrix-valued Markov functions are the only known example beyond this one  [30].

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 Factas 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  [45]. 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, now Factas, 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 ([8],  [45], [38]), 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  [41]. 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 6.1. The technique is implemented in the software FindSources3D, see Section 3.4.3.