Stability Analysis of a Variant of the Prony Method

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 of 𝑇₂-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.

1. Introduction

In this paper we consider the problem of recovering the parameters 𝑏, 𝜆_𝑗, and 𝐶_𝑗  𝑗=1,…,𝑘, from a set of measurements 𝑦_𝑖, 𝑖=1,…,𝑛 with 𝑛≥2𝑘+1 and a given exponential model 𝜇(𝑡)=𝑏+^𝑘_𝑗=1𝐶_𝑗𝑒^{−𝜆𝑗𝑡}. (1.1) The measurement 𝑦_𝑖 corresponds to the value of expression (1.1) for 𝑡=𝑖Δ𝑡, Δ𝑡>0. Hence we get the set of equations 𝑦_𝑖=𝑏+^𝑘_𝑗=1𝐶_𝑗𝑒^{−𝑖𝜆𝑗Δ𝑡}𝑖=1,…,𝑛. (1.2) This kind of exponential fit appears in the case of 𝑇₂-weighted magnetic resonance images of the brain. 𝑇₂ 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; 𝑇₂ weighted refers to a signal produced with long TR and TE.

The exponents 𝜆_𝑗 are the nonlinear parameters and correspond to the relaxation rates associated with the different tissues present in the images. The coefficients 𝑏,𝐶₁,…,𝐶_𝑘 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 𝑘=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 𝑏=𝐶₀ and 𝜆₀=0, and write 𝜇(𝑡)=^𝑘_𝑗=0𝐶_𝑗𝑒^{−𝜆𝑗𝑡}. (2.1) The Prony type methods, also known as polynomial methods, use the fact that 𝜇(𝑡) satisfies a difference equation of the form 𝛿_𝑘+2𝐸^𝑘+1+⋯+𝛿₂𝐸+𝛿₁𝜇(𝑡)=0, (2.2) where 𝐸 is the forward shift operator 𝐸𝜇(𝑡)=𝜇(𝑡+Δ𝑡), (2.3) and the values 𝛽_𝑗=𝑒^{−𝜆𝑗Δ𝑡} are the roots of the polynomial 𝑃(𝑧)=𝛿_𝑘+2𝑧^𝑘+1+⋯+𝛿₂𝑧+𝛿₁=0, (2.4) which is called the characteristic polynomial of the associated difference equation (2.2). Evaluating (2.2) for 𝑡_𝑖=𝑖Δ𝑡,𝑖=1,…,𝑛, we get the following set of linear equations: 𝛿_𝑘+2𝜇𝑡_𝑘+2+⋯+𝛿₁𝜇𝑡₁𝛿=0⋮⋮_𝑘+2𝜇𝑡_𝑛+⋯+𝛿₁𝜇𝑡_{𝑛−𝑘−1}=0. (2.5) Let 𝛿=[𝛿₁,…𝛿_𝑘+2]^𝑇, 𝑣=[𝑣₁,…,𝑣_𝑛]^𝑇=[𝜇(𝑡₁),…,𝜇(𝑡_𝑛)]^𝑇 and 𝑋_𝛿=𝛿₁𝛿⋮⋱₁𝛿_𝑘+2𝛿⋱⋮_𝑘+2. (2.6) Then the above system of linear equations can be rewritten in matrix form as 𝑋^𝑇_𝛿𝑣=0. (2.7) Alternatively, if 𝑊_𝑣 is the (𝑛−𝑘−1)×(𝑘+2) Hankel matrix 𝑊_𝑣=𝑣₁𝑣₂⋯𝑣_𝑘+2𝑣⋅⋅⋯⋅⋅⋅⋯⋅⋅⋅⋯⋅_{𝑛−𝑘−1}⋅⋯𝑣_𝑛 (2.8) then 𝑊_𝑣𝛿=𝑋^𝑇_𝛿𝑣. (2.9) If there is no noise in the observation data, then 𝑦_𝑖=𝜇(𝑡_𝑖) and the coefficients 𝛿_𝑖 can be determined from the equivalent system of equations 𝑋^𝑇_𝛿𝑦=0, (2.10) where 𝑦=[𝑦₁,…𝑦_𝑛]^𝑇. Then the 𝛽_𝑗 are computed as the roots of 𝑃(𝑧) and finally 𝜆_𝑗=−log𝛽_𝑗Δ𝑡𝑗=0,…,𝑘. (2.11) In the presence of noisy data instead of solving the system (2.10), we consider the nonlinear optimization problem min_𝛿𝑦^𝑇𝑋_𝛿𝑋^𝑇_𝛿𝑦. (2.12) In order to obtain a nontrivial solution, it is necessary to impose restrictions over the parameters 𝛿_𝑗; 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 min_𝛿𝑦^𝑇𝑋_𝛿𝑋⁺_𝛿𝑦subjectto‖𝛿‖²₂=1, (2.13) where 𝑋⁺_𝛿 is the Moore Penrose generalized inverse of 𝑋_𝛿. Osborne et al. considered 𝑋_𝛿 to be full rank, hence 𝑋⁺_𝛿=(𝑋^𝑇_𝛿𝑋_𝛿)⁻¹𝑋^𝑇_𝛿.

