# Stability Analysis of a Variant of the Prony Method

Authors: Rodney Jaramillo; Marianela Lentini
Journal: Mathematical Problems in Engineering
 (2012)

---

## Abstract

Prony type methods are used in many engineering applications to determine the exponential fit corresponding to a dataset. In this paper we study a variant of Prony's method that was used by Martín-Landrove et al., in a process of segmentation ofT2-weighted MRI brain images. We show the equivalence between that method and the classical Prony method and study the stability of the computed solutions with respect to noise in the data set. In particular, we show that the relative error in the calculation of the exponential fit parameters is linear with respect to noise in the data. Our analysis is based on classical results from linear algebra, matrix computation theory, and the theory of stability for roots of polynomials.


---

## Body

## 1. Introduction

In this paper we consider the problem of recovering the parametersb, λj, and Cj  j=1,…,k, from a set of measurements yi, i=1,…,n with n≥2k+1 and a given exponential model
(1.1)μ(t)=b+∑j=1k‍Cje-λjt.
The measurement yi corresponds to the value of expression (1.1) for t=iΔt, Δt>0. Hence we get the set of equations
(1.2)yi=b+∑j=1k‍Cje-iλjΔti=1,…,n.
This kind of exponential fit appears in the case of T2-weighted magnetic resonance images of the brain. T2 is the transverse relaxation time in the process of magnetic resonance; it measures for how long the transverse magnetization lasts in a uniform external magnetic field. Two parameters, TR (repetition time) and TE (echo time) characterize the acquisition of the signal; T2 weighted refers to a signal produced with long TR and TE.The exponentsλj are the nonlinear parameters and correspond to the relaxation rates associated with the different tissues present in the images. The coefficients b,C1,…,Ck are the linear parameters and are related to the noise, and the proportions of those tissues in the given images. In previous works Martín-Landrove et al. [1, 2], it is considered a variant of the Prony method to recover the relaxation rates, the noise and the proportions. Paluszny et al. [3] compared that variant of the Prony method with the Variable Projection Method of Golub and Pereyra, [4–6]. The method proved to be reliable for low noise levels and k=1,2,3. The Variable Projection Method proved to be more robust in the presence of larger noise levels but required ten times more computational time to get results comparable to those of the variant of the Prony method. In this paper we first show that the method proposed by Martín-Landrove is equivalent to the Prony method described by Kahn et al. in [7], and then we study the behavior of the linear systems and the conditioning of the polynomial roots that have to be computed to obtain the model parameters using the method.


## 2. Prony Method

In the model formulation (1.1) set b=C0 and λ0=0, and write
(2.1)μ(t)=∑j=0k‍Cje-λjt.
The Prony type methods, also known as polynomial methods, use the fact that μ(t) satisfies a difference equation of the form
(2.2)(δk+2Ek+1+⋯+δ2E+δ1)μ(t)=0,
where E is the forward shift operator
(2.3)(Eμ)(t)=μ(t+Δt),
and the values βj=e-λjΔt are the roots of the polynomial
(2.4)P(z)=δk+2zk+1+⋯+δ2z+δ1=0,
which is called the characteristic polynomial of the associated difference equation (2.2). Evaluating (2.2) for ti=iΔt,i=1,…,n, we get the following set of linear equations:
(2.5)δk+2μ(tk+2)+⋯+δ1μ(t1)=0⋮⋮δk+2μ(tn)+⋯+δ1μ(tn-k-1)=0.
Let δ=[δ1,…δk+2]T, v=[v1,…,vn]T=[μ(t1),…,μ(tn)]T and
(2.6)Xδ=[δ1⋮⋱δ1δk+2⋱⋮δk+2].
Then the above system of linear equations can be rewritten in matrix form as
(2.7)XδTv=0.
Alternatively, if Wv is the (n-k-1)×(k+2) Hankel matrix
(2.8)Wv=[v1v2⋯vk+2··⋯···⋯···⋯·vn-k-1·⋯vn]
then
(2.9)Wvδ=XδTv.
If there is no noise in the observation data, then yi=μ(ti) and the coefficients δi can be determined from the equivalent system of equations
(2.10)XδTy=0,
where y=[y1,…yn]T. Then the βj are computed as the roots of P(z) and finally
(2.11)λj=-log(βj)Δtj=0,…,k.
In the presence of noisy data instead of solving the system (2.10), we consider the nonlinear optimization problem
(2.12)minδyTXδXδTy.
In order to obtain a nontrivial solution, it is necessary to impose restrictions over the parameters δj; each choice of restrictions characterizes a particular version of the Prony method. For example the modified Prony method described in Osborne and Smyth [8, 9], uses
(2.13)minδyTXδXδ+ysubjectto∥δ∥22=1,
where Xδ+ is the Moore Penrose generalized inverse of Xδ. Osborne et al. considered Xδ to be full rank, hence Xδ+=(XδTXδ)-1XδT.We consider a different optimization problem, which appears in the literature, to compare with the Martín-Landrove method:(2.14)minδyTXδXδTysubjecttoδk+2=1.
The above methods and others, such as classical Prony method, Pisarenko's method, and the linear predictor method, are described in [7–11]. Once the nonlinear parameters have been found, the linear ones are computed as the least squares solution of the linear system obtained by replacing the nonlinear parameters in (1.2), a separation of variables approach.


