Geometry and Singularities of the Prony mapping

Dmitry Batenkov Department of Mathematics
Weizmann Institute of Science
Rehovot 76100
Israel dima.batenkov@weizmann.ac.il http://www.wisdom.weizmann.ac.il/ dmitryb and Yosef Yomdin yosef.yomdin@weizmann.ac.il http://www.wisdom.weizmann.ac.il/ yomdin

Abstract.

Prony mapping provides the global solution of the Prony system of equations

\Sigma_{i=1}^{n}A_{i}x_{i}^{k}=m_{k},\ k=0,1,\dots,2n-1.

This system appears in numerous theoretical and applied problems arising in Signal Reconstruction. The simplest example is the problem of reconstruction of linear combination of $\delta$ -functions of the form $g(x)=\sum_{i=1}^{n}a_{i}\delta(x-x_{i})$ , with the unknown parameters $a_{i},\ x_{i},\ i=1,\dots,n,$ from the “moment measurements” $m_{k}=\int x^{k}g(x)dx.$

Global solution of the Prony system, i.e. inversion of the Prony mapping, encounters several types of singularities. One of the most important ones is a collision of some of the points $x_{i}.$ The investigation of this type of singularities has been started in [21] where the role of finite differences was demonstrated.

In the present paper we study this and other types of singularities of the Prony mapping, and describe its global geometry. We show, in particular, close connections of the Prony mapping with the “Vieta mapping” expressing the coefficients of a polynomial through its roots, and with hyperbolic polynomials and “Vandermonde mapping” studied by V. Arnold.

Key words and phrases:

Singularities, Signal acquisition, Non-linear models, Moments inversion.

2000 Mathematics Subject Classification:

94A12 62J02, 14P10, 42C99

This research is supported by the Adams Fellowship Program of the Israel Academy of Sciences and Humanities, ISF grant 264/09 and the Minerva Foundation.

1. Introduction

Prony system appears as we try to solve a very simple “algebraic signal reconstruction” problem of the following form: assume that the signal $F(x)$ is known to be a linear combination of shifted $\delta$ -functions:

F\left(x\right)=\sum_{j=1}^{d}a_{j}\delta\left(x-x_{j}\right).

(1.1)

We shall use as measurements the polynomial moments:

m_{k}=m_{k}\left(F\right)=\int x^{k}F\left(x\right)\operatorname{d}x.

(1.2)

After substituting $F$ into the integral defining $m_{k}$ we get

m_{k}(F)=\int x^{k}\sum_{j=1}^{d}a_{j}\delta(x-x_{j})\operatorname{d}x=\sum_{j=1}^{d}a_{j}x_{j}^{k}.

Considering $a_{j}$ and $x_{j}$ as unknowns, we obtain equations

m_{k}\left(F\right)=\sum_{j=1}^{d}a_{j}x_{j}^{k},\;k=0,1,\dots.

(1.3)

This infinite set of equations (or its part, for $k=0,1,\dots,2d-1$ ), is called Prony system. It can be traced at least to R. de Prony (1795, [19]) and it is used in a wide variety of theoretical and applied fields. See [2] for an extensive bibligoraphy on the Prony method.

In writing Prony system (1.3) we have assumed that all the nodes $x_{1},\dots,x_{d}$ are pairwise different. However, as a right-hand side $\mu=(m_{0},\dots,m_{2d-1})$ of (1.3) is provided by the actual measurements of the signal $F$ , we cannot guarantee a priori, that this condition is satisfied for the solution. Moreover, we shall see below that multiple nodes may naturally appear in the solution process. In order to incorporate possible collisions of the nodes, we consider “confluent Prony systems”.

Assume that the signal $F(x)$ is a linear combination of shifted $\delta$ -functions and their derivatives:

F\left(x\right)=\sum_{j=1}^{s}\sum_{\ell=0}^{d_{j}-1}a_{j,\ell}\delta^{\left(\ell\right)}\left(x-x_{j}\right).

(1.4)

Definition 1.1.

For $F\left(x\right)$ as above, the vector $D\left(F\right)\stackrel{{\scriptstyle\text{def}}}{{=}}(d_{1},\dots,d_{s})$ is the multiplicity vector of $F$ , $s=s\left(F\right)$ is the size of its support, $T\left(F\right)\stackrel{{\scriptstyle\text{def}}}{{=}}\left(x_{1},\dots,x_{s}\right)$ , and $\operatorname{rank}\left(F\right)\stackrel{{\scriptstyle\text{def}}}{{=}}\sum_{j=1}^{s}d_{j}$ is its rank. For avoiding ambiguity in these definitions, it is always understood that $a_{j,d_{j}-1}\neq 0$ for all $j=1,\dots,s$ (i.e. $d_{j}$ is the maximal index for which $a_{j,d_{j}-1}\neq 0$ ).

For the moments $m_{k}=m_{k}(F)=\int x^{k}F(x)\operatorname{d}x$ we now get

m_{k}=\sum_{j=1}^{s}\sum_{\ell=0}^{d_{j}-1}a_{j,\ell}\frac{{k!}}{{(k-\ell)!}}x_{j}^{k-\ell}.

Considering $x_{i}$ and $a_{j,\ell}$ as unknowns, we obtain a system of equations

\sum_{j=1}^{s}\sum_{\ell=0}^{d_{j}-1}\frac{k!}{\left(k-\ell\right)!}a_{j,\ell}x_{j}^{k-\ell}=m_{k},\quad k=0,1,\dots,2d-1,

(1.5)

which is called a confluent Prony system of order $d$ with the multiplicity vector $D=\left(d_{1},\dots,d_{s}\right)$ . The original Prony system (1.3) is a special case of the confluent one, with $D$ being the vector $(1,\dots,1)$ of length $d$ .

The system (1.5) arises also in the problem of reconstructing a planar polygon $P$ (or even an arbitrary semi-analytic quadrature domain) from its moments

m_{k}(\chi_{P})=\iint_{\mathbb{R}^{2}}z^{k}\chi_{P}\operatorname{d}x\operatorname{d}y,\;z=x+\imath y,

where $\chi_{P}$ is the characteristic function of the domain $P\subset\mathbb{R}^{2}$ . This problem is important in many areas of science and engineering [11]. The above yields the confluent Prony system

m_{k}=\sum_{j=1}^{s}\sum_{i=0}^{d_{j}-1}c_{i,j}k(k-1)\cdots(k-i+1)z_{j}^{k-i},\qquad c_{i,j}\in\mathbb{C},\;z_{j}\in\mathbb{C}\setminus\left\{0\right\}.

Definition 1.2.

For a given multiplicity vector $D=\left(d_{1},\dots,d_{s}\right)$ , its order is $\sum_{j=1}^{s}d_{j}$ .

As we shall see below, if we start with the measurements $\mu(F)=\mu=(m_{0},\dots,m_{2d-1})$ , then a natural setting of the problem of solving the Prony system is the following:

Problem 1.3 (Prony problem of order $d$ ).

Given the measurements

\mu=(m_{0},\dots,m_{2d-1})\in\mathbb{C}^{2d}

in the right hand side of (1.5), find the multiplicity vector $D=(d_{1},\dots,d_{s})$ of order $r=\sum_{j=1}^{s}d_{j}\leq d$ , and find the unknowns $x_{j}$ and $a_{j,\ell},$ which solve the corresponding confluent Prony system (1.5) with the multiplicity vector $D$ (hence, with solution of rank $r$ ).

It is extremely important in practice to have a stable method of inversion. Many research efforts are devoted to this task (see e.g. [3, 7, 10, 17, 18, 20] and references therein). A basic question here is the following.

Problem 1.4 (Noisy Prony problem).

Given the noisy measurements

\tilde{\mu}=(\tilde{m_{0}},\dots,\tilde{m}_{2d-1})\in\mathbb{C}^{2d}

and an estimate of the error $\left|\tilde{m}_{k}-m_{k}\right|\leq\varepsilon_{k}$ , solve Problem 1.3 so as to minimize the reconstruction error.

In this paper we study the global setting of the Prony problem, stressing its algebraic structure. In Section 2 the space where the solution is to be found (Prony space) is described. It turns out to be a vector bundle over the space of the nodes $x_{1},\dots,x_{d}$ . We define also three mappings: “Prony”, “Taylor”, and “Stieltjes” ones, which capture the essential features of the Prony problem and of its solution process.

In Section 3 we investigate solvability conditions for the Prony problem. The answer leads naturally to a stratification of the space of the right-hand sides, according to the rank of the associated Hankel-type matrix and its minors. The behavior of the solutions near various strata turns out to be highly nontrivial, and we present some initial results in the description of the corresponding singularities.

In Section 4, we study the multiplicity-restricted Prony problem, fixing the collision pattern of the solution, and derive simple bounds for the stability of the solution via factorization of the Jacobian determinant of the corresponding Prony map.

In Section 5 we consider the rank-restricted Prony problem, effectively reducing the dimension to $2r$ instead of $2d$ , where $r$ is precisely the rank of the associated Hankel-type matrix. In this formulation, the Prony problem is solvable in a small neighborhood of the exact measurement vector.

In Section 6 we study one of the most important singularities in the Prony problem: collision of some of the points $x_{i}.$ The investigation of this type of singularities has been started in [21] where the role of finite differences was demonstrated. In the present paper we introduce global bases of finite differences, study their properties, and prove that using such bases we can resolve in a robust way at least the linear part of the Prony problem at and near colliding configurations of the nodes.

In Section 7 we discuss close connections of the Prony problem with hyperbolic polynomials and “Vandermonde mapping” studied by V.I.Arnold in [1] and by V.P.Kostov in [13, 14, 15], and with “Vieta mapping” expressing the coefficients of a polynomial through its roots. We believe that questions arising in theoretical study of Prony problem and in its practical applications justify further investigation of these connections, as well as further applications of Singularity Theory.

Finally, in Appendix A we describe a solution method for the Prony system based on Padé approximation.

2. Prony, Stieltjes and Taylor Mappings

In this section we define “Prony”, “Taylor”, and “Stieltjes” mappings, which capture some essential features of the Prony problem and of its solution process. The main idea behind the spaces and mappings introduced in this section is the following: associate to the signal $F(x)=\sum_{i=1}^{d}a_{i}\delta(x-x_{i})$ the rational function $R(z)=\sum_{i=1}^{d}\frac{{a_{i}}}{{z-x_{i}}}$ . (In fact, $R$ is the Stieltjes integral transform of $F$ ). The functions $R$ obtained in this way can be written as $R(z)=\frac{{P(z)}}{{Q(z)}}$ with $\deg P\leq\deg Q-1,$ and they satisfy $R(\infty)=0$ . Write $R$ as $R(z)=\sum_{i=1}^{d}\frac{{za_{i}}}{{1-x_{i}/z}}.$ Developing the summands into geometric progressions we conclude that $R(z)=\sum_{k=0}^{\infty}m_{k}(\frac{1}{z})^{k+1},$ with $m_{k}=\sum_{i=1}^{d}a_{i}x_{i}^{k}$ , so the moment measurements $m_{k}$ in the right hand side of the Prony system (1.3) are exactly the Taylor coefficients of $R(z)$ . We shall see below that this correspondence reduces solution of the Prony system to an appropriate Padé approximation problem.

Definition 2.1.

For each $w=\left(x_{1},\dots,x_{d}\right)\in\mathbb{C}^{d}$ , let $s=s\left(w\right)$ be the number of distinct coordinates $\tau_{j}$ , $j=1,\dots,s$ , and denote $T\left(w\right)=\left(\tau_{1},\dots,\tau_{s}\right)$ . The multiplicity vector is $D=D\left(w\right)=\left(d_{1},\dots,d_{s}\right)$ , where $d_{j}$ is the number of times the value $\tau_{j}$ appears in $\left\{x_{1},\dots,x_{d}\right\}.$ The order of the values in $T\left(w\right)$ is defined by their order of appearance in $w$ .

Example 2.2.

For $w=\left(3,1,2,1,0,3,2\right)$ we have $s\left(w\right)=4$ , $T\left(w\right)=\left(3,1,2,0\right)$ and $D\left(w\right)=\left(2,2,2,1\right)$ .

Remark 2.3.

Note the slight abuse of notations between Definition 1.1 and Definition 2.1. Note also that the order of $D\left(w\right)$ equals to $d$ for all $w\in\mathbb{C}^{d}$ .

Definition 2.4.

For each $w\in\mathbb{C}^{d}$ , let $s=s\left(w\right),\;T\left(w\right)=\left(\tau_{1},\dots,\tau_{s}\right)$ and $D\left(w\right)=\left(d_{1},\dots,d_{s}\right)$ be as in Definition 2.1.

(1)

$V_{w}$ is the vector space of dimension $d$ containing the linear combinations

$g=\sum_{j=1}^{s}\sum_{\ell=0}^{d_{j}-1}\gamma_{j,\ell}\delta^{\left(\ell\right)}\left(x-\tau_{j}\right)$ (2.1)

of $\delta$ -functions and their derivatives at the points of $T\left(w\right)$ . The “standard basis” of $V_{w}$ is given by the distributions

$\delta_{j,\ell}=\delta^{\left(\ell\right)}\left(x-\tau_{j}\right),\qquad j=1,\dots,s\left(w\right);\;\ell=0,\dots,d_{j}-1.$ (2.2)

(2)

$W_{w}$ is the vector space of dimension $d$ of all the rational functions with poles $T\left(w\right)$ and multiplicities $D\left(w\right)$ , vanishing at $\infty:$

R\left(z\right)=\frac{P\left(z\right)}{Q\left(z\right)},\qquad Q\left(z\right)=\prod_{j=1}^{s}\left(z-\tau_{j}\right)^{d_{j}},\;\deg P\left(z\right)<\deg Q\leqslant d.

The “standard basis” of $W_{w}$ is given by the elementary fractions

R_{j,\ell}=\frac{1}{\left(z-\tau_{j}\right)^{\ell}},\qquad j=1,\dots,s;\;\ell=1,\dots,d_{j}.

