Structured low rank decomposition of multivariate Hankel matrices

TL;DR

This paper introduces a new algorithm for decomposing multivariate Hankel matrices into low-rank components, leveraging algebraic structures and eigenvector methods, with improved numerical stability and error correction techniques.

Contribution

It presents a novel multivariate generalization of the Pencil method using algebraic properties of Artinian Gorenstein algebras for low-rank Hankel matrix decomposition.

Findings

The algorithm effectively decomposes multivariate Hankel matrices into polynomial-exponential series.

The method is robust to noise and includes a rescaling technique for high-amplitude frequencies.

A new Newton iteration improves accuracy by correcting errors in input moments.

Abstract

We study the decomposition of a multivariate Hankel matrix H\_$\sigma$ as a sum of Hankel matrices of small rank in correlation with the decomposition of its symbol $\sigma$ as a sum of polynomial-exponential series. We present a new algorithm to compute the low rank decomposition of the Hankel operator and the decomposition of its symbol exploiting the properties of the associated Artinian Gorenstein quotient algebra A\_$\sigma$. A basis of A\_$\sigma$ is computed from the Singular Value Decomposition of a sub-matrix of the Hankel matrix H\_$\sigma$. The frequencies and the weights are deduced from the generalized eigenvectors of pencils of shifted sub-matrices of H $\sigma$. Explicit formula for the weights in terms of the eigenvectors avoid us to solve a Vandermonde system. This new method is a multivariate generalization of the so-called Pencil method for solving Prony-type decomposition problems. We analyse its numerical behaviour in the presence of noisy input moments, and describe a rescaling technique which improves the numerical quality of the reconstruction for frequencies of high amplitudes. We also present a new Newton iteration, which converges locally to the closest multivariate Hankel matrix of low rank and show its impact for correcting errors on input moments.

Figures

Equations

$$H=[{\sigma }_{i+j}{]}_{0\le i\le l,0\le j\le m}$$

$$H_{\sigma }:L_0({\mathbb{C}}^{\mathbb{N}})$$

$$(p_k{)}_k$$

$$H=[{\sigma }_{\alpha +\beta }{]}_{\alpha \in A,\beta \in B}$$

$$H_{\sigma }:L_0({\mathbb{C}}^{{\mathbb{N}}^n})$$

$$(p_{\alpha }{)}_{\alpha }$$

$$\sigma (\mathbf{y})=\sum _{\alpha \in {\mathbb{N}}^n}{\sigma }_{\alpha }\frac{{\mathbf{y}}^{\alpha }}{\alpha !}\in \mathbb{C}[[\mathbf{y}]]$$

$$\begin{array}[]{rcl}\sigma:\mathbb{C}[\mathbf{x}]&\rightarrow&\mathbb{C}\\ p=\sum_{\alpha\in A\subset{\mathbb{N}}^{n}}p_{\alpha}\mathbf{x}^{\alpha}&\mapsto&\langle\sigma\mid p\rangle=\sum_{\alpha\in A\subset{\mathbb{N}}^{n}}p_{\alpha}\sigma_{\alpha}.\end{array}$$

$$\forall \sigma \in \mathbb{C}[[\mathbf{y}]],\forall p,q\in \mathbb{C}[\mathbf{x}],⟨p\star \sigma \mid q⟩=⟨\sigma \mid pq⟩.$$

$$p\star \sigma =\sum _{\beta \in {\mathbb{N}}^n}(\sum _{\alpha \in {\mathbb{N}}^n}{p}_{\alpha }{\sigma }_{\alpha +\beta })\frac{{\mathbf{y}}^{\alpha }}{\alpha !}.$$

$$H_{\sigma }:\mathbb{C}[\mathbf{x}]$$

$$p$$

$$\begin{array}[]{rcl}H^{A,B}_{\sigma}:\langle\mathbf{x}^{B}\rangle&\rightarrow&\langle\mathbf{y}^{A}\rangle\\ p=\sum_{\beta\in B}p_{\beta}\mathbf{x}^{\beta}&\mapsto&\sum_{\alpha\in A}(\sum_{\beta\in B}p_{\alpha}\sigma_{\alpha+\beta})\frac{\mathbf{y}^{\alpha}}{\alpha!}=p\star\sigma_{|\langle\mathbf{x}^{A}\rangle}\end{array}$$

$$H_{\sigma }^{A,B}=[{\sigma }_{\alpha +\beta }{]}_{\alpha \in A,\beta \in B}.$$

$$H_{\sigma}^{A,B}=\left[\begin{array}[]{cccc}1&2&3&4\\ 2&4&5&7\\ 3&5&6&8\\ \end{array}\right]$$

$$H_{\sigma }:(p_k)\in L_0({\mathbb{C}}^{\mathbb{N}})\mapsto (\sum _k{p}_k{\sigma }_{k+l}{)}_{l\in \mathbb{N}}\in {\mathbb{C}}^{\mathbb{N}}$$

$${\sigma }_n=\sum _{i=1}^{r^{\prime }}{\omega }_i(n){\xi }_i^n$$

$$\sigma (y)=\sum _{n\in \mathbb{N}}{\sigma }_n\frac{y^n}{n!}=\sum _{i=1}^{r^{\prime }}{\tilde{\omega }}_i(y)e^{{\xi }_{i}y}$$

$$\forall p\in R,⟨{\mathbf{e}}_{\xi }|p⟩=\sum _{\alpha \in {\mathbb{N}}^n}{p}_{\alpha }{\xi }^{\alpha }=p(\xi ).$$

$$H_{y}^{A,A}=\left[\begin{array}[]{ccc}0&1&0\\ 1&0&0\\ 0&0&0\\ \end{array}\right]\neq\lambda_{1}\,\left[\begin{array}[]{ccc}1&\xi_{1}&\xi_{1}^{2}\\ \xi_{1}&\xi_{1}^{2}&\xi_{1}^{3}\\ \xi_{1}^{2}&\xi_{1}^{3}&\xi_{1}^{4}\\ \end{array}\right]+\lambda_{2}\,\left[\begin{array}[]{ccc}1&\xi_{2}&\xi_{2}^{2}\\ \xi_{2}&\xi_{2}^{2}&\xi_{2}^{3}\\ \xi_{2}^{2}&\xi_{2}^{3}&\xi_{2}^{4}\\ \end{array}\right]$$

$$\omega {\mathbf{e}}_{\xi }=\sum _{i=1}^{r_1}{\omega }_{1,i}{\mathbf{e}}_{{\xi }_{1,i}}+\sum _{i=1}^{r_2}{\omega }_{2,i}{\mathbf{e}}_{{\xi }_{2,i}}.$$

$$\mathcal{A}=\mathbb{C}[\mathbf{x}]/I={\mathcal{A}}_1\oplus \cdots \oplus {\mathcal{A}}_{r^{\prime }}$$

$${\mathbf{u}}_{{\xi }_i}^2(\mathbf{x})\equiv {\mathbf{u}}_{{\xi }_i}(\mathbf{x}),\sum _{i=1}^r{\mathbf{u}}_{{\xi }_i}(\mathbf{x})\equiv 1.$$

$$\begin{array}[]{cccl}\mathcal{M}_{g}:&\mathcal{A}&\to&\mathcal{A}\\ &h&\mapsto&\mathcal{M}_{g}(h)=g\,h.\end{array}$$

$$\begin{array}[]{cccl}{\mathcal{M}}_{g}^{\scriptscriptstyle\mathsf{T}}:&\mathcal{A}^{\ast}&\to&\mathcal{A}^{\ast}\\ &\Lambda&\mapsto&{\mathcal{M}}_{g}^{\scriptscriptstyle\mathsf{T}}(\Lambda)=\Lambda\circ{\mathcal{M}}_{g}=g\star\Lambda.\end{array}$$

$$M_g^{\mathsf{T}}B({\xi }_i)=g({\xi }_i)B({\xi }_i).$$

$$\sum _{i=1}^r⟨\sigma |b_i{b}_j^{\prime }⟩{\lambda }_i=0.$$

$$H_{g\star \sigma }^{B,B^{\prime }}=(M_g^B{)}^{\mathsf{T}}{H}_{\sigma }^{B,B^{\prime }}=H_{\sigma }^{B,B^{\prime }}{M}_g^{B^{\prime }}.$$

$$(H_{g\star \sigma }^{B,B^{\prime }}{)}_{[i,j]}=⟨\sigma \mid g{b}_i{b}_j^{\prime }⟩$$

$$(H_{g\star \sigma }^{B,B^{\prime }}{)}_{[i,j]}=⟨\sigma \mid b_{i}g{b}_j^{\prime }⟩$$

$$\begin{array}[]{ccccc}H_{\sigma}^{d_{1},d_{2}}&:&{\mathbb{C}[\mathbf{x}]}_{d_{2}}&\to&({\mathbb{C}[\mathbf{x}]}_{d_{1}})^{*}\\ &&p&\mapsto&p\star\sigma\end{array}$$

$$H_{\sigma }^{d_1,d_2}=US{V}^{\mathsf{T}}$$

Full text

AMS classification: 14Q20, 68W30, 47B35, 15B05

Keywords: Hankel; polynomial; exponential series; low rank decomposition; eigenvector; Singular Value Decomposition.

1 Introduction

Structured matrices such as Toeplitz or Hankel matrices appear in many problems. They are naturally associated to operations on polynomials or series [10]. The correlation with polynomial algebra can be exploited to accelerate matrix computations [8]. The associated algebraic model provides useful information on the problem to be solved or the phenomena to be analysed. Understanding its structure often yields a better insight on the problem and its solution. In many cases, an efficient way to analyze the structure of the underlying models is to decompose the structured matrix into a sum of low rank matrices of the same structure. This low rank decomposition has applications in many domains [17] and appears under different formulations [16, 4, 3].

In this paper, we study specifically the class of Hankel matrices. We investigate the problem of decomposing a Hankel matrix as a sum of indecomposable Hankel matrices of low rank. Natural questions arise. What are the indecomposable Hankel matrices? Are they necessarily of rank $1$ ? How to compute a decomposition of a Hankel matrix as a sum of indecomposable Hankel matrices? Is the structured low rank decomposition of a Hankel matrix unique ?

