Comparing two Fitting Algorithms for Determining the Cole–Cole Parameters in Blood Glucose Problems

: This paper addresses the nonlinear inverse problem of estimating the parameters of the Cole–Cole model used to describe the behavior of the complex permittivity of blood samples. Such a model provides an efﬁcient and accurate representation of biological tissues in the entire frequency band considered and reduces the complexity of the experimental data to a few parameters. In this way, it is possible to extract a “synthetic view” of the dielectric properties of tissues in such a way that more information on the glucose concentration can be derived, in addition to the resonance peak or phase shift. In order to perform the ﬁtting of the Cole–Cole model, two different algorithms were used and compared: the Levenberg–Marquardt and the variable projection algorithms. Synthetic data present in the literature were used to evaluate the performances obtainable with these methods. In particular, Monte Carlo analysis was used in order to evaluate the accuracy and the precision that these two methods provide in the process of estimating the parameters involved, with respect to the starting points of the parameters. The results obtained showed that the variable projection algorithm always outperformed the Levenberg–Marquardt one, although the former has a greater computational burden than the latter.


Introduction
Diabetes is a metabolic disorder that afflicts millions of people in the world [1]. It degrades the cell's ability to absorb glucose from the bloodstream because of the improper regulation of the insulin hormone. For this reason, great efforts have been dedicated to the development of non-invasive glucose monitoring devices, which may considerably improve the quality of life for diabetics [2].
This work is particularly concerned with microwave sensor technology that relies on the change in the dielectric and conductivity properties of blood plasma as a function of the glucose concentration in order to track such a change.
In this framework, developing accurate and precise fitting methods for blood models, at different glucose concentrations, is essential for the development of robust electromagnetic (EM)-based techniques that could be employed for non-invasive, continuous glucose monitoring. Indeed, accurate electromagnetic tissue modeling is of paramount importance since it affects the simulation stage required for sensor design [3]. Moreover, extracting a "synthetic view" (in terms of a few parameters) of the sensor response data is essential for analyzing patterns and possibly extracting more information, besides the resonance peak or phase shift, about glucose concentration.
In this paper, the aim was just to address the fitting problem. More in detail, starting from the dielectric spectrum, which is assumed known over a certain number of frequencies, we aimed at estimating the parameters of a single-pole Cole-Cole model. As is well known, this entails solving a nonlinear inverse problem, which here was addressed by two different methods: the classical Levenberg-Marquardt method [4,5] and the variable projection algorithm [6]. We evaluated how sensitive the two methods are with respect to the starting points of the parameters and with what accuracy and precision these parameters can be estimated. In order to check the two methods, we first generated synthetic relative permittivities by employing a single-pole Cole-Cole model, using data from the literature [7] as true values for its parameters, and then solved an inverse problem in order to trace these values by resorting to the two aforesaid methods.
Although higher-order models could perform better, we considered a first-order model to perform the comparison in the simplest possible case.

Cole-Cole Model
The Cole-Cole model [8] is widely used to describe the complex relative permittivity of biological tissues, ε r (ω) = ε(ω)/ε 0 , and its equation is: (1) in which N is the number of poles and thus the order of the model, ε ∞ = lim ω→∞ ε r (ω) is the permittivity at high frequencies, σ is the static ionic conductivity and ε sn = lim ω→0 ε r (ω), and τ n , and α n are the static permittivity, the relaxation time constant, and the so-called distribution parameter of the n-th addend of the summation, respectively. Such a model incorporates the Debye model [9]. Indeed, the main difference between the Debye and the Cole-Cole models is that the latter includes the exponent 1 − α, with 0 ≤ α ≤ 1. When the exponent becomes smaller, the relaxation time distribution becomes broader, i.e., the transition between low-and high-frequency values becomes wider and the peak on the imaginary part of the spectrum also becomes wider. The complexity of both the structure and composition of biological material is such that the dispersion region of each pole may be broadened by multiple contributions to it. The broadening of the dispersion could be empirically accounted for by using the Cole-Cole model [10], which is expected to give more accurate dielectric spectrum curve fitting.