We consider a different optimization problem, which appears in the literature, to compare with the Martín-Landrove method: min_𝛿𝑦^𝑇𝑋_𝛿𝑋^𝑇_𝛿𝑦subjectto𝛿_𝑘+2=1. (2.14) 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 case 𝑘=3 and 𝑛≥7. Then, the system of equations (1.2) can be written as 𝑦₁=𝑏+𝐶₁𝛽₁+𝐶₂𝛽₂+𝐶₃𝛽₃𝑦₂=𝑏+𝐶₁𝛽²₁+𝐶₂𝛽²₂+𝐶₃𝛽²₃⋮𝑦_𝑛=𝑏+𝐶₁𝛽^𝑛₁+𝐶₂𝛽^𝑛₂+𝐶₃𝛽^𝑛₃. (3.1) From (3.1) and defining 𝑞_𝑖=𝑦_𝑖−𝑦_𝑖+1 we get 𝑞₁=𝐶₁1−𝛽₁𝛽₁+𝐶₂1−𝛽₂𝛽₂+𝐶₃1−𝛽₃𝛽₃⋮𝑞_𝑛−1=𝐶₁1−𝛽₁𝛽₁^𝑛−1+𝐶₂1−𝛽₂𝛽₂^𝑛−1+𝐶₃1−𝛽₃𝛽₃^𝑛−1. (3.2) The dimension of the system can be reduced by using the transformation 𝑞_𝑗+1−𝑞_𝑗𝛽₁, to get 𝑞₂−𝑞₁𝛽₁=𝐶₂1−𝛽₂𝛽₂𝛽₂−𝛽₁+𝐶₃1−𝛽₃𝛽₃𝛽₃−𝛽₁⋮𝑞_𝑛−1−𝑞_𝑛−2𝛽₁=𝐶₂1−𝛽₂𝛽₂^𝑛−2𝛽₂−𝛽₁+𝐶₃1−𝛽₃𝛽₃^𝑛−2𝛽₃−𝛽₁. (3.3) Then 𝑞₂=𝑞₁𝛽₁+𝐶₂1−𝛽₂𝛽₂𝛽₂−𝛽₁+𝐶₃1−𝛽₃𝛽₃𝛽₃−𝛽₁⋮𝑞_𝑛−1=𝑞_𝑛−2𝛽₁+𝐶₂1−𝛽₂𝛽₂^𝑛−2𝛽₂−𝛽₁+𝐶₃1−𝛽₃𝛽₃^𝑛−2𝛽₃−𝛽₁. (3.4) We now apply the transformation 𝑞_𝑗+1−𝑞_𝑗𝛽₂, and the new set can be written as 𝑞₃−𝑞₂𝛽₂=𝑞₂𝛽₁−𝑞₁𝛽₁𝛽₂+𝐶₃1−𝛽₃𝛽₃𝛽₃−𝛽₁𝛽₃−𝛽₂⋮𝑞_𝑛−1−𝑞_𝑛−2𝛽₂=𝑞_𝑛−2𝛽₁−𝑞_𝑛−3𝛽₁𝛽₂+𝐶₃1−𝛽₃𝛽₃^𝑛−3𝛽₃−𝛽₁𝛽₃−𝛽₂. (3.5) The equations above are equivalent to 𝑞₃=𝑞₂𝛽₁+𝛽₂−𝑞₁𝛽₁𝛽₂+𝐶₃1−𝛽₃𝛽₃𝛽₃−𝛽₁𝛽₃−𝛽₂𝑞₄=𝑞₃𝛽₁+𝛽₂−𝑞₂𝛽₁𝛽₂+𝐶₃1−𝛽₃𝛽²₃𝛽₃−𝛽₁𝛽₃−𝛽₂⋮𝑞_𝑛−1=𝑞_𝑛−2𝛽₁+𝛽₂−𝑞_𝑛−3𝛽₁𝛽₂+𝐶₃1−𝛽₃𝛽₃^𝑛−3𝛽₃−𝛽₁𝛽₃−𝛽₂. (3.6) Finally we use the transformation 𝑞_𝑗+1−𝑞_𝑗𝛽₃. At this final step we get the following system: 𝑞₄=𝑞₃𝛽₁+𝛽₂+𝛽₃−𝑞₂𝛽₁𝛽₂+𝛽₁𝛽₃+𝛽₂𝛽₃+𝑞₁𝛽₁𝛽₂𝛽₃⋮𝑞_𝑛−1=𝑞_𝑛−2𝛽₁+𝛽₂+𝛽₃−𝑞_𝑛−3𝛽₁𝛽₂+𝛽₁𝛽₃+𝛽₂𝛽₃+𝑞_𝑛−4𝛽₁𝛽₂𝛽₃. (3.7) In matrix notations we have 𝑀𝑤=𝑄, (3.8) where [𝑤₁,𝑤₂,𝑤₃]=[𝛽₁+𝛽₂+𝛽₃,−𝛽₁𝛽₂−𝛽₁𝛽₃−𝛽₂𝛽₃,𝛽₁𝛽₂𝛽₃], 𝑞𝑀=³𝑞₂𝑞₁𝑞₄𝑞₃𝑞₂𝑞⋮⋮⋮_𝑛−2𝑞_𝑛−3𝑞_𝑛−4𝑞,𝑄=⁴𝑞₅⋮𝑞_𝑛−1. (3.9) If the data values 𝑦_𝑖 are noiseless, the values 𝑤₁,  𝑤₂, and 𝑤₃ may be obtained from the previous system of equations and the 𝛽_𝑗 are the roots of the polynomial 𝛼(𝑧)=𝑧³−𝑤₁𝑧²−𝑤₂𝑧−𝑤₃. (3.10) 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 arbitrary 𝑘, let 𝑀^(𝑘) be the (𝑛−𝑘−1)×𝑘 matrix 𝑀^(𝑘)=𝑞_𝑘⋅⋅𝑞₁𝑞_𝑘+1⋅⋅𝑞₂𝑞⋮⋮⋮_𝑛−2⋅⋅𝑞_{𝑛−1−𝑘},𝑄^(𝑘)=𝑞_𝑘+1𝑞_𝑘+2⋮𝑞_𝑛−1,𝑤^(𝑘)=𝑤₁^(𝑘)⋮𝑤_𝑘^(𝑘). (3.11) In this case the coefficients of the polynomial 𝛼(𝑧) are the symmetric functions of 𝛽₁,…,𝛽_𝑘 defined as 𝑤₁^(𝑘)=𝛽₁+⋯+𝛽_𝑘,𝑤₂^(𝑘)=−_𝑙≠𝑟𝛽_𝑙𝛽_𝑟,𝑤₃^(𝑘)=_{𝑙≠𝑟,𝑙≠𝑠,𝑠≠𝑟}𝛽_𝑙𝛽_𝑟𝛽_𝑠,⋮𝑤^(𝑘)_𝑘−1=(−1)^𝑘𝑘_𝑗=1^𝑘_𝑙≠𝑗𝛽_𝑙,𝑤_𝑘^(𝑘)=(−1)^𝑘𝑘+1_𝑗=1𝛽_𝑗. (3.12) These coefficients are determinated by the solution of the system 𝑀^(𝑘)𝑤^(𝑘)=𝑄^(𝑘). (3.13) Finally, the 𝛽_𝑗 are the roots of the polynomial 𝛼^(𝑘)(𝑧)=𝑧^𝑘−^𝑘_𝑗=1𝑤_𝑗^(𝑘)𝑧^𝑘−𝑗. (3.14) 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. Let 𝑅 be the 𝑘×𝑘 matrix defined as follows: for 𝑘=1 set 𝑅=1 and for 𝑘>1𝑅𝑖,𝑗=1,𝑖𝑓𝑖=𝑗,−1,𝑖𝑓𝑗=𝑖+1,0,𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒. (3.15) In addition, let 𝑃(𝑧) and 𝛼^(𝑘)(𝑧) be the polynomials defined in (2.4) and (3.14), respectively. If 𝛿=[𝛿₁,…,𝛿_𝑘+1,1] is the solution of the optimization problem (2.14), then 𝑤^(𝑘)=𝑅⁻¹[𝛿_𝑘,…,𝛿₁]^𝑇 satisfies 𝑀^(𝑘)𝑤^(𝑘)−𝑄^(𝑘)=𝑋^𝑇_𝛿𝑦. (3.16) Moreover, 𝑃(𝑧)=(𝑧−1)𝛼^(𝑘)(𝑧). (3.17)

