A Critical Evaluation and Modification of the Padé–Laplace Method for Deconvolution of Viscoelastic Spectra

This paper shows that using the Padé–Laplace (PL) method for deconvolution of multi-exponential functions (stress relaxation of polymers) can produce ill-conditioned systems of equations. Analysis of different sets of generated data points from known multi-exponential functions indicates that by increasing the level of Padé approximants, the condition number of a matrix whose entries are coefficients of a Taylor series in the Laplace space grows rapidly. When higher levels of Padé approximants need to be computed to achieve stable modes for separation of exponentials, the problem of generating matrices with large condition numbers becomes more pronounced. The analysis in this paper discusses the origin of ill-posedness of the PL method and it was shown that ill-posedness may be regularized by reconstructing the system of equations and using singular value decomposition (SVD) for computation of the Padé table. Moreover, it is shown that after regularization, the PL method can deconvolute the exponential decays even when the input parameter of the method is out of its optimal range.


Introduction
The representation of a function as the sum of exponential decays can be observed in many areas of science and engineering such as solid-state physics, chemical kinetics, biology, and rheology of polymers [1]. For example, in relaxation phenomenon of polymers, the experimental data can be mathematically described as a linear combination of multiple exponential functions. The parameters involved in the linear combination of exponential functions (amplitudes and decay constants) have physical significance, and to extract these constants, an inverse problem that is inherently ill-posed [2] must be solved. In the relaxation phenomena of polymers, the decay constants are related to the relaxation times and the amplitudes are the corresponding weights. The set of relaxation times and their corresponding weights represents the discrete relaxation spectrum that is considered as the fingerprint of each polymer and can be used in formulation of constitutive equations and in prediction of rheological properties of polymers such as zero shear viscosity, dynamic moduli, and dynamic viscosities. The relaxation spectra of polymers are related to the dynamics of polymer chains, which highly depends on the molecular characteristics of polymer chains such as their chemical structure, polymer chain architecture, molecular weight, and molecular weight distribution. It is important to note that the number of exponential modes is not known a priori and must be determined. In [1] and references therein, there are several proposed approaches for solving the problem of separation of exponentials. Additionally, a detailed review of some other numerical techniques of exponential analysis can be found in [3]. Among the numerical procedures presented for multi-exponential analysis, the Padé-Laplace (PL) method developed by Yeramian and Claverie [4] can deconvolute exponential decays without using initial guesses for the constants [3]. In addition to that, the number of exponential modes is an outcome of the numerical procedure and only one input parameter is required to extract the exponential modes [4]. The PL method combines the Laplace transform and Padé approximation to address the ill-posed problem of separation of exponentials, and its key step in deconvolution of exponentials is the computation of Padé approximants [4], for which Yeramian and Claverie [4]. used the algorithm proposed by Longman [5]. However, as described by Tang and Norris [6], the Longman algorithm can be unstable, especially when the first few data points are close to zero. In a response to Tang and Norris [6], Yeramian [7] indicated that the Longman method is just an efficient numerical tool that allows one to calculate the coefficients of the Padé approximants. More importantly, Yeramian [7] mentioned that "To solve the linear system equations corresponding to each Padé approximant we may use simple determinants".
The PL method has been used in deconvolution of the relaxation spectra of polymers [8][9][10][11][12], but the points raised by Tang and Norris [6] regarding the instability of the Longman method and in general the problematic computation of Padé approximants have been, so far, overlooked. Knowing the fact that the computation of Padé approximants is a crucial step in PL numerical procedure, an analysis of the PL method, specifically the computation of Padé approximants and the possibility of ill-posedness caused by this computation, is still lacking. Therefore, the objective of this paper was to analyze the PL method in separating exponential functions with an emphasis on the computation of Padé approximants in PL theory.
The paper is organized as follows: The second section gives a brief overview of the Padé-Laplace theory, discusses the numerical implementation, and presents some indications of the propensity for ill-posedness. The third section analyzes the computation of Padé approximants from discrete data points generated by known multi-exponential functions. In this section, we show that the PL method can create a system of equations with large condition numbers. In section four, a method for regularization of the PL method is used and we demonstrate its success in several examples. The paper closes with the conclusions of the results.