Curve Fitting Algorithms
Let be x the vector of model parameters and P its length and M the number of frequency points for which the measures are taken. We define the data vector ( stands for transposition) y = y(ω 1 ) . . . y(ω m ) . . . y(ω M ) in which the m-th component of the vector y is the observed value y(ω m ). Let also ε r = ε r (ω 1 ; x) . . . ε r (ω m ; x) . . . ε r (ω M ; x) be the model vector, here given by Equation (1)), with ε r (ω m ; x) being the estimation at ω m .
Solving the least-squares problem means findingx such that: in which the function to minimize, Ψ = 1 2 ε r (x) − y 2 2 , is the 2 quadratic norm of the misfit r = ε r (x) − y, which is a nonlinear function such that r : R P → R M with P M. We addressed the nonlinear fitting problem with two methods: the Levenberg-Marquardt algorithm (LMA) and the variable projection algorithm (VPA).
The Levenberg-Marquardt algorithm [4,5] acts more as a gradient-descent method when the parameters are far from their optimal value and acts more as the Gauss-Newton method when the parameters are close to their optimal value [11]. The equation for the step h at the kth iteration is: where J is the Jacobian of f and λ k is the damping parameter. It controls both the magnitude and direction of h, and it is chosen at each iteration. It can be shown [5] that, at each iteration, Equation (3) solves the minimization problem over a reduced set of admissible solutions, i.e., those that satisfy h ≤ R(λ), limiting the correction step to within a region near x k . The radius of the trust region R = R(λ) is a strictly decreasing function with lim λ→∞ R(λ) = 0. When λ k = 0, the step h is identical to that of the Gauss-Newton method, i.e., the same direction and maximum magnitude. As λ → ∞, h tends towards the steepest descent direction, with the magnitude tending towards 0. Based on the above, we inferred the qualitative update rule for λ k+1 : If Ψ(x k + h) < Ψ(x k ), then the quadratic approximation works well, and we can extend the trust region, i.e., it will be λ k+1 < λ k . Otherwise, the step is unsuccessful and we reduce the trust region, i.e., it will be λ k+1 > λ k ; in this way, the next step tends towards the negative gradient method, and a lower value of Ψ is more likely to be found.
The MATLAB implementation was used, in particular the lscurvefit function with the Levenberg-Marquardt option.
The variable projection algorithm [6] is a method used to solve separable nonlinear least-squares problems. The least-squares problem is said to be separable when the model parameters can be separated into two sets of parameters, one that enters linearly into the model, c = [c 1 , . . . , c k ], and another set of parameters that enter the model nonlinearly, a = [a 1 , . . . , a l ], so that x = [c, a]. For each observation y m of a separable nonlinear leastsquares problem, the model is a linear combination of nonlinear functions that depend on nonlinear parameters, and the model function can be written as ε r (ω) = ∑ k j=1 c j φ j (ω; a). The functional Ψ is written in terms of residual vector r as: in which the columns of the matrix Φ are the nonlinear functions φ j (ω; a). The linear parameters c could be obtained from the knowledge of a, by solving the linear least-squares problem: which stands for the minimum-norm solution of the linear least-squares problem for fixed a, where Φ(a) † is the Moore-Penrose generalized inverse of Φ(a). By replacing this in Equation (4), we obtain the variable projection functional: The variable projection algorithm consists of two steps: first minimizing Equation (6) with an iterative nonlinear method and then using the optimal value found for a to solve for c in Equation (5) [12]. The principal advantage is that the iterative nonlinear algorithm used to solve the first minimization problem works in a reduced space and less initial guesses are necessary. A robust implementation in MATLAB, called VARPRO [13], was adapted and used to deal with complex-valued problems, choosing the Levenberg-Marquardt option for the solution of Equation (6).

