Prony analysis [Technical report]

Consider the following system of signals:

\[S = \{ (a_j, x_j) \}, \quad j = 1, \ldots , N\]

where \(a_j\) are signal amplitudes and \(x_j\) are signal positions. \(\{a_j\}\) and \(\{x_j\}\) may be represent by real or complex numbers.

Define the Prony map as a map \(\{ (a_j, x_j) \} \mapsto \{m_k\}\) of the following form: \[\label{eq:prony} m_k = \sum_{j=1}^N a_j x_j^k.\]

Numbers \(\{m_k\}\) are moments of the system, equation \eqref{eq:prony} is the Prony equation system.

There is a special case of the Prony system:

\[S_t = \{ (a_j, \mu_j) \}; \quad m_k = \sum_{j=1}^N a_j e^{-2\pi i \mu_j k \delta}\]

where \(\delta \in \mathbb{R}_+\) is a parameter. We say that it is the trigonometric Prony map. If \(|x_j| = 1\), we can reduce general Prony map to the trigonometric case by \(x_j = e^{-2\pi i \mu_j \delta}\).

In the practical application, we usually know only moments of a source system and we want to reconstruct original amplitudes and positions of signals. This means that we should build the inverse Prony map. Unfortunately, measurements of moments frequently include some noise. A feature of the inverse Prony map is the following: even small noise can produce a significant reconstruction error. Therefore, it is very important to understand geometry of the Prony map and the inverse Prony map.

In this technical report, we consider some examples of the direct and inverse Prony map and build corresponding plots. You can also view the source code on the R language for each plot.

See also: (Akinshin 2015), (Azaïs in press), (Batenkov 2014), (Batenkov 2014a), (Batenkov 2014b), (Batenkov 2014c), (Batenkov 2013), (Batenkov 2013a), (Candès 2014), (Candès 2013), (Demanet 2013), (Donoho 1992), (Donoho 1989), (Duval 2013), (Fernandez-Granda 2013), (Heckel 2014), (Levy 1981), (Liao 2014), (McCutchen 1967), (Minami 1985), (Moitra 2014), (Odendaal 1994), (Slepian 1978), (Yomdin 2010), (Yomdin 2005), (Demanet 2014), (Morgenshtern 2014).

Let’s consider the following system:

\[m_k = \sum_{j=1}^N a_j x_j^k, \quad k \in \{0, 1, \ldots, 2N-1\}.\]

Introduce variables \(\{u_j\}\):

\[\begin{cases} u_1 = x_1 + x_2 + \ldots + x_N, \\ u_2 = x_1x_2 + x_1x_3 + \ldots + x_{N-1}x_N, \\ \ldots \\ u_N = x_1x_2\ldots x_N. \end{cases}\]

From the Viet’s theorem, we have:

\[\label{eq3} m_{l} = m_{l-1}u_1 - m_{l-2}u_2 + \ldots + (-1)^{N-1}m_{l-N} u_N.\]

Thus, \(\{u_j\}\) are solution of the following system of equations:

\[\begin{pmatrix} m_{N-1} & -m_{N-2} & \ldots & (-1)^{N-1} m_0 \\ m_{N} & -m_{N-1} & \ldots & (-1)^{N-1} m_1 \\ \vdots & \vdots & \ddots & \vdots \\ m_{2N-2} & -m_{2N-3} & \ldots & (-1)^{N-1} m_{N-1} \end{pmatrix} \times \begin{pmatrix} u_1 \\ u_2 \\ \vdots \\ u_N \end{pmatrix} = \begin{pmatrix} m_N \\ m_{N+1} \\ \vdots \\ m_{2N-1} \\ \end{pmatrix}.\]

We can find \(\{x_j\}\) from the following equation:

\[x^N - u_1 x^{N-1} + \ldots + (-1)^N u_N = 0.\]

We can find \(\{a_j\}\) as a solution of the following system:

\[\begin{pmatrix} 1 & 1 & \ldots & 1 \\ x_1 & x_2 & \ldots & x_d \\ x_1^2 & x_2^2 & \ldots & x_N^2 \\ \vdots & \vdots & \ddots & \vdots \\ x_1^{N-1} & x_2^{N-1} & \ldots & x_N^{N-1} \\ \end{pmatrix} \times \begin{pmatrix} a_1 \\ a_2 \\ a_3\\ \vdots \\ a_N \end{pmatrix} = \begin{pmatrix} m_0 \\ m_1 \\ m_2 \\ \vdots \\ m_{N-1} \end{pmatrix}.\]

Consider the following system of signals: \(S = \{ a, b, \mu, \nu \}\). We assume that we know the following \(\{m_k\}\):