These questions have simple answers for non-structured or dense matrices: The indecomposable dense matrices are the matrices of rank one, which are the tensor product of two vectors. The Singular Value Decomposition of a dense matrix yields a decomposition as a minimal sum of rank one matrices, but this decomposition is not unique.

It turns out that for the Hankel structure, the answers to these questions are not so direct and involve the analysis of the so-called symbol associated to the Hankel matrix. The symbol is a formal power series defined from the coefficients of the Hankel matrix. As we will see, the structured decomposition of an Hankel matrix is closely related to the decomposition of the symbol as a sum of polynomial-exponential series.

The decomposition of the symbol of a Hankel matrix is a problem, which has a long history. The first work on this problem is probably due to Gaspard-Clair-François-Marie Riche de Prony [2]. He proposed a method to reconstruct a sum of exponentials from the values at equally spaced data points, by computing a polynomial in the kernel of a Hankel matrix, and deducing the decomposition from the roots of this polynomial. Since then, many works have been developed to address the decomposition problem in the univariate case, using linear algebra tools on Hankel matrices such as Pencil method [13], ESPRIT method [26] or MUSIC method [28]. Other methods such as [12], approximate Prony Method [7], [24] use minimization techniques, to recover the frequencies or the weights in the sum of exponential functions. See [23][chap. 1] for a survey on some of these approaches.

The numerical behavior of these methods has also been studied. The condition number of univariate Hankel matrices, which decomposition involves real points has been investigated in [29], [5]. It is shown that it grows exponentially with the dimension of the matrix. The condition number of Vandermonde matrices of general complex points has been studied recently in [21]. In [6], the numerical sensitivity of the generalized eigenvalues of pencils of Hankel matrices appearing in Prony’s method has been analysed.

The development of multivariate decomposition methods is more recent. Extension of the univariate approaches have been considered e.g. in [1], [25], [22]. These methods project the problem in one dimension and solve several univariate decomposition problems to recover the multivariate decomposition by least square minimization from a grid of frequencies. In [22], [14], the decomposition problem is transformed into the solution of an overdetermined polynomial system associated to the kernel of these Hankel matrices, and the frequencies of the exponential terms are found by general polynomial solver. These methods involves Hankel matrices of size exponential in the number of variables of the problem or moments of order at least twice the number of terms in the decomposition. In [27], an $H$-basis of the ideal defining the frequencies is computed from Hankel matrices built from moments of big enough order. Tables of multiplications are deduced from the $H$-basis and their eigenvalues yield the frequencies of the exponential terms. The weights are computed as the solution of a Vandermonde linear system. Moments of order bigger than twice the degree of an H-basis are involved in this construction.

Contributions

We study the decomposition of a multivariate Hankel matrix as a sum of Hankel matrices of small rank in correlation with the decomposition of its symbol $\sigma$ as a sum of polynomial-exponential series. We show how to recover efficiently this decomposition from the structure of the quotient algebra ${\mathcal{A}}_{\sigma}$ of polynomials modulo the kernel of the corresponding Hankel operator $H_{\sigma}$. In particular, a basis of ${\mathcal{A}}_{\sigma}$ can be extracted from any maximal non-zero minor of the matrix of $H_{\sigma}$. We also show how to compute the matrix of multiplication by a variable in the basis of ${\mathcal{A}}_{\sigma}$ from sub-matrices of the matrix of $H_{\sigma}$. We describe how the frequencies of the polynomial-exponential decomposition of the symbol can be deduced from the eigenvectors of these matrices. Exploiting properties of these multiplication operators, we show that the weights of the decomposition can be recovered directly from these eigenvectors, avoiding the solution of a Vandermonde system. We present a new algorithm to compute the low rank decomposition of $H_{\sigma}$ and the decomposition of its symbol as a sum of polynomial-exponential series from sub-matrices of the matrix of $H_{\sigma}$. A basis of ${\mathcal{A}}_{\sigma}$ is computed from the Singular Value Decomposition of a sub-matrix. The frequencies and the weights are deduced from the generalized eigenvectors of pencils of sub-matrices of $H_{\sigma}$. This new method is a multivariate generalization of the so-called Pencil method for solving Prony-type decomposition problems. It can be used to decompose series as sums of polynomial-exponential functions from moments and provides structured low rank decomposition of multivariate Hankel matrices. We analyse its numerical behaviour in the presence of noisy input moments, for different numbers of variables, of exponential terms of the symbol and different amplitudes of the frequencies. We present a rescaling technique, which improves the numerical quality of the reconstruction for frequencies of high amplitudes. We also present a new Newton iteration, which converges locally to the multivariate Hankel matrix of a given rank the closest to a given input Hankel matrix. Numerical experimentations show that the Newton iteration combined the decomposition method allows to compute accurately and efficiently the polynomial-exponential decomposition of the symbol, even for noisy input moments.

Structure of the paper

The next section describes multivariate Hankel operators, their symbol and the generalization of Kronecker theorem, which establishes a correspondence between Hankel operators of finite rank and polynomial-exponential series. In Section 3, we recall techniques exploiting the properties of multiplication operators for solving polynomial systems and show how they can be used for the Artinian Gorenstein algebra associated to Hankel operators of finite rank. In Section 4, we describe in details the decomposition algorithm. Finally in section 5, we present numerical experimentations, showing the numerical behaviour of the decomposition method for noisy input moments and the improvements obtained by rescaling and by an iterative projection method.

2 Hankel matrices and operators

Hankel matrices are structured matrices of the form

[TABLE]

where the entry $\sigma_{i+j}$ of the $i^{\mathrm{th}}$ row and the $j^{\mathrm{th}}$ columns depends only on the sum $i+j$. By reversing the order of the columns or the rows, we obtain Toeplitz matrices, which entries depend on the difference of the row and column indices. Exploiting their structure leads to superfast methods for many linear algebra operations such as matrix-vector product, solution of linear systems, …(see e.g. [8]).

A Hankel matrix is a sub-matrix of the matrix of the Hankel operator associated to a sequence $\sigma=(\sigma_{k})\in\mathbb{C}^{{\mathbb{N}}}$:

[TABLE]

where $L_{0}(\mathbb{C}^{{\mathbb{N}}})$ is the set of sequences of $\mathbb{C}^{{\mathbb{N}}}$ with a finite support.

Multivariate Hankel matrices have a similar structure of the form

[TABLE]

where $A,B\subset{\mathbb{N}}^{n}$ are subsets of multi-indices $\alpha=(\alpha_{1},\ldots,\alpha_{n})\in{\mathbb{N}}^{n}$ $\beta=(\beta_{1},\ldots,\beta_{n})\in{\mathbb{N}}^{n}$ indexing respectively the rows and columns. Multivariate Hankel operators are associated to multi-index sequences $\sigma=(\sigma_{\alpha})_{\alpha\in{\mathbb{N}}^{n}}\in\mathbb{C}^{{\mathbb{N}}^{n}}$:

[TABLE]

where $L_{0}(\mathbb{C}^{{\mathbb{N}}^{n}})$ is the set of sequences of $\mathbb{C}^{{\mathbb{N}}^{n}}$ with a finite support. In order to describe the algebraic properties of Hankel operators, we will identify hereafter the space $L_{0}(\mathbb{C}^{{\mathbb{N}}^{n}})$ with the ring $\mathbb{C}[\mathbf{x}]=\mathbb{C}[x_{1},\ldots,x_{n}]$ of polynomials in the variables $\mathbf{x}=(x_{1},\ldots,x_{n})$ with coefficients in $\mathbb{C}$.

The set of multi-index sequences $\mathbb{C}^{{\mathbb{N}}^{n}}$ can be identified with the ring of formal power series $\mathbb{C}[[y_{1},\ldots,y_{n}]]=\mathbb{C}[[\mathbf{y}]]$. A sequence $\sigma=(\sigma_{\alpha})_{\alpha}$ is identified with the series

[TABLE]

where $\mathbf{y}^{\alpha}=y_{1}^{\alpha_{1}}\cdots y_{n}^{\alpha_{n}}$, $\alpha!=\prod_{i=1}^{n}\alpha_{i}!$ for $\alpha=(\alpha_{1},\ldots,\alpha_{n})\in{\mathbb{N}}^{n}$. It can also be interpreted as a linear functional on polynomials as follows:

[TABLE]

The identification of $\sigma$ with an element of $\mathbb{C}[\mathbf{x}]^{\ast}=\operatorname{Hom}_{\mathbb{C}}(\mathbb{C}[\mathbf{x}],\mathbb{C})$ is uniquely defined by its coefficients $\sigma_{\alpha}=\langle\sigma\mid\mathbf{x}^{\alpha}\rangle$ for $\alpha\in{\mathbb{N}}^{n}$, which are called the moments of $\sigma$. This allows us to identify the dual $\mathbb{C}[\mathbf{x}]^{\ast}$ with $\mathbb{C}[[\mathbf{y}]]$ or with the set of multi-index sequences $\mathbb{C}^{{\mathbb{N}}^{n}}$.

The dual space $\mathbb{C}[\mathbf{x}]^{\ast}\equiv\mathbb{C}[[\mathbf{y}]]$ has a natural structure of $\mathbb{C}[\mathbf{x}]$-module, defined as follows:

[TABLE]

For a polynomial $p=\sum_{\alpha\in{\mathbb{N}}^{n}}p_{\alpha}\mathbf{x}^{\alpha}$ with $p_{\alpha}=0$ for almost all $\alpha\in{\mathbb{N}}^{n}$, we have

[TABLE]

We check that $p\star\sigma=p(\partial_{1},\ldots,\partial_{n})(\sigma)=p({\boldsymbol{\partial}})(\sigma)$ where ${\boldsymbol{\partial}}=(\partial_{1},\ldots,\partial_{n})$ and $\partial_{i}$ is the derivation with respect to the variable $y_{i}$.