Now we are ready to formally define the Prony space ${\mathcal{}P}_{d}$ and the Stieltjes space ${\mathcal{}S}_{d}$ .

Definition 2.5.

The Prony space ${\mathcal{}P}_{d}$ is the vector bundle over $\mathbb{C}^{d}$ , consisting of all the pairs

\left(w,g\right):\quad w\in\mathbb{C}^{d},\;g\in V_{w}.

The topology on ${\mathcal{}P}_{d}$ is induced by the natural embedding ${\mathcal{}P}_{d}\subset\mathbb{C}^{d}\times{\mathcal{}D},$ where ${\mathcal{}D}$ is the space of distributions on $\mathbb{C}$ with its standard topology.

Definition 2.6.

The Stieltjes space ${\mathcal{}S}_{d}$ is the vector bundle over $\mathbb{C}^{d}$ , consisting of all the pairs

\left(w,\gamma\right):\qquad w\in\mathbb{C}^{d},\;\gamma\in W_{w}.

The topology on ${\mathcal{}S}_{d}$ is induced by the natural embedding ${\mathcal{}S}_{d}\subset\mathbb{C}^{d}\times{\mathcal{}R}$ , where ${\mathcal{}R}$ is the space of complex rational functions with its standard topology.

Definition 2.7.

The Stieltjes mapping ${\mathcal{}SM}:{\mathcal{}P}_{d}\rightarrow{\mathcal{}S}_{d}$ is defined by the Stieltjes integral transform: for $(w,g)\in{\mathcal{}P}_{d}$

{\mathcal{}SM}\left(\left(w,g\right)\right)=\left(w,\gamma\right),\qquad\gamma\left(z\right)=\int_{-\infty}^{\infty}\frac{g\left(x\right)\operatorname{d}x}{z-x}.

Sometimes we abuse notation and write for short ${\mathcal{}SM}\left(g\right)=\gamma$ , with the understanding that ${\mathcal{}SM}$ is also a map ${\mathcal{}SM}:V_{w}\to W_{w}$ for each $w\in\mathbb{C}^{d}$ .

The following fact is immediate consequence of the above definitions.

Proposition 2.8.

${\mathcal{}SM}$ is a linear isomorphism of the bundles ${\mathcal{}P}_{d}$ and ${\mathcal{}S}_{d}$ (for each $w\in\mathbb{C}^{d}$ , ${\mathcal{}SM}$ is a linear isomorphism of the vector spaces $V_{w}$ and $W_{w}$ ). In the standard bases of $V_{w}$ and $W_{w}$ , the map ${\mathcal{}SM}$ is diagonal, satisfying

{\mathcal{}SM}\left(\delta_{j,\ell}\right)=\left(-1\right)^{\ell}\ell!R_{j,\ell}\left(z\right).

Furthermore, for any $\left(w,g\right)\in{\mathcal{}P}_{d}$

{\mathcal{}SM}\left(g\right)=\underbrace{\frac{P\left(z\right)}{Q\left(z\right)}}_{\text{irreducible}},\qquad\deg P<\deg Q=\operatorname{rank}\left(g\right)\leqslant d.

(2.3)

Definition 2.9.

The Taylor space ${\mathcal{}T}_{d}$ is the space of complex Taylor polynomials at infinity of degree $2d-1$ of the form $\sum_{k=0}^{2d-1}m_{k}(\frac{1}{z})^{k+1}$ . We shall identify ${\mathcal{}T}_{d}$ with the complex space $\mathbb{C}^{2d}$ with the coordinates $m_{0},\dots,m_{2d-1}$ .

Definition 2.10.

The Taylor mapping ${\mathcal{}TM}:{\mathcal{}S}_{d}\to{\mathcal{}T}_{d}$ is defined by the truncated Taylor development at infinity:

{\mathcal{}TM}\left(\left(w,\gamma\right)\right)=\sum_{k=0}^{2d-1}\alpha_{k}\left(\frac{1}{z}\right)^{k+1},\qquad\text{ where }\gamma\left(z\right)=\sum_{k=0}^{\infty}\alpha_{k}\left(\frac{1}{z}\right)^{k+1}.

We identify ${\mathcal{}TM}\left(\left(w,\gamma\right)\right)$ as above with $\left(\alpha_{0},\dots,\alpha_{2d-1}\right)\in\mathbb{C}^{2d}.$ Sometimes we write for short ${\mathcal{}TM}\left(\gamma\right)=\left(\alpha_{0},\dots,\alpha_{2d-1}\right)$ .

Finally, we define the Prony mapping ${\mathcal{}PM}$ which encodes the Prony problem.

Definition 2.11.

The Prony mapping ${\mathcal{}PM}:{\mathcal{}P}_{d}\to\mathbb{C}^{2d}$ for $\left(w,g\right)\in{\mathcal{}P}_{d}$ is defined as follows:

{\mathcal{}{\mathcal{}PM}}\left(\left(w,g\right)\right)=\left(m_{0},\dots,m_{2d-1}\right)\in\mathbb{C}^{2d},\qquad m_{k}=m_{k}\left(g\right)=\int x^{k}g\left(x\right)\operatorname{d}x.

By the above definitions, we have

{\mathcal{}PM}={\mathcal{}TM}\circ{\mathcal{}SM}.

(2.4)

Solving the Prony problem for a given right-hand side $(m_{0},\dots,m_{2d-1})$ is therefore equivalent to inverting the Prony mapping ${\mathcal{}PM}$ . As we shall elaborate in the subsequent section, the identity (2.4) allows us to split this problem into two parts: inversion of ${\mathcal{}TM}$ , which is, essentially, the Padé approximation problem, and inversion of ${\mathcal{}SM}$ , which is, essentially, the decomposition of a given rational function into the sum of elementary fractions.

3. Solvability of the Prony problem

3.1. General condition for solvability

In this section we provde a necessary and sufficient condition for the Prony problem to have a solution (which is unique, as it turns out by Proposition 3.2). As mentioned in the end of the previous section, our method is based on inverting (2.4) and thus relies on the solution of the corresponding (diagonal) Padé approximation problem [4].

Problem 3.1 (Diagonal Padé approximation problem).

Given $\mu=\left(m_{0},\dots,m_{2d-1}\right)\in\mathbb{C}^{2d}$ , find a rational function $R_{d}(z)=\frac{P\left(z\right)}{Q\left(z\right)}\in{\mathcal{}S}_{d}$ with $\deg P<\deg Q\leqslant d$ , such that the first $2d$ Taylor coefficients at infinity of $R_{d}(z)$ are $\left\{m_{k}\right\}_{k=0}^{2d-1}$ .

Proposition 3.2.

A solution to Problem 3.1, if exists, is unique.

Proof.

Writing $R\left(z\right)=\frac{P\left(z\right)}{Q\left(z\right)},\;R_{1}\left(z\right)=\frac{P_{1}\left(z\right)}{Q_{1}\left(z\right)}$ , with $\deg P<\deg Q\leqslant d$ and $\deg P_{1}<\deg Q_{1}\leqslant d$ , we get

R-R_{1}=\frac{PQ_{1}-P_{1}Q}{QQ_{1}},

and this function, if nonzero, can have a zero of order at most $2d-1$ at infinity. ∎

Let us summarize the above discussion with the following statement.

Proposition 3.3.

The tuple

\left\{s,\;D=(d_{1},\dots,d_{s}),\;r=\sum_{j=1}^{s}d_{j}\leq d,\;X=\left\{x_{j}\right\}_{j=1}^{s},\;A=\left\{a_{j,\ell}\right\}_{j=1,\dots,s;\;\ell=0,\dots,d_{j}-1}\right\}

is a (unique, up to a permutation of the nodes $\left\{x_{j}\right\}$ ) solution to Problem 1.3 with right-hand side

\mu=(m_{0},\dots,m_{2d-1})\in\mathbb{C}^{2d}

if and only if the rational function

R_{D,X,A}\left(z\right)=\sum_{j=1}^{s}\sum_{\ell=1}^{d_{j}}\left(-1\right)^{\ell-1}\left(\ell-1\right)!\frac{a_{j,\ell-1}}{\left(z-x_{j}\right)^{\ell}}=\sum_{k=0}^{2d-1}\frac{m_{k}}{z^{k+1}}+O\left(z^{-2d-1}\right)

is a (unique) solution to Problem 3.1 with input $\mu$ . In that case,

R_{D,X,A}\left(z\right)=\int_{-\infty}^{\infty}\frac{g\left(x\right)\operatorname{d}x}{z-x}\qquad\text{ where }\;g\left(x\right)=\sum_{j=1}^{s}\sum_{\ell=0}^{d_{j}-1}a_{j,\ell}\delta^{\left(\ell\right)}\left(x-x_{j}\right),

i.e. $R_{D,X,A}\left(z\right)$ is the Stieltjes transform of $g\left(x\right)$ .

Proof.

This follows from the definitions of Section 2, (2.4), Proposition 3.2 and the fact that the problem of representing a given rational function as a sum of elementary fractions of the specified form (i.e. inverting ${\mathcal{}SM}$ ) is always uniquely solvable up to a permutation of the poles. ∎

The next result provides necessary and sufficient conditions for the solvability of Problem 3.1. It summarizes some well-known facts in the theory of Padé approximation, related to “normal indices” (see, for instance, [4]). However, these facts are not usually formulated in the literature on Padé approximation in the form we need in relation to the Prony problem. Consequently, we give a detailed proof of this result in Appendix A. This proof contains, in particular, some facts which are important for understanding the solvability issues of the Prony problem.

Definition 3.4.

Given a vector $\mu=\left(m_{0},\dots,m_{2d-1}\right)$ , let $\tilde{M}_{d}$ denote the $d\times\left(d+1\right)$ Hankel matrix

\tilde{M}_{d}=\begin{bmatrix}m_{0}&m_{1}&m_{2}&\dots&m_{d}\\ m_{1}&m_{2}&m_{3}&\dots&m_{d+1}\\ \adots&\adots&\adots&\adots&\adots\\ m_{d-1}&m_{d}&m_{d+1}&\dots&m_{2d-1}\end{bmatrix}.

(3.1)

For each $e\leqslant d$ , denote by $\tilde{M}_{e}$ the $e\times\left(e+1\right)$ submatrix of $\tilde{M}_{d}$ formed by the first $e$ rows and $e+1$ columns, and let $M_{e}$ denote the corresponding square matrix.

Theorem 3.5.

Let $\mu=(m_{0},\dots,m_{2d-1})$ be given, and let $r\leqslant d$ be the rank of the Hankel matrix $\tilde{M}_{d}$ as in (3.1). Then Problem 3.1 is solvable for the input $\mu$ if and only if the upper left minor $\left|M_{r}\right|$ of $\tilde{M}_{d}$ is non-zero.

As an immediate consequence of Theorem 3.5 and Proposition 3.3, we obtain the following result.

Theorem 3.6.

Let $\mu=(m_{0},\dots,m_{2d-1})$ be given, and let $r\leqslant d$ be the rank of the Hankel matrix $\tilde{M}_{d}$ as in (3.1). Then Problem 1.3 with input $\mu$ is solvable if and only if the upper left minor $\left|M_{r}\right|$ of $\tilde{M}_{d}$ is non-zero. The solution, if exists, is unique, up to a permutation of the nodes $\left\{x_{j}\right\}$ . The multiplicity vector $D=\left(d_{1},\dots,d_{s}\right)$ , of order $\sum_{j=1}^{s}d_{j}=r$ , of the resulting confluent Prony system of rank $r$ is the multiplicity vector of the poles of the rational function $R_{D,X,A}\left(z\right)$ , solving the corresponding Padé problem.

As a corollary we get a complete description of the right-hand side data $\mu\in\mathbb{C}^{2d}$ for which the Prony problem is solvable (unsolvable). Define for $r=1,\dots,d$ sets $\Sigma_{r}\subset\mathbb{C}^{2d}$ (respectively, $\Sigma^{\prime}_{r}\subset\mathbb{C}^{2d}$ ) consisting of $\mu\in\mathbb{C}^{2d}$ for which the rank of $\tilde{M}_{d}=r$ and $|M_{r}|\neq 0$ (respectively, $|M_{r}|=0)$ . The set $\Sigma_{r}$ is a difference $\Sigma_{r}=\Sigma_{r}^{1}\setminus\Sigma_{r}^{2}$ of two algebraic sets: $\Sigma_{r}^{1}$ is defined by vanishing of all the $s\times s$ minors of $\tilde{M}_{d},\ r<s\leq d,$ while $\Sigma_{r}^{2}$ is defined by vanishing of $|M_{r}|.$ In turn, $\Sigma^{\prime}_{r}=\Sigma_{r}^{{}^{\prime}1}\setminus\Sigma_{r}^{{}^{\prime}2},$ with $\Sigma_{r}^{{}^{\prime}1}=\Sigma_{r}^{1}\cap\Sigma_{r}^{2}$ and $\Sigma_{r}^{{}^{\prime}2}$ defined by vanishing of all the $r\times r$ minors of $\tilde{M}_{d}.$ The union $\Sigma_{r}\cup\Sigma^{\prime}_{r}$ consists of all $\mu$ for which the rank of $\tilde{M}_{d}=r,$ which is $\Sigma_{r}^{1}\setminus\Sigma_{r}^{{}^{\prime}2}.$

Corollary 3.7.

The set $\Sigma$ (respectively, $\Sigma^{\prime}$ ) of $\mu\in\mathbb{C}^{2d}$ for which the Prony problem is solvable (respectively, unsolvable) is the union $\Sigma=\cup_{r=1}^{d}\Sigma_{r}$ (respectively, $\Sigma^{\prime}=\cup_{r=1}^{d}\Sigma^{\prime}_{r}$ ). In particular, $\Sigma^{\prime}\subset\{\mu\in\mathbb{C}^{2d},\det M_{d}=0\}.$

So for a generic right hand side $\mu$ we have $|M_{d}|\neq 0$ , and the Prony problem is solvable. On the algebraic hypersurface of $\mu$ for which $|M_{d}|=0,$ the Prony problem is solvable if $M_{d-1}\neq 0$ , etc.