Proof. The solution 𝛿=[𝛿₁,…,𝛿_𝑘+1,𝛿_𝑘+2] of (2.14) satisfies 𝛿_𝑘+2=1. In the case we are considering 𝛽₀=1 is a root of 𝑃(𝑧), which implies that −^𝑘+1_𝑗=1𝛿_𝑗=1. (3.18) Then we have 𝑀^(𝑘)𝑤^(𝑘)−𝑄^(𝑘)=𝑀^(𝑘)𝑅⁻¹𝑅𝑤^(𝑘)−𝑄^(𝑘)=𝑀^(𝑘)𝑅⁻¹𝛿_𝑘⋮𝛿₁−𝑦_𝑘+1−𝑦_𝑘+2⋮𝑦_𝑛−1−𝑦_𝑛=𝑀^(𝑘)𝑅⁻¹𝛿_𝑘⋮𝛿₁−−∑^𝑘+1_𝑗=1𝛿_𝑗𝑦_𝑘+1−𝑦_𝑘+2⋮−∑^𝑘+1_𝑗=1𝛿_𝑗𝑦_𝑛−1−𝑦_𝑛=𝑀^(𝑘)𝑅⁻¹𝛿_𝑘⋮𝛿₁+𝑦_𝑘+1⋅⋅𝑦_𝑘+1𝑦⋮⋮⋮⋮_𝑛−1⋅⋅𝑦_𝑛−1𝛿_𝑘⋮𝛿₁+⋯𝛿_𝑘+2𝑦_𝑘+2⋮𝑦_𝑛+𝛿_𝑘+1𝑦_𝑘+1⋮𝑦_𝑛−1=𝑦₁𝑦₂⋯𝑦_𝑘+2𝑦₂𝑦₃⋯𝑦_𝑘+3𝑦₃𝑦₄⋯𝑦_𝑘+4𝑦⋅⋅⋯._{𝑛−𝑘−1}𝑦_𝑛−𝑘⋯𝑦_𝑛𝛿₁𝛿₂⋮𝛿_𝑘+1𝛿_𝑘+2=𝑊_𝑦𝛿. (3.19) Then, using (2.9), it follows that 𝑀^(𝑘)𝑤^(𝑘)−𝑄^(𝑘)=𝑋^𝑇_𝛿𝑦. (3.20) Now, for the polynomial 𝑃(𝑧) we have 𝑃(𝑧)=𝛿_𝑘+2𝑧^𝑘+1+⋯+𝛿₂𝑧+𝛿₁=(𝑧−1)𝑧^𝑘−^𝑘_𝑗=1𝛿_𝑗𝑧^𝑘−1−^𝑘−1_𝑗=1𝛿_𝑗𝑧^𝑘−2−⋯−𝛿₁+𝛿₂𝑧−𝛿₁=(𝑧−1)𝑧^𝑘−𝑤₁^(𝑘)𝑧^𝑘−1−𝑤₂^(𝑘)𝑧^𝑘−2⋯−𝑤^(𝑘)_𝑘−1𝑧−𝑤_𝑘^(𝑘)=(𝑧−1)𝑧^𝑘−^𝑘_𝑗=1𝑤_𝑗^(𝑘)𝑧^𝑘−𝑗=(𝑧−1)𝛼^(𝑘)(𝑧). (3.21)