## 3. An Alternative Formulation of the Prony Method

To simplify the explanation let us consider the casek=3 and n≥7. Then, the system of equations (1.2) can be written as
(3.1)y1=b+C1β1+C2β2+C3β3y2=b+C1β12+C2β22+C3β32⋮yn=b+C1β1n+C2β2n+C3β3n.
From (3.1) and defining qi=yi-yi+1 we get
(3.2)q1=C1(1-β1)β1+C2(1-β2)β2+C3(1-β3)β3⋮qn-1=C1(1-β1)β1n-1+C2(1-β2)β2n-1+C3(1-β3)β3n-1.
The dimension of the system can be reduced by using the transformation qj+1-qjβ1, to get
(3.3)q2-q1β1=C2(1-β2)β2(β2-β1)+C3(1-β3)β3(β3-β1)⋮qn-1-qn-2β1=C2(1-β2)β2n-2(β2-β1)+C3(1-β3)β3n-2(β3-β1).
Then
(3.4)q2=q1β1+C2(1-β2)β2(β2-β1)+C3(1-β3)β3(β3-β1)⋮qn-1=qn-2β1+C2(1-β2)β2n-2(β2-β1)+C3(1-β3)β3n-2(β3-β1).
We now apply the transformation qj+1-qjβ2, and the new set can be written as
(3.5)q3-q2β2=q2β1-q1β1β2+C3(1-β3)β3(β3-β1)(β3-β2)⋮qn-1-qn-2β2=qn-2β1-qn-3β1β2+C3(1-β3)β3n-3(β3-β1)(β3-β2).
The equations above are equivalent to
(3.6)q3=q2(β1+β2)-q1β1β2+C3(1-β3)β3(β3-β1)(β3-β2)q4=q3(β1+β2)-q2β1β2+C3(1-β3)β32(β3-β1)(β3-β2)⋮qn-1=qn-2(β1+β2)-qn-3β1β2+C3(1-β3)β3n-3(β3-β1)(β3-β2).
Finally we use the transformation qj+1-qjβ3. At this final step we get the following system:
(3.7)q4=q3(β1+β2+β3)-q2(β1β2+β1β3+β2β3)+q1β1β2β3⋮qn-1=qn-2(β1+β2+β3)-qn-3(β1β2+β1β3+β2β3)+qn-4β1β2β3.
In matrix notations we have
(3.8)Mw=Q,
where [w1,w2,w3]=[β1+β2+β3,-β1β2-β1β3-β2β3,β1β2β3],
(3.9)M=[q3q2q1q4q3q2⋮⋮⋮qn-2qn-3qn-4],Q=[q4q5⋮qn-1].
If the data values yi are noiseless, the values w1,  w2, and w3 may be obtained from the previous system of equations and the βj are the roots of the polynomial
(3.10)α(z)=z3-w1z2-w2z-w3.
Once the roots have been computed the nonlinear parameters can be calculated using (2.11), and likewise the linear parameters, as stated before.In general, for an arbitraryk, let M(k) be the (n-k-1)×k matrix
(3.11)M(k)=[qk··q1qk+1··q2⋮⋮⋮qn-2··qn-1-k],Q(k)=[qk+1qk+2⋮qn-1],w(k)=[w1(k)⋮wk(k)].
In this case the coefficients of the polynomial α(z) are the symmetric functions of β1,…,βk defined as
(3.12)w1(k)=β1+⋯+βk,w2(k)=-∑l≠r‍βlβr,w3(k)=∑l≠r,l≠s,s≠r‍βlβrβs,⋮wk-1(k)=(-1)k∑j=1k‍(∏l≠jk‍βl),wk(k)=(-1)k+1∏j=1k‍βj.
These coefficients are determinated by the solution of the system
(3.13)M(k)w(k)=Q(k).
Finally, the βj are the roots of the polynomial
(3.14)α(k)(z)=zk-∑j=1k‍wj(k)zk-j.
Next we will study the relationship between the solution obtained by the procedure described above and the Prony method described in Section 2.Theorem 3.1.
LetR be the k×k matrix defined as follows: for k=1 set R=1 and for k>1(3.15)R(i,j)={1,ifi=j,-1,ifj=i+1,0,otherwise.
In addition, let P(z) and α(k)(z) be the polynomials defined in (2.4) and (3.14), respectively. If δ=[δ1,…,δk+1,1] is the solution of the optimization problem (2.14), then w(k)=R-1[δk,…,δ1]T satisfies
(3.16)M(k)w(k)-Q(k)=XδTy.
Moreover,
(3.17)P(z)=(z-1)α(k)(z).Proof.
The solutionδ=[δ1,…,δk+1,δk+2] of (2.14) satisfies δk+2=1. In the case we are considering β0=1 is a root of P(z), which implies that
(3.18)-∑j=1k+1‍δj=1.
Then we have
(3.19)M(k)w(k)-Q(k)=M(k)R-1Rw(k)-Q(k)=M(k)R-1[δk⋮δ1]-[yk+1-yk+2⋮yn-1-yn]=M(k)R-1[δk⋮δ1]-[-[∑j=1k+1‍δj]yk+1-yk+2⋮-[∑j=1k+1‍δj]yn-1-yn]=M(k)R-1[δk⋮δ1]+[yk+1··yk+1⋮⋮⋮⋮yn-1··yn-1][δk⋮δ1]+⋯δk+2[yk+2⋮yn]+δk+1[yk+1⋮yn-1]=[y1y2⋯yk+2y2y3⋯yk+3y3y4⋯yk+4··⋯.yn-k-1yn-k⋯yn][δ1δ2⋮δk+1δk+2]=Wyδ.
Then, using (2.9), it follows that
(3.20)M(k)w(k)-Q(k)=XδTy.
Now, for the polynomial P(z) we have
(3.21)P(z)=δk+2zk+1+⋯+δ2z+δ1=(z-1)(zk-(∑j=1k‍δj)zk-1-(∑j=1k-1‍δj)zk-2-⋯-(δ1+δ2)z-δ1)=(z-1)(zk-w1(k)zk-1-w2(k)zk-2⋯-wk-1(k)z-wk(k))=(z-1)(zk-∑j=1k‍wj(k)zk-j)=(z-1)α(k)(z).Theorem 3.2.
Let us suppose that there exists a unique solution to the optimization problem (2.14). A vector δ∈ℝk+2 is the solution of the optimization problem (2.14) if and only if R-1[δk,…,δ1]T is the least squares solution of the linear equation (3.13).Proof.
Letδ∈ℝk+2 be the solution of the optimization problem (2.14), ζ=R-1[δk,…,δ1]T, and let ψ be the least squares solution of the linear system (3.13). From Theorem 3.1 it follows that
(3.22)∥XδTy∥=∥M(k)ζ-Q(k)∥≥minz∥M(k)z-Q(k)∥=∥M(k)ψ-Q(k)∥.
Let us consider ξ∈ℝk given by
(3.23)[ξk,…,ξ1]T=Rψ
and γ∈ℝk+2 defined as
(3.24)γ=[ξ1,…,ξk,-1-∑j=1k‍ξj,1]T.
By Theorem 3.1 we have
(3.25)∥M(k)ψ-Q(k)∥=∥XγTy∥,
and therefore ∥XδTy∥≥∥XγTy∥. By hypothesis δ is the only solution of the optimization problem, thus δ=γ, which implies that [δk,…,δ1]=[ξk,…,ξ1]. Now we have Rζ=Rψ and thus ζ=ψ.
In a similar way we can prove thatγ in (3.24) is the solution of the optimization problem (2.14), whenever ψ is the least squares solution of the linear equation (3.13).


## 4. Stability of the Alternative Prony Method

In (3.13), arrays M(k) and Q(k) depend on the vector of measurements y. In the presence of noisy data the equation that really holds is
(4.1)(M(k)+ϵ1)w(k)=Q(k)+ϵ2,
where ϵ1 and ϵ2 depend on the level of noise. Then it is necessary to determine the condition number of the matrix in order to see the accuracy in the coefficients given by w(k). Let us consider (1.2) for n=2k+1. In this situation M(k) is a nonsingular square matrix of dimension k×k. In the following lemma we establish a factorization of the matrix M(k) that will be used in the stability analysis of the system (3.13). Our analysis is similar to the one described in [12].Lemma 4.1.
Let us considern=2k+1. Let G and H be the k×k matrices defined by
(4.2)gi,j=βji-1,hi,j=βik-j,
and let D be the k×k diagonal matrix given by di,i=Ci(1-βi)βi. Then
(4.3)M(k)=GDH.Proof.
By definitionqi=yi-yi+1, i=1,…,2k. Then
(4.4)qi=∑j=1k‍Cjβji(1-βj).
Thus
(4.5)M(k)=[qk⋯q1·⋯··⋯·q2k-1⋯qk]=[∑j=1k‍Cjβjk(1-βj)⋯∑j=1k‍Cjβj(1-βj)·⋯··⋯·∑j=1k‍Cjβj2k-1(1-βj)⋯∑j=1k‍Cjβjk(1-βj)]=[1⋯1β1⋯βk·⋯··⋯·β1k-1⋯βkk-1][d1,1⋯·.⋯·.⋯·.⋯dk,k][β1k-1⋯β11β2k-1⋯β21·⋯··βkk-1⋯βk1]=GDH.Lemma 4.2.
Let us considern=2k+1. Let ηM(k), ηG and ηD be the condition numbers in the infinite norm of the matrices M(k), G, and D, respectively. Then
(4.6)ηM(k)≤k2ηG2ηD.Proof.
LetηH be the condition number in the infinite norm of the matrix H. Using Lemma 4.1 it follows that
(4.7)ηM(k)≤ηGηDηH.
To get (4.6) it is enough to observe that
(4.8)ηH=∥H∥∞∥H-1∥∞=∥HT∥1∥(H-1)T∥1=∥HT∥1∥(HT)-1∥1=∥G∥1∥G-1∥1≤k2∥G∥∞∥G-1∥∞=k2ηG.The following theorem establishes an estimate for the variationΔw(k) in the vector components w(k) in (3.13), as it depends on the noise level in the vector of the measurements y=[y1,…,yn]T.Theorem 4.3.
Letn=2k+1, y=[y1,…,yn]T, Δy=[Δy1,…,Δyn]T, and
(4.9)Λ1=min{∥M(k)∥∞,∥q(k)∥∞}.
Let δ>0 such that
(4.10)2k3ηG2ηDδ<Λ1.
If ∥Δy∥∞<δ, then the variation Δw(k) in the solution of the perturbed system
(4.11)(M(k)+ΔM(k))(w(k)+Δw(k))=q(k)+Δq(k)
satisfies
(4.12)∥Δw(k)∥∞∥w(k)∥∞≤4k3ηG2ηDδΛ1+O(δ2).Proof.
Letϵ=2δk/Λ1. If ∥Δy∥∞<δ, then ∥ΔM(k)∥∞≤2kδ and ∥Δq(k)∥∞≤2kδ. Therefore
(4.13)∥Δq(k)∥∞≤2kδ∥q(k)∥∞Λ1=ϵ∥q(k)∥∞.
In a similar way we get
(4.14)∥ΔM(k)∥∞≤ϵ∥M(k)∥∞.
By using (4.6) we see that
(4.15)ϵηM(k)≤ϵk2ηG2ηD<1.
From (4.13), (4.14), and (4.15), and using Theorem 2.7.2 in [13], it follows that
(4.16)∥Δw(k)∥∞∥w(k)∥∞≤2ϵηM(k)1-ϵηM(k)≤2ϵηM(k)+O(ϵ2ηM(k)2)≤4k3ηG2ηDδΛ1+O(δ2).
Once we have estimated the vector w~(k)=w(k)+Δw(k), it remains to consider the corresponding impact in the calculation of the roots of the polynomial
(4.17)α~(k)(z)=zk-∑j=1k‍w~j(k)zk-j.
Let βj be one of the roots of the polynomial α(k)(z) and suppose that β~j=βj+Δβj is the root of the polynomial α~(k)(z) closest to βj. The following is an estimate for Δβj which follows from the theory of stability for simple roots of polynomials, see [14].Lemma 4.4.
Let0<θ≪1 such that
(4.18)∥Δw(k)∥∞∥w(k)∥∞<θ.
Then there exits Λ2>0 such that
(4.19)|Δβj|≤Λ2∥Δw(k)∥∞∥w(k)∥∞.From (4.12) we see that there is a constant Λ3 such that
(4.20)∥Δw(k)∥∞∥w(k)∥∞≤(4k3ηG2ηDΛ1+Λ3δ)δ.
Thus, if δ in Theorem 4.3 is sufficiently small such that
(4.21)δ(4k3ηG2ηDΛ1+Λ3δ)=θ<1,
then the conditions of Lemma 4.4 are satisfied and therefore
(4.22)∥Δβ∥∞≤Λ2∥Δw(k)∥∞∥w(k)∥∞.
With the computed β~j,  1≤j≤k, we should analyze the least squares solution of the perturbed system
(4.23)[1β~1⋯β~k··⋯···⋯·1β~12k+1⋯β~k2k+1][b~C~1·C~k]=y+Δy,
which will be written as B~C~=y~. The following theorem provides an estimate of ∥ΔC∥2, being C~=C+ΔC the solution of the perturbed system. Note that the following result is the first one in this paper where we require the nonlinear parameters in the model (1.1) to be positive.Theorem 4.5.
Letn, y, Δy, and δ be as in the statement of Theorem 4.3, with ∥Δy∥∞<δ. Let η^β be the condition number of matrix B in the Euclidean norm and consider
(4.24)Λ4=min{∥B∥2,∥y∥2}.
Let us suppose that δ satisfies (4.10) and (4.21) and let
(4.25)ρδ=δ2k+1max{k(2k+1)(4k3ηG2ηDΛ1+Λ3δ)Λ2,1}<min{1,Λ4η^β}.
Then
(4.26)∥ΔC∥2∥C∥2≤ρδΛ4(η^βcosσ+η^β2tanσ)+O(δ2),
where
(4.27)sinσ=∥y-y*∥2∥y∥2,
and y* is the orthogonal projection of the vector y on the subspace spanned by the columns of B.Proof.
Letϵ=ρδ/Λ4. Using the mean value theorem we see that
(4.28)|(βl+Δβl)j-βlj|≤j(βl*)j-1Δβj,
with βl* between βj and βj+Δβj and 0<βj<1. Then
(4.29)∥ΔB∥∞=max1≤i≤2k+1(∑l=1k‍|(βl+Δβl)j-βlj|)≤k(2k+1)∥Δβ∥∞.
From (4.20) and (4.22) it follows that
(4.30)∥ΔB∥2≤2k+1∥ΔB∥∞≤k(2k+1)3/2∥Δβ∥∞≤k(2k+1)3/2Λ2δ(4k3ηG2ηDΛ1+Λ3δ)≤ρδΛ4∥B∥2.
On the other hand
(4.31)∥Δy∥2≤2k+1∥Δy∥∞≤2k+1δ≤ρδΛ4∥y∥2.
By definition
(4.32)ρδΛ4<1η^β,
and (4.26) follows from (4.30), (4.31), and (4.32), by Theorem 5.3.1 in [13].If the value ofρδ in (4.25) is given by
(4.33)ρδ=k(2k+1)3/2(4k3ηG2ηDΛ1+Λ3δ)Λ2δ
and sinσ=0, then
(4.34)∥ΔC∥2∥C∥2≤4k4(2k+1)3/2Λ2Λ1Λ4ηG2ηDη^βδ+O(δ2).
This estimate is similar to the is studied one obtained in reference [12], in which it is studied the stability of a confluent Prony system. Let us suppose, without loss of generality, that β1<β2⋯<βk. From [15] we see that
(4.35)ηG=max{k,g}max1≤j≤kgj,
where
(4.36)g=∑j=1k‍βjk-1,gj=∏l≠jk‍1+βl|βj-βl|.D is a diagonal matrix, therefore. We have
(4.37)ηD=maxj|Cjβj(1-βj)|minj|Cjβj(1-βj)|.
We do not know a reasonable upper bound for the condition number η^B.It follows that the method is well conditioned if each one of the following statements is satisfied: there isσ>0 such that βj≥σ for 1≤j≤k, the powers βj are widely separated, and the number of nonlinear parameters remains small. The proximity between two βj values leads to a large value of the condition number of G, ηG, and hence to the ill-conditioning of the algorithm. This characteristic is consistent with the observations of Varah, [16], who remarked the ill conditioning behaviour for any exponential fitting algorithm, provided that there are two parameters βj very close.


## 5. Numerical Results

In this section we present some numerical experiments to see the actual behavior of our implementation of the alternative Prony method described above. We have results for models with two and three exponentials.Since the application which motivated our work is related toT2 magnetic resonance brain images, we will take the linear parameters Cj to be positive and such that ∑j=1k‍Cj=255, and the nonlinear ones to lie in the interval [0.65,26], see [17].Given a particular modely(t), the data values yi are generated evaluating y(ti) with Δt=0.04, ti=iΔt and adding noise. We have used two different kinds of noise: Gaussian and Rician.In MRI, data are intrinsically complex valued and corrupted with zero mean Gaussian distributed noise with equal variance [18]. MR magnitude images are formed by simply taking the square root of the two independent Gaussian random variables (real and imaginary parts) pixel by pixel. For an MR magnitude image, the Rician probability density function of the measure pixel intensity x is given by
(5.1)p(x∣y)=xδ2e-(x2+y2)/2δ2I0(yxδ2),
where I0 is the modified zeroth-order Bessel function of the first kind, y is the underlying noise-free signal amplitude, and δ denotes the standard deviation of the Gaussian noise in the real and imaginary parts.The figures are log-log graphs, the horizontal axis corresponds to the noise, and the vertical axis corresponds to the relative errors in the computations of the linear and nonlinear parameters(5.2)db=|Δbb|,dCj=|ΔCjCj|,dλj=|Δλjλj|.
For each example, we first show the results for one run corresponding to a particular value of the noise and then we show the average value of the relative errors for one hundred runs. The Gaussian noise has zero mean, and the variance, δ, varies between 10-10 and 1. For the Rice noise we use (yi,δ), with δ varying between 10-10 and 1 as the parameters for each simulation.For Figures1, 2, 3 and 4 we used the model
(5.3)y(t)=10+80e-11t+165e-20t.
Note the linear dependence of the errors with respect to level of noise. We also noticed that, in general, the errors in the nonlinear parameters are slightly smaller than the errors in the linear parameters. In Figures 1 and 2 we have the results using Gaussian noise. Figures 3 and 4 show the results for the same model and Rice noise; as the behaviour is similar for Gaussian and Rice noises, we will only consider Rice noise for the rest of the experiments.Figure 1
Relative errors for the first model corresponding to Gaussian noise.Figure 2
Relative error average, after 100 runs, for the first model corresponding to Gaussian noise.Figure 3
Relative errors for the first model corresponding to Rice noise.Figure 4
Relative error average, after 100 runs, for the first model corresponding to Rice noise.Next we present the results for the model(5.4)y(t)=5+100e-11t+150e-12t.
In Figures 5 and 6 the graphs show a deterioration, as the theory suggests, because the nonlinear parameters are closer than in the previous example. In this case the errors in the linear parameters are greater, for more than one order of magnitude, than the errors on the nonlinear parameters.Figure 5
Relative errors for the second model corresponding to Rice noise.Figure 6
Relative error average, after 100 runs, for the second model corresponding to Rice noise.Finally, we use the model(5.5)y(t)=5+70e-2t+85e-10t+95e-18t.
The results are shown in Figures 7 and 8. Again, there is a deterioration, with respect to the first example, which is caused by the addition of one term to the exponential model.Figure 7
Relative errors for the third model corresponding to Rice noise.Figure 8
Relative error average, after 100 runs, for the third model corresponding to Rice noise.Now we present Table1 with the computed condition numbers for the different relevant matrices for each one of the three examples.Table 1
Condition numbers.
η
G
η
D
η
^
B
Example 1
16.88724660783
2.22604845378
53.04339947209
Example 2
130.20508573840
1.54342570227
397.12140292138
Example 3
184.8409648621
4.7771637988
269.0372904019


## 6. Conclusions

In this work we have developed a stability analysis for the alternative formulation of the Prony method presented by Martín-Landrove et al. [1, 2]. The analysis shows the linear dependence between the perturbation of the data and the relative errors in the computed values for the model linear and nonlinear parameters. It is also shown that the errors in the linear parameters depend upon both the number and the closeness of the nonlinear parameters.



---

*Source: 390645-2012-11-19.xml*