 # Prony analysis [Technical report]

•  Andrey Akinshin
•  Awaiting Activation

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

# 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.  ## 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]$$) 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. ## Real Prony, Geometry (WIP)           The hyperbola $$S_F^x$$.