Theorem 3.2. Let us suppose that there exists a unique solution to the optimization problem (2.14). A vector 𝛿∈ℝ^𝑘+2 is the solution of the optimization problem (2.14) if and only if 𝑅⁻¹[𝛿_𝑘,…,𝛿₁]^𝑇 is the least squares solution of the linear equation (3.13).

Proof. Let 𝛿∈ℝ^𝑘+2 be the solution of the optimization problem (2.14), 𝜁=𝑅⁻¹[𝛿_𝑘,…,𝛿₁]^𝑇, and let 𝜓 be the least squares solution of the linear system (3.13). From Theorem 3.1 it follows that 𝑋^𝑇_𝛿𝑦=𝑀^(𝑘)𝜁−𝑄^(𝑘)≥min_𝑧𝑀^(𝑘)𝑧−𝑄^(𝑘)=𝑀^(𝑘)𝜓−𝑄^(𝑘). (3.22) Let us consider 𝜉∈ℝ^𝑘 given by 𝜉_𝑘,…,𝜉₁^𝑇=𝑅𝜓 (3.23) and 𝛾∈ℝ^𝑘+2 defined as 𝛾=𝜉₁,…,𝜉_𝑘,−1−^𝑘_𝑗=1𝜉_𝑗,1^𝑇. (3.24) By Theorem 3.1 we have 𝑀^(𝑘)𝜓−𝑄^(𝑘)=𝑋^𝑇_𝛾𝑦, (3.25) and therefore ‖𝑋^𝑇_𝛿𝑦‖≥‖𝑋^𝑇_𝛾𝑦‖. By hypothesis 𝛿 is the only solution of the optimization problem, thus 𝛿=𝛾, which implies that [𝛿_𝑘,…,𝛿₁]=[𝜉_𝑘,…,𝜉₁]. Now we have 𝑅𝜁=𝑅𝜓 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 𝑀^(𝑘) and 𝑄^(𝑘) depend on the vector of measurements 𝑦. In the presence of noisy data the equation that really holds is 𝑀^(𝑘)+𝜖₁𝑤^(𝑘)=𝑄^(𝑘)+𝜖₂, (4.1) where 𝜖₁ and 𝜖₂ 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 𝑤^(𝑘). Let us consider (1.2) for 𝑛=2𝑘+1. In this situation 𝑀^(𝑘) is a nonsingular square matrix of dimension 𝑘×𝑘. In the following lemma we establish a factorization of the matrix 𝑀^(𝑘) 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 consider 𝑛=2𝑘+1. Let 𝐺 and 𝐻 be the 𝑘×𝑘 matrices defined by 𝑔_𝑖,𝑗=𝛽_𝑗^𝑖−1,ℎ_𝑖,𝑗=𝛽_𝑖^𝑘−𝑗, (4.2) and let 𝐷 be the 𝑘×𝑘 diagonal matrix given by 𝑑_𝑖,𝑖=𝐶_𝑖(1−𝛽_𝑖)𝛽_𝑖. Then 𝑀^(𝑘)=𝐺𝐷𝐻. (4.3)