Let us now consider some examples.

Example 3.8.

Let us fix $d=1,2,\dots$ . Consider $\mu=(m_{0},\dots,m_{2d-1})\in\mathbb{C}^{2d}$ , the right hand sides of the Prony problem, to be of the form $\mu=\mu_{\ell}=\left(\delta_{k\ell}\right)=(0,\dots,0,\underbrace{1}_{\text{ position $\ell$+1}},0,\dots,0)$ , with all the $m_{k}=0$ besides $m_{\ell}=1,\ \ell=0,\dots,2d-1,$ and let $\tilde{M}_{d}^{\ell}$ be the corresponding matrix.

Proposition 3.9.

The rank of $\tilde{M}_{d}^{\ell}$ is equal to $\ell+1$ for $\ell\leq d-1$ , and it is equal to $2d-\ell$ for $\ell\geq d$ . The corresponding Prony problem is solvable for $\ell\leq d-1$ , and it is unsolvable for $\ell\geq d$ .

Proof.

For $d=5$ and $\ell=2,4,5,9$ , the corresponding matrices $\tilde{M}_{\ell}^{d}$ are as follows.

	$\displaystyle\tilde{M}_{5}^{2}$	$\displaystyle=$	$\displaystyle\begin{bmatrix}0&0&1&0&0&0\\ 0&1&0&0&0&0\\ 1&0&0&0&0&0\\ 0&0&0&0&0&0\\ 0&0&0&0&0&0\end{bmatrix},\;\tilde{M}_{5}^{4}=\begin{bmatrix}0&0&0&0&1&0\\ 0&0&0&1&0&0\\ 0&0&1&0&0&0\\ 0&1&0&0&0&0\\ 1&0&0&0&0&0\end{bmatrix},\qquad\text{(solvable)}$
	$\displaystyle\tilde{M}_{5}^{5}$	$\displaystyle=$	$\displaystyle\begin{bmatrix}0&0&0&0&0&1\\ 0&0&0&0&1&0\\ 0&0&0&1&0&0\\ 0&0&1&0&0&0\\ 0&1&0&0&0&0\end{bmatrix},\;\tilde{M}_{5}^{9}=\begin{bmatrix}0&0&0&0&0&0\\ 0&0&0&0&0&0\\ 0&0&0&0&0&0\\ 0&0&0&0&0&0\\ 0&0&0&0&0&1\end{bmatrix}.\qquad\text{(unsolvable)}$

In general, the matrices $\tilde{M}_{d}^{\ell}$ have the same pattern as in the special cases above, so their rank is $\ell+1$ for $\ell\leqslant d-1$ , and $2d-\ell$ for $\ell\geqslant d$ , as stated above. Application of Theorem 3.6 completes the proof. ∎

In fact, $\mu_{\ell}$ is a moment sequence of

F\left(x\right)=\frac{1}{\ell!}\delta^{\left(\ell\right)}\left(x\right),

and this signal belongs to ${\mathcal{}P}_{d}$ if and only if $\ell\leqslant d-1$ . In notations of Corollary 3.7 we have

\displaystyle\begin{aligned} \mu_{\ell}&\in&\Sigma_{\ell+1},&&\ell\leqslant d-1,\\ \mu_{\ell}&\in&\Sigma^{\prime}_{2d-\ell},&&\ell\geqslant d.\end{aligned}

It is easy to provide various modifications of the above example. In particular, for $\mu=\tilde{\mu}_{\ell}=\left(0,\dots,0,1,1,\dots,1\right)$ , the result of Proposition 3.9 remains verbally true.

Example 3.10.

Another example is provided by $\mu_{\ell_{1},\ell_{2}}$ , with all the $m_{k}=0$ besides $m_{\ell_{1}}=1,\;m_{\ell_{2}}=1,\ 0\leq\ell_{1}<d\leq\ell_{2}\leq 2d-1.$ For $\ell_{1}<\ell_{2}-d+1$ the rank of the correspondent matrix $\tilde{M}_{d}$ is $r=2d+\ell_{1}-\ell_{2}+1$ while $|M_{r}|=0$ , so the Prony problem for such $\mu_{\ell_{1},\ell_{2}}$ is unsolvable. For $d=5$ and $\ell_{1}=2,\;\ell_{2}=8$ the matrix is as follows:

\tilde{M}_{5}^{\left(2,8\right)}=\begin{bmatrix}0&0&1&0&0&0\\ 0&1&0&0&0&0\\ 1&0&0&0&0&0\\ 0&0&0&0&0&1\\ 0&0&0&0&1&0\end{bmatrix}.

3.2. Near-singular inversion

The behavior of the inversion of the Prony mapping near the unsolvability stratum $\Sigma^{\prime}$ and near the strata where the rank of $\tilde{M}_{d}$ drops, turns out to be pretty complicated. In particular, in the first case at least one of the nodes tends to infinity. In the second case, depending on the way the right-hand side $\mu$ approaches the lower rank strata, the nodes may remain bounded, or some of them may tend to infinity. In this section we provide one initial result in this direction, as well as some examples. A comprehensive description of the inversion of the Prony mapping near $\Sigma^{\prime}$ and near the lower rank strata is important both in theoretical study and in applications of Prony-like systems, and we plan to provide further results in this direction separately.

Theorem 3.11.

As the right-hand side $\mu\in\mathbb{C}^{2d}\setminus\Sigma^{\prime}$ approaches a finite point $\mu_{0}\in\Sigma^{\prime},$ at least one of the nodes $x_{1},\dots,x_{d}$ in the solution tends to infinity.

Proof.

By assumptions, the components $m_{0},\dots,m_{2d-1}$ of the right-hand side $\mu=(m_{0},\dots,m_{2d-1})\in\mathbb{C}^{2d}$ remain bounded as $\mu\rightarrow\mu_{0}$ . By Theorem 6.14, the finite differences coordinates of the solution ${\mathcal{}PM}^{-1}(\mu)$ remain bounded as well. Now, if all the nodes are also bounded, by compactness we conclude that ${\mathcal{}PM}^{-1}(\mu)\rightarrow\omega\in{\mathcal{}P}_{d}.$ By continuity in the distribution space (Lemma 6.6) we have ${\mathcal{}PM}(\omega)=\mu_{0}$ . Hence the Prony problem with the right-hand side $\mu_{0}$ has a solution $\omega\in{\mathcal{}P}_{d},$ in contradiction with the assumption that $\mu_{0}\in\Sigma^{\prime}$ .∎

Example 3.12.

Let us consider an example: $d=2$ and $\mu_{0}=(0,0,1,0)$ . Here the rank $\ell$ of $\tilde{M}_{2}$ is $2$ , and $|M_{2}|=0$ , so by Theorem 3.6 we have $\mu_{0}\in\Sigma^{\prime}_{2}\subset\Sigma^{\prime}$ . Consider now a perturbation $\mu(\epsilon)=(0,\epsilon,1,0)$ of $\mu_{0}$ . For $\epsilon\neq 0$ we have $\mu(\epsilon)\in\Sigma_{2}\subset\Sigma,$ and the Prony system is solvable for $\mu_{\epsilon}$ . Let us write an explicit solution: the coefficients $c_{0},c_{1}$ of the polynomial $Q(z)=c_{0}+c_{1}z+z^{2}$ we find from the system (A. $\star\star$ ):

\begin{bmatrix}0&\epsilon\\ \epsilon&1\end{bmatrix}\begin{bmatrix}c_{0}\\ c_{1}\end{bmatrix}=\begin{bmatrix}-1\\ 0\end{bmatrix},

whose solution is $c_{1}=-\frac{1}{\epsilon},\ c_{0}=\frac{1}{{\epsilon^{2}}}.$ Hence the denominator $Q(z)$ of $R(z)$ is $Q(z)=\frac{1}{{\epsilon^{2}}}-\frac{1}{\epsilon}z+z^{2}$ , and its roots are $x_{1}=\frac{{1+\imath\sqrt{3}}}{{2\epsilon}},\ x_{2}=\frac{{1-\imath\sqrt{3}}}{{2\epsilon}}$ . The coefficients $b_{0},b_{1}$ of the numerator $P(z)=b_{0}+b_{1}z$ we find from (A. $\star$ ):

\begin{bmatrix}0&0\\ 0&\epsilon\end{bmatrix}\begin{bmatrix}-\frac{1}{\epsilon}\\ 1\end{bmatrix}=\begin{bmatrix}b_{1}\\ b_{0}\end{bmatrix},

i.e. $b_{1}=0,\ b_{0}=\epsilon$ . Thus the solution of the associated Padé problem is

R(z)=\frac{{P(z)}}{{Q(z)}}=\frac{\epsilon}{{(z-x_{1})(z-x_{2})}}=\frac{\epsilon^{2}}{\imath\sqrt{3}}\frac{1}{(z-x_{1})}-\frac{\epsilon^{2}}{\imath\sqrt{3}}\frac{1}{(z-x_{2})}.

Finally, the (unique up to a permutation) solution of the Prony problem for $\mu_{\epsilon}$ is

a_{1}=\frac{\epsilon^{2}}{\imath\sqrt{3}},\ a_{2}=-\frac{\epsilon^{2}}{\imath\sqrt{3}},\quad x_{1}=\frac{{1+\imath\sqrt{3}}}{{2\epsilon}},\ x_{2}=\frac{{1-\imath\sqrt{3}}}{{2\epsilon}}.

As $\epsilon$ tends to zero, the nodes $x_{1},x_{2}$ tend to infinity while the coefficients $a_{1},a_{2}$ tend to zero.

As it was shown above, for a given $\mu\in\Sigma$ (say, with pairwise different nodes) the rank of the matrix $\tilde{M}_{d}$ is equal to the number of the nodes in the solution for which the corresponding $\delta$ -function enters with a non-zero coefficients. So $\mu$ approaches a certain $\mu_{0}$ belonging to a stratum of a lower rank of $\tilde{M}_{d}$ if and only if some of the coefficients $a_{j}$ in the solution tend to zero. We do not analyze all the possible scenarios of such a degeneration, noticing just that if $\mu_{0}\in\Sigma^{\prime},$ i.e., the Prony problem is unsolvable for $\mu_{0}$ , then Theorem 3.11 remains true, with essentially the same proof. So at least one of the nodes, say, $x_{j},$ escapes to infinity. Moreover, one can show that $a_{j}x_{j}^{2d-1}$ cannot tend to zero - otherwise the remaining linear combination of $\delta$ -functions would provide a solution for $\mu_{0}$ .

If $\mu_{0}\in\Sigma,$ i.e., the Prony problem is solvable for $\mu_{0},$ all the nodes may remain bounded, or some $x_{j}$ may escape to infinity, but in such a way that $a_{j}x_{j}^{2d-1}$ tends to zero.

4. Multiplicity-restricted Prony problem

Consider Problem 1.4 at some point $\mu_{0}\in\Sigma$ . By definition, $\mu_{0}\in\Sigma_{r_{0}}$ for some $r_{0}\leq d$ . Let $\mu_{0}={\mathcal{}PM}\left(\left(w_{0},g_{0}\right)\right)$ for some $\left(w_{0},g_{0}\right)\in{\mathcal{}P}_{d}$ . Assume for a moment that the multiplicity vector $D_{0}=D\left(g_{0}\right)=\left(d_{1},\dots d_{s_{0}}\right)$ , $\sum_{j=1}^{s_{0}}d_{j}=r_{0}$ , has a non-trivial collision pattern, i.e. $d_{j}>1$ for at least one $j=1,\dots,s_{0}$ . It means, in turn, that the function $R_{D_{0},X,A}\left(z\right)$ has a pole of multiplicity $d_{j}$ . Evidently, there exists an arbitrarily small perturbation $\tilde{\mu}$ of $\mu_{0}$ for which this multiple pole becomes a cluster of single poles, thereby changing the multiplicity vector to some $D^{\prime}\neq D_{0}$ . While we address this problem in Section 6 via the bases of divided differences, in this section we consider a “multiplicity-restricted” Prony problem.

Definition 4.1.

Let $\mathbf{x}=\left(x_{1},\dots,x_{s}\right)\in\mathbb{C}^{s}$ and $D=\left(d_{1},\dots,d_{s}\right)$ with $d=\sum_{j=1}^{s}d_{j}$ be given. The $d\times d$ confluent Vandermonde matrix is

V=V\left(\mathbf{x},D\right)=V\left(x_{1},d_{1},\dots,x_{s},d_{s}\right)=\left[\begin{array}[]{cccc}\mathbf{v_{1,0}}&\mathbf{v_{2,0}}&\dotsc&\mathbf{v_{s,0}}\\ \mathbf{v_{1,1}}&\mathbf{v_{2,1}}&\dotsc&\mathbf{v_{s,1}}\\ &&\dotsc\\ \mathbf{v_{1,d-1}}&\mathbf{v_{2,d-1}}&\dotsc&\mathbf{v_{s,d-1}}\end{array}\right]

(4.1)

where the symbol $\mathbf{v_{j,k}}$ denotes the following $1\times d_{j}$ row vector

\mathbf{v_{j,k}}\stackrel{{\scriptstyle\text{def}}}{{=}}\left[\begin{array}[]{cccc}x_{j}^{k},&kx_{j}^{k-1},&\dots&,k\left(k-1\right)\cdots\left(k-d_{j}\right)x_{j}^{k-d_{j}+1}\end{array}\right].

Proposition 4.2.

The matrix $V$ defines the linear part of the confluent Prony system (1.5) in the standard basis for $V_{w}$ , namely,

V\left(x_{1},d_{1},\dots,x_{s},d_{s}\right)\begin{bmatrix}a_{1,0}\\ \vdots\\ a_{1,d_{1}-1}\\ \vdots\\ \\ a_{s,d_{s}-1}\end{bmatrix}=\begin{bmatrix}m_{0}\\ m_{1}\\ \vdots\\ \\ \\ m_{d-1}\end{bmatrix}.