The Padé-Laplace Theory
This section gives a brief overview of the PL method and presents important issues related to this numerical procedure. The complete details of the PL theory can be found in [4].
As mentioned earlier, the PL method uses a combination of Laplace transform and Padé approximants to extract the number of exponential modes, amplitudes, and decay constants associated with each mode, whose summation results in given experimental data obtained in discrete time intervals. In other words, the experimental data f (t) may be expressed as where n, α k , and β k are number of modes, amplitudes, and decay constants of mode k, respectively. In general, the exponents can be complex numbers. Applying the Laplace transform to Equation (1) gives where Re(p) > sup k [Re(β k )] and F(p) is In Equation (2), amplitudes appear as the residues of F(p) and exponents are the poles. The first step of the PL method is to express F(p) as a polynomial function using a Taylor series expansion about some point p 0 truncated to some order K such as where (d k F(p)/dp k )(p) is given by Since the values of f (t) are known at discrete time intervals, c k can be calculated by numerical integration.
The second step of the PL method is to construct the Padé approximant of the polynomial found by the Taylor expansion, Equation (4). In other words, the polynomial function should be expressed by the division of two polynomials as All the papers that implemented the PL method [4,[8][9][10][11][12] considered the condition b 0 = 1 for the construction of the Padé approximants, Equation (6), which can result in ill-posedness of the PL method. In effect, as will be demonstrated in this section, after applying b 0 = 1, a system of equations will be extracted from Equation (6) that, by performing analysis and numerical computations later in the paper, will show that these equations are ill-conditioned.
Upon constructing the Padé approximants, α k and β k can be extracted by finding the poles and residues of the Padé approximants in a comparison between the right-hand sides of Equations (6) and (2).
The only input parameter of this method is p 0 , that is, the point over which the Taylor series is expanded. Theoretically, the results of the computation must be independent of p 0 ; however, this is not the case in reality. According to previous papers [4,[8][9][10][11][12], the problem with some values of p 0 is attributed to the round-off errors. Therefore, it was suggested that p 0 must be chosen in an optimal range. Aubard et al. [13] proposed an optimal range as a rule of thumb in the interval between the largest and smallest values of absolute values of β k , i.e., [inf k |β k |, sup k |β k |]. Hereafter, by optimal range for a function, we suggest the optimal range proposed by Aubard et al. [13] As a practical supposition, Hellen suggested that a good choice for p 0 can be the inverse of the time that it takes for the data points to decay to the half of the initial value [14]. However, later we will show that it is possible to achieve the results even when p 0 is out of the optimal range. After Taylor expansion in Laplace space, by increasing the level of the Padé approximant, poles and residues are calculated at each Padé level and then the stable modes are identified. The stable modes are the modes that appear in Padé table and remain in the subsequent Padé levels. The number of stable modes determines the number of exponential modes. Therefore, the number of modes is an outcome of the PL numerical procedure.
To construct the Padé approximant from Equation (6), a system of linear equations should be solved. Considering b 0 = 1, and expanding the summations in Equation (6), one arrives at Multiplying both sides by the denominator of the right-hand side and comparing the terms with identical powers of (p − p 0 ), we find (see Appendix A) The third part of Equation (8) can be expanded as and this system of equations must be solved for b values. After doing so, the a k values will be calculated from the second part of Equation (8). Solving the system of equation given by Equation (9) is the major part of the PL method. The coefficient matrix in this system of equations, Equation (9), is a Toeplitz matrix whose entries are the coefficients of the Taylor expansion calculated in the Laplace space. In the next section, we will show that the Toeplitz matrices that appear in computation of Padé approximants are close to singular. In effect, it will be shown that for different sets of synthesized data points from known exponential functions the condition numbers of the Toeplitz matrices become quite large, which results in ill-conditioned systems of equations.
The level of the Padé approximants determines the size of the Toeplitz matrix. For example, when n = 3, the Padé approximant, which is shown by [2/3], is and is related to a 3 × 3 Toeplitz matrix Therefore, to achieve higher levels of Padé approximants, the size of the Toeplitz matrix must increase.

Ill-Posedness of the PL Method
The condition number of coefficient matrix in a system of linear equations is the most important indicator in analyzing the stability of computations and numerical sensitivity [15]. In Section 2, we explained that to extract the Padé coefficients, a system of linear equations, whose coefficient matrix is a Toeplitz matrix, must be solved. To show the structure of the Toeplitz matrix in the deconvolution process of the sum of exponential functions using the PL method, we consider different sets of data points generated from known exponential functions and use the procedure explained in Section 2 to construct the Toeplitz matrix for each Padé level. For all the calculations in this paper, we used the trapezoidal rule for the numerical integration of Equation (5). Table 1 shows the Toeplitz matrices for the different Padé levels calculated for a three-component function f (t), which is 2-norm condition number = 6.85 × 10 10 Infinity-norm condition number = 8.05 × 10 10 In Table 1, the Toeplitz matrices corresponding to each Padé level and their 2-norm and infinity-norm condition numbers are given. The coefficient matrices, given by Equation (9), associated with each Padé level in Table 1 are given by As one can observe, by increasing the Padé levels, the condition numbers grow rapidly. Therefore, the systems of equations become ill-conditioned. For computations in Table 1, p 0 was chosen in the optimal range. The very high condition number is a manifestation of the unstable and inaccurate numerical calculations. Table 2 shows the Toeplitz matrices for the same function, Equation (12), calculated by a different value of p 0 still in the optimal range. However, as the variation of p 0 changes the condition numbers, the condition numbers are still quite large, which is the sign of ill-conditioned systems of equations.  For the next example, we consider f (t) that is made of two exponential decays Aubard et al. [13] considered that p 0 = 0.0043 for this two-component function is in the optimal range. Toeplitz matrices for different Padé levels of Equation (14) are given in Table 3. Like the computation results for the three-component function, Equation (12), although p 0 was chosen in the optimal range, the condition numbers of the Toeplitz matrices are still quite large.  Table 4 shows the Toeplitz matrices and their condition numbers for the same function at different Padé levels, where p 0 is out of the optimal range. Interestingly, the condition numbers of coefficient matrices for the p 0 value that is out of the optimal range are still large, while they are smaller than condition numbers calculated with p 0 in the optimal range. By increasing the number of exponential modes, higher levels of Padé approximants should be calculated to reach the stable modes and, thus, by increasing the size of Toeplitz matrix, the condition numbers become very large. For example, consider the function f (t) that consists of five exponential decays   As mentioned earlier, the objective of the PL theory is to deconvolute the exponential decays from discrete experimental data f (t) measured at different time intervals. Thus, the integration by Equation (5) must be conducted numerically. To further analyze the PL theory, the entries of the coefficient matrix will be expressed in terms of the parameters of the problem. In other words, the coefficient matrix will be constructed for a general case where the function f (t) is considered in the form given by Equation (1) and expresses the entries of the coefficient matrix in terms of α i and β i . Using Equations (4) and (5), the entries of the coefficient matrix in Equation (9) are given by

Padé Approximant 2-Norm Condition Number Infinity-Norm Condition Number
Now, by plugging f (t) from Equation (1), in the form of Equation (16), we arrive at where the conditions Re(p 0 ) > Re(β i ) for all values of β i and Re(k) > −1 must be satisfied. Thus, for [1/2], [2/3], and [3/4] Padé levels the coefficient matrix will be presented as respectively. The determinants of the coefficient matrices for Padé levels [1/2], [2/3], and [3/4] when the number of exponential decays in f (t) are 2, 3, and 4, respectively, are given by As one can observe, the coefficient matrix becomes singular if the number of exponential decays in f (t) is less than the number of the Padé level. In effect, we have This analysis indicates that achieving the higher levels of Padé levels can result in singular Toeplitz matrices. Moreover, the calculations show that the coefficient matrices become rank-deficient, and, in fact, the rank of the coefficient matrix equals the number of exponential decays in function f (t). For example, we have It means that the coefficient matrix for each Padé level higher than the number of exponential decays is rank-deficient, which is consistent with the determinants of the coefficient matrices.
The results presented herein indicate that PL is an ill-posed method for the separation of exponential functions. Therefore, in contrast to the previous thought that the method is taking advantage of the properties of analyticity of Laplace transform to deal with the ill-posed problem of separation of exponentials, [9] the computation of Padé approximants generates ill-conditioned systems of equations. Although the PL method proposes a powerful numerical procedure for deconvolution of the exponential modes, it is likely to encounter ill-conditioned problems attributable to the ill-posedness of Padé table computations. In the next section, we will exploit a numerical algorithm to regularize the PL method and demonstrate that this algorithm can successfully resolve the ill-posedness of the PL numerical procedure.

Regularization of the PL Method
To resolve the ill-posedness of the PL numerical procedure, one must regularize this method. As mentioned in Section 2, after considering the condition b 0 = 1, Equation (6) will be expressed as an ill-conditioned system of linear equations. It is important to note that the Longman algorithm [5] also implements the same condition b 0 = 1 (see Equation (7) in [5]). Therefore, the instability reported by Tang and Norris [6] for using the Longman algorithm might be attributed to using the same coefficient conditions.
Knowing the origin of ill-posedness, one may regularize the PL method by changing the coefficient condition b 0 = 1, and reconstruct the system of equations in a way that eliminates the ill-conditioning.
Recently, Gonnet et al. [16] proposed a numerical algorithm for computation of the Padé table using singular value decomposition (SVD). In this numerical algorithm, instead of using the coefficient condition b 0 = 1, they considered where b 0 b 1 · · · b n T is a vector whose components are the coefficients of the polynomial in the denominator of the Padé approximant and · 2 is the vector 2-norm operator. After computations, the output of the numerical algorithm presents a polynomial in the denominator of the Padé approximants in the form of Using the normalization condition, Equation (30), in computation of the Padé table helps to eliminate the ill-conditioning problem. In effect, the numerical algorithm proposed by Gonnet et al. [16] can regularize the ill-posed problem of Padé computations. Hereafter, we call the regularized PL method RPL, which represents the PL method where Padé table is computed using the SVD solver developed by Gonnet et al. [16]. In the following, we show the capability of RPL in separation of exponentials. Table 6 shows the results of the deconvolution of generated data points of the function f (t) in Equation (12) using RPL. As shown in Table 6, three stable modes appear in [4/5] and [5/6] Padé approximants. It is expected that these stable modes remain in the computations after increasing the Padé levels. On the other hand, considering the rapid growth of condition number of Toeplitz matrices, as shown earlier, achieving high levels of Padé approximants is not possible without regularization. Table 6. Deconvolution of three-component function f (t) = 150e −0.004t + 19e −0.04t + 30e −0.06t using RPL when p 0 = 0.04. 2.0717494 × 10 −1 −2.0717494 × 10 −1 0.0000000 × 10 0 0.0000000 × 10 0 0.0000000 × 10 0 0.0000000 × 10 0 Table 7 shows the results of deconvolution of generated data points of the function f (t) in Equation (12) using RPL in [10/11] and [11/12] Padé approximants. Table 7. High levels of Padé approximants in deconvolution of three-component function f (t) = 150e −0.004t + 19e −0.04t + 30e −0.06t using RPL when p 0 = 0.04.

Deconvolution of Noisy Data by RPL
To show the capability of regularization in tackling the deconvolution process, consider a two-component function, f (t) in Equation (14), and analyze it with RPL after adding white Gaussian noise with the signal-to-noise ratio of SNR = 10 dB. Figure 1 shows the plots of function f (t) before and after adding the white Gaussian noise.  Table 12 shows the result of deconvolution of the noisy data shown in Figure 1 by RPL. The deconvolution results indicate that the RPL can find the exponential decays from the data points perturbed by a noise with SNR = 10. However, by increasing the Padé levels, the results start to deviate from the actual values. This deviation is attributable to the cumulative errors involved in the numerical integration of the noisy data.
As explained in Section 2, to calculate the Taylor expansion coefficients, numerical integration must be performed on the discrete data points. The noise results in cumulative error in numerical integration. For low levels of Padé approximants, where a small number of coefficients need to be calculated, the RPL is capable of deconvolution of the data points; however, by increasing the number of levels, the numerical error results in deviation from the actual values. Figure 2 shows the plots of function f (t) before and after adding white Gaussian noise of SNR = 10 dB.   Tables 13 and 14 give the deconvolution results for Equation (32) before and after adding the noise, respectively. Like the results we found in the case of Equation (14), the RPL is capable of finding the exponential modes after disturbing the data points. The deviation observed in the results shown in Table 14 is attributable to perturbation of data with the noise. Table 13. Deconvolution of three-component function f (t) = 40e −0.002t + 35e −0.009t − 60e −0.05t , using RPL when p 0 = 0.05.

Conclusions
The Padé-Laplace (PL) method is a powerful numerical scheme for deconvolution of Maxwellian modes from stress relaxation data of polymers obtained in discrete time intervals. The PL method needs only one parameter to perform the computations and does not require any initial guesses for the number of modes and parameters (amplitude and exponents) associated with each mode. The amplitudes and their corresponding exponents convey important information relating to the rheological behavior of polymers. In effect, the relaxation spectrum is the fingerprint of any polymer that is necessary to formulate a constitutive equation and can be used to predict its rheological behavior. A crucial step in this numerical procedure is constructing the Padé approximants that is an ill-posed problem. Since 1987, when the PL method was developed for separation of exponentials, the potential problem of ill-posedness attributable to the computation of the Padé table has been overlooked. In this paper, it was shown that the computation of the Padé approximants can result in ill-conditioned systems of equations. Therefore, it was shown that, apart from its elegant mathematical structure, the PL method that was believed to be able to solve the ill-posed problem of separation of exponentials using the properties of Laplace transform of an analytic function [9] can produce ill-conditioned systems of equations. As numerical computations demonstrate, the condition number of a matrix whose entries are the coefficient of Taylor expansion grows rapidly. A regularization of this method is possible by reconstructing the system of equations and using singular value decomposition (SVD) for computation of the Padé table. After regularization, the numerical computation indicates that the PL method can deconvolute data points even when p 0 , the only input parameter of the method, is chosen out of its optimal range. However, this occurs at the expense of calculating more levels of Padé approximants to achieve the stable modes. The analysis shown in this paper recommends applying the same regularization method in cases where the extended version of the PL method [9] was used to deconvolute experimental data. Although the focus of this paper in terms of application was on the deconvolution of viscoelastic spectrum of polymers, the results of this work will be fruitful in other areas such as analysis of chemical relaxation signals [13], voltage decays [14], fluorescence intensity decay [17,18], NMR relaxation data [19][20][21], and transient electric birefringence decay [22]. c n + b 1 c n−1 + b 2 c n−2 + · · · + b n c 0 = 0 c n+1 + b 1 c n + b 2 c n−1 + · · · + b n c 1 = 0 c n+2 + b 1 c n+1 + b 2 c n + · · · + b n c 2 = 0 . . . c 2n−1 + b 1 c 2n−2 + b 2 c 2n−3 + · · · + b n c n−1 = 0 (A5)