Proof. By definition 𝑞_𝑖=𝑦_𝑖−𝑦_𝑖+1, 𝑖=1,…,2𝑘. Then 𝑞_𝑖=^𝑘_𝑗=1𝐶_𝑗𝛽^𝑖_𝑗1−𝛽_𝑗. (4.4) Thus 𝑀^(𝑘)=𝑞_𝑘⋯𝑞₁𝑞⋅⋯⋅⋅⋯⋅_2𝑘−1⋯𝑞_𝑘=^𝑘_𝑗∑₌₁𝐶_𝑗𝛽^𝑘_𝑗1−𝛽_𝑗⋯∑^𝑘_𝑗=1𝐶_𝑗𝛽_𝑗1−𝛽_𝑗⋅⋯⋅⋅⋯⋅_𝑘𝑗∑₌₁𝐶_𝑗𝛽_𝑗^2𝑘−11−𝛽_𝑗⋯∑^𝑘_𝑗=1𝐶_𝑗𝛽^𝑘_𝑗1−𝛽_𝑗=𝛽1⋯1₁⋯𝛽_𝑘𝛽⋅⋯⋅⋅⋯⋅_1𝑘−1⋯𝛽_𝑘^𝑘−1𝑑_1,1⋯⋅.⋯⋅.⋯⋅.⋯𝑑_𝑘,𝑘𝛽₁^𝑘−1⋯𝛽₁1𝛽₂^𝑘−1⋯𝛽₂1𝛽⋅⋯⋅⋅_𝑘𝑘−1⋯𝛽_𝑘1=𝐺𝐷𝐻. (4.5)