(4.2)

Definition 4.3.

Let ${\mathcal{}PM}\left(w_{0},g_{0}\right)=\mu_{0}\in\Sigma_{r_{0}}$ with $D\left(g_{0}\right)=D_{0}$ and $s\left(g_{0}\right)=s_{0}$ . Let ${\mathcal{}P}_{D_{0}}$ denote the following subbundle of ${\mathcal{}P}_{d}$ of dimension $s_{0}+r_{0}$ :

{\mathcal{}P}_{D_{0}}=\left\{\left(w,g\right)\in{\mathcal{}P}_{d}:\quad D\left(g\right)=D_{0}\right\}.

The multiplicity-restricted Prony mapping ${\mathcal{}{\mathcal{}PM}}_{D_{0}}^{*}:{\mathcal{}P}_{D_{0}}\to\mathbb{C}^{s_{0}+r_{0}}$ is the composition

{\mathcal{}{\mathcal{}PM}}_{D_{0}}^{*}=\pi\circ{\mathcal{}PM}\restriction_{{\mathcal{}P}_{D_{0}}},

where $\pi:\mathbb{C}^{2d}\to\mathbb{C}^{s_{0}+r_{0}}$ is the projection map on the first $s_{0}+r_{0}$ coordinates.

Inverting this ${\mathcal{}{\mathcal{}PM}}_{D_{0}}^{*}$ represents the solution of the confluent Prony system (1.5) with fixed structure $D_{0}$ from the first $k=0,1,\dots,s_{0}+r_{0}-1$ measurements.

Theorem 4.4 ([7]).

Let $\mu_{0}^{*}={\mathcal{}{\mathcal{}PM}}_{D_{0}}^{*}\left(\left(w_{0},g_{0}\right)\right)\in\mathbb{C}^{s_{0}+r_{0}}$ with the unperturbed solution $g_{0}=\sum_{j=1}^{s_{0}}\sum_{\ell=0}^{d_{j}-1}a_{j,\ell}\delta^{\left(\ell\right)}\left(x-\tau_{j}\right)$ . In a small neighborhood of $\left(w_{0},g_{0}\right)\in{\mathcal{}P}_{D_{0}}$ , the map ${\mathcal{}{\mathcal{}PM}}_{D_{0}}^{*}$ is invertible. Consequently, for small enough $\varepsilon$ , the multiplicity-restricted Prony problem with input data $\tilde{\mu}^{*}\in\mathbb{C}^{r_{0}+s_{0}}$ satisfying $\|\tilde{\mu}^{*}-\mu_{0}^{*}\|\leq\varepsilon$ has a unique solution. The error in this solution satisfies

	$\displaystyle\left\|\Delta a_{j,\ell}\right\|$	$\displaystyle\leq$	$\displaystyle\frac{2}{\ell!}\left(\frac{2}{\delta}\right)^{s_{0}+r_{0}}\left(\frac{1}{2}+\frac{s_{0}+r_{0}}{\delta}\right)^{d_{j}-\ell}\left(1+\frac{\left\|a_{j,\ell-1}\right\|}{\left\|a_{j,d_{j}-1}\right\|}\right)\varepsilon,$
	$\displaystyle\left\|\Delta\tau_{j}\right\|$	$\displaystyle\leq$	$\displaystyle\frac{2}{d_{j}!}\left(\frac{2}{\delta}\right)^{s_{0}+r_{0}}\frac{1}{\left\|a_{j,d_{j}-1}\right\|}\varepsilon,$

where $\delta\stackrel{{\scriptstyle\text{def}}}{{=}}\min_{i\neq j}\left|\tau_{i}-\tau_{j}\right|$ (for consistency we take $a_{j,-1}=0$ in the above formula).

Proof outline.

The Jacobian of ${\mathcal{}{\mathcal{}PM}}_{D_{0}}^{*}$ can be easily computed, and it turns out to be equal to the product

{\mathcal{}J}_{{\mathcal{}{\mathcal{}PM}}_{D_{0}}^{*}}=V\left(\tau_{1},d_{1}+1,\dots,\tau_{s_{0}},d_{s_{0}}+1\right)\operatorname{diag}\left\{E_{j}\right\}

where $V$ is the confluent Vandermonde matrix (4.1) on the nodes $\left(\tau_{1},\dots,\tau_{s_{0}}\right)$ , with multiplicity vector

\tilde{D}_{0}=\left(d_{1}+1,\dots,d_{s_{0}}+1\right),

while $E$ is the $\left(d_{j}+1\right)\times\left(d_{j}+1\right)$ block

E_{j}=\begin{bmatrix}1&0&0&\cdots&0\\ 0&1&0&\cdots&a_{j,0}\\ \vdots&\vdots&\vdots&\ \ddots&\vdots\\ 0&0&0&\cdots&a_{j,d_{j}-1}\end{bmatrix}.

Since $\mu_{0}\in\Sigma_{r}$ , the highest order coefficients $a_{j,d_{j}-1}$ are nonzero. Furthermore, since all the $\tau_{j}$ are distinct, the matrix $V$ is nonsingular. Local invertability follows. To estimate the norm of the inverse, use bounds from [6].∎

Remark 4.5.

Note that as two nodes collide ( $\delta\to 0$ ), the inversion of the multiplicity-restricted Prony mapping ${\mathcal{}{\mathcal{}PM}}_{D_{0}}^{*}$ becomes ill-conditioned proportionally to $\delta^{-\left(s_{0}+r_{0}\right)}$ .

Let us stress that we are not aware of any general method of inverting ${\mathcal{}{\mathcal{}PM}}_{D_{0}}^{*}$ , i.e. solving the multiplicity-restricted confluent Prony problem with the smallest possible number of measurements. As we demonstrate in [5], such a method exists for a very special case of a single point, i.e. $s=1$ .

5. Rank-restricted Prony problem

Recall that the Prony problem consists in inverting the Prony mapping ${\mathcal{}PM}:{\mathcal{}P}_{d}\rightarrow{\mathcal{}T}_{d}$ . So, given $\mu=(m_{0},\dots,m_{2d-1})\in{\mathcal{}T}_{d}$ we are looking for $(w,g)\in{\mathcal{}P}_{d}$ such that $m_{k}(g)=\int x^{k}g(x)dx=m_{k}$ , with $k=0,1,\dots,2d-1$ . If $\mu\in\Sigma_{r}$ with $r<d$ , then in fact any neighborhood of $\mu$ will contain points from the non-solvability set $\Sigma^{\prime}$ . Indeed, consider the following example.

Example 5.1.

Slightly modifying the construction of Example 3.10, consider $\mu_{\ell_{1},\ell_{2},\epsilon}\in\mathbb{C}^{2d}$ with all the $m_{k}=0$ besides $m_{\ell_{1}}=1$ and $m_{\ell_{2}}=\epsilon$ , such that $\ell_{2}>\ell_{1}+d-1$ . For example, if $d=5$ and $\ell_{1}=2,\;\ell_{2}=8$ , the corresponding matrix is

\tilde{M}_{5}^{\left(2,8,\epsilon\right)}=\begin{bmatrix}0&0&1&0&0&0\\ 0&1&0&0&0&0\\ 1&0&0&0&0&0\\ 0&0&0&0&0&\epsilon\\ 0&0&0&0&\epsilon&0\end{bmatrix}.

For $\epsilon=0$ the Prony problem is solvable, while for any small perturbation $\epsilon\neq 0$ it becomes unsolvable. However, if we restrict the whole problem just to $d=3$ , it remains solvable for any small perturbation of the input.

We therefore propose to consider the rank-restricted Prony problem analogous to the construction of Section 4, but instead of fixing the multiplicity $D\left(g\right)$ we now fix the rank $r$ (recall Definition 1.1).

Definition 5.2.

Denote by ${\mathcal{}P}_{r}$ the following vector bundle:

{\mathcal{}P}_{r}=\left\{\left(w,g\right):\quad w\in\mathbb{C}^{r},\;g\in V_{w}\right\},

where $V_{w}$ is defined exactly as in Definition 2.4, replacing $d$ with $r$ .

Likewise, we define the Stieltjes bundle of order $r$ as follows.

Definition 5.3.

Denote by ${\mathcal{}S}_{r}$ the following vector bundle:

{\mathcal{}S}_{r}=\left\{\left(w,\gamma\right):\qquad w\in\mathbb{C}^{r},\;\gamma\in W_{w}\right\},

where $W_{w}$ is defined exactly as in Definition 2.4, replacing $d$ with $r$ .

The Stieltjes mapping acts naturally as a map ${\mathcal{}SM}:{\mathcal{}P}_{r}\to{\mathcal{}S}_{r}$ with exactly the same definition as Definition 2.7.

The restricted Taylor mapping ${\mathcal{}TM}_{r}:{\mathcal{}S}_{r}\to\mathbb{C}^{2r}$ is, as before, given by the truncated development at infinity to the first $2r$ Taylor coefficients.

Definition 5.4.

Let $\pi:\mathbb{C}^{2d}\to\mathbb{C}^{2r}$ denote the projection operator onto the first $2r$ coordinates. Denote $\Sigma_{r}^{*}\stackrel{{\scriptstyle\text{def}}}{{=}}\pi\left(\Sigma_{r}\right)$ . The rank-restricted Prony mapping ${\mathcal{}{\mathcal{}PM}}_{r}^{*}:{\mathcal{}P}_{r}\to\Sigma_{r}^{*}$ is given by by

{\mathcal{}{\mathcal{}PM}}_{r}^{*}\left(\left(w,g\right)\right)=\left(m_{0},\dots,m_{2r-1}\right),\qquad m_{k}=m_{k}\left(g\right)=\int x^{k}g\left(x\right)\operatorname{d}x.

Remark 5.5.

${\mathcal{}P}_{r}$ can be embedded in ${\mathcal{}P}_{d}$ , for example by the map $\Xi_{r}:{\mathcal{}P}_{r}\to{\mathcal{}P}_{d}$

\Xi_{r}:\;\left(w,g\right)\in{\mathcal{}P}_{r}\longmapsto\left(w^{\prime},g^{\prime}\right)\in{\mathcal{}P}_{d}:\qquad w^{\prime}=\left(x_{1},\dots,x_{r},\underbrace{0,\dots 0}_{\times\left(d-r\right)}\right),\;g^{\prime}=g.

With this definition, ${\mathcal{}{\mathcal{}PM}}_{r}^{*}$ can be represented also as the composition

{\mathcal{}{\mathcal{}PM}}_{r}^{*}=\pi\circ{\mathcal{}PM}\circ\Xi_{r}.

Proposition 5.6.

The rank-restricted Prony mapping satisfies

{\mathcal{}{\mathcal{}PM}}_{r}^{*}={\mathcal{}TM}_{r}\circ{\mathcal{}SM}.

Inverting ${\mathcal{}{\mathcal{}PM}}_{r}^{*}$ represents the solution of the rank-restricted Prony problem. Unlike in the multiplicity-restricted setting of Section 4, here we allow two or more nodes to collide (thereby changing the multiplicty vector $D\left(g\right)$ of the solution).

The basic fact which makes this formulation useful is the following result.

Theorem 5.7.

Let $\mu_{0}^{*}\in\Sigma_{r}^{*}$ . Then in a small neighborhood of $\mu_{0}^{*}\in\mathbb{C}^{2r}$ , the Taylor mapping ${\mathcal{}TM}_{r}$ is continuously invertible.

Proof.

This is a direct consequence of the solution method to the Padé approximation problem described in Appendix A. Indeed, if the rank of $\tilde{M}_{r}$ is full, then it remains so in a small neighborhood of the entire space $\mathbb{C}^{2r}$ . Therefore, the system (A. $\star\star$ ) remains continuously invertible, producing the coefficients of the denominator $Q\left(z\right)$ . Consequently, the right-hand side of (A. $\star$ ) depends continuously on the moment vector $\mu^{*}=\left(m_{0},\dots,m_{2r-1}\right)\in\mathbb{C}^{2r}$ . Again, since the rank always remains full, the polynomials $P\left(z\right)$ and $Q\left(z\right)$ cannot have common roots, and thereby the solution $R=\frac{P}{Q}={\mathcal{}TM}_{r}^{-1}\left(\mu^{*}\right)$ depends continuously on $\mu^{*}$ (in the topology of the space of rational functions). ∎

In the next section, we consider the remaining problem: how to invert ${\mathcal{}SM}$ in this setting.

6. Collision singularities and bases of finite differences

6.1. Introduction

Collision singularities occur in Prony systems as some of the nodes $x_{i}$ in the signal $F(x)=\sum_{i=1}^{d}a_{i}\delta(x-x_{i})$ approach one another. This happens for $\mu$ near the discriminant stratum $\Delta\subset\mathbb{C}^{2d}$ consisting of those $(m_{0},\dots,m_{2d-1})$ for which some of the coordinates $\left\{x_{j}\right\}$ in the solution collide, i.e. the function $R_{D,X,A}\left(z\right)$ has multiple poles (or, nontrivial multiplicity vector $D$ ). As we shall see below, typically, as $\mu$ approaches $\mu_{0}\in\Delta$ , i.e. some of the nodes $x_{i}$ collide, the corresponding coefficients $a_{i}$ tend to infinity. Notice, that all the moments $m_{k}=m_{k}(F)$ remain bounded. This behavior creates serious difficulties in solving “near-colliding” Prony systems, both in theoretical and practical settings. Especially demanding problems arise in the presence of noise. The problem of improvement of resolution in reconstruction of colliding nodes from noisy measurements appears in a wide range of applications. It is usually called a “super-resolution problem” and a lot of recent publications are devoted to its investigation in various mathematical and applied settings. See [8] and references therein for a very partial sample.

Here we continue our study of collision singularities in Prony systems, started in [21]. Our approach uses bases of finite differences in the Prony space ${\mathcal{}P}_{r}$ in order to “resolve” the linear part of collision singularities. In these bases the coefficients do not blow up any more, even as some of the nodes collide.

Example 6.1.