Identifying $L_{0}(\mathbb{C}^{{\mathbb{N}}^{n}})$ with $\mathbb{C}[\mathbf{x}]$, the Hankel operator (1) is nothing else than the operator of multiplication by $\sigma$:

[TABLE]

Truncated Hankel operators are obtained by restriction of Hankel operators. For $A,B\subset{\mathbb{N}}^{n}$, let $\langle\mathbf{x}^{B}\rangle\subset\mathbb{C}[\mathbf{x}]$, $\langle\mathbf{y}^{A}\rangle\subset\mathbb{C}[[\mathbf{y}]]$ be the vector spaces spanned respectively by the monomials $\mathbf{x}^{\beta}$ for $\beta\in B$ and $\mathbf{y}^{\alpha}$ for $\alpha\in A$. The truncated Hankel operator of $\sigma$ on $A,B$ is

[TABLE]

The matrix of $H_{\sigma}^{A,B}$ in the bases $(\mathbf{x}^{\beta})_{\beta\in B}$ and $(\frac{\mathbf{y}^{\alpha}}{\alpha!})_{\alpha\in A}$ is of the form:

[TABLE]

It is also called the moment matrix of $\sigma$. Multivariate Hankel matrices have a structure, which can be exploited to accelerate linear algebra operations (see e.g. [20] for more details).

Example

Consider the series $\sigma=1+2y_{1}+3y_{2}+4\frac{y_{1}^{2}}{2}+5y_{1}y_{2}+6\frac{y_{2}^{2}}{2}+7\frac{y_{1}^{3}}{6}+8\frac{y_{1}^{2}y_{2}}{2}+\cdots\in\mathbb{C}[[y_{1},y_{2}]]$. Its truncated Hankel matrix on $A=[(0,0),(1,0),(0,1)]$ (corresponding to the monomials $1,x_{1},x_{2}$), $B=[(0,0),(1,0),(0,1)$, $(2,0)]$ (corresponding to the monomials $1,x_{1},x_{2},x_{1}^{2}$) is

[TABLE]

For $d\in{\mathbb{N}}$, we denote by $\mathbb{C}[\mathbf{x}]_{d}$ the vector space of polynomials of degree $\leq d$. Its dimension is $s_{d}={n+d\choose n}$. For $d,d^{\prime}\in{\mathbb{N}}$, we denote by $H_{\sigma}^{d,d^{\prime}}$ the Hankel matrix of $\sigma$ on the subset of monomials in $\mathbf{x}$ respectively of degree $\leq d$ and $\leq d^{\prime}$. We also denote by $H_{\sigma}^{d,d^{\prime}}$ the corresponding truncated Hankel operator of $H_{\sigma}$ from $\mathbb{C}[\mathbf{x}]_{d^{\prime}}$ to $(\mathbb{C}[\mathbf{x}]_{d})^{\ast}$. More generally, for $U=\{u_{1},\ldots,u_{l}\}\subset\mathbb{C}[\mathbf{x}]$, $V=\{v_{1},\ldots,v_{m}\}\subset\mathbb{C}[\mathbf{x}]$, the Hankel matrix of $\sigma$ on $U$, $V$ is $H_{\sigma}^{U,V}=(\langle\sigma\mid u_{i}\,v_{j}\rangle)_{1\leq i\leq l,1\leq j\leq m}$. We use the same notation $H_{\sigma}^{U,V}$ for the truncated Hankel operator from $\langle V\rangle$ to $\langle U\rangle^{\ast}$.

2.1 Hankel operator of finite rank

We are interested in structured decompositions of Hankel matrices (resp. operators) as sums of Hankel matrices (resp. operators) of low rank. This raises the question of describing the Hankel operators of finite rank and leads to the problem of decomposing them into indecomposable Hankel operators of low rank.

A first answer is given by the celebrated theorem of Kronecker [15].

Theorem 2.1** (Kronecker Theorem).**

The Hankel operator

[TABLE]

is of finite rank $r$, if and only if, there exist polynomials $\omega_{1},\ldots,\omega_{r^{\prime}}\in\mathbb{C}[y]$ and $\xi_{1},\ldots,\xi_{r^{\prime}}\in\mathbb{C}$ distinct s.t.

[TABLE]

with $\sum_{i=1}^{r^{\prime}}(\deg(\omega_{i})+1)=r$.

This results says that the Hankel operator $H_{\sigma}$ is of finite rank, if and only if, its symbol $\sigma$ is of the form

[TABLE]

for some univariate polynomials $\tilde{\omega}_{i}(y)\in\mathbb{C}[y]$ and distinct complex numbers $\xi_{i}$ $i=1,\ldots,r^{\prime}$. Moreover, the rank of $H_{\sigma}$ is $r=\sum_{i=1}^{r^{\prime}}(\deg(\tilde{\omega}_{i})+1)$.

The previous result admits a direct generalization to multivariate Hankel operators, using polynomial-exponential series.

Definition 2.2**.**

For $\xi=(\xi_{1},\ldots,\xi_{n})\in\mathbb{C}^{n}$, we denote $\mathbf{e}_{\xi}(\mathbf{y})=e^{y_{1}\xi_{1}+\cdots+y_{n}\xi_{n}}=\sum_{\alpha\in{\mathbb{N}}^{n}}\xi^{\alpha}\,\frac{\mathbf{y}^{\alpha}}{\alpha!}\in\mathbb{C}[[\mathbf{y}]]$ where $\xi^{\alpha}=\xi_{1}^{\alpha_{1}}\cdots\xi_{n}^{\alpha_{n}}$ for $\alpha=(\alpha_{1},\ldots,\alpha_{n})\in{\mathbb{N}}^{n}$.

Let $\mathcal{POLYEXP}(\mathbf{y})=\left\{\sigma=\sum_{i=1}^{r}\omega_{i}(\mathbf{y})\mathbf{e}_{\xi_{i}}(\mathbf{y})\in\mathbb{C}[[\mathbf{y}]]\mid\xi_{i}\in\mathbb{C}^{n},\omega_{i}(\mathbf{y})\in\mathbb{C}[\mathbf{y}]\right\}$ be the set of polynomial-exponential series. The polynomials $\omega_{i}(\mathbf{y})$ are called the weights of $\sigma$ and $\xi_{i}$ the frequencies.

For $\omega(\mathbf{y})\in\mathbb{C}[\mathbf{y}]$, we denote by $\mu(\omega)$ the dimension of the vector space spanned by $\omega$ and its derivatives $\partial^{\alpha}\omega(\mathbf{y})=\partial_{1}^{\alpha_{1}}\cdots\partial_{n}^{\alpha_{n}}\omega(\mathbf{y})$ of any order for $\alpha=(\alpha_{1},\ldots,\alpha_{n})\in{\mathbb{N}}^{n}$.

The next theorem characterizes the multivariate Hankel operators of finite rank in terms of their symbol [19]:

Theorem 2.3** (Generalized Kronecker Theorem).**

Let $\sigma(\mathbf{y})\in\mathbb{C}[[\mathbf{y}]]$. Then $\operatorname{rank}H_{\sigma}=r<\infty$, if and only if, $\sigma(\mathbf{y})=\sum_{i=1}^{r^{\prime}}\omega_{i}(\mathbf{y})\mathbf{e}_{{\xi_{i}}}(\mathbf{y})\in\mathcal{POLYEXP}(\mathbf{y})$ with $\omega_{i}(\mathbf{y})\in\mathbb{C}[\mathbf{y}]\setminus\{0\}$ and $\xi_{i}\in\mathbb{C}^{n}$ pairwise distinct, with $r=\sum_{i=1}^{r^{\prime}}{\mu}(\omega_{i})$ where ${\mu}(w_{i})$ is the dimension of the inverse system spanned by $\omega_{i}(\mathbf{y})$ and all its derivatives $\mathbf{\partial}^{\alpha}\omega_{i}(\mathbf{y})$ for $\alpha=(\alpha_{1},\ldots,\alpha_{n})\in{\mathbb{N}}^{n}$.

Example 2.4**.**

For $\xi\in\mathbb{C}^{n}$, the series $\mathbf{e}_{\xi}(\mathbf{y})=\sum_{\alpha\in{\mathbb{N}}^{n}}\xi^{\alpha}\,\frac{\mathbf{y}^{\alpha}}{\alpha!}=e^{\mathbf{y}\cdot\xi}$ represents the linear functional corresponding to the evaluation at $\xi$:

[TABLE]

The Hankel operator $H_{\mathbf{e}_{\xi}}:p\mapsto p\star\mathbf{e}_{\xi}=p(\xi)\mathbf{e}_{\xi}$ is of rank $1$, since its image is spanned by $\mathbf{e}_{\xi}$. For $A,B\subset{\mathbb{N}}^{n}$, the Hankel matrix of $\mathbf{e}_{\xi}$ is $H_{\xi}^{A,B}=[\xi^{\beta+\alpha}]_{\beta\in B,\alpha\in A}.$ If $H_{\xi}^{A,B}\neq 0$, it is a matrix of rank $1$.

Hankel operators associated to evaluations $\mathbf{e}_{\xi}$ are of rank $1$. As shown in the next example, a Hankel operator of rank $>1$ is not necessarily the sum of such Hankel operators of rank $1$.

Example 2.5**.**

For $n=1$ and $\sigma=y$, we check that $H_{y}$ is of rank $2$, but it cannot be decomposed as a sum of two rank-one Hankel operators. If $A=\{1,x,x^{2}\}$, we have

[TABLE]

for $\lambda_{1},\lambda_{2},\xi_{1},\xi_{2}\in\mathbb{C}$. This shows that the symbol $y$ is indecomposable as a sum of polynomial-exponential series, though it defines an Hankel operator of rank $2$.

Definition 2.6**.**

For $\sigma\in\mathbb{C}[[\mathbf{y}]]$, we say that $\sigma$ is indecomposable if ${\sigma}$ cannot be written as a sum ${\sigma}={\sigma_{1}}+{\sigma_{2}}$ with $\operatorname{img}H_{\sigma}=\operatorname{img}H_{\sigma_{1}}\oplus\operatorname{img}H_{\sigma_{2}}$.

