Prony analysis [Technical report]

Introduction

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).

Solution approaches

General complex Prony, \(N\)D

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}.\]

Trigonometric real Prony, 2D

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.\]

Numerical examples

General real Prony, 1D

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.

Prony-1D map

Prony-1D inverse map

Trigonometric real Prony, 2D

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]\))

Prony-2D inverse map (contour3d)

In fact, here you can see the inverse Prony map. Each surface corresponds to a sphere in the four-dimensional moment space with given radius (\(r=\{1.5, 3.5, 5.5, 7.5\}\) for red, green, blue, and purple surfaces respectively) and center (\((1, 1, \mu_U)\)). This means that if \(S_U\) and \(S_V\) belong to inner region of a surface of radius \(r\), the corresponded distance in the moment space will less that \(2r\).

Trigonometric noisy real Prony, 2D

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.

Prony-2D Noisy problem

Real Prony, Geometry (WIP)

The hyperbola \( S_F^x \).

References

  1. A. Akinshin, D. Batenkov, Y. Yomdin. Accuracy of spike-train Fourier reconstruction for colliding nodes. arXiv:1502.06932 [math.CA] (2015). Link

  2. Jean-Marc Azaïs, Yohann de Castro, Fabrice Gamboa. Spike detection from inaccurate samplings. Applied and Computational Harmonic Analysis (in press). Link

  3. Dmitry Batenkov. Numerical stability bounds for algebraic systems of Prony type and their accurate solution by decimation. arXiv preprint arXiv:1409.3137 (2014). Link

  4. Dmitry Batenkov. Accurate solution of near-colliding Prony systems via decimation and homotopy continuation. arXiv:1501.00160 [cs, math] (2014). Link

  5. D. 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).

  6. D. Batenkov, Y. Yomdin. Geometry and Singularities of the Prony mapping. Journal of Singularities 10, 1–25 (2014).

  7. 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).

  8. D. Batenkov, Y. Yomdin. On the accuracy of solving confluent Prony systems. SIAM J.Appl.Math. 73, 134–154 (2013).

  9. Emmanuel J. Candès, Carlos Fernandez-Granda. Towards a Mathematical Theory of Super-resolution. Communications on Pure and Applied Mathematics 67, 906–956 (2014). Link

  10. Emmanuel J. Candès, Carlos Fernandez-Granda. Super-Resolution from Noisy Data. Journal of Fourier Analysis and Applications 19, 1229–1254 (2013). Link

  11. Laurent 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).

  12. L. Demanet, N. Nguyen. The recoverability limit for superresolution via sparsity. Preprint (2014).

  13. D.L. Donoho. Superresolution via sparsity constraints. SIAM Journal on Mathematical Analysis 23, 1309–1331 (1992).

  14. D. L. Donoho, P. B. Stark. Uncertainty principles and signal recovery. SIAM J. Appl.Math. 49, 906–931 (1989).

  15. Vincent Duval, Gabriel Peyré. Exact support recovery for sparse spikes deconvolution. arXiv preprint arXiv:1306.6909 (2013). Link

  16. Carlos Fernandez-Granda. Support detection in super-resolution. 145–148 In Proc. of 10th Sampling Theory and Applications (SAMPTA). (2013). Link

  17. Reinhard Heckel, Veniamin I. Morgenshtern, Mahdi Soltanolkotabi. Super-Resolution Radar. arXiv:1411.6272 [cs, math] (2014). Link

  18. S. 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).

  19. Wenjing Liao, Albert Fannjiang. MUSIC for Single-Snapshot Spectral Estimation: Stability and Super-resolution. arXiv:1404.1484 [cs, math] (2014). Link

  20. C. W. McCutchen. Superresolution in microscopy and the Abbe resolution limit. J. Opt. Soc. Am. 57, 1190–1190 (1967).

  21. 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).

  22. Ankur Moitra. The Threshold for Super-resolution via Extremal Functions. arXiv:1408.1681 [cs, math, stat] (2014). Link

  23. V. I. Morgenshtern, E. J. Candes. Stable super-resolution of positive sources: the discrete setup. Preprint (2014).

  24. 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).

  25. D. Slepian. Prolate spheroidal wave functions. Fourier Analysis and uncertainty, V. - The discrete case. Bell System Technical Journal 57, 1371–1430 (1978).

  26. Y. Yomdin. Singularities in algebraic data acquisition. Real and complex singularities, London Math. Soc. Lecture Note Ser. 380, 378–396 (2010).

  27. Y. Yomdin. Some quantitative results in singularity theory. Ann. Polon. Math. 87, 277–-299 (2005).

[Someone else is editing this]

You are editing this file