Let $r=2$ , and consider the signal $F=a_{1}\delta\left(x-x_{1}\right)+a_{2}\delta\left(x-x_{2}\right)$ with

	$\displaystyle x_{1}$	$\displaystyle=$	$\displaystyle t,\;x_{2}=t+\epsilon,$
	$\displaystyle a_{1}$	$\displaystyle=$	$\displaystyle-\epsilon^{-1},\;a_{2}=\epsilon^{-1}.$

The corresponding Prony system is

\left(a_{1}x_{1}^{k}+a_{2}x_{2}^{k}=\right)m_{k}=kt^{k-1}+\underbrace{\sum_{j=2}^{k}{k\choose j}t^{k-j}\epsilon^{j-1}}_{\stackrel{{\scriptstyle\text{def}}}{{=}}\rho_{k}\left(t,\epsilon\right)},\qquad k=0,1,2,3.

As $\epsilon\to 0$ , the Prony system as above becomes ill-conditioned and the coefficients $\left\{a_{j}\right\}$ blow up, while the measurements remain bounded. Note that

\tilde{M}_{2}=\begin{bmatrix}0&1&2t+\rho_{2}\left(t,\epsilon\right)\\ 1&2t+\rho_{2}\left(t,\epsilon\right)&3t^{2}+\rho_{3}\left(t,\epsilon\right)\end{bmatrix},

therefore $\operatorname{rank}\tilde{M}_{2}=2$ and $\left|M_{2}\right|=1\neq 0$ , i.e. the Prony problem with input $\left(m_{0},\dots,m_{3}\right)$ remains solvable for all $\epsilon$ . However, the standard basis $\left\{\delta\left(x-x_{1}\right),\;\delta\left(x-x_{2}\right)\right\}$ degenerates, and in the limit it is no more a basis. If we represent the solution

F_{\epsilon}\left(x\right)=-\frac{1}{\epsilon}\delta\left(x-t\right)+\frac{1}{\epsilon}\delta\left(x-t-\epsilon\right)

in the basis

	$\displaystyle\Delta_{1}\left(x_{1},x_{2}\right)$	$\displaystyle=$	$\displaystyle\delta\left(x-x_{1}\right),$
	$\displaystyle\Delta_{2}\left(x_{1},x_{2}\right)$	$\displaystyle=$	$\displaystyle\frac{1}{x_{1}-x_{2}}\delta\left(x-x_{1}\right)+\frac{1}{x_{2}-x_{1}}\delta\left(x-x_{2}\right),$

then we have

F_{\epsilon}\left(x\right)=1\cdot\Delta_{2}\left(t,t+\epsilon\right),

i.e. the coefficients in this new basis are just $\left\{b_{1}=0,\;b_{1}=1\right\}$ . As $\epsilon\to 0$ , in fact we have

\Delta_{2}\left(t,t+\epsilon\right)\to\delta^{\prime}\left(x-t\right),

where the convergence is in the topology of the bundle ${\mathcal{}P}_{r}$ .

Our goal in this section is to generalize the construction of Example 6.1 and [21] to handle the general case of colliding configurations.

6.2. Divided finite differences

For modern treatment of divided differences, see e.g. [9, 12, 16]. We follow [9] and adopt what has become by now the standard definition.

Definition 6.2.

Let an arbitrary sequence of points $w=\left(x_{1},x_{2},\dots,\right)$ be given (repetitions are allowed). The (n-1)-st divided difference $\mathord{\kern 4.29993pt\vrule width=0.6pt,height=5.6pt,depth=-0.28pt\kern-4.29993pt\Delta}^{n-1}\left(w\right):\Pi\to\mathbb{C}$ is the linear functional on the space $\Pi$ of polynomials, associating to each $p\in\Pi$ its (uniquely defined) $n$ -th coefficient in the Newton form

p\left(x\right)=\sum_{j=1}^{\infty}\left\{\mathord{\kern 4.29993pt\vrule width=0.6pt,height=5.6pt,depth=-0.28pt\kern-4.29993pt\Delta}^{j-1}\left(x_{1},\dots,x_{j}\right)p\right\}\cdot q_{j-1,w}\left(x\right),\qquad q_{i,w}\left(x\right)\stackrel{{\scriptstyle\text{def}}}{{=}}\prod_{j=1}^{i}\left(x-x_{j}\right).

(6.1)

It turns out that this definition can be extended to all sufficiently smooth functions for which the interpolation problem is well-defined.

Definition 6.3 ([9]).

For any smooth enough function $f$ , defined at least on $x_{1},\dots,x_{n}$ , the divided finite difference $\mathord{\kern 4.29993pt\vrule width=0.6pt,height=5.6pt,depth=-0.28pt\kern-4.29993pt\Delta}^{n-1}\left(x_{1},\dots,x_{n}\right)f$ is the $n$ -th coefficient in the Newton form (6.1) of the Hermite interpolation polynomial $P_{n}$ , which agrees with $f$ and its derivatives of appropriate order on $x_{1},\dots,x_{n}:$

f^{\left(\ell\right)}\left(x_{j}\right)=P_{n}^{\left(\ell\right)}\left(x_{j}\right):\qquad 1\leqslant j\leqslant n,\;0\leqslant\ell<d_{j}\stackrel{{\scriptstyle\text{def}}}{{=}}\#\left\{i:\quad x_{i}=x_{j}\right\}.

(6.2)

Therefore, each divided difference can be naturally associated with an element of the Prony space (see Item 5 in Proposition 6.4 and Definition 6.5 below for an accurate statement).

Let us now summarize relevant properties of the functional $\mathord{\kern 4.29993pt\vrule width=0.6pt,height=5.6pt,depth=-0.28pt\kern-4.29993pt\Delta}$ which we shall use later on.

Proposition 6.4.

For $w=\left(x_{1},\dots,x_{n}\right)\in\mathbb{C}^{n}$ , let $s\left(w\right),\;T\left(w\right)$ and $D\left(w\right)$ be defined according to Definition 2.1. Let $q_{n,w}\left(z\right)=\prod_{j=1}^{s}\left(z-\tau_{j}\right)^{d_{j}}$ be defined as in (6.1).

(1)

The functional $\mathord{\kern 4.29993pt\vrule width=0.6pt,height=5.6pt,depth=-0.28pt\kern-4.29993pt\Delta}^{n-1}\left(x_{1},\dots,x_{n}\right)$ is a symmetric function of its arguments, i.e. it depends only on the set $\left\{x_{1},\dots,x_{n}\right\}$ but not on its ordering.

(2)

$\mathord{\kern 4.29993pt\vrule width=0.6pt,height=5.6pt,depth=-0.28pt\kern-4.29993pt\Delta}^{n-1}\left(x_{1},\dots,x_{n}\right)$ is a continuous function of the vector $\left(x_{1},\dots,x_{n}\right)$ . In particular, for any test function $f$

\lim_{\left(x_{1},\dots,x_{n}\right)\to\left(t_{1},\dots,t_{n}\right)}\mathord{\kern 4.29993pt\vrule width=0.6pt,height=5.6pt,depth=-0.28pt\kern-4.29993pt\Delta}^{n-1}\left(x_{1},\dots,x_{n}\right)f=\mathord{\kern 4.29993pt\vrule width=0.6pt,height=5.6pt,depth=-0.28pt\kern-4.29993pt\Delta}^{n-1}\left(t_{1},\dots,t_{n}\right)f.

(3)

$\mathord{\kern 4.29993pt\vrule width=0.6pt,height=5.6pt,depth=-0.28pt\kern-4.29993pt\Delta}$ may be computed by the recursive rule

\mathord{\kern 4.29993pt\vrule width=0.6pt,height=5.6pt,depth=-0.28pt\kern-4.29993pt\Delta}^{n-1}\left(x_{1},\dots,x_{n}\right)f=\begin{cases}\frac{\mathord{\kern 3.52356pt\vrule width=0.6pt,height=5.6pt,depth=-0.28pt\kern-3.52356pt\Delta}^{n-2}\left(x_{2},\dots,x_{n}\right)f-\mathord{\kern 3.52356pt\vrule width=0.6pt,height=5.6pt,depth=-0.28pt\kern-3.52356pt\Delta}^{n-2}\left(x_{1},\dots,x_{n-1}\right)f}{x_{n}-x_{1}}&x_{1}\neq x_{n},\\ \left\{\frac{\operatorname{d}}{\operatorname{d}\xi}\mathord{\kern 4.29993pt\vrule width=0.6pt,height=5.6pt,depth=-0.28pt\kern-4.29993pt\Delta}^{n-2}\left(\xi,x_{2},\dots,x_{n-1}\right)f\right\}|_{\xi=x_{n}},&x_{1}=x_{n},\end{cases}

(6.3)

where $\mathord{\kern 4.29993pt\vrule width=0.6pt,height=5.6pt,depth=-0.28pt\kern-4.29993pt\Delta}^{0}\left(x_{1}\right)f=f\left(x_{1}\right).$

(4)

Let $f_{z}\left(x\right)=\left(z-x\right)^{-1}$ . Then for all $z\notin\left\{x_{1},\dots,x_{n}\right\}$

\mathord{\kern 4.29993pt\vrule width=0.6pt,height=5.6pt,depth=-0.28pt\kern-4.29993pt\Delta}^{n-1}\left(x_{1},\dots,x_{n}\right)f_{z}=\frac{1}{q_{n,w}\left(z\right)}.

(6.4)

(5)

By (6.2), $\mathord{\kern 4.29993pt\vrule width=0.6pt,height=5.6pt,depth=-0.28pt\kern-4.29993pt\Delta}^{n-1}\left(x_{1},\dots,x_{n}\right)$ is a linear combination of the functionals

\delta^{\left(\ell\right)}\left(x-\tau_{j}\right),\qquad 1\leqslant j\leqslant s,\;0\leqslant\ell<d_{j}.

In fact, using (6.4) we obtain the Chakalov’s expansion (see [9])

\mathord{\kern 4.29993pt\vrule width=0.6pt,height=5.6pt,depth=-0.28pt\kern-4.29993pt\Delta}^{n-1}\left(x_{1},\dots,x_{n}\right)=\sum_{j=1}^{s}\sum_{\ell=0}^{d_{j}-1}a_{j,\ell}\delta^{\left(\ell\right)}\left(x-\tau_{j}\right),

(6.5)

where the coefficients $\left\{a_{j,\ell}\right\}$ are defined by the partial fraction decomposition¹¹1The coefficients $\left\{a_{j,\ell}\right\}$ may be readily obtained by the Cauchy residue formula $a_{j,\ell}=\frac{1}{\left(d_{j}-1-\ell\right)!}\lim_{z\to\tau_{j}}\left(\frac{\operatorname{d}}{\operatorname{d}z}\right)^{d_{j}-1-\ell}\left\{\frac{\left(z-\tau_{j}\right)^{\ell+1}}{q_{n,w}\left(z\right)}\right\}.$

\frac{1}{q_{n,w}\left(z\right)}=\sum_{j=1}^{s}\sum_{\ell=0}^{d_{j}-1}\frac{\ell!a_{j,\ell}}{\left(z-\tau_{j}\right)^{\ell+1}}.

(6.6)

(6)

By (6.5) and (6.6)

\mathord{\kern 4.29993pt\vrule width=0.6pt,height=5.6pt,depth=-0.28pt\kern-4.29993pt\Delta}^{n-1}\left(\underbrace{t,\dots,t}_{\times n}\right)=\frac{1}{\left(n-1\right)!}\delta^{\left(n-1\right)}\left(x-t\right).

(6.7)

(7)

Popoviciu’s refinement lemma [9, Proposition 23]: for every index subsequence

1\leqslant\sigma\left(1\right)<\sigma\left(2\right)<\dots<\sigma\left(k\right)\leqslant n,

there exist coefficients $\alpha\left(j\right)$ such that

\mathord{\kern 4.29993pt\vrule width=0.6pt,height=5.6pt,depth=-0.28pt\kern-4.29993pt\Delta}^{k-1}\left(x_{\sigma\left(1\right)},\dots,x_{\sigma\left(k\right)}\right)=\sum_{j=\sigma\left(1\right)-1}^{\sigma\left(k\right)-k}\alpha\left(j\right)\mathord{\kern 4.29993pt\vrule width=0.6pt,height=5.6pt,depth=-0.28pt\kern-4.29993pt\Delta}^{k-1}\left(x_{j+1},x_{j+2},\dots,x_{j+k}\right).

(6.8)

Based on the above, we may now identify $\mathord{\kern 4.29993pt\vrule width=0.6pt,height=5.6pt,depth=-0.28pt\kern-4.29993pt\Delta}$ with elements of the bundle ${\mathcal{}P}_{r}$ .

Definition 6.5.

Let $w=\left(x_{1},\dots,x_{r}\right)\in\mathbb{C}^{r}$ , and $X=\left\{n_{1},n_{2},\dots,n_{\alpha}\right\}\subseteq\left\{1,2,\dots,r\right\}$ of size $\left|X\right|=\alpha$ be given. Let the elements of $X$ be enumerated in increasing order, i.e.

1\leqslant n_{1}<n_{2}<\dots<n_{\alpha}\leqslant r.

Denote by $w_{X}$ the vector

w_{X}\stackrel{{\scriptstyle\text{def}}}{{=}}\left(x_{n_{1}},x_{n_{2}},\dots,x_{n_{\alpha}}\right)\in\mathbb{C}^{\alpha}.

Then we denote

\Delta_{X}\left(w\right)\stackrel{{\scriptstyle\text{def}}}{{=}}\mathord{\kern 4.29993pt\vrule width=0.6pt,height=5.6pt,depth=-0.28pt\kern-4.29993pt\Delta}^{\alpha-1}\left(w_{X}\right).

We immediately obtain the following result.

Lemma 6.6.

For all $w\in\mathbb{C}^{r}$ and $X\subseteq\left\{1,2,\dots,r\right\}$ , we have $\Delta_{X}\left(w\right)\in V_{w}$ . Moreover, letting $\alpha=\left|X\right|$ we have