Proposition 2.7**.**

The series $\omega(\mathbf{y})\,\mathbf{e}_{\xi}(\mathbf{y})$ with $\omega(\mathbf{y})\in\mathbb{C}[\mathbf{y}]\setminus\{0\}$ and $\xi\in\mathbb{C}^{n}$ is indecomposable.

Proof.

Let $\sigma=\omega\,\mathbf{e}_{\xi}$ and $r=\mu(\omega)$ the rank of $H_{\sigma}$. Suppose that $\sigma=\sigma_{1}+\sigma_{2}$ with $\operatorname{img}H_{\sigma}=\operatorname{img}H_{\sigma_{1}}\oplus\operatorname{img}H_{\sigma_{2}}$. We assume that the rank of $H_{\sigma_{1}}$ is minimal. By the Generalized Kronecker Theorem, $\sigma_{1}=\sum_{i=1}^{r_{1}}\omega_{1,i}\,\mathbf{e}_{\xi_{1,i}}$, $\sigma_{2}=\sum_{i=1}^{r_{2}}\omega_{2,i}\,\mathbf{e}_{\xi_{2,i}}$ with $\omega_{l,i}\in\mathbb{C}[\mathbf{x}]$, $\xi_{l,i}\in\mathbb{C}^{n}$ and

[TABLE]

By the independence of the polynomial-exponential series [19][Lem. 2.7], we can assume that $\xi_{1,1}=\xi_{2,1}=\xi$ and that $\omega=\omega_{1,1}+\omega_{2,1}$ (possibly with $\omega_{2,1}=0$) and that $\omega_{1,i}=-\omega_{2,i}$ for $i=2,\ldots,r_{1}=r_{2}$ . As $\operatorname{rank}H_{\sigma_{1}}=\sum_{i=1}^{r_{1}}\mu(\omega_{1,i})$ is minimal, we can assume moreover that $\omega_{1,i}=0$ for $i=2,\ldots,r_{1}$, that is, $r_{1}=r_{2}=1$. Then, we have $\sigma=\omega\,\mathbf{e}_{\xi}$, $\sigma_{1}=\omega_{1}\,\mathbf{e}_{\xi}$ $\sigma_{2}=\omega_{2}\,\mathbf{e}_{\xi}$ with $\omega=\omega_{1}+\omega_{2}$. As $\operatorname{img}H_{\sigma_{i}}=\langle\partial^{\alpha}(\omega_{i})\,\mathbf{e}_{\xi}\rangle$, $i=1,2$, we have $\operatorname{img}H_{\sigma_{1}}\cap\operatorname{img}H_{\sigma_{2}}\ni\mathbf{e}_{\xi}$ and $\operatorname{img}H_{\sigma}$ is not the direct sum of $\operatorname{img}H_{\sigma_{1}}$ and $\operatorname{img}H_{\sigma_{2}}$. This shows that $\sigma$ is indecomposable. ∎

The goal of this paper is to present a method to decompose the symbol of a Hankel operator as a sum of indecomposable polynomial-exponential series from truncated Hankel matrices.

3 Structured decomposition of Hankel matrices

In this section, we show how the decomposition of the symbol $\sigma$ of a Hankel operator $H_{\sigma}$ as a sum of polynomial-exponential series reduces to the solution of polynomial equations. This corresponds to the decomposition of $H_{\sigma}$ as a sum of Hankel matrices of low rank. We first recall classical techniques for solving polynomial systems and show how these methods can be applied on the Hankel matrix $H_{\sigma}$, to compute the decomposition.

3.1 Solving polynomial equations by eigenvector computation

A quotient algebra ${\mathcal{A}}=\mathbb{C}[\mathbf{x}]/I$ is Artinian if it is of finite dimension over $\mathbb{C}$. In this case, the ideal $I$ defines a finite number of roots ${\mathcal{V}}(I)=\{\xi_{1},\ldots,\xi_{r^{\prime}}\}=\{\xi\in\mathbb{C}^{n}\mid\forall p\in I,p(\xi)=0\}$ and we have a decomposition of ${\mathcal{A}}$ as a sum of sub-algebras:

[TABLE]

where ${\mathcal{A}}_{i}=\mathbf{u}_{\xi_{i}}{\mathcal{A}}\sim\mathbb{C}[\mathbf{x}]/Q_{i}$ and $Q_{i}$ is the primary component of $I$ associated to the root $\xi_{i}\in\mathbb{C}^{n}$. The elements $\mathbf{u}_{1},\ldots,\mathbf{u}_{r^{\prime}}$ satisfy the relations

[TABLE]

The polynomials $\mathbf{u}_{\xi_{1}},\ldots,\mathbf{u}_{\xi_{r^{\prime}}}$ are called *idempotents * of ${\mathcal{A}}$. The dimension of ${\mathcal{A}}_{i}$ is the multiplicity of the point $\xi_{i}$. For more details, see [9][Chap. 4].

For $g\in\mathbb{C}[\mathbf{x}]$, the multiplication operator ${\mathcal{M}}_{g}$ is defined by

[TABLE]

The transpose ${\mathcal{M}}_{g}^{\scriptscriptstyle\mathsf{T}}$ of the multiplication operator ${\mathcal{M}}_{g}$ is

[TABLE]

The main property that we will use to recover the roots is the following [9][Thm. 4.23]:

Proposition 3.1**.**

Let $I$ be an ideal of $\mathbb{C}[\mathbf{x}]$ and suppose that ${\mathcal{V}}(I)=\{{\xi}_{1},{\xi}_{2},\ldots,{\xi}_{r^{\prime}}\}$. Then

for all $g\in{\mathcal{A}}$, the eigenvalues of ${\mathcal{M}}_{g}$ and ${\mathcal{M}}_{g}^{{\scriptscriptstyle\mathsf{T}}}$ are the values $g(\mathbf{\xi}_{1}),\ldots,g({\xi}_{r^{\prime}})$ of the polynomial $g$ at the roots with multiplicities ${\mu}_{i}=\dim{\mathcal{A}}_{i}$.

• 2.

The eigenvectors common to all ${\mathcal{M}}_{g}^{{\scriptscriptstyle\mathsf{T}}}$ with $g\in{\mathcal{A}}$ are - up to a scalar - the evaluations ${\mathbf{e}}_{{\xi}_{1}},\ldots,{\mathbf{e}}_{{\xi}_{r^{\prime}}}$.

If $B=\{b_{1},\ldots,b_{r}\}$ is a basis of ${\mathcal{A}}$, then the coefficient vector of the evaluation $\mathbf{e}_{{\xi}_{i}}$ in the dual basis of $B$ is $\left[\left\langle\mathbf{e}_{{\xi_{i}}}|b_{j}\right\rangle\right]_{\beta\in B}=[b_{j}({\xi}_{i})]_{i=1\ldots r}=B(\xi_{i})$. The previous proposition says that if $M_{g}$ is the matrix of $\mathcal{M}_{g}$ in the basis $B$ of ${\mathcal{A}}$, then

[TABLE]

If moreover the basis $B$ contains the monomials $1,x_{1},x_{2},\ldots,x_{n}$, then the common eigenvectors of $M_{g}^{{\scriptscriptstyle\mathsf{T}}}$ are of the form ${\mathbf{v}}_{i}=c\,[1,\xi_{i,1},\ldots,\xi_{i,n},\ldots]$ and the root $\xi_{i}$ can be computed from the coefficients of ${\mathbf{v}}_{i}$ by taking the ratio of the coefficients of the monomials $x_{1},\ldots,x_{n}$ by the coefficient of $1$: $\xi_{i,k}=\frac{{\mathbf{v}}_{i,k+1}}{{\mathbf{v}}_{i,1}}$. Thus computing the common eigenvectors of all the matrices $M_{g}^{{\scriptscriptstyle\mathsf{T}}}$ for $g\in{\mathcal{A}}$ yield the roots ${\xi}_{i}$ ($i=1,\ldots,r$).

In practice, it is enough to compute the common eigenvectors of $M_{x_{1}}^{{\scriptscriptstyle\mathsf{T}}},\ldots,M_{x_{n}}^{{\scriptscriptstyle\mathsf{T}}}$, since $\forall g\in\mathbb{C}[\mathbf{x}],M_{g}^{{\scriptscriptstyle\mathsf{T}}}=g(M_{x_{1}}^{{\scriptscriptstyle\mathsf{T}}},\ldots,M_{x_{n}}^{{\scriptscriptstyle\mathsf{T}}})$. Therefore, the common eigenvectors $M_{x_{1}}^{{\scriptscriptstyle\mathsf{T}}},\ldots,M_{x_{n}}^{{\scriptscriptstyle\mathsf{T}}}$ are also eigenvectors of any $M_{g}^{{\scriptscriptstyle\mathsf{T}}}$.

The multiplicity structure, that is the dual $Q_{i}^{\perp}$ of each primary component $Q_{i}$ of $I$, also called the inverse system of the point $\xi_{i}$ can be deduced by linear algebra tools (see e.g. [18]).

In the case of simple roots, we have the following property [9][Chap. 4]:

Proposition 3.2**.**

If the roots $\{{\xi}_{1},{\xi}_{2},\ldots,{\xi}_{r}\}$ of $I$ are simple (i.e. $\mu_{i}=\dim{\mathcal{A}}_{i}=1$) then we have the following:

$\mathbbm{u}=\{\mathbf{u}_{\xi_{1}},\ldots,\mathbf{u}_{\xi_{r}}\}$* is a basis of ${\mathcal{A}}$.*

• 2.

The polynomials ${\mathbf{u}}_{\xi_{1}},\ldots,\mathbf{u}_{\xi_{r}}$ are interpolation polynomials at the roots $\xi_{i}$: $\mathbf{u}_{\xi_{i}}(\xi_{j})=1$ if $i=j$ and [math] otherwise.

• 3.

