Curve_fit for a function that returns a numpy array

Question

I know the library curve_fit of scipy and its power to fitting curves. I have read many examples here and in the documentation, but I cannot solve my problem. For example, I have 10 files (chemical structers but it does not matter) and ten experimental energy values. I have a function inside a class that calculates for each structure the theoretical energy for some parameters and it returns a numpy array with the theoretical energy values.

I want to find the best parameters to have the theoretical values nearest to the experimental ones. I will furnish here the minimum exemple of my code

This is the class function that reads the experimental energy files, extracts the correct substring and returns the values as a numpy array. The self.path is just the directory and self.nPoints = 10. It is not so important, but I furnish for the sake of completeness

def experimentalValues(self):
        os.chdir(self.path)
        energy = np.zeros(self.nPoints)
        for i in range(1, self.nPoints):
            f = open("p_" + str(i + 1) + ".xyz", "r")
            energy[i] = float(f.readlines()[1].split()[1])
            f.close()
        os.chdir('..')
        return energy

I calculate the theoretical value with this class function that takes two numpy arrays as arguments, lets say

sigma = np.full(nSubstrate, 2.)
epsilon = np.full(nSubstrate, 0.15)

where nSubstrate = 9

Here there is the class function. It reads files and does two nested loops to calculate for each file the theoretical value and return it to a numpy array.

def theoreticalEnergy(self, epsilon, sigma):
        os.chdir(self.path)
        cE = np.zeros(self.nPoints)
        for n in range(0, self.nPoints):
            filenameXYZ = "p_" + str(n + 1) + "_extended.xyz"

            allCoordinates = np.loadtxt(filenameXYZ, skiprows = 0, usecols = (1, 2, 3))
            substrate = allCoordinates[0:self.nSubstrate]
            surface = allCoordinates[self.nSubstrate:]
            for i in range(0, substrate.shape[0]):
                positionAtomI = np.array(substrate[i][:])
                for j in range(0, surface.shape[0]):
                    positionAtomJ = np.array(surface[j][:])
                    distanceIJ = self.distance(positionAtomI, positionAtomJ)
                    cE[n] += self.LennardJones(distanceIJ, epsilon[i], sigma[i])
                
        os.chdir('..')
        return cE

Again, for the sake of completeness the Lennard Jones class function is defined as

def LennardJones(self, distance, epsilon, sigma):
        repulsive = (sigma/distance) ** 12.
        attractive = (sigma/distance) ** 6.
        potential = 4. * epsilon* (repulsive - attractive)
        return potential

where in this case all the arguments are scalar as the return value. To conclude the problem presentation I have 3 ingredients:

• a numpy array with the experimental data
• two numpy arrays with a guess for the parameters sigma and epsilon
• a function that takes the last parameters and returns a numpy vector with the values to be fitted.

How can I solve this problem like the approach described in the documentation https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html?

Answer 1

Curve fitting

The curve_fit fits a function f(w, x[i]) to points y[i] by finding w that minimizes sum((f(w, x[i] - y[i])**2 for i in range(n)). As you will read in the first line after the function definition

[It uses] non-linear least squares to fit a function, f, to data.

It refers to least_squares where it states

Given the residuals f(x) (an m-D real function of n real variables) and the loss function rho(s) (a scalar function), least_squares finds a local minimum of the cost function F(x):

Curve fitting is a kind of convex-cost multi-objective optimization. Since the each individual cost is convex, you can add all of them and that will still be a convex function. Notice that the decision variables (the parameters to be optimized) are the same in every point.

Your problem

In my understanding for each energy level you have a different set of parameters, if you write it as a curve fitting problem, the objective function could be expressed as sum((f(w[i], x[i]) - y[i])**2 ...), where y[i]is determined by the energy level. Since each of the terms in the sum is independent on the other terms, this is equivalent to finding each group of parametersw[i]separately minimizing(f(w[i], x[i]) - y[i])**2`.

Convexity

Convexity is a very convenient property for optimization because it ensures that you will have only one minimum in the parameter space. I am not doing a detailed analysis but have reasonable doubts about the convexity of your energy function.

1. The Lennard Jones function has the difference of a repulsive and an attractive force both with negative even exponent on the distance this alone is very unlikely to be convex.

2. The sum of multiple local functions centered at different positions has no defined convexity.

3. Molecular energy, or crystal energy, or protein folding are well known to be non-convex.

A few days ago (on a bike ride) I was thinking about this, how the molecules will be configured in a global minimum energy, and I was wondering if it finds that configuration so rapidly because of quantum tunneling effects.

Non-convex optimization

The non-convex (global) optimization is different from (non-linear) least-squares, in the sense that when a local minimum is found the process don't return immediately, it start making new attempts in different regions of the search spaces. If the function is smooth you can still take advantage of a gradient based local optimization method, but the complexity is still NP.

A classic global optimization method is the Simulated annenaling, if you have a chemical background I think you will have some insights reading about it. Once upon a time, simulated annealing was provided in scipy.optimize.

You will find a few global optimization methods in scipy.optimize. I would encourage you to try Basin hopping, since it was successfully applied to similar problems, as you can read in the references.

I hope this drop you on the right way to your solution. But, be aware that you will probably need to spend, learning how to use the function and will need to make some decisions. You will need to find a balance of accuracy, simplicity, efficiency.

If you want better solution take the time to derive the gradient of the cost function (you can return two values f, and df, where df is the gradient of f with respect to the decision variables).