Lemma 4.2. Let us consider 𝑛=2𝑘+1. Let 𝜂_𝑀(𝑘), 𝜂_𝐺 and 𝜂_𝐷 be the condition numbers in the infinite norm of the matrices 𝑀^(𝑘), 𝐺, and 𝐷, respectively. Then 𝜂_𝑀(𝑘)≤𝑘²𝜂²_𝐺𝜂_𝐷. (4.6)

Proof. Let 𝜂_𝐻 be the condition number in the infinite norm of the matrix 𝐻. Using Lemma 4.1 it follows that 𝜂_𝑀(𝑘)≤𝜂_𝐺𝜂_𝐷𝜂_𝐻. (4.7) To get (4.6) it is enough to observe that 𝜂_𝐻=‖𝐻‖_∞𝐻⁻¹_∞=𝐻^𝑇₁(𝐻⁻¹)^𝑇₁=𝐻^𝑇₁(𝐻^𝑇)⁻¹₁=‖𝐺‖₁𝐺⁻¹₁≤𝑘²‖𝐺‖_∞𝐺⁻¹_∞=𝑘²𝜂_𝐺. (4.8)

The following theorem establishes an estimate for the variation Δ𝑤^(𝑘) in the vector components 𝑤^(𝑘) in (3.13), as it depends on the noise level in the vector of the measurements 𝑦=[𝑦₁,…,𝑦_𝑛]^𝑇.

Theorem 4.3. Let 𝑛=2𝑘+1, 𝑦=[𝑦₁,…,𝑦_𝑛]^𝑇, Δ𝑦=[Δ𝑦₁,…,Δ𝑦_𝑛]^𝑇, and Λ₁=min𝑀^(𝑘)_∞,𝑞^(𝑘)_∞. (4.9) Let 𝛿>0 such that 2𝑘³𝜂²_𝐺𝜂_𝐷𝛿<Λ₁. (4.10) If ‖Δ𝑦‖_∞<𝛿, then the variation Δ𝑤^(𝑘) in the solution of the perturbed system 𝑀^(𝑘)+Δ𝑀^(𝑘)𝑤^(𝑘)+Δ𝑤^(𝑘)=𝑞^(𝑘)+Δ𝑞^(𝑘) (4.11) satisfies Δ𝑤^(𝑘)_∞𝑤^(𝑘)_∞≤4𝑘³𝜂²_𝐺𝜂_𝐷𝛿Λ₁+𝑂𝛿². (4.12)

Proof. Let 𝜖=2𝛿𝑘/Λ₁. If ‖Δ𝑦‖_∞<𝛿, then ‖Δ𝑀^(𝑘)‖_∞≤2𝑘𝛿 and ‖Δ𝑞^(𝑘)‖_∞≤2𝑘𝛿. Therefore Δ𝑞^(𝑘)_∞𝑞≤2𝑘𝛿^(𝑘)∞Λ₁=𝜖𝑞^(𝑘)_∞. (4.13) In a similar way we get Δ𝑀^(𝑘)_∞≤𝜖𝑀^(𝑘)_∞. (4.14) By using (4.6) we see that 𝜖𝜂_𝑀(𝑘)≤𝜖𝑘²𝜂²_𝐺𝜂_𝐷<1. (4.15) From (4.13), (4.14), and (4.15), and using Theorem 2.7.2 in [13], it follows that Δ𝑤^(𝑘)_∞𝑤^(𝑘)_∞≤2𝜖𝜂_𝑀(𝑘)1−𝜖𝜂_𝑀(𝑘)≤2𝜖𝜂_𝑀(𝑘)+𝑂𝜖²𝜂²_𝑀(𝑘)≤4𝑘³𝜂²_𝐺𝜂_𝐷𝛿Λ₁+𝑂𝛿². (4.16) Once we have estimated the vector 𝑤^(𝑘)=𝑤^(𝑘)+Δ𝑤^(𝑘), it remains to consider the corresponding impact in the calculation of the roots of the polynomial 𝛼^(𝑘)(𝑧)=𝑧^𝑘−^𝑘_𝑗=1𝑤_𝑗^(𝑘)𝑧^𝑘−𝑗. (4.17) Let 𝛽_𝑗 be one of the roots of the polynomial 𝛼^(𝑘)(𝑧) and suppose that 𝛽_𝑗=𝛽_𝑗+Δ𝛽_𝑗 is the root of the polynomial 𝛼^(𝑘)(𝑧) closest to 𝛽_𝑗. The following is an estimate for Δ𝛽_𝑗 which follows from the theory of stability for simple roots of polynomials, see [14].