\[\label{eq:TProny2D-mk} m_k = U(k\delta) = a e^{-2\pi i \mu k \delta} + b e^{-2\pi i \nu k \delta} = a (e^{-2\pi i \mu \delta})^k + b (e^{-2\pi i \nu \delta})^k, \quad k = 0, 1, 2.\]

Introduce

\[x = e^{-2\pi i \mu \delta}, \quad y = e^{-2\pi i \nu \delta}.\]

Then the source system can be rewritten in the following form:

\[m_k = ax^k + by^k.\]

The solution can be obtained from the following formulas:

\[M_k = |m_k|, \quad \phi = -2\pi \mu \delta, \quad \theta = -2\pi \nu \delta, \quad \Delta = \phi - \theta,\]

\[\Delta = \arccos \Bigg( \dfrac{2 M_1^2 - M_0^2 - M_2^2}{2(M_0^2 - M_1^2)} \Bigg)\]

\[D = M_0^2 - 4\dfrac{(M_1^2 - M_0^2)^2}{M_2^2 + 3M_0^2 - 4M_1^2}\]

\[a = \dfrac{M_0 \pm \sqrt{D}}{2}, \quad b = \dfrac{M_0 \mp \sqrt{D}}{2}\]

\[\theta = -\arccos \Big( \dfrac {\Re m_1 (a\cos\Delta + b) + \Im m_1 (a\sin\Delta)} {(a\cos\Delta + b)^2 + (a\sin\Delta)^2} \Big), \quad \phi = \theta + \Delta.\]

In the one-dimensional case, we can write the following system for the first two moments:

\[\begin{cases} m_0 = a_0, \\ m_1 = a_0 x_0. \end{cases}\]

Below you can see the direct and inverse Prony map for this case.

Consider the two following sets of two-dimensional systems:

\[\begin{array}{l} S_U(\mu_U) = (a = \{1, 1\}, \mu = \{ -\mu_U/2, \mu_U/2 \}), \\ S_V(a_1, a_2, \mu_V) = (a = \{a_1, a_2\}, \mu = \{ -\mu_V/2, \mu_V/2 \}). \end{array}\]

Define the trigonometric Prony map for each system:

\[\begin{cases} m_{U,k}(\mu_V) = e^{-2\pi i (-\mu_U/2) k \delta} + e^{-2\pi i (\mu_U/2) k \delta}, \\ m_{V,k}(a_1, a_2, \mu_V) = a_1 e^{-2\pi i (-\mu_V/2) k \delta} + a_2 e^{-2\pi i (\mu_V/2) k \delta}. \end{cases}\]

Define a scalar field in the \((a_1, a_2, \mu_V)\) space:

\[F(a_1, a_2, \mu_V) = \sqrt{\Delta m_0^2 + \Delta m_1^2 + \Delta m_2^2 + \Delta m_3^2}, \quad \Delta m_k = m_{U,k} - m_{V,k}.\]

You can see equipotential surfaces of the field on the following figure (\(\delta = 0.2\), \(\mu_U \in [0.1, 5]\))

General solving approach for general complex Prony problems doesn’t work well for noisy system. Consider the following numerical example: we have two-dimensional trigonometry Prony system (\(S_1\)). Let’s add some noise to the moments and try to reconstruct original system by the general approach (\(S_2\)) and the trigonometric approach (\(S_3\)). The numerical results are in the table \ref{tbl:noisy}.

\label{tbl:noisy}

\(x_1\) | \(x_2\) | \(a_1\) | \(a_2\) | \(m_0\) | \(m_1\) | \(m_2\) | \(m_3\) | |
---|---|---|---|---|---|---|---|---|

\(S_1\) | 0.891-0.454i | 0.951-0.309i | 1+0i | 1+0i | 2+0i | 1.84-0.76i | 1.4-1.4i | 0.74-1.8i |

\(S_2\) | 0.919-0.385i | 1.93-6.15i | 2.04+0.09i | 0.0005-0.0002i | 2.04+0.09i | 1.91-0.71i | 1.46-1.39i | 0.76-1.7i |

\(S_3\) | 0.953-0.304i | 0.916-0.402i | 1.02+0i | 1.02+0i | 2.04+0.09i | 1.91-0.71i | 1.46-1.39i | 0.76-1.7i |

\(S_2-S_1\) | 0.0277+0.0691i | 0.98-5.84i | 1.04+0.09i | -0.999-0i | 0.0425+0.0877i | 0.067+0.0542i | 0.0624+0.0105i | 0.019+0.0951i |

\(S_3-S_1\) | 0.062+0.15i | -0.0355-0.0932i | 0.0222+0i | 0.0222+0i | 0.0425+0.0877i | 0.067+0.0542i | 0.0624+0.0105i | 0.019+0.0951i |

We can see that the general approach gives a big error for \(x_2\) and \(a_2\). The trigonometric approach solves the noisy problem well.

