Curve fitting with constraints

[jean-patrick francoia]

Hi,

I'm trying to use cvxpy to fit experimental data (they 're obeying a parametrized model). I tried lmfit and other stuff included in scipy, but it didn't work because it is impossible to define complex constraints with these modules.

Here is my model, defined in the scipy way:

 def model(t, y0, h, A1, x1, s1, A2, x2, s2, A3, x3, s3, A4, x4, s4):

            y = y0 + h * (
                          A1 * (0.5 + 0.5 * erf((t - x1) / (s1 * 2**0.5))) +
                          A2 * (0.5 + 0.5 * erf((t - x2) / (s2 * 2**0.5))) +
                          A3 * (0.5 + 0.5 * erf((t - x3) / (s3 * 2**0.5))) +
                          A4 * (0.5 + 0.5 * erf((t - x4) / (s4 * 2**0.5)))
                         )

            return y

Basically, "t" and "y" are array of values obtained experimentally. And y = f(t), very basic.
I'm trying to fit the experimental data, and get values for y0, h, A1, A2, A3, A4, x1, x2, x3, x4, s1, s2, s3 and s4. I would like to try a basic approach, a least square fitting.

For now I'm just getting started. I tried this:

     

        t = cvx.Variable(self.x)
        y = cvx.Variable(self.y)
        y0 = cvx.Variable()
        h = cvx.Variable()
        A1 = cvx.Variable()
        A2 = cvx.Variable()
        A3 = cvx.Variable()
        A4 = cvx.Variable()
        x1 = cvx.Variable()
        x2 = cvx.Variable()
        x3 = cvx.Variable()
        x4 = cvx.Variable()
        s1 = cvx.Variable()
        s2 = cvx.Variable()
        s3 = cvx.Variable()
        s4 = cvx.Variable()

        obj = cvx.Minimize(y0 + h * (
                          A1 * (0.5 + 0.5 * erf((t - x1) / (s1 * 2**0.5))) +
                          A2 * (0.5 + 0.5 * erf((t - x2) / (s2 * 2**0.5))) +
                          A3 * (0.5 + 0.5 * erf((t - x3) / (s3 * 2**0.5))) +
                          A4 * (0.5 + 0.5 * erf((t - x4) / (s4 * 2**0.5)))
                         ))
        prob = cvx.Problem(obj)
        prob.solve()

I'm well aware I didn't define any constraint for now, I'm just starting. I'm trying to undesrtand how the module works for now.
The previous code gives:

ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

Could you give me a hand please, and help me to start ?

[Steven Diamond]

CVXPY only solves convex optimization problems. Moreover, you have to format your problem so it satisfies the DCP rules (see dcp.stanford.edu and the tutorial at cvxpy.org).

The main upshot of these limitations is you can only use CVXPY functions to formulate your problem. You can't combine the erf function from scipy with CVXPY variables.

Also, the problem you're trying to solve is not convex since erf is not convex. There are convex optimization based methods for solving problems of this form, but they require a solid background in convex optimization to use. If you think you're up for it check out these slides on sequential convex programming: http://web.stanford.edu/class/ee364b/lectures/seq_slides.pdf

The basic idea is to solve a series of convex approximations to your problem. In this case you would take the first order Taylor expansion of your function and use that in the least-squares objective. 

You might be interest in generalized least-squares solvers like CERES. These solvers won't give you the same expressivity as CVXPY in formulating constraints, but they'll automate more of the work involved in using sequential convex programming.