{\mathcal{}SM}\left(\Delta_{X}\left(w\right)\right)=\mathord{\kern 4.29993pt\vrule width=0.6pt,height=5.6pt,depth=-0.28pt\kern-4.29993pt\Delta}^{\alpha-1}\left(w_{X}\right)\frac{1}{z-x}=\frac{1}{q_{\alpha,w_{X}}\left(z\right)}.

(6.9)

Finally, $\left(w,\Delta_{X}\left(w\right)\right)$ is a continuous section of ${\mathcal{}P}_{r}$ .

6.3. Constructing a basis

The following result is well-known, see e.g. [9, Proposition 35].

Theorem 6.7.

Denote $N_{j}=\left\{1,2,\dots,j\right\}$ for $j=1,2,\dots,r$ . Then for every $w\in\mathbb{C}^{r}$ , the collection

\left\{\Delta_{N_{j}}\left(w\right)\right\}_{j=1}^{r}

is a basis for $V_{w}$ .

There are various proofs of this statement. Below we show how to construct sets which do not necessarily remain basis for all $w\in\mathbb{C}^{r}$ , but only for $w$ in a small neighborhood of a given $w_{0}\in\mathbb{C}^{r}.$ Theorem 6.7 will then follow as a special case of this construction.

Informally, if two coordinates $x_{i}$ and $x_{j}$ can collide, then it is necessary to allow them to be glued by some element of the basis, i.e. we will need $\Delta_{X}\left(w\right)$ where $i,j\in X$ (in Theorem 6.7 all coordinates might be eventually glued into a single point because $w$ is unrestricted.) In order to make this statement formal, let us introduce a notion of configuration, which is essentially a partition of the set of indices.

Definition 6.8.

A configuration ${\mathcal{}C}$ is a partition of the set $N_{r}=\left\{1,2,\dots,r\right\}$ into $s=s\left({\mathcal{}C}\right)$ disjoint nonempty subsets

\sqcup_{i=1}^{s}X_{i}=N_{r},\qquad\left|X_{i}\right|=d_{i}>0.

The multiplicity vector of ${\mathcal{}C}$ is

T\left({\mathcal{}C}\right)=\left(d_{1},\dots,d_{s}\right).

Every configuration defines a continuous family of divided differences as follows.

Definition 6.9.

Let a configuration ${\mathcal{}C}=\left\{X_{j}\right\}_{j=1}^{s\left({\mathcal{}C}\right)}$ . Enumerate each $X_{j}$ in increasing order of its elements

X_{j}=\left\{n_{1}^{j}<n_{2}^{j}<\dots n_{d_{j}}^{j}\right\}

and denote for every $m=1,2,\dots,d_{j}$

X_{j,m}\stackrel{{\scriptstyle\text{def}}}{{=}}\left\{n_{k}^{j}:\;k=1,2,\dots,m\right\}.

For every $w\in\mathbb{C}^{r}$ , the collection ${\mathcal{}B_{C}}\left(w\right)\subset V_{w}$ is defined as follows:

{\mathcal{}B_{C}}\left(w\right)\stackrel{{\scriptstyle\text{def}}}{{=}}\left\{\Delta_{X_{j,m}}\left(w\right)\right\}_{j=1,\dots,s\left({\mathcal{}C}\right)}^{m=1,\dots,d_{j}}.

Now we formally define when a partition is “good” with respect to a point $w\in\mathbb{C}^{r}$ .

Definition 6.10.

The point $w=\left(x_{1},\dots,x_{r}\right)\in\mathbb{C}^{r}$ is subordinated to the configuration ${\mathcal{}C}=\left\{X_{j}\right\}_{j=1}^{s\left({\mathcal{}C}\right)}$ if whenever $x_{k}=x_{\ell}$ for a pair of indices $k\neq\ell$ , then necessarily $k,\ell\in X_{j}$ for some $X_{j}$ .

Now we are ready to formulate the main result of this section.

Theorem 6.11.

For a given $w_{0}\in\mathbb{C}^{r}$ and a configuration ${\mathcal{}C}$ , the collection ${\mathcal{}B_{C}}\left(w_{0}\right)$ is a basis for $V_{w_{0}}$ if and only if $w_{0}$ is subordinated to ${\mathcal{}C}$ . In this case, ${\mathcal{}B_{C}}\left(w\right)$ is a continuous family of bases for $V_{w}$ in a sufficiently small neighborhood of $w_{0}$ .

Let us first make a technical computation.

Lemma 6.12.

For a configuration ${\mathcal{}C}$ and a point $w\in\mathbb{C}^{r}$ , consider for every fixed $j=1,\dots,s\left({\mathcal{}C}\right)$ the set

S_{j}\stackrel{{\scriptstyle\text{def}}}{{=}}\left\{\Delta_{X_{j,m}}\left(w\right)\right\}_{m=1}^{d_{j}}.

(6.10)

(1)

Define for any pair of indices $1\leqslant k\leqslant\ell\leqslant d_{j}$ the index set

X_{j,k:\ell}\stackrel{{\scriptstyle\text{def}}}{{=}}\left\{n_{k}^{j}<n_{k+1}^{j}<\dots<n_{\ell}^{j}\right\}\subseteq X_{j}=X_{j,1:d_{j}}=X_{j,d_{j}}.

Then

\Delta_{X_{j,k:\ell}}\left(w\right)\in\operatorname{span}S_{j}.

(2)

For an arbitrary subset $Y\subseteq X_{j}$ (and not necessarily containing segments of consecutive indices), we also have

$\Delta_{Y}\left(w\right)\in\operatorname{span}S_{j}.$

Proof.

For clarity, we denote $y_{i}=x_{n_{i}^{j}}$ and $\left[k:\ell\right]=\Delta_{X_{j,k:\ell}}\left(w\right)$ . By (6.3) we have in all cases (including repeated nodes)

\left(y_{\ell}-y_{k}\right)\left[k:\ell\right]=\left[k+1:\ell\right]-\left[k:\ell-1\right].

(6.11)

The proof of the first statement is by backward induction on $n=\ell-k$ . We start from $n=d_{j}$ , and obviously $\left[1:d_{j}\right]\in S_{j}$ . In addition, by definition of $S_{j}$ we have $\left[1:m\right]\in S_{j}$ for all $m=1,\dots,d_{j}$ . Therefore, in order to obtain all $\left[k:\ell\right]$ with $\ell-k=n-1$ , we apply (6.11) several times as follows.

$\displaystyle\left[2:n\right]$	$\displaystyle=$	$\displaystyle\left(y_{n}-y_{1}\right)\left[1:n\right]+\left[1:n-1\right]$
$\displaystyle\left[3:n+1\right]$	$\displaystyle=$	$\displaystyle\left(y_{n+1}-y_{2}\right)\underleftrightarrow{\left[2:n+1\right]}+\underline{\left[2:n\right]}$
	$\displaystyle\dots$
$\displaystyle\left[d_{j}-n+2:d_{j}\right]$	$\displaystyle=$	$\displaystyle\left(y_{d_{j}}-y_{d_{j}-n+1}\right)\underleftrightarrow{\left[d_{j}-n+1:d_{j}\right]}+\underline{\left[d_{j}-n+1:d_{j}-1\right]}$

Here the symbol $\underline{\cdots}$ under a term means that the term is taken directly from the previous line, while $\underleftrightarrow{\cdots}$ indicates that the induction hypothesis is used. In the end, the left-hand side terms are shown to belong to $\operatorname{span}S_{j}$ .

In order to prove the second statement, we employ the first statement, (6.8) and Proposition 6.4, Item 1. ∎

Proof of Theorem 6.11.

In one direction, assume that $w_{0}=\left(x_{1},\dots,x_{r}\right)$ is subordinated to ${\mathcal{}C}$ . It is sufficient to show that every element of the standard basis (2.2) belongs to $\operatorname{span}\left\{{\mathcal{}B_{C}}\left(w_{0}\right)\right\}$ .

Let $\tau_{j}\in T\left(w_{0}\right)$ , let $d_{j}$ be the corresponding multiplicity, and let $Y_{j}\subseteq N_{r}$ denote the index set of size $d_{j}$

Y_{j}\stackrel{{\scriptstyle\text{def}}}{{=}}\left\{i:\quad x_{i}=\tau_{j}\right\}.

By the definition of subordination, there exists an element in the partition of ${\mathcal{}C}$ , say $X_{k}$ , for which $Y_{j}\subseteq X_{k}$ . By Lemma 6.12 we conclude that for all subsets $Z\subseteq Y_{j}$ ,

\Delta_{Z}\left(w_{0}\right)\in\operatorname{span}\left\{\Delta_{X_{k,m}}\left(w_{0}\right)\right\}_{m=1}^{\left|X_{k}\right|}\subseteq\operatorname{span}\left\{{\mathcal{}B_{C}}\left(w_{0}\right)\right\}.

By (6.7), $\Delta_{Z}\left(w_{0}\right)$ is nothing else but

\Delta_{Z}\left(w_{0}\right)=\mathord{\kern 4.29993pt\vrule width=0.6pt,height=5.6pt,depth=-0.28pt\kern-4.29993pt\Delta}^{\left|Z\right|-1}\left(\underbrace{\tau_{j},\dots,\tau_{j}}_{\times\left|Z\right|}\right)=\frac{1}{\left(\left|Z\right|-1\right)!}\delta^{\left(\left|Z\right|-1\right)}\left(x-\tau_{j}\right).

This completes the proof of the necessity. In the other direction, assume by contradiction that $x_{k}=x_{\ell}=\tau$ but nevertheless there exist two distinct elements of the partition ${\mathcal{}C}$ , say $X_{\alpha}$ and $X_{\beta}$ such that $k\in X_{\alpha}$ and $\ell\in X_{\beta}.$ Let the sets $\left\{S_{j}\right\}_{j=1}^{s\left({\mathcal{}C}\right)}$ be defined by (6.10). Again, by Lemma 6.12 and (6.7) we conclude that

\delta\left(x-\tau\right)\in\operatorname{span}S_{\alpha}\cap\operatorname{span}S_{\beta}.

But notice that ${\mathcal{}B_{C}}\left(w_{0}\right)=\bigcup_{j=1}^{s\left({\mathcal{}C}\right)}S_{j}$ and $\sum_{j=1}^{s}\left|S_{j}\right|=d$ , therefore by counting dimensions we conclude that

\dim\operatorname{span}\left\{{\mathcal{}B_{C}}\left(w_{0}\right)\right\}<d,

in contradiction to the assumption that ${\mathcal{}B_{C}}\left(w_{0}\right)$ is a basis.

Finally, one can evidently choose a sufficiently small neighborhood $U\subset\mathbb{C}^{r}$ of $w_{0}$ such that for all $w\in U$ , no new collisions are introduced, i.e. $w$ is still subordinated to ${\mathcal{}C}$ . The continuity argument (Lemma 6.6) finishes the proof.∎

Remark 6.13.

Another possible method of proof is to consider the algebra of elementary fractions in the Stieltjes space ${\mathcal{}S}_{r}$ , and use the correspondence (6.9).

As we mentioned, Theorem 6.7 follows as a corollary of Theorem 6.11 for the configuration ${\mathcal{}C}$ consisting of a single partition set $N_{r}$ .

6.4. Resolution of collision singularities

Let $\mu_{0}^{*}\in\Sigma_{r}^{*}\subset\mathbb{C}^{2r}$ be given, and let $\left(w_{0},g_{0}\right)\in{\mathcal{}P}_{r}$ be a solution to the (rank-restricted) Prony problem. The point $w_{0}$ is uniquely defined up to a permutation of the coordinates, so we just fix a particular permutation. Let $T\left(w_{0}\right)=\left(\tau_{1},\dots,\tau_{s}\right)$ .

Our goal is to solve the rank-restricted Prony problem for every input $\mu^{*}\in\mathbb{C}^{2r}$ in a small neighborhood of $\mu_{0}^{*}$ . According to Theorem 5.7, this amounts to a continuous representation of the solution $R_{\mu^{*}}\left(z\right)=\frac{P_{\mu^{*}}\left(z\right)}{Q_{\mu^{*}}\left(z\right)}={\mathcal{}TM}_{r}^{-1}\left(\mu^{*}\right)$ to the corresponding diagonal Padé approximation problem as an element of the bundle ${\mathcal{}P}_{r}$ .

Define $\delta=\min_{i\neq j}\left|\tau_{i}-\tau_{j}\right|$ to be the “separation distance” between the clusters. Since the roots of $Q_{\mu^{*}}$ depend continuously on $\mu^{*}$ and the degree of $Q_{\mu^{*}}$ does not drop, we can choose some $\mu_{1}^{*}$ sufficiently close to $\mu_{0}^{*}$ , for which

(1)

all the roots of $Q_{\mu_{1}^{*}}\left(z\right)$ are distinct, and
(2)

these roots can be grouped into $s$ clusters, such that each of the elements of the $j$ -th cluster is at most $\delta/3$ away from $\tau_{j}$ .

Enumerate the roots of $Q_{\mu_{1}^{*}}$ within each cluster in an arbitrary manner. This choice enables us to define locally (in a neighborhood of $\mu_{1}^{*}$ ) $r$ algebraic functions $x_{1}\left(\mu^{*}\right),\dots,x_{r}\left(\mu^{*}\right)$ , satisfying

Q_{\mu^{*}}\left(z\right)=\prod_{j=1}^{s}\left(z-x_{j}\left(\mu^{*}\right)\right).

Then we extend these functions by analytic continuation according to the above formula into the entire neighborhood of $\mu_{0}^{*}$ . Consequently,

w\left(\mu^{*}\right)\stackrel{{\scriptstyle\text{def}}}{{=}}\left(x_{1}\left(\mu^{*}\right),\dots,x_{r}\left(\mu^{*}\right)\right)

is a continuous (multivalued) algebraic function in a neighborhood of $\mu_{0}^{*}$ , satisfying

w\left(\mu_{0}^{*}\right)=w_{0}.