The matrix of ${\mathcal{M}}_{g}$ in the basis $\mathbbm{u}$ is the diagonal matrix $\operatorname{diag}(g(\xi_{1}),\ldots,g(\xi_{r}))$.

This proposition tells us that if $g$ is separating the roots, i.e. $g(\xi_{i})\neq g(\xi_{j})$ for $i\neq j$, then the eigenvectors of ${\mathcal{M}}_{g}$ are, up to a scalar, interpolation polynomials at the roots.

3.2 Artinian Gorenstein algebra of a multivariate Hankel operator

We associate to a Hankel operator $H_{\sigma}$, the quotient ${\mathcal{A}}_{\sigma}=\mathbb{C}[\mathbf{x}]/I_{\sigma}$ of the polynomial ring $\mathbb{C}[\mathbf{x}]$ modulo the kernel $I_{\sigma}=\{p\in\mathbb{C}[\mathbf{x}]\mid\forall q\in R,\langle\sigma\mid pq\rangle=0\}$ of $H_{\sigma}$. We check that $I_{\sigma}$ is an ideal of $\mathbb{C}[\mathbf{x}]$, so that ${\mathcal{A}}_{\sigma}$ is an algebra.

As ${\mathcal{A}}_{\sigma}=\mathbb{C}[\mathbf{x}]/I_{\sigma}\sim\operatorname{img}H_{\sigma}$, the operator $H_{\sigma}$ is of finite rank $r$, if and only if, ${\mathcal{A}}_{\sigma}$ is Artinian of dimension $\dim_{\mathbb{C}}{\mathcal{A}}_{\sigma}=r$ .

A quotient algebra ${\mathcal{A}}$ is called Gorenstein if its dual ${\mathcal{A}}^{*}=\operatorname{Hom}_{\mathbb{C}}({\mathcal{A}},\mathbb{C})$ is a free ${\mathcal{A}}$-module generated by one element.

In our context, we have the following equivalent properties [19]:

$\sigma=\sum_{i=1}^{r^{\prime}}\omega_{i}(\mathbf{y})\mathbf{e}_{\xi_{i}}(\mathbf{y})$ with $\omega_{i}\in\mathbb{C}[\mathbf{y}]$, $\xi_{i}\in\mathbb{C}^{n}$ and $\sum_{i=1}^{r^{\prime}}\mu(\omega_{i})=r$,

• 2.

$H_{\sigma}$ is of rank $r$,

• 3.

${\mathcal{A}}_{\sigma}$ is an Artinian Gorenstein algebra of dimension $r$.

The following proposition shows that the frequencies $\xi_{i}$ and the weights $\omega_{i}$ can be recovered from the ideal $I_{\sigma}$ (see [19] for more details):

Proposition 3.3**.**

If $\sigma(\mathbf{y})=\sum_{i=1}^{r^{\prime}}\omega_{i}(\mathbf{y})\mathbf{e}_{{\xi_{i}}}(\mathbf{y})$ with $\omega_{i}(\mathbf{y})\in\mathbb{C}[\mathbf{y}]\setminus\{0\}$ and $\xi_{i}\in\mathbb{C}^{n}$ pairwise distinct, then we have the following properties:

The points ${\xi}_{1},{\xi}_{2},\ldots,{\xi}_{r^{\prime}}\in\mathbb{C}^{n}$ are the common roots of the polynomials in $I_{\sigma}=\ker H_{\sigma}=\{p\in\mathbb{C}[\mathbf{x}]\mid\forall q\in\mathbb{C}[\mathbf{x}],\langle\sigma|pq\rangle=0\}$.

• 2.

The series $\omega_{i}(\mathbf{y})\mathbf{e}_{\xi_{i}}$ is a generator of the inverse system of $Q_{i}^{\bot}$, where $Q_{i}$ is the primary component of $I_{\sigma}$ associated to $\xi_{i}$ such that $\dim\mathbb{C}[\mathbf{x}]/Q_{i}=\mu(\omega_{i})$.

This result tells us that the problem of decomposing $\sigma$ as a sum of polynomial-exponential series reduces to the solution of the polynomial equations $p=0$ for $p$ in the kernel $I_{\sigma}$ of $H_{\sigma}$.

Another property that will be helpful to determine a basis of ${\mathcal{A}}_{\sigma}$ is the following:

Lemma 3.4**.**

Let $B=\{b_{1},\ldots,b_{r}\}$, $B^{\prime}=\{b_{1}^{\prime},\ldots,b_{r}^{\prime}\}\subset\mathbb{C}[\mathbf{x}]$. If the matrix $H_{\sigma}^{B,B^{\prime}}=(\langle\sigma|b_{i}b_{j}^{\prime}\rangle)_{1\leq i,j\leq r}$ is invertible, then B and $B^{\prime}$ are linearly independent in ${\mathcal{A}}_{\sigma}$.

Proof.

Suppose that $H_{\sigma}^{B,B^{\prime}}$ is invertible. If there exists $p=\sum_{i}\lambda_{i}b_{i}$ ($\lambda_{i}\in\mathbb{C}$) such that $p\equiv 0$ in ${\mathcal{A}}_{\sigma}$. Then $p\star\sigma=0$ and $\forall q\in R$, $\langle\sigma|pq\rangle=0$. In particular, for $j=1,\ldots,r$ we have

[TABLE]

As $H_{\sigma}^{B,B^{\prime}}$ is invertible, $\lambda_{i}=0$ for $i=1,\ldots,r$ and $B$ is a family of linearly independent elements in ${\mathcal{A}}_{\sigma}$. Since we have $(H_{\sigma}^{B,B^{\prime}})^{{\scriptscriptstyle\mathsf{T}}}=H_{\sigma}^{B^{\prime},B}$,  we prove by a similar argument that $H_{\sigma}^{B,B^{\prime}}$ invertible also implies that $B^{\prime}$ is linearly independent in ${\mathcal{A}}_{\sigma}$. ∎

By this Lemma, bases of ${\mathcal{A}}_{\sigma}$ can be computed by identifying non-zero minors of maximal size of the matrix of $H_{\sigma}$.

Proposition 3.5**.**

Let $B,B^{\prime}$ be basis of ${\mathcal{A}}_{\sigma}$ and $g\in\mathbb{C}[\mathbf{x}]$. We have

[TABLE]

where $M_{g}^{B}$ (resp. $M_{g}^{B^{\prime}}$) is the matrix of the multiplication by $g$ in the basis $B$ (resp. $B^{\prime}$) of ${\mathcal{A}}_{\sigma}$.

Proof.

Let $B=\{b_{1},\ldots,b_{r}\},B^{\prime}=\{b^{\prime}_{1},\ldots,b^{\prime}_{r}\}$ be two bases of ${\mathcal{A}}_{\sigma}$. We have $g\,{b}_{j}=\sum_{i=1}^{r}m_{i,j}{b}^{\prime}_{i}+\kappa$ where $m_{i,j}$ is the $(i,j)$ entry of the matrix $M_{g}^{B}$ of multiplication by $g$ in the basis $B$ and $\kappa\in I_{\sigma}$. Then,

[TABLE]

Similarly, we have $g\,b^{\prime}_{j}=\sum_{i=1}^{r}m^{\prime}_{i,j}b^{\prime}_{i}+\kappa^{\prime}$ where $m^{\prime}_{i,j}$ is the $(i,j)$ entry of the matrix $M_{g}^{B^{\prime}}$ of multiplication by $g$ in the basis $B^{\prime}$ and $\kappa^{\prime}\in I_{\sigma}$. For ${1\leq i,j\leq r}$, the entry $(i,j)$ of ${H}_{g\star\sigma}^{B,B^{\prime}}$ is

[TABLE]

This concludes the proof of the relations (2). ∎

We deduce the following property:

Proposition 3.6**.**

Let $\sigma(\mathbf{y})=\sum_{i=1}^{r}\omega_{i}(\mathbf{y})\mathbf{e}_{\xi_{i}}(\mathbf{y})$ with $\omega_{i}\in\mathbb{C}[\mathbf{y}]\setminus\{0\}$ and $\xi_{i}\in\mathbb{C}^{n}$ distinct and let $B,B^{\prime}$ be bases of ${\mathcal{A}}_{\sigma}$. We have the following properties:

For $g\in\mathbb{C}[\mathbf{x}]$, $M_{g}^{B^{\prime}}=(H_{\sigma}^{B,B^{\prime}})^{-1}H_{g\star\sigma}^{B,B^{\prime}}$, $(M_{g}^{B})^{{\scriptscriptstyle\mathsf{T}}}=H_{g\star\sigma}^{B,B^{\prime}}(H_{\sigma}^{B,B^{\prime}})^{-1}$.

• 2.

For $g\in\mathbb{C}[\mathbf{x}]$, the generalized eigenvalues of $(H_{g\star\sigma}^{B,B^{\prime}},H_{\sigma}^{B,B^{\prime}})$ are $g({\xi}_{i})$ with multiplicity ${\mu}_{i}={\mu}(\omega_{i})$, $i=1,\ldots,r$.

• 3.

The generalized eigenvectors common to all $(H_{g\star\sigma}^{B,B^{\prime}},H_{\sigma}^{B,B^{\prime}})$ for $g\in\mathbb{C}[\mathbf{x}]$ are - up to a scalar - $(H_{\sigma}^{B,B^{\prime}})^{-1}\,B(\xi_{i})$, $i=1,\ldots,r$.

Proof.

The two first points are direct consequences of Propositions 3.5 and 3.1. The third point is also a consequence of Proposition 3.1, since the coordinate vector of the evaluation $\mathbf{e}_{\xi_{i}}$ in the dual basis of $B$ is $B(\xi_{i})$ for $i=1,\ldots,r$. ∎

