Fitting Curves by Simulated Annealing

T.J.P.Penna

Instituto de Física, Universidade Federal Fluminense,

Av. Litorânea, s/n 24210-340 Niterói, RJ, Brazil e-mail: gfitjpp@vmhpo.uff.br

October 20, 1996

Abstract

Simulated annealing is shown to be a powerful method to fitting curves. This method is particularly appropriate to fit curves with discontinuous derivatives where the traditional methods fail. Recent generalizations of the simulated annealing method enable easy and efficient implementation on microcomputers.

to appear in Computers in Physics

Fitting observational data (by experiments or numerical simulations) to a given function is a common task in science, sometimes to test a theoretical model or as interpolating purposes, etc. The most popular method consists in minimizing the squared deviation function

ε^{2} = ∑_{i=1}^{N} [y_{i} − y(x_{i}; a_{1} … a_{M})]^{2}  (1)

where y_{i} (i = 1 … N) are the data, y(x) ≡ y(x, a_{1} … a_{M}) is the function to be fit, with the M free parameters (a_{1} … a_{M}). Let us consider the data without error bars, for the sake of simplicity. The errors σ_{i} (i = 1 … N) can be introduced dividing each term on the right side of ( 1 ) by σ_{i}. The specific methods to fitting curves frequently demand a good estimate of the parameters to be fit[1]. In addition, derivatives of the function have to been evaluated. This fact restricts the application of those method to fit well-behaved functions. In this note, we show that a simulated annealing approach[2] can be useful in situations where the traditional methods fails as, for example, for discontinuous derivative functions. Moreover, the implementation is easier and, due to recent developments, the performance is also good.

The simulated annealing is widely used in the minimization of nonconvex (multiple minima) functions – as that presented in eq.(1) – as alternative to gradient descent methods. As far as we know, it has not been applied in fitting curves because other methods are more efficient than the original version of the simulated annealing. Here we present some improvements in this approach that enable it to work as an efficient and simple algorithm. The method is implemented as follows: the variables are allowed to change (a jump in the configuration space), according to a visiting distribution that gives the probability of occurrence of a jump of a given size (we will return to this point in more detail later). If this jump decreases the value of a cost function then the changes are always accepted, just like gradient descent. However, if the jump is in direction of increasing the cost function, the changes will be accepted according to an acceptance probability, depending on an artificial temperature. This step avoids the trapping in local minima. Then, the temperature is gradually cooled, according to a recipe that depends on the choice of the visiting distribution, and the process is repeat until the desired precision is reached.