A. Akinshin, D. Batenkov, Y. Yomdin. Accuracy of spike-train Fourier reconstruction for colliding nodes.

*arXiv:1502.06932 [math.CA]*(2015). LinkJean-Marc Azaïs, Yohann de Castro, Fabrice Gamboa. Spike detection from inaccurate samplings.

*Applied and Computational Harmonic Analysis*(in press). LinkDmitry Batenkov. Numerical stability bounds for algebraic systems of Prony type and their accurate solution by decimation.

*arXiv preprint arXiv:1409.3137*(2014). LinkDmitry Batenkov. Accurate solution of near-colliding Prony systems via decimation and homotopy continuation.

*arXiv:1501.00160 [cs, math]*(2014). LinkD. Batenkov, N. Sarig, Y. Yomdin. Accuracy of algebraic Fourier reconstruction for shifts of several signals.

*Sampling Theory in Signal and Image Processing***13**, 151–173 (2014).D. Batenkov, Y. Yomdin. Geometry and Singularities of the Prony mapping.

*Journal of Singularities***10**, 1–25 (2014).D. Batenkov, Y. Yomdin. Algebraic signal sampling, Gibbs phenomenon and Prony-type systems.. In

*Proceedings of the 10th International Conference on Sampling Theory and Applications (SAMPTA)*. (2013).D. Batenkov, Y. Yomdin. On the accuracy of solving confluent Prony systems.

*SIAM J.Appl.Math.***73**, 134–154 (2013).Emmanuel J. Candès, Carlos Fernandez-Granda. Towards a Mathematical Theory of Super-resolution.

*Communications on Pure and Applied Mathematics***67**, 906–956 (2014). LinkEmmanuel J. Candès, Carlos Fernandez-Granda. Super-Resolution from Noisy Data.

*Journal of Fourier Analysis and Applications***19**, 1229–1254 (2013). LinkLaurent Demanet, Deanna Needell, Nam Nguyen. Super-resolution via superset selection and pruning. In

*Proceedings of the 10th International Conference on Sampling Theory and Applications (SAMPTA)*. (2013).L. Demanet, N. Nguyen. The recoverability limit for superresolution via sparsity.

*Preprint*(2014).D.L. Donoho. Superresolution via sparsity constraints.

*SIAM Journal on Mathematical Analysis***23**, 1309–1331 (1992).D. L. Donoho, P. B. Stark. Uncertainty principles and signal recovery.

*SIAM J. Appl.Math.***49**, 906–931 (1989).Vincent Duval, Gabriel Peyré. Exact support recovery for sparse spikes deconvolution.

*arXiv preprint arXiv:1306.6909*(2013). LinkCarlos Fernandez-Granda. Support detection in super-resolution. 145–148 In

*Proc. of 10th Sampling Theory and Applications (SAMPTA)*. (2013). LinkReinhard Heckel, Veniamin I. Morgenshtern, Mahdi Soltanolkotabi. Super-Resolution Radar.

*arXiv:1411.6272 [cs, math]*(2014). LinkS. Levy, P. K. Fullagar. Reconstruction of a sparse spike train from a portion of its spectrum and application to high-resolution deconvolution.

*Geophysics***46**, 1235–1243 (1981).Wenjing Liao, Albert Fannjiang. MUSIC for Single-Snapshot Spectral Estimation: Stability and Super-resolution.

*arXiv:1404.1484 [cs, math]*(2014). LinkC. W. McCutchen. Superresolution in microscopy and the Abbe resolution limit.

*J. Opt. Soc. Am.***57**, 1190–1190 (1967).K. Minami, S. Kawata, S. Minami. Superresolution of Fourier Transform spectra by autoregressive model fitting with singular value decomposition.

*Appl. Optics***24**, 162–167 (1985).Ankur Moitra. The Threshold for Super-resolution via Extremal Functions.

*arXiv:1408.1681 [cs, math, stat]*(2014). LinkV. I. Morgenshtern, E. J. Candes. Stable super-resolution of positive sources: the discrete setup.

*Preprint*(2014).J. Odendaal, E. Barnard, C. W. I. Pistorius. Two-dimensional superresolution radar imaging using the MUSIC algorithm.

*IEEE Transactions on Antennas and Propagation***42**, 1386–1391 (1994).D. Slepian. Prolate spheroidal wave functions.

*Fourier Analysis and uncertainty, V. - The discrete case. Bell System Technical Journal***57**, 1371–1430 (1978).Y. Yomdin. Singularities in algebraic data acquisition.

*Real and complex singularities, London Math. Soc. Lecture Note Ser.***380**, 378–396 (2010).Y. Yomdin. Some quantitative results in singularity theory.

*Ann. Polon. Math.***87**, 277–-299 (2005).