After this “pre-processing” step, we can solve the rank-restricted Prony problem in this neighborhood of $\mu_{0}^{*}$ , as follows.

The input to the problem is a measurement vector $\mu^{*}=\left(m_{0},\dots,m_{2r-1}\right)\in\mathbb{C}^{2r}$ , which is in a small neighborhood of $\mu_{0}^{*}$ .

(1)

Construct the function $w=w\left(\mu^{*}\right)$ as described above.
(2)

Build the basis ${\mathcal{}B_{C}}\left(w\right)=\left\{\Delta_{X_{j,\ell}}\left(w\right)\right\}_{j=1,\dots,s\left({\mathcal{}C}\right)}^{\ell=1,\dots,d_{j}}$ for $V_{w}$ .

(3)

Find the coefficients $\left\{\beta_{j,\ell}\right\}_{j=1,\dots,s\left({\mathcal{}C}\right)}^{\ell=1,\dots,d_{j}}$ such that

{\mathcal{}SM}\left(\sum_{j,\ell}\beta_{j,\ell}\Delta_{X_{j,\ell}}\left(w\right)\right)=R\left(z\right),

by solving the linear system

\underbrace{\sum_{j,\ell}\beta_{j,\ell}\left(w\right)\Delta_{X_{j,\ell}}\left(w\right)}_{=g\left(w\right)}\left(x^{k}\right)=m_{k}\left(=\int x^{k}g\left(w\right)\left(x\right)\operatorname{d}x\right),\qquad k=0,1,\dots,2r-1.

(6.12)

Algorithm 1 Solving rank-restricted Prony problem with collisions.

Theorem 6.14.

The coordinates $\left\{\beta_{j,\ell}\right\}$ of the solution to the rank-restricted Prony problem, given by Algorithm 1, are (multivalued) algebraic functions, continuous in a neighborhood of the point $\mu_{0}^{*}$ .

Proof.

Since the divided differences $\Delta_{j,\ell}\left(w\right)$ are continuous in $w$ , then clearly for each $k=0,1,\dots,2r-1$ the functions

\nu_{j,\ell,k}\left(w\right)=\Delta_{j,\ell}\left(w\right)\left(x^{k}\right)=\mathord{\kern 4.29993pt\vrule width=0.6pt,height=5.6pt,depth=-0.28pt\kern-4.29993pt\Delta}^{\ell-1}\left(w_{X_{j,\ell}}\right)\left(x^{k}\right)

are continuous²²2In fact, $\nu_{j,\ell,k}\left(w\right)$ are symmetric polynomials in some of the coordinates of $w$ . in $w$ , and hence continuous, as multivalued functions, in a neighborhood of $\mu_{0}^{*}$ . Since ${\mathcal{}B_{C}}\left(w\left(\mu^{*}\right)\right)$ remains a basis in a (possibly smaller) neighborhood of $\mu_{0}^{*}$ , the system (6.12), taking the form

\sum_{j,\ell}\nu_{j,\ell,k}\left(w\right)\beta_{j,\ell}\left(w\right)=m_{k},\qquad k=0,1,\dots,2r-1,

remains non-degenerate in this neighborhood. We conclude that the coefficients $\left\{\beta_{j,\ell}\left(w\left(\mu^{*}\right)\right)\right\}$ are multivalued algebraic functions, continuous in a neighborhood of $\mu_{0}^{*}$ . ∎

7. Real Prony space and hyperbolic polynomials

In this section we shall restrict ourselves to the real case. Notice that in many applications only real Prony systems are used. On the other hand, considering the Prony problem over the real numbers significantly simplifies some constructions. In particular, we can easily avoid topological problems, related with the choice of the ordering of the points $x_{1},\dots,x_{d}\in\mathbb{C}.$ So in a definition of the real Prony space $R{\mathcal{}P}_{d}$ we assume that the coordinates $x_{1},\ldots,x_{d}$ are taken with their natural ordering $x_{1}\leq x_{2}\leq\dots\leq x_{d}$ . Accordingly, the real Prony space $R{\mathcal{}P}_{d}$ is defined as the bundle $(w,g),\ w\in\prod_{d}\subset{\mathbb{R}}^{d},g\in RV_{w}.$ Here $\prod_{d}$ is the prism in ${\mathbb{R}}^{d}$ defined by the inequalities $x_{1}\leq x_{2}\leq\dots\leq x_{d}$ , and $RV_{w}$ is the space of linear combinations with real coefficients of $\delta$ -functions and their derivatives with the support $\{x_{1},\ldots,x_{d}\},$ as in Definition 2.4. The Prony, Stieltjes and Taylor maps are the restrictions to the real case of the complex maps defined above.

In this paper we just point out a remarkable connection of the real Prony space and mapping with hyperbolic polynomials, and Vieta and Vandermonde mappings studied in Singularity Theory (see [1, 13, 14, 15] and references therein).

Hyperbolic polynomials (in one variable) are real polynomials $Q(z)=z^{d}+\sum_{j=1}^{d}\lambda_{j}z^{d-j},$ with all $d$ of their roots real. We denote by $\Gamma_{d}$ the space of the coefficients $\Lambda=(\lambda_{1},\ldots,\lambda_{d})\subset{\mathbb{R}}^{d}$ of all the hyperbolic polynomials, and by $\hat{\Gamma}_{d}$ the set of $\Lambda\in\Gamma_{d}$ with $\lambda_{1}=0,\ |\lambda_{2}|\leq 1.$ Recalling (2.3), it is evident that all hyperbolic polynomials appear as the denominators of the irreducible fractions in the image of $R{\mathcal{}P}_{d}$ by ${\mathcal{}SM}$ . This shows, in particular, that the geometry of the boundary $\partial\Gamma$ of the hyperbolicity domain $\Gamma$ is important in the study of the real Prony map ${\mathcal{}PM}$ : it is mapped by ${\mathcal{}PM}$ to the boundary of the solvability domain of the real Prony problem. This geometry has been studied in a number of publications, from the middle of 1980s. In [13] V. P. Kostov has shown that $\hat{\Gamma}$ possesses the Whitney property: there is a constant $C$ such that any two points $\lambda_{1},\lambda_{2}\in\hat{\Gamma}$ can be connected by a curve inside $\hat{\Gamma}$ of the length at most $C\|\lambda_{2}-\lambda_{1}\|$ . “Vieta mapping” which associates to the nodes $x_{1}\leq x_{2}\leq\dots\leq x_{d}$ the coefficients of $Q(z)$ having these nodes as the roots, is also studied in [13]. In our notations, Vieta mapping is the composition of the Stieltjes mapping ${\mathcal{}SM}$ with the projection to the coefficients of the denominator.

In [1] V.I.Arnold introduced and studied the notion of maximal hyperbolic polynomial, relevant in description of $\hat{\Gamma}$ . Furthermore, the Vandermonde mapping ${\mathcal{}V}:{\mathbb{R}}^{d}\rightarrow{\mathbb{R}}^{d}$ was defined there by

\displaystyle\begin{cases}y_{1}=a_{1}x_{1}+\ldots+a_{d}x_{d},\\ \dots\\ y_{d}=a_{1}x_{1}^{d}+\ldots+a_{d}x_{d}^{d},\end{cases}

with $a_{1},\ldots,a_{d}$ fixed. In our notations ${\mathcal{}V}$ is the restriction of the Prony mapping to the pairs $(w,g)\in R{\mathcal{}P}_{d}$ with the coefficients of $g$ in the standard basis of $RV_{w}$ fixed. It was shown in [1] that for $a_{1},\ldots,a_{d}>0$ ${\mathcal{}V}$ is a one-to-one mapping of $\prod_{d}$ to its image. In other words, the first $d$ moments uniquely define the nodes $x_{1}\leq x_{2}\leq\dots\leq x_{d}$ . For $a_{1},\ldots,a_{d}$ with varying signs, this is no longer true in general. This result is applied in [1] to the study of the colliding configurations.

Next, the “Vandermonde varieties” are studied in [1], which are defined by the equations

\displaystyle\begin{cases}a_{1}x_{1}+\ldots+a_{d}x_{d}&=\alpha_{1},\\ &\dots\\ a_{1}x_{1}^{\ell}+\ldots+a_{d}x_{d}^{\ell}&=\alpha_{\ell}.\end{cases}

\displaystyle\ell\leqslant d.

It is shown that for $a_{1},\ldots,a_{d}>0$ the intersections of such varieties with $\prod_{d}$ are either contractible or empty. Finally, the critical points of the next Vandermonde equation on the Vandermond variety are studied in detail, and on this base a new proof of Kostov’s theorem is given.

We believe that the results of [1, 13] and their continuation in [14, 15] and other publications are important for the study of the Prony problem over the reals, and we plan to present some results in this direction separately.

Appendix A Proof of Theorem 3.5

Recall that we are interested in finding conditions for which the Taylor mapping ${\mathcal{}TM}:\;{\mathcal{}S}_{d}\to{\mathcal{}T}_{d}$ is invertible. In other words, given

S\left(z\right)=\sum_{k=0}^{2d-1}m_{k}\left(\frac{1}{z}\right)^{k+1},

we are looking for a rational function $R\left(z\right)\in{\mathcal{}S}_{d}$ such that

S\left(z\right)-R\left(z\right)=\frac{d_{1}}{z^{2d+1}}+\frac{d_{2}}{z^{2d+2}}+\dots.

(A.1)

Write $R\left(z\right)=\frac{P\left(z\right)}{Q\left(z\right)}$ with $Q\left(z\right)=\sum_{j=0}^{d}c_{j}z^{j}$ and $P\left(z\right)=\sum_{i=0}^{d-1}b_{i}z^{i}$ . Multiplying (A.1) by $Q\left(z\right)$ , we obtain

Q\left(z\right)S\left(z\right)-P\left(z\right)=\frac{e_{1}}{z^{d+1}}+\frac{e_{2}}{z^{d+2}}+\dots.

(A.2)

Proposition A.1.

The identity (A.2), considered as an equation on $P$ and $Q$ with $\deg P<\deg Q\leq d$ , always has a solution.

Proof.

Substituting the expressions for $S,\ P$ and $Q$ into (A.2) we get

\left(c_{0}+c_{1}z+\dots+c_{d}z^{d}\right)\left(\frac{m_{0}}{z}+\frac{m_{1}}{z^{2}}+\dots\right)-b_{0}-\dots-b_{{d-1}}z^{{d-1}}=\frac{e_{1}}{z^{{d+1}}}+\dots.

(A.3)

The highest degree of $z$ in the left hand side of (A.3) is $d-1$ . So equating to zero the coefficients of $z^{s}$ in (A.3) for $s=d-1,\dots,-d$ we get the following systems of equations:

\begin{bmatrix}0&0&0&m_{0}\\ 0&0&m_{0}&m_{1}\\ \adots&\adots\\ m_{0}&m_{1}&\dots&m_{d-1}\end{bmatrix}\begin{bmatrix}c_{1}\\ c_{2}\\ \vdots\\ c_{d}\end{bmatrix}=\begin{bmatrix}b_{d-1}\\ b_{d-2}\\ \vdots\\ b_{0}\end{bmatrix}.

(A.

\star

)

From this point on, the equations become homogeneous:

\begin{bmatrix}m_{0}&m_{1}&\dots&m_{d}\\ m_{1}&m_{2}&\dots&m_{d+1}\\ \adots&\adots\\ m_{d-1}&m_{d}&\dots&m_{2d-1}\end{bmatrix}\begin{bmatrix}c_{0}\\ c_{1}\\ \vdots\\ c_{d}\end{bmatrix}=\begin{bmatrix}0\\ 0\\ \vdots\\ 0\end{bmatrix}.

(A.

\star\star

)

The homogeneous system (A. $\star\star$ ) has the Hankel-type $d\times\left(d+1\right)$ matrix $\tilde{M}_{d}=\left(m_{i+j}\right)$ with $0\leqslant i\leqslant d-1$ and $0\leqslant j\leqslant d$ . This system has $d$ equations and $d+1$ unknowns $c_{0},\dots,c_{d}$ . Consequently, it always has a nonzero solution $c_{0},\dots,c_{d}$ . Now substituting these coefficients $c_{0},\dots,c_{d}$ of $Q$ into the equations (A. $\star$ ) we find the coefficients $b_{0},\dots,b_{d-1}$ of the polynomial $P$ , satisfying (A. $\star$ ). Notice that if $c_{j}=0$ for $j\geqslant\ell+1$ then it follows from the structure of the equations (A. $\star$ ) that $b_{j}=0$ for $j\geq\ell$ . Hence these $P,Q$ provide a solution of (A.2), satisfying $\deg P<\deg Q\leq d,$ and hence belonging to ${\mathcal{}S}_{d}.$ ∎

However, in general (A.2) does not imply (A.1). This implication holds only if $\deg Q=d$ . The following proposition describes a possible “lost of accuracy” as we return from (A.2) to (A.1) and $\deg Q<d$ :

Proposition A.2.

Let (A.2) be satisfied with the highest nonzero coefficient of $Q$ being $c_{\ell},\ \ell\leq d$ . Then

S(z)-\frac{{P(z)}}{{Q(z)}}=\frac{d_{1}}{z^{{d+\ell+1}}}+\frac{d_{2}}{z^{{d+\ell+2}}}+\dots.

(A.4)

Proof.

We notice that if the leading nonzero coefficient of $Q$ is $c_{\ell}$ then we have

\frac{1}{Q}=\frac{1}{{z^{\ell}}}(\frac{1}{{c_{\ell}+\frac{c_{\ell-1}}{z}+\dots}})=\frac{1}{{z^{\ell}}}(f_{0}+f_{1}\frac{1}{z}+\dots).

So multiplying (A.2) by $\frac{1}{Q}$ we get (A.4). ∎

Proof of Theorem 3.5.