Geman and Geman[3] showed that, in order to reach a global minimum, if the visiting distribution is a Gaussian one, the temperature must decrease logarithmically with time. This dependence strongly discourages the use of simulated annealing as modeling data algorithm. Szu and Hartley[4] proposed to use a Cauchy-Lorentz distribution (semi-local distribution) instead of the Gaussian one (local distribution). In this case, the temperature decreases like the inverse of time, therefore the algorithm is more efficient. Recently, Tsallis and Stariolo[5] proposed a new visiting distribution that recovers both cited before as special cases. A new parameter q_{v} ∈ [1, 3[ is introduced, such q_{v} = 1 recovers the classical simulated annealing (with Gaussian distribution of jumps, sometimes called Boltzmann machine), q_{v} = 2 recovers the fast simulated annealing (Cauchy machine). The probability of trial value for a parameter x changing from x_{t} to x_{t+1} is[5]

g_{q_{V}} (Δx_{t}) ∝ \frac{T (t)^{−D/(3−q_{v})}}{\{1 + (q_{v} − 1)(Δx_{t})^{2}/T (t)^{2/(3−q_{v})}\}^{[1/(q_{v}−1)]+[(D−1)/2]}} (2)

where T (t) is the artificial temperature at time t, D is the dimension of the function to be minimized and Δx_{t} = x_{t+1} − x_{t} is the size of the jump. For arbitrary D, the temperature should be cooled as

T (t) = T (1)\frac{2^{q_{v}−1} − 1}{(1 + t)^{q_{v}−1} − 1}.  (3)

Therefore the temperature decreases as a (q_{v}−1)-power of time, consequently in a much more efficient way than the earlier strategies.

Now, let us focus the discussion on the acceptance probability as function of temperature. The most frequent used is the so-called Metropolis algorithm[6]. Since the temperature is an artificial parameter the only condition to be respected is whether the acceptance probability leads to correct equilibrium states. Within the framework of a generalized statistical mechanics introduced by Tsallis[7], we can define a generalized version of the Metropolis algorithm[8]: the new value of the function to be minimized is accepted with a probability given by

p = min (1, (1 − (1 − q_{a}) (y(x_{t+1}) − y(x_{t})) /T (t))^{(1/(1−q_{a})})  (4)

where q_{a} is a new parameter. We recover the Metropolis algorithm in the q_{a} → 1 limit. The temperature for the acceptance probability can be different from the temperature for the visiting distribution. Here, we made them equal, for simplicity. Eq.(4) also leads to equilibrium states for all values of q_{a}[5, 8]. The most important feature of the eq.(4) – for our purposes – is that approximate solutions are found even at higher temperatures (for negative values of q_{a}) than had we used the Metropolis algorithm. Therefore, starting from a high temperature we do not need good initial values.

In fig.1, we present a complete code written in C where we try to fit a Gaussian function centered at origin. The variance and the amplitude are the parameters to be fitted. In order to show this method is useful to fitting curves, we present a situation where there are local minima. In the fourth line of the code, we define the function to be fitted (it will used also to generating the input data): two combinated Gaussian functions, with dispersion A for negative values of X, and dispersion C for the positive ones. The amplitude is given by parameter B. For input data we have chosen A = 2.0, B = 1.0 and C = 2.15. In this test case, we will try to fit a single Gaussian function (frequently in modeling data we do not know, a priori, the actual distribution of data, so this example should not be considered as a biased one). The error landscape is strongly dependent of the distribution of the input data, i.e., whether we have more data either in the positive or in the negative section of axis X. In this test, the input data were equally spaced in the X-axis in the interval [−1.7, 2.3]. In fig.2, one can see how the variance evolves in time when we started from a low temperature (30 times smaller than initial value of the error). The system is initially trapped in a local minimum (A = 2.0, in this case) but, after about 100 iterations, it goes to the global minimum (A = 2.1137 ± 0.0005). In this we can also see how accurate is the value of A, checking the fluctuations. In this example the existence of multiple minima turns evident but the temporary trapping can be avoided starting from higher temperatures. In this test, we follow a different approach from the one suggested by Stariolo and Tsallis[5]: we have allowed to change only one parameter by step. In our computer tests, this approach has been seen as more efficient than the original one in finding the global minimum. It is straightforward to extend this version for more than two parameters, just changing the value of D. Finally, we show the best fit for the test case, in fig.3.

In summary, we have applied successfully the generalized simulated annealing to fitting curves. The fast cooling scheme and a generalization of the Metropolis algorithm make this method competitive with the traditional others, with additional advantages as it avoids derivatives evaluation and good initial values. This method was used[9] to fit more complex function as a Lévy-distribution

P (I, ψ, γ) = \frac{1}{π} ∫_{0}^{∞} exp(−γq^{ψ}) cos(qI)dq  (5)

to experimental data from leaky faucets, with γ and ψ as free parameters and the integral was evaluated numerically[1]. Acknowledgments

Financial support from the Brazilian Agencies CNPq and FINEP is gratefully acknowledged.

References

[1] W.H. Press, S.A. Teukolsky, W.T. Vetterling and B.P. Flannery, Numerical Recipes in C, Cambridge University Press, NY (1992).

[2] S. Kirkpatrick, C.D. Gellat and M.P. Vecchi, Science 220, 671 (1983). S. Kirkpatrick, J.Stat.Phys 34, 975 (1984).

[3] S.Geman and D. Geman, IEEE Trans. Patt. Anan. Mach. Int., PAMI-6 6, 721 (1984).

[4] H.Szu and R. Hartley, Phys.Lett. A122, 157 (1987).

[5] D.A. Stariolo and C.Tsallis, in Annual Review of Computational Physics II, ed. D.Stauffer, World Scientific Po., Singapore (1994).

[6] N. Metropolis, A.W. Rosenbluth, M.N. Rosenbluth, A.H. Teller and E. Teller, J.Chem.Phys. 21, 1087 (1953).

[7] C. Tsallis, J.Stat.Phys. 52, 479 (1988).

[8] T.J.P. Penna, Phys.Rev. E (1994).

[9] T.J.P.Penna, P.M.C. de Oliveira, J.C. Sartorelli, W.M.Gonçalves and R.D. Pinto, “Long Range Anticorrelations: Is the Healthy Heart a Leaky Faucet ?” , preprint (1994).

Figure Captions

Fig.1 A code written in C, for fitting a gaussian function. Fig.2 Time evolution of the variance. Fig.3 Input data (squares) fitted by a Gaussian function (solid line).