Lemma 4.4. Let 0<𝜃≪1 such that Δ𝑤^(𝑘)_∞𝑤^(𝑘)_∞<𝜃. (4.18) Then there exits Λ₂>0 such that Δ𝛽_𝑗≤Λ₂Δ𝑤^(𝑘)_∞𝑤^(𝑘)_∞. (4.19)

From (4.12) we see that there is a constant Λ₃ such that Δ𝑤^(𝑘)_∞𝑤^(𝑘)_∞≤4𝑘³𝜂²_𝐺𝜂_𝐷Λ₁+Λ₃𝛿𝛿. (4.20) Thus, if 𝛿 in Theorem 4.3 is sufficiently small such that 𝛿4𝑘³𝜂²_𝐺𝜂_𝐷Λ₁+Λ₃𝛿=𝜃<1, (4.21) then the conditions of Lemma 4.4 are satisfied and therefore Δ𝛽_∞≤Λ₂Δ𝑤^(𝑘)_∞𝑤^(𝑘)_∞. (4.22) With the computed 𝛽_𝑗,  1≤𝑗≤𝑘, we should analyze the least squares solution of the perturbed system 1𝛽₁⋯𝛽_𝑘1𝛽⋅⋅⋯⋅⋅⋅⋯⋅_12𝑘+1⋯𝛽_𝑘^2𝑘+1𝑏̃𝐶₁⋅𝐶_𝑘=𝑦+Δ𝑦, (4.23) which will be written as 𝐵𝐶=𝑦. The following theorem provides an estimate of ‖Δ𝐶‖₂, being 𝐶=𝐶+Δ𝐶 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. Let 𝑛, 𝑦, Δ𝑦, and 𝛿 be as in the statement of Theorem 4.3, with ‖Δ𝑦‖_∞<𝛿. Let 𝜂_𝛽 be the condition number of matrix 𝐵 in the Euclidean norm and consider Λ₄=min‖𝐵‖₂,𝑦₂. (4.24) Let us suppose that 𝛿 satisfies (4.10) and (4.21) and let 𝜌_𝛿=𝛿√2𝑘+1max𝑘(2𝑘+1)4𝑘³𝜂²_𝐺𝜂_𝐷Λ₁+Λ₃𝛿Λ₂Λ,1<min1,⁴𝜂_𝛽. (4.25) Then ‖Δ𝐶‖₂‖𝐶‖₂≤𝜌_𝛿Λ₄𝜂_𝛽+𝜂cos𝜎^2𝛽tan𝜎+𝑂𝛿², (4.26) where sin𝜎=𝑦−𝑦^∗₂𝑦₂, (4.27) and 𝑦^∗ is the orthogonal projection of the vector 𝑦 on the subspace spanned by the columns of 𝐵.