Assume that the rank of $\tilde{M}_{d}$ is $r\leq d,$ and that $|M_{r}|\neq 0.$ Let us find a polynomial $Q(z)$ of degree $r$ of the form $Q(z)=z^{r}+\sum_{j=0}^{r-1}c_{j}z^{j},$ whose coefficients satisfy system (A. $\star\star$ ). Put $\mathbf{c}_{r}=(c_{0},\dots,c_{r-1},1)^{T}$ and consider a linear system $\tilde{M}_{r}\mathbf{c}_{r}=0$ . Since by assumptions $|M_{r}|\neq 0,$ this system has a unique solution. Extend this solution by zeroes, i.e. put $\mathbf{c}_{d}=(c_{0},\dots,c_{r-1},1,0,\dots,0)^{T}.$ We want $\mathbf{c}_{d}$ to satisfy (A. $\star\star$ ), which is $\tilde{M}_{d}\mathbf{c}_{d}=0$ . This fact is immediate for the first $r$ rows of $\tilde{M}_{d}$ . But since the rank of $\tilde{M}_{d}$ is $r$ by the assumption, its other rows are linear combinations of the first $r$ ones. Hence $\mathbf{c}_{d}$ satisfies (A. $\star\star$ ).

Now the equations (A. $\star$ ) produce a polynomial $P(z)$ of degree at most $r-1$ . So we get a rational function $R(z)=\frac{{P(z)}}{{Q(z)}}\in{\mathcal{}S}_{r}\subseteq{\mathcal{}S}_{d}$ which solves the Padé problem (A.2), with $\deg Q(z)=r$ . Write $R(z)=\sum_{k=0}^{\infty}\alpha_{k}(\frac{1}{z})^{k+1}$ . By Proposition A.2 we have $m_{k}=\alpha_{k}$ till $k=d+r-1$ .

Now, the Taylor coefficients $\alpha_{k}$ of $R(z)$ satisfy a linear recurrence relation

m_{k}=-\sum_{{s=1}}^{r}c_{s}m_{{k-s}},\qquad k=r,r+1,\dots.

(A.5)

Considering the rows of the system $\tilde{M}_{d}\mathbf{c}_{d}=0$ we see that $m_{k}$ satisfy the same recurrence relation (A.5) till $k=d+r-1$ (we already know that $m_{k}=\alpha_{k}$ till $k=d+r-1$ ). We shall show that in fact $m_{k}$ satisfy (A.5) till $k=2d-1.$

Consider a $d\times r$ matrix $\bar{M}_{d}$ formed by the first $r$ columns of $M_{d}$ , and denote its row vectors by $\mathbf{v}_{i}=(m_{i,0},\dots,m_{i,r-1}),\ i=1,\dots,d-1$ . The vectors $\mathbf{v}_{i}$ satisfy

\mathbf{v}_{i}=-\sum_{{s=1}}^{r}c_{s}\mathbf{v}_{{i-s}},\quad i=r,\dots,d-1,

(A.6)

since their coordinates satisfy (A.5) till $k=d+r-1$ . Now $\mathbf{v}_{0},\dots,\mathbf{v}_{r-1}$ are linearly independent, and hence each $\mathbf{v}_{i},\ i=r,\dots,d-1,$ can be expressed as

\mathbf{v}_{i}=\sum_{{s=0}}^{{r-1}}\gamma_{{i,s}}\mathbf{v}_{s}.

(A.7)

Denote by $\mathbf{\tilde{v}}_{i}=(m_{i,0},\dots,m_{i,d}),\ i=1,\dots,d-1$ the row vectors of $\tilde{M}_{d}$ . Since by assumptions the rank of $\tilde{M}_{d}$ is $r$ , the vectors $\mathbf{\tilde{v}}_{i}$ can be expressed through the first $r$ of them exactly in the same form as $\mathbf{v}_{i}$ :

\mathbf{\tilde{v}}_{i}=\sum_{{s=0}}^{{r-1}}\gamma_{{i,s}}\mathbf{\tilde{v}}_{s},\quad i=r,\dots,d-1.

(A.8)

Now the property of a system of vectors to satisfy the linear recurrence relation (A.6) depends only on the coefficients $\gamma_{i,s}$ in their representation (A.7) or (A.8). Hence from (A.6) we conclude that the full rows $\mathbf{\tilde{v}}_{i}$ of $\tilde{M}_{d}$ satisfy the same recurrence relation. Coordinate-wise this implies that $m_{k}$ satisfy (A.5) till $k=2d-1,$ and hence $m_{k}=\alpha_{k}$ till $k=2d-1.$ So $R(z)$ solves the original Problem 3.1.

In the opposite direction, assume that $R(z)$ solves Problem 3.1, and that the representation $R(z)=\frac{{P(z)}}{{Q(z)}}\in{\mathcal{}S}_{r}\subset{\mathcal{}S}_{d}$ is irreducible, i.e. $\deg Q=r$ . Write $Q(z)=z^{r}+\sum_{j=0}^{r-1}c_{j}z^{j}$ . Then $m_{k}$ , being the Taylor coefficients of $R(z)$ till $k=2d-1$ , satisfy a linear recurrence relation (A.5): $m_{k}=-\sum_{s=1}^{r}c_{s}m_{k-s},\ k=r,r+1,\dots,2d-1.$ Applying this relation coordinate-wise to the rows of $\tilde{M}_{d}$ we conclude that all the rows can be linearly expressed through the first $r$ ones. So the rank of $\tilde{M}_{d}$ is at most $r$ .

It remains to show that the left upper minor $|M_{r}|$ is non-zero, and hence the rank of $\tilde{M}_{d}$ is exactly $r$ .

By Proposition 3.3, if the decomposition of $R\left(z\right)$ in the standard basis is

R\left(z\right)=\sum_{j=1}^{s}\sum_{\ell=1}^{d_{j}}a_{j,\ell-1}\frac{\left(-1\right)^{\ell-1}\left(\ell-1\right)!}{\left(z-x_{j}\right)^{\ell}},

where $\sum_{j=1}^{s}d_{j}=r$ and $\left\{x_{j}\right\}$ are pairwise distinct, then the Taylor coefficients of $R\left(z\right)$ are given by (1.5). Clearly, we must have $a_{j,d_{j}-1}\neq 0$ for all $j=1,\dots,s$ , otherwise $\deg Q<r$ , a contradiction. Now consider the following well-known representation of $M_{r}$ as a product of three matrices (see e.g. [7]):

M_{r}=V\left(x_{1},d_{1},\dots,x_{s},d_{s}\right)\times\operatorname{diag}\left\{A_{j}\right\}_{j=1}^{s}\times V\left(x_{1},d_{1},\dots,x_{s},d_{s}\right)^{T},

(A.9)

where $V\left(\dots\right)$ is the confluent Vandermonde matrix (4.1) and each $A_{j}$ is the following $d_{j}\times d_{j}$ block:

A_{j}\stackrel{{\scriptstyle\text{def}}}{{=}}\begin{bmatrix}a_{j,0}&a_{j,1}&\cdots&\cdots&a_{j,d_{j}-1}\\ a_{j,1}&&&{d_{j}-1\choose d_{j}-2}a_{j,d_{j}-1}&0\\ \cdots&&&\cdots&0\\ &{d_{j}-1\choose 2}a_{j,d_{j}-1}&0&\cdots&0\\ a_{j,d_{j}-1}&0&\cdots&\cdots&0\end{bmatrix}.

The formula (A.9) can be checked by direct computation. Since $\left\{x_{j}\right\}$ are pairwise distinct and $a_{j,d_{j}-1}\neq 0$ for all $j=1,\dots,s$ , we immediately conclude that $\left|M_{r}\right|\neq 0$ .

This finishes the proof of Theorem 3.5. ∎

References

[1] V.I. Arnol’d. Hyperbolic polynomials and Vandermonde mappings. Functional Analysis and Its Applications, 20(2):125–127, 1986.
[2] J.R. Auton. Investigation of Procedures for Automatic Resonance Extraction from Noisy Transient Electromagnetics Data. Volume III. Translation of Prony’s Original Paper and Bibliography of Prony’s Method. Technical report, Effects Technology Inc., Santa Barbara, CA, 1981.
[3] R. Badeau, B. David, and G. Richard. Performance of ESPRIT for estimating mixtures of complex exponentials modulated by polynomials. Signal Processing, IEEE Transactions on, 56(2):492–504, 2008.
[4] G.A. Baker and P. Graves-Morris. Pade approximants. Part 1: Basic theory. Encyclopedia of Mathematics and its applications, 1981.
[5] D. Batenkov. Complete Algebraic Reconstruction of Piecewise-Smooth Functions from Fourier Data. arXiv preprint arXiv:1211.0680.
[6] D. Batenkov. On the norm of inverses of confluent Vandermonde matrices. arXiv preprint arXiv:1212.0172.
[7] D. Batenkov and Y. Yomdin. On the accuracy of solving confluent Prony systems. To appear in SIAM J.Appl.Math. Arxiv preprint arXiv:1106.1137.
[8] E. Candes and C. Fernandez-Granda. Towards a mathematical theory of super-resolution. To appear in Communications on Pure and Applied Mathematics, 2012.
[9] C. de Boor. Divided differences. Surv. Approx. Theory, 1:46–69, 2005.
[10] D.L. Donoho, M. Elad, and V.N. Temlyakov. Stable recovery of sparse overcomplete representations in the presence of noise. Information Theory, IEEE Transactions on, 52(1):6–18, 2006.
[11] B. Gustafsson, C. He, P. Milanfar, and M. Putinar. Reconstructing planar domains from their moments. Inverse Problems, 16(4):1053–1070, 2000.
[12] W. Kahan and R.J. Fateman. Symbolic computation of divided differences. ACM SIGSAM Bulletin, 33(2):7–28, 1999.
[13] V.P. Kostov. On the geometric properties of Vandermonde’s mapping and on the problem of moments. Proceedings of the Royal Society of Edinburgh: Section A Mathematics, 112(3-4):203–211, 1989.
[14] V.P. Kostov. Root arrangements of hyperbolic polynomial-like functions. Revista matemática complutense, 19(1):197–225, 2006.
[15] V.P. Kostov. On root arrangements for hyperbolic polynomial-like functions and their derivatives. Bulletin des sciences mathematiques, 131(5):477–492, 2007.
[16] L.M. Milne-Thompson. The calculus of finite differences. Macmillan, 1933.
[17] T. Peter, D. Potts, and M. Tasche. Nonlinear approximation by sums of exponentials and translates. SIAM Journal on Scientific Computing, 33(4):1920, 2011.
[18] D. Potts and M. Tasche. Parameter estimation for exponential sums by approximate Prony method. Signal Processing, 90(5):1631–1642, 2010.
[19] R. Prony. Essai experimental et analytique. J. Ec. Polytech.(Paris), 2:24–76, 1795.
[20] P. Stoica and N. Arye. MUSIC, maximum likelihood, and Cramer-Rao bound. IEEE Transactions on Acoustics, Speech and Signal Processing, 37(5):720–741, 1989.
[21] Y. Yomdin. Singularities in Algebraic Data Acquisition. In M. Manoel, M.C.R. Fuster, and C.T.C. Wall, editors, Real and Complex Singularities. Cambridge University Press, 2010.

	$\displaystyle\left\|\Delta a_{j,\ell}\right\|$	$\displaystyle\leq$	$\displaystyle\frac{2}{\ell!}\left(\frac{2}{\delta}\right)^{s_{0}+r_{0}}\left(\frac{1}{2}+\frac{s_{0}+r_{0}}{\delta}\right)^{d_{j}-\ell}\left(1+\frac{\left\|a_{j,\ell-1}\right\|}{\left\|a_{j,d_{j}-1}\right\|}\right)\varepsilon,$
	$\displaystyle\left\|\Delta\tau_{j}\right\|$	$\displaystyle\leq$	$\displaystyle\frac{2}{d_{j}!}\left(\frac{2}{\delta}\right)^{s_{0}+r_{0}}\frac{1}{\left\|a_{j,d_{j}-1}\right\|}\varepsilon,$

Geometry and Singularities of the Prony mapping

Abstract.

Key words and phrases:

2000 Mathematics Subject Classification:

1. Introduction

Definition 1.1.

Definition 1.2.

Problem 1.3 (Prony problem of order dd).

Problem 1.4 (Noisy Prony problem).

2. Prony, Stieltjes and Taylor Mappings

Definition 2.1.

Example 2.2.

Remark 2.3.

Definition 2.4.

Definition 2.5.

Definition 2.6.

Definition 2.7.

Proposition 2.8.

Definition 2.9.

Definition 2.10.

Definition 2.11.

3. Solvability of the Prony problem

3.1. General condition for solvability

Problem 3.1 (Diagonal Padé approximation problem).

Proposition 3.2.

Proof.

Proposition 3.3.

Proof.

Definition 3.4.

Theorem 3.5.

Theorem 3.6.

Corollary 3.7.

Example 3.8.

Proposition 3.9.

Proof.

Example 3.10.

3.2. Near-singular inversion

Theorem 3.11.

Proof.

Example 3.12.

4. Multiplicity-restricted Prony problem

Definition 4.1.

Proposition 4.2.

Definition 4.3.

Theorem 4.4 ([7]).

Proof outline.

Remark 4.5.

5. Rank-restricted Prony problem

Example 5.1.

Definition 5.2.

Definition 5.3.

Definition 5.4.

Remark 5.5.

Proposition 5.6.

Theorem 5.7.

Proof.

6. Collision singularities and bases of finite differences

6.1. Introduction

Example 6.1.

6.2. Divided finite differences

Definition 6.2.

Definition 6.3 ([9]).

Proposition 6.4.

Definition 6.5.

Lemma 6.6.

6.3. Constructing a basis

Theorem 6.7.

Definition 6.8.

Definition 6.9.

Definition 6.10.

Theorem 6.11.

Lemma 6.12.

Proof.

Proof of Theorem 6.11.

Remark 6.13.

6.4. Resolution of collision singularities

Theorem 6.14.

Proof.

7. Real Prony space and hyperbolic polynomials

Appendix A Proof of Theorem 3.5

Problem 1.3 (Prony problem of order $d$ ).