Numerical Simulations
The generation of the synthetic complex relative permittivity of blood plasma relied on quadratic fits to glucose-dependent Cole-Cole parameters reported in [7]; in particular, we considered two different concentrations, 100 mg/dL and 250 mg/dL, one normal and the other of severe diabetes, respectively, in accordance with the diagnostic criteria in [1]. The data vector consisted of M = 1000 points in the frequency range 500 MHz-20 GHz.
In gradient-like algorithms (such as those used in this paper), the choice of the initial point is a crucial factor for the convergence of the procedure. For the single-pole model case, it is fairly easy to exploit the physical meaning of the parameters to infer an initial estimate. However, since the noise can invalidate the initial estimate, we propose to study the robustness of the two algorithms with respect to the initial point. To this end, a Monte Carlo analysis was performed, iteratively evaluating the deterministic algorithms using a set of N = 1000 uniformly distributed random initial points arranged in a 5D hypercube of the parameter space in order to statistically characterize the results. Each side of the hypercube represents an interval containing the range of variation of each parameter for the glucose concentrations considered.
We ran simulations on a machine with Intel i9-10850K (10 physical cores), 32 GB RAM, and Ubuntu 21.04, and we took advantage of the Parallel Computing Toolbox, using the parfor loop for running the 1000 simulations.
The intervals for generating the random initial value for each parameter (of the Cole-Cole model) were chosen from the data tabulated in [7]. In particular, the widths of these intervals were the same for each glucose concentration and were: [1,5] for ε ∞ , [1,150] for ε s , [1 × 10 −14 , 1 × 10 −11 ] for τ, [0.1 − 1 × 10 −9 , 0.1 + 1 × 10 −9 ] for α, and [0, 5] for σ. These intervals are relatively large compared to the values taken from [7] in order to test the two algorithms in sufficiently stressful situations. Only the range of variation of α is extremely small because the model used in [7] practically fixes it a priori to 0.1. Obviously, it must be taken into account that the VPA requires only generating the values for τ and α.
For an initial bland qualitative assessment, we established evaluation intervals (the same for generations) for the estimated parameters so that we could assert that a reconstruction is "good" if it falls within these ranges, and "wrong" otherwise.
be the vector of the parameter estimates returned by the two algorithms at the i-th simulation, and letx (i) denote one of its five ele- be the sample mean and standard deviation, respectively, calculated for each parameter.
For a quantitative evaluation of the performance of the two algorithms, we then define multiple figures of merit for statistically characterizing the results of Monte Carlo analysis. For each parameter, Equations (7) and (8) define measures of accuracy and precision, respectively, defined over the entire set of reconstructions. However, such measures can be greatly affected by estimates that are very far from the true value, x true , where the latter represents one of the five elements of the vector of reference valueŝ x true = [ ε ∞ true , ε s true , τ true , α true , σ true ]. For this reason, we also introduced A cut and P cut in Equations (9) and (10) in order to define the accuracy and precision measures, respectively, which instead dampen the effect of the above isolated events. The subscript cut means that they were calculated on a subset obtained by eliminating ζ% of reconstructions with lower values and ζ% of reconstructions with higher values, in which 0 < ζ < 50.

Results
We conducted many numerical simulations by widening more and more the generation intervals. In this paper, we report the case where the generation intervals are very large except for the α interval, due to the above explanation. In all these experiments, following the qualitative criterion mentioned above, the VPA always provided good estimations, while the same was not true for the LMA. In particular, in the case considered, the LMA provided the wrong estimates in about 260 simulations out of 1000, for each of the two glucose concentrations, while no wrong estimate was returned by the VPA. Here, by the "wrong" estimate, we mean that at least one component of the parameter vector x had a value that was outside its generation range. Graphical representations of those qualitative results are provided in Figure 1. Consistent with the qualitative results, considering the whole set of 1000 estimates, the VPA exhibited excellent accuracy and precision, while this was not the case for LMA, as can be observed in Tables 1 and 2. For this reason, for the VPA only, the results calculated by means of Equations (7) and (8) are reported, whilst for the LMA, the results derived from Equations (9) and (10), obtained by cutting 25% of the lowest values and 25% of the highest values, were also considered.

Parameter estimates
Regarding execution times, the LMA took 7 s to finish 1000 simulations, while the VPA took around 170 s.

Discussion
In this paper, we faced the problem of fitting the dielectric spectrum of a blood sample in order to estimate the parameters of the single-pole Cole-Cole model. In particular, we compared the performance of two different algorithms, the LMA and VPA, in terms of the accuracy and precision with respect to the starting points of the parameters.
For the parameter range considered, the VPA outperformed the LMA in robustness with respect to the initial point of the algorithm. However, analyzing the figures of merit related to the LMA, it became clear that there were erroneous reconstructions so far from the true value such that they heavily deteriorated the (standard) accuracy and precision, while the cut versions did not suffer from this problem. In fact, once the erroneous ones were removed, in all other simulations, the algorithm converged to the true values.
On the other hand, the VPA gave these good results because less initial guesses were necessary and because the iterative nonlinear algorithm used to solve the first minimization problem works in a reduced space. The big disadvantage of the VPA, or at least of the implementation used in this paper, was the execution time. This was certainly due to the numerous singular-value decompositions (SVDs) that the algorithm calculates in its runtime.
The results are promising, and the research will continue by evaluating the algorithms in increasingly realistic scenarios, including adding noise to the synthetic data.

Data Availability Statement:
The data presented in this study are available upon request from the corresponding author.

Conflicts of Interest:
The authors declare no conflict of interest.