Proof. Let 𝜖=𝜌_𝛿/Λ₄. Using the mean value theorem we see that 𝛽_𝑙+Δ𝛽_𝑙^𝑗−𝛽^𝑗_𝑙≤𝑗𝛽^∗_𝑙^𝑗−1Δ𝛽_𝑗, (4.28) with 𝛽^∗_𝑙 between 𝛽_𝑗 and 𝛽_𝑗+Δ𝛽_𝑗 and 0<𝛽_𝑗<1. Then ‖Δ𝐵‖_∞=max_{1≤𝑖≤2𝑘+1}^𝑘_𝑙=1𝛽_𝑙+Δ𝛽_𝑙^𝑗−𝛽^𝑗_𝑙≤𝑘(2𝑘+1)Δ𝛽_∞. (4.29) From (4.20) and (4.22) it follows that ‖Δ𝐵‖₂≤√2𝑘+1‖Δ𝐵‖_∞≤𝑘(2𝑘+1)^3/2Δ𝛽_∞≤𝑘(2𝑘+1)^3/2Λ₂𝛿4𝑘³𝜂²_𝐺𝜂_𝐷Λ₁+Λ₃𝛿≤𝜌_𝛿Λ₄‖𝐵‖₂. (4.30) On the other hand Δ𝑦₂≤√2𝑘+1Δ𝑦_∞≤√𝜌2𝑘+1𝛿≤^𝛿Λ₄𝑦₂. (4.31) By definition 𝜌_𝛿Λ₄<1𝜂_𝛽, (4.32) 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 𝜌_𝛿=𝑘(2𝑘+1)^3/24𝑘³𝜂²_𝐺𝜂_𝐷Λ₁+Λ₃𝛿Λ₂𝛿 (4.33) and sin𝜎=0, then ‖Δ𝐶‖₂‖𝐶‖₂≤4𝑘⁴(2𝑘+1)^3/2Λ₂Λ₁Λ₄𝜂²_𝐺𝜂_𝐷𝜂_𝛽𝛿+𝑂𝛿². (4.34) 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 𝛽₁<𝛽₂⋯<𝛽_𝑘. From [15] we see that 𝜂_𝐺=max𝑘,𝑔max_{1≤𝑗≤𝑘}𝑔_𝑗, (4.35) where 𝑔=^𝑘_𝑗=1𝛽_𝑗^𝑘−1,𝑔_𝑗=^𝑘_𝑙≠𝑗1+𝛽_𝑙𝛽_𝑗−𝛽_𝑙. (4.36)𝐷 is a diagonal matrix, therefore. We have 𝜂_𝐷=max_𝑗𝐶_𝑗𝛽_𝑗1−𝛽_𝑗min_𝑗𝐶_𝑗𝛽_𝑗1−𝛽_𝑗. (4.37) We do not know a reasonable upper bound for the condition number 𝜂_𝐵.

It follows that the method is well conditioned if each one of the following statements is satisfied: there is 𝜎>0 such that 𝛽_𝑗≥𝜎 for 1≤𝑗≤𝑘, the powers 𝛽_𝑗 are widely separated, and the number of nonlinear parameters remains small. The proximity between two 𝛽_𝑗 values leads to a large value of the condition number of 𝐺, 𝜂_𝐺, 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 𝛽_𝑗 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 to 𝑇₂ magnetic resonance brain images, we will take the linear parameters 𝐶_𝑗 to be positive and such that ∑^𝑘_𝑗=1𝐶_𝑗=255, and the nonlinear ones to lie in the interval [0.65,26], see [17].

Given a particular model 𝑦(𝑡), the data values 𝑦_𝑖 are generated evaluating 𝑦(𝑡_𝑖) with Δ𝑡=0.04, 𝑡_𝑖=𝑖Δ𝑡 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 𝑥 is given by 𝑝=𝑥𝑥∣𝑦𝛿²𝑒^{−(𝑥2+𝑦2)/2𝛿2}𝐼₀𝑦𝑥𝛿², (5.1) where 𝐼₀ is the modified zeroth-order Bessel function of the first kind, 𝑦 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) 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⁻¹⁰ and 1. For the Rice noise we use (𝑦_𝑖,𝛿), with 𝛿 varying between 10⁻¹⁰ and 1 as the parameters for each simulation.

For Figures 1, 2, 3 and 4 we used the model 𝑦(𝑡)=10+80𝑒^−11𝑡+165𝑒^−20𝑡. (5.3) 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+100𝑒^−11𝑡+150𝑒^−12𝑡. (5.4) 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+70𝑒^−2𝑡+85𝑒^−10𝑡+95𝑒^−18𝑡. (5.5) 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 Table 1 with the computed condition numbers for the different relevant matrices for each one of the three examples.

Table 1 

Condition numbers.

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.

Acknowledgment

This research was partially supported by Vicerrectoría de Investigación de la Universidad Nacional de Colombia, grant numbers 12427 and 16084.