This proposition shows that the matrices of multiplication by an element $g$ in ${\mathcal{A}}$, and thus the roots $\{\xi_{1},\ldots,\xi_{r}\}$ and their multiplicity structure, can be computed from truncated Hankel matrices, provided we can determine bases $B$, $B^{\prime}$ of ${\mathcal{A}}_{\sigma}$. In practice, it is enough to compute the generalized eigenvectors common to $(H_{x_{i}\star\sigma}^{B,B^{\prime}},H_{\sigma}^{B,B^{\prime}})$ for $i=1,\ldots,n$ to recover the roots. As $H_{x_{i}\star\sigma}^{B,B^{\prime}}=H_{\sigma}^{x_{i}B,B^{\prime}}=H_{\sigma}^{B,x_{i}\,B^{\prime}}$, the decomposition can be computed from sub-matrices of $H_{\sigma}^{B,B^{\prime+}}$ or $H_{\sigma}^{B^{+},B^{\prime}}$ where $B^{+}=B\cup x_{1}B\cup\cdots\cup x_{n}B$, $B^{\prime+}=B^{\prime}\cup x_{1}B^{\prime}\cup\cdots\cup x_{n}B^{\prime}$.

4 Decomposition algorithm

We are given the first moments $\sigma_{\alpha},|\alpha|\leq d$ of the series $\sigma(\mathbf{y})=\sum_{i=1}^{r}\omega_{i}\mathbf{e}_{\xi_{i}}(\mathbf{y})$ with $\omega_{i}\in\mathbb{C}\backslash{(0)}$ and $\xi_{i}\in\mathbb{C}^{n}.$ The goal is to recover the number of terms r, the constant weights $\omega_{i}$ and the frequencies $\xi_{i}$ of the series $\sigma(\mathbf{y})$.

4.1 Computation of a basis

The first problem is to find automatically bases $B_{1}$ and $B_{2}$ of the quotient algebra ${\mathcal{A}}_{\sigma}$, of maximal sizes such that $H_{\sigma}^{B_{1},B_{2}}$ is invertible. Using Proposition 3.6, we will compute the multiplication matrices $M_{g}^{B_{2}}$ for $g=x_{i},i=1,\ldots,n$. The frequencies $\xi_{j}$ and the weights $\omega_{j}$, $j=1,\ldots,r$ will be deduced from their eigenvectors, as described in section 4.2.

Given the set of moments $(\sigma_{\alpha})_{|\alpha|\leq d}$, we create two sets $A_{1}=(\mathbf{x}^{\alpha})_{|\alpha|\leq d_{1}}$ and $A_{2}=(\mathbf{x}^{\beta})_{|\beta|\leq d_{2}}$ of monomials such that $\alpha$ and $\beta$ are multi-indices in ${\mathbb{N}}^{n}$ with $|\alpha|\leq d_{1}$ and $|\beta|\leq d_{2}$. The degrees $d_{1}$ and $d_{2}$ are chosen such that $d_{1}+d_{2}<d$. Let $N_{1}=|A_{1}|$ and $N_{2}=|A_{2}|$. The truncated Hankel operator associated to $\sigma$ is:

[TABLE]

The Hankel matrix in these two monomial sets $A_{1}$ and $A_{2}$ is defined by $H_{\sigma}^{d_{1},d_{2}}=[\sigma_{(\alpha+\beta)}]_{\begin{subarray}{c}\mid\alpha\mid\leq d_{1}\\ \mid\beta\mid\leq d_{2}\end{subarray}}$.

Computing the singular value decomposition of $H_{\sigma}^{d_{1},d_{2}}$, we obtain

[TABLE]

where $S$ is the diagonal matrix of all singular values of $H_{\sigma}^{d_{1},d_{2}}$ arranged in a decreasing order, $U$ is an unitary matrix whose columns are the left singular vectors of $H_{\sigma}^{d_{1},d_{2}}$, $V$ is an unitary matrix whose columns are the right singular vectors of $H_{\sigma}^{d_{1},d_{2}}$. We denote by $U^{{\scriptscriptstyle\mathsf{H}}}$ the hermitian transpose of $U$ and $\overline{V}$ the conjugate of $V$.

Let $u_{i}=[u_{\alpha,i}]_{\alpha\in A_{1}}$ and $v_{j}=[v_{\beta,j}]_{\beta\in A_{2}}$ be respectively the $i^{\textrm{th}}$ and $j^{\textrm{th}}$ columns of $U^{{\scriptscriptstyle\mathsf{H}}}$ and $\overline{V}$. They are vectors respectively in $\mathbb{C}^{N_{1}}$ and $\mathbb{C}^{N_{2}}$. We denote by $u_{i}(\mathbf{x})=u_{i}^{{\scriptscriptstyle\mathsf{T}}}\,A_{1}=\sum_{|\alpha|\leq d_{1}}u_{\alpha,i}\mathbf{x}^{\alpha}$ and $v_{j}(\mathbf{x})=v_{j}^{{\scriptscriptstyle\mathsf{T}}}\,A_{2}=\sum_{|\beta|\leq d_{2}}v_{\beta,j}\mathbf{x}^{\beta}$ the corresponding polynomials. The bases formed by these first $r$ polynomials are denoted $U_{r}^{{\scriptscriptstyle\mathsf{H}}}:=(u_{i}(\mathbf{x}))_{i=1,\ldots,r}$ and $\overline{V}_{r}:=(v_{j}(\mathbf{x}))_{j=1,\ldots,r}$. We will also denote by $U_{r}^{{\scriptscriptstyle\mathsf{H}}}$ (resp. $\overline{V}_{r}$) the corresponding coefficient matrix, formed by the first rows (resp. columns) of $U^{\scriptscriptstyle\mathsf{H}}$ (resp. $\overline{V}$). We denote by $S_{r}$ the diagonal matrix of the first $r$ rows and columns of $S$, formed by the first $r$ singular values.

Proposition 4.1**.**

Let $\sigma=\sum_{i=1}^{r^{\prime}}\omega_{i}(\mathbf{y})\mathbf{e}_{\xi_{i}}$ with $\omega_{i}\in\mathbb{C}[\mathbf{y}]$, $\xi_{i}\in\mathbb{C}^{n}$ and $\sum_{i=1}^{r^{\prime}}\mu(\omega_{i})=r$. If $\operatorname{rank}H_{\sigma}^{d_{1},d_{2}}=r$, then the sets of polynomials $U_{r}^{{\scriptscriptstyle\mathsf{H}}}$ and $\overline{V}_{r}$ are bases of ${\mathcal{A}}_{\sigma}$. The matrix $M_{x_{i}}^{\overline{V}_{r}}$ associated to the multiplication operator by $x_{i}$ in the basis $\overline{V}_{r}$ of ${\mathcal{A}}_{\sigma}$ is $M_{x_{i}}^{\overline{V}_{r}}=S_{r}^{-1}\,U_{r}^{{\scriptscriptstyle\mathsf{H}}}\,H_{x_{i}\star\sigma}^{d_{1},d_{2}}\,\overline{V}_{r}$ $i=1,\ldots,n$.

Proof.

The $(i,j)$ entry of the matrix $H_{\sigma}^{U_{r}^{{\scriptscriptstyle\mathsf{H}}},\overline{V}_{r}}$ of the truncated Hankel operator of $\sigma$ with respect to $U_{r}^{{\scriptscriptstyle\mathsf{H}}}$ and $\overline{V}_{r}$ is equal to:

[TABLE]

Using the SVD decomposition of $H_{\sigma}^{d_{1},d_{2}}$, we have

[TABLE]

since $U^{H}\,U=\mathrm{Id}_{N_{1}}$, $V^{T}\,\overline{V}=\mathrm{Id}_{N_{2}}$. As $r=\operatorname{rank}H_{\sigma}^{d_{1},d_{2}}$, $S_{r}$ is invertible and by Lemma 3.4, $U_{r}^{{\scriptscriptstyle\mathsf{H}}}$ and $\overline{V}_{r}$ are linearly independent in ${\mathcal{A}}_{\sigma}$, which is a vector space of dimension $r$. Thus they are bases of ${\mathcal{A}}_{\sigma}$.

Let $H_{{x_{i}}\star\sigma}^{U_{r}^{{\scriptscriptstyle\mathsf{H}}},\overline{V}_{r}}$ be the matrix of the truncated Hankel operator of ${x_{i}}\star\sigma$ on the two bases $U_{r}^{{\scriptscriptstyle\mathsf{H}}}$ and $\overline{V}_{r}$. A similar computation yields $H_{{x_{i}}\star\sigma}^{U_{r}^{{\scriptscriptstyle\mathsf{H}}},\overline{V}_{r}}=U_{r}^{{\scriptscriptstyle\mathsf{H}}}H_{x_{i}\star\sigma}^{d_{1},d_{2}}\overline{V}_{r}$, where $H_{x_{i}\star\sigma}^{d_{1},d_{2}}$ is the matrix of the truncated Hankel operator of ${x}_{i}\star\sigma$ in the bases $A_{1}$ and $A_{2}$ for all $i=1,\dots,n$. Since $S_{r}$ is an invertible matrix, by Proposition 3.6 we obtain $M_{\mathbf{x}_{i}}^{\overline{V}_{r}}=(H_{\sigma}^{U_{r}^{{\scriptscriptstyle\mathsf{H}}},\overline{V}_{r}})^{-1}H_{x_{i}\star\sigma}^{U_{r}^{{\scriptscriptstyle\mathsf{H}}},\overline{V}_{r}}=S_{r}^{-1}U_{r}^{{\scriptscriptstyle\mathsf{H}}}H_{x_{i}\star\sigma}^{d_{1},d_{2}}\overline{V}_{r}$. ∎

By this proposition $U^{{\scriptscriptstyle\mathsf{H}}}_{r}$ and $\overline{V}_{r}$) are bases of ${\mathcal{A}}_{\sigma}$. By Proposition 3.6, the eigenvalues of $M_{x_{i}}^{\overline{V}_{r}}$ are the $i^{\text{th}}$ coordinates $x_{i}(\xi_{j})=\xi_{j,i}$ of the roots $\xi_{j}$ for $i=1,\ldots,n,j=1,\ldots,r$.

4.2 Computation of the weights

The weight $\omega_{i},i=1,\ldots,r$ of the decomposition of $\sigma$ can be easily computed using the eigenvectors of all $M_{x_{j}}^{\overline{V}_{r}},j=1,\ldots,n$ as follows.

Proposition 4.2**.**

Let $\sigma=\sum_{i=1}^{r}\omega_{i}\,\mathbf{e}_{\xi_{i}}$ with $\omega_{i}\in\mathbb{C}\setminus\{0\}$, $\xi_{i}=(\xi_{i,1},\ldots,\xi_{i,n})\in\mathbb{C}^{n}$ distinct and let $M_{x_{j}}^{\overline{V}_{r}}$ be the matrix of multiplication by $x_{j}$ in the basis $\overline{V}_{r}$. Let $\mathbf{v}_{i}$ be a common eigenvector of $M_{x_{j}}^{\overline{V}_{r}}$, $j=1,..,n$ for the eigenvalues $\xi_{i,j}$. Then the weight of $\mathbf{e}_{\xi_{i}}$ in the decomposition of $\sigma$ is

[TABLE]

Proof.

According to Proposition 3.2, the eigenvectors of the multiplication operator ${\mathcal{M}}_{x_{i}}$ are, up to scalar, the interpolation polynomials $\mathbf{u}_{i}(\mathbf{x})$ at the roots. Let $\mathbf{u}_{\xi_{i}}$ be the coefficient vector associated to $\mathbf{u}_{\xi_{i}}(\mathbf{x})$ in the basis $\overline{V}_{r}$ of ${\mathcal{A}}_{\sigma}$. Let $\mathbf{v}_{i}=\lambda\mathbf{u}_{\xi_{i}}$ be the eigenvector of $M_{x_{i}}^{\overline{V}_{r}}$ associated to the eigenvalue $\xi_{j,i}$ for $j=1,\ldots,r,i=1,\ldots,n$ such that $\mathbf{v}_{i}(\mathbf{x})={A_{2}}^{{\scriptscriptstyle\mathsf{T}}}\,\tilde{\mathbf{v}}_{i}=\sum_{|\beta|\leq d_{2}}\tilde{\mathbf{v}}_{i\beta}\mathbf{x}^{\beta}$ where $\tilde{\mathbf{v}}_{i}=\overline{V}_{r}\,\mathbf{v}_{i}$. Applying the series on all the idempotents, we obtain

[TABLE]

Therefore, we have $\omega_{i}=\frac{\langle\sigma\mid\lambda\mathbf{u}_{\xi_{i}}(\mathbf{x})\rangle}{\lambda}=\frac{\langle\sigma\mid\mathbf{v}_{i}(\mathbf{x})\rangle}{\lambda}=\frac{\langle\sigma\mid\mathbf{v}_{i}(\mathbf{x})\rangle}{\mathbf{v}_{i}(\xi_{i})}$ because of $\mathbf{v}_{i}(\xi_{i})=(\lambda\mathbf{u}_{\xi_{i}})(\xi_{i})=\lambda$. Then

[TABLE]

where $[1]$ is the vector of coefficients of the polynomial $1$ in the monomial basis $A_{1}=(\mathbf{x}^{\alpha})_{|\alpha|\leq d_{1}}$ and

[TABLE]

We deduce that $\omega_{i}=\frac{[1]^{{\scriptscriptstyle\mathsf{T}}}\,H_{\sigma}^{d_{1},d_{2}}\,\overline{V}_{r}\,\mathbf{v}_{i}}{[\xi_{i}^{\alpha}]_{\alpha\in A_{2}}^{{\scriptscriptstyle\mathsf{T}}}\overline{V}_{r}\,\mathbf{v}_{i}}.$ ∎

4.3 Algorithm

We describe now the algorithm to recover the sum $\displaystyle\sigma(\mathbf{y})=\sum_{j=1}^{r}\omega_{j}\,\mathbf{e}_{\xi_{j}}(\mathbf{y})$, $\omega_{j}\in\mathbb{C}\setminus\{0\},\xi_{j}\in\mathbb{C}^{n}$, from the first coefficients $(\sigma_{\alpha})_{|\alpha|\leq d}$ of the formal power series $\displaystyle\sigma=\sum_{\alpha}\sigma_{\alpha}\frac{\mathbf{y}^{\alpha}}{\alpha!}$.

5 Experimentation

In this section, we present numerical experimentations for the decomposition of $\sigma=\sum_{\alpha\in{\mathbb{N}}}\sigma_{\alpha}\frac{\mathbf{y}^{\alpha}}{\alpha!}$ from its moments $\sigma_{\alpha}$. For a given number of variables $n$ and a fixed degree $d$, we compute the coefficients $\sigma_{\alpha}=\sigma(\mathbf{x}^{\alpha})=\sum_{i=1}^{r}\omega_{i}\xi_{i}^{\alpha}$ such that $|\sigma_{\alpha}|\leq d$ where $\omega_{j}\in\mathbb{C}\backslash{(0)}$ and $\xi_{i}=(\xi_{i,1},\ldots,\xi_{i,n}),i=1,\ldots,r$ have random uniform distributions such that $0.5\,M\leq|\xi_{i,j}|\leq 1.5\,M$, $-\pi\leq arg({\xi_{i,j}})\leq\pi$, $0.5\leq|\omega_{i}|\leq 1$ and $-\pi\leq arg({\omega_{i}})\leq\pi$, for $M\geq 1$. To analyse the numerical behaviour of our method, we compare the results with the known frequencies and weights used to compute the moments of $\sigma$.

We use Maple 16 to implement the algorithms. The arithmetic operations are performed on complex numbers, with a numerical precision fixed to $Digits=15$.

5.1 Numerical behavior against perturbation

We apply random perturbations on the moments of the form $\sigma_{\alpha}+\epsilon(p_{\alpha}+\boldsymbol{i}\,q_{\alpha})$ where $p_{\alpha}$ and $q_{\alpha}$ are two random numbers in $[-1,1]$ with a uniform distribution, and $\varepsilon=10^{-e}$ where $e$ is a fixed positive integer.

To measure the consistency of our algorithm, we compute the maximum error between the input frequencies $\xi_{i}$ and the output frequencies $\xi_{i}^{{}^{\prime}}$, and between the input weights $\omega_{i}$ and the output weights $\omega_{i}^{{}^{\prime}}$:

[TABLE]

In each computation, we compute the average of the maximum errors resulting from $10$ tests.

In Figures 1a and 1b, we study the evolution of the error in terms of the perturbation $\varepsilon=10^{(-e)}$, for a fixed degree $d=10$, a number of variables $n=3$, different ranks $r=5,10,20,30$ and for two different amplitudes of the frequencies $M=1$ and $M=100$.

In Figure 1a for $M=1$, the lower error is for the lower rank $r=5$. Between $\varepsilon\approx 10^{-12}$ and $\varepsilon=1$, the error err increases in terms of the perturbation as err$=exp(t\,e)$ for some slope $t\approx 1$. The slope $t$ remains approximately constant but the error increases slightly with the rank $r$. Before $\varepsilon=10^{-13}$, it is approximately constant (approximately $10^{-12}$ for $r=5$) This is due to the fact, in this range, the perturbation is lower than the numerical precision.

In Figure 1b for $M=100$, the lower error is also for the lower rank. The error has almost a constant value when $e$ varies. It is bigger than for $M=1$ for small perturbations. For $r=5,10$ the error slightly increases between $e=-2$ and $e=0$, with a similar slope. This figure clearly shows that the error degrades significantly from $M=1$ to $M=100$ and that the degradation increases rapidly with the rank $r$.

In Figure 2a, we fix the number of variables $n=3$, the rank $r=20$ and we change the degree $d$ which induces a change in the dimensions of the Hankel matrices. For $e\in[-19,0]$, the error decreases when we increase the degree from $d=8$ to $d=10$. It is slightly lower when $d=12$ than when $d=10$, and error is similar for $d=10$ and $d=16$. This increase of the precision with the degree can be related to ratio of number of moments by the number of values to recover in the decomposition.

In Figure 2b, we fix the degree $d=10$ and the rank $r=15$ and we change the number of variables $n=2,3,4,5$. The dimension of the matrices increases polynomially with $n$. We observe that the error decreases quickly with $n$. It shows that the precision improves significantly with the dimension.

5.2 Numerical Rank

To compute the number $r$ of terms in the decomposition of $\sigma$, we arrange the diagonal entries in the decreasing order $s_{1}\geq s_{2}\geq\cdots\geq s_{r}>s_{r+1}\geq\cdots\geq{s_{N_{2}}}$ and we determine the numerical rank of $H_{\sigma}^{d_{1},d_{2}}$ by fixing the largest integer $r$ such that $s_{r}/{s_{1}}\geq\epsilon$.

It is known that the ill-conditioning of the Hankel matrix associated to Prony’s method is in the origin of a numerical instability with respect to perturbed measurements

[TABLE]

In our algorithm the computation of the numerical rank can be affected by this instability. We can explain this instability, using a reasoning close to [27], as follows.

We denote by $s_{j}$ (resp. $\tilde{s}_{j}$) the $j^{\mathrm{th}}$ largest singular value of $H:=H_{\sigma}^{d_{1},d_{2}}$ (resp. $\tilde{H}:=H_{\tilde{\sigma}}^{d_{1},d_{2}}$). The perturbation result for singular values satisfies the estimate (see [11])

[TABLE]

Then, as long as the perturbation is small relative to the conditioning of the problem, that is

[TABLE]

then $\|s_{j}-\tilde{s}_{j}\|\leq\frac{1}{2}s_{r}\;\forall j$ and therefore $\tilde{s}_{r}\geq\frac{1}{2}s_{r}$ and $\tilde{s}_{r+1}\leq\frac{1}{2}s_{r}$. Hence by taking $\varepsilon\leq\frac{1}{2}s_{r}$ as a threshold level we will be sure that the rank is calculated correctly.

But the problem may be badly ill-conditioned and then such a level will not be reasonable. In fact

[TABLE]

where $v_{i,d_{1}}=(\xi_{i}^{|\alpha|})_{|\alpha|\leq d_{1}}$ (resp. $v_{i,d_{2}}=(\xi_{i}^{|\beta|})_{|\beta|\leq d_{2}}$) is the $i^{\mathrm{th}}$ column of the Vandermonde matrix $V_{d_{1}}=(\xi_{i}^{\alpha})_{\begin{subarray}{c}1\leq i\leq r\\ |\alpha|\leq d_{1}\end{subarray}}$ (resp. $V_{d_{2}}=(\xi_{i}^{\beta})_{\begin{subarray}{c}1\leq i\leq r\\ |\beta|\leq d_{2}\end{subarray}}$). Then

[TABLE]

where $C=\text{diag}\left((\omega_{i})_{1\leq i\leq r}\right)$ is the diagonal matrix with $\omega_{i}$ on the diagonal.

Now, using the fact that

[TABLE]

we remark that if $V_{d_{2}}$ (resp. $V_{d_{1}}$) is ill-conditioned then $\|V_{d_{2}}x\|_{2}$ (resp. $\|V_{d_{1}}CV_{d_{2}}^{T}x\|$) may be very small and $s_{r}(H)$ is small as well. This situation can also be produced if $\max_{1\leq i\leq r}\omega_{i}$ is very small. In our numerical experiments, the $\omega_{i}$ are chosen randomly in $[0.5,1]$ and then they don’t cause any numerical instability.

On the other hand, the $\xi_{i}$ vary in such a way that their amplitude can be large, which can generate very ill-conditioned Vandermonde matrices. In fact, it is known (see [21]), that for a nonsingular univariate Vandermonde matrix $V=(a_{i}^{j})_{0\leq i,j\leq n-1}$, where $(a_{i})_{0\leq i\leq n-1}$ denotes a vector of $n$ distinct knots, the condition number of $V$ is exponential in $n$ if $\max_{0\leq i\leq n-1}|a_{i}|>1$ or in $k$ if $|a_{i}|<1$ for at least $k$ knots $a_{i}$. Therefore an $n\times n$ Vandermonde matrix is badly ill-conditioned unless all knots lie in or near the disc $D(0,1)=\{x:\;|x|\leq 1\}$ and unless they lie mostly on or near its boundary $C(0,1)$.

In the multivariate case, it appears that the condition number of multivariate Vandermonde matrices has the same behavior as in univariate case. That is, it is exponential in the highest degree of the entries.

According to the foregoing, when the amplitude $M$ of the frequencies increases (even for moderate values of $M$) the numerical rank calculated by truncating the singular values of $H$ will be different from the exact rank of $H$. An idea to remedy this problem is to rescale the frequencies $\xi_{i}$ in order to obtain points with coordinates close to the unitary circle $C(0,1)$.

5.3 Rescaling

As we have seen in Figures 1a and 1b, the error increases significantly with the amplitude $M$. To remedy this issue, we present a rescaling technique and its numerical impact. It’s done like this:

For a chosen non-zero constant $\lambda$, we transform the input moments of the series as follows:

[TABLE]

which corresponds to the scaling on the frequencies $e_{\xi}(\lambda\mathbf{y})=e_{\lambda\xi}(\mathbf{y})$.

• 2.

We compute decomposition of $\tilde{\sigma}(\mathbf{y})=\sigma(\lambda\mathbf{y})$ from the moments $\tilde{\sigma}_{\alpha}=\lambda^{|\alpha|}\sigma_{\alpha}$.

• 3.

We apply the inverse scaling on the computed frequencies $\tilde{\xi}_{i}$ which gives \xi_{i}=\frac{\tilde{\xi}_{i}}{\lambda}=\big{(}\frac{\tilde{\xi}_{i,1}}{\lambda},\ldots,\frac{\tilde{\xi}_{i,n}}{\lambda}\big{)}.

To determine the scaling factor $\lambda$, we use $\lambda:=\frac{1}{m}$ where $m=\frac{max_{|\alpha|=d}|\sigma_{\alpha}|}{max_{|\alpha|={d-1}}|\sigma_{\alpha}|}$. This is justified as follows: If $|\omega_{j}|\leq 1,j=1,\dots,r$, then $|\sigma_{\alpha}|=|\sum_{j=1}^{r}\omega_{j}{\xi_{j}}^{\alpha}|\simeq M^{d}$ for $|\alpha|=d$ big and for $M$ is the highest modulus of frequencies. Similarly $|\sigma_{\alpha^{{}^{\prime}}}|\simeq M^{d-1}$ for $|\alpha^{{}^{\prime}}|=d-1$. Then we have $m=\frac{max_{|\alpha|=d}|\sigma_{\alpha}|}{max_{|\alpha|={d-1}}|\sigma_{\alpha}|}\approx M$.

To study the numerical influence of the rescaling, we compute the maximum relative error between the input frequencies $\xi_{i}$ and the output frequencies $\tilde{\xi}_{i}$, and the maximum error between the input weights $\omega_{i}$ and the output weights $\tilde{\omega}_{i}$, and we take their maximum:

[TABLE]

$\mathrm{where}\hskip 2.84526pt\mathrm{err}(\omega_{i},\tilde{\omega}_{i})=\max_{1\leq i\leq r}|\omega_{i}-\tilde{\omega}_{i}|\hskip 2.84526pt\mathrm{and}\hskip 2.84526pt\mathrm{rel.err}(\xi_{i},\tilde{\xi}_{i})=\max_{1\leq i\leq r}\frac{{\|\xi_{i}-\tilde{\xi}_{i}\|}_{2}}{{\|\xi_{i}\|}_{2}}$.

In Figure 3, we see the influence of the rescaling on the maximum relative error. The perturbation on the moments is of the order $\varepsilon=10^{-6}$. Each curve for $r=5,10,20,30,$ has almost a constant evolution with the increasing values of $M$ between $10^{2}$ and $10^{10}$. The maximum relative error is lower when $M=100$ than when $M=1$ which is confirmed with the results shown in Figures 1a and 1b. When we increase $r$ the maximum relative error decreases slightly.

In conclusion, the rescaling has an important influence on the computation of the maximum relative error when the modulus $M$ of points is quite big.

The scaling of moments by some computed factor $\lambda$ also enhances the computation of the numerical rank $r$ and leads to a better decomposition as we have seen in 5.2.

5.4 Newton iteration

Given a perturbation $\tilde{\sigma}=\sum_{\alpha}\tilde{\sigma}_{\alpha}\frac{\mathbf{y}^{\alpha}}{\alpha!}$ of a polynomial-exponential series $\sigma=\sum_{i=1}^{r}\omega_{i}\mathbf{e}_{\xi_{i}}(\mathbf{y})$, we want to remove the perturbation on $\tilde{\sigma}$ by computing the polynomial-exponential series of rank $r$, which is the closest to the perturbed series $\tilde{\sigma}$. Starting from an approximate decomposition, using the previous method on the perturbed data, we apply a Newton-type method to minimize the distance between the input series and a weighted sum of $r$ exponential terms.

To evaluate the distance between the series, we use the first moments $\tilde{\sigma}_{\alpha}$ for $\alpha\in A$, where $A$ is a finite subset of ${\mathbb{N}}^{n}$. For $\alpha\in A$, let $F_{\alpha}(\Xi)=\sum_{i=1}\omega_{i}\xi_{i}^{\alpha}-\tilde{\sigma}_{\alpha}$ be the error function for the moment $\tilde{\sigma}_{\alpha}$, where $\omega_{i},\xi_{i,j}$ are variables. We denote by $\Xi=(\xi_{i,j})_{1\leq i\leq r,0\leq j\leq n}$ this set of variables, with the convention that $\xi_{i,0}=\omega_{i}$ for $i=1,\ldots,r$. Let $I=[1,r]\times[0,n]=\{(i,j)\mid 1\leq i\leq r,0\leq j\leq n\}$ be the indices of the variables and $N=(n+1)\,r=|I|$. We denote by $F(\Xi)=(F_{\alpha}(\Xi))_{\alpha\in A}$ the vector of these error functions.

We want to minimize the distance

[TABLE]

Let $M(\Xi_{i})=[w_{i}\xi_{i}^{\alpha}]_{\alpha\in A}$. We denote by $V(\Xi)=\left(\partial_{(i,j)}M(\Xi_{i})\right)_{(i,j)\in I}$ the $|A|\times N$ Vandermonde-like matrix, which columns are the vectors $\partial_{(i,j)}M(\Xi_{i})$. The gradient of $E(\Xi)$ is

[TABLE]

where $\partial_{(i,j)}$ is the derivation with respect to $\Xi_{i,j}$ for $(i,j)\in I$. We denote by $V(\Xi)=\left(\partial_{(i,j)}M(\Xi_{i})\right)_{(i,j)\in I}$ the $|A|\times N$ Vandermonde-like matrix, which columns are the vectors $\partial_{(i,j)}M(\Xi_{i}),(i,j)\in I$.

To find a local minimizer of $E(\Xi)$, we compute a solution of the system $\nabla E(\Xi)=0$, by Newton method. The Jacobian of $\nabla E(\Xi)$ with respect to the variables $\Xi$ is

[TABLE]

Notice that $\partial_{{(i,j)}}\partial_{(i^{\prime},j^{\prime})}M(\Xi_{i})=0$ if $i\neq i^{\prime}$ so that the first matrix is a block diagonal matrix. Then, Newton iteration takes the form:

[TABLE]

To study the numerical influence of Newton method, we compute the maximum absolute error between the input frequencies $\tilde{\xi}_{i}$ and the output frequencies $\tilde{\xi}_{i}$, and the maximum error between the input weights $\omega_{i}$ and the output weights $\tilde{\omega}_{i}$ as in (5).

Figures 4a and 4b show that Newton iterations improve the error. The error decreases by a factor of $\approx 10^{2}$ compared to the computation without Newton iterations. In Figure 4b for $M=100$ the error is smaller than without Newton iterations by a similar order of magnitude (see in Figure 1b).

References