Surface fitting with constraints

I have a 3D-surface Fitting Problem. I want to find the best fit function for my lab-data with some constraints I already tried the Curve Fitting Toolbox and the Optimization Toolbox but I am not sure if they can help me with my problem.

I got 3 vectors of lab-data

xdata, ydata, zdata

I want to surface-fit the model with some polynomial equation z = f(x,y) = a+b*x + c*y + d*x^2... so that the R^2 value with respect to the lab-datapoints is as close to 1 as possible.

This works wonderful with the Curve Fitting Toolbox. But now I need to add some constraints:

z(1) < z(2) for y(1) > y(2) (if y or is increasing, z has to decrease monotonically)

z = 1 for x = 0 OR y = 0 (the 'edges' of the surface plot have to be 1 in the z-dimension)

Maybe it is more a mathematical/modelling problem than a programming problem.

Is there any way to solve this in matlab?

Answer 1

[Torsten]

The first constraint gives dz/dx < 0, thus b+2*d*x(i) < 0 for all i.

The second constraint gives a=1 in your model equation.

You can use MATLAB's "lsqlin" to solve this kind of problem.

Best wishes

Torsten.

[Andreas Harter]

Thank you. I'm struggling to implemente this into the lsqlin formula. Can you show me the syntax for the two constraints that you mentioned?

I would have to use this form: x = lsqlin(C,d,A,b,Aeq,beq,lb,ub). C are my xdata and ydata labpoints, and d is my zdata-vector. But how do I define (A,b) (nonequality constraints) and (Aeq, beq) (linear equality constraints)?

Edit: The second constraint is: z = 1 if x OR y = 0. I didn't write it correctly above. So a=1 does not fulfill this constraint.

[Torsten]

n=length(x);
C=[x y x.^2];
d=z-1;
A=[ones(n,1) zeros(n,1) 2*x];
b=zeros(n,1);
p=lsqlin(C,d,A,b);

x,y and z are column vectors.

Best wishes

Torsten.

[John D'Errico]

Note that the constraint Torsten describes is only a necessary constraint for monotonicity. Thus if one enforces that constraint at a set of points x(i), then the function will have a negative derivative at those points, but it MAY easily have a positive derivative at other intermediate points.

So the constraint described is NECESSARY but not SUFFICIENT to enforce monotonicity.

As well, be careful.

Consider the model with unit constant term

   z(x,y) = 1 + a*x + b*y + c*x*y + d*x*y^2 + ...

Note that z(0,0) is 1. HOWEVER, z(x,0) is NOT always 1 for all x.

In order for the model to be zero for x==0 OR y==0, so along each axis, then you may have terms only of the form

c*x^n*y^m

where n and m are all at least 1. Any term with n=0 OR m=0 will fail the edge constraint. I explain all of this in my answer too.

Answer 2

[John D'Errico]

Ugh. This could become a long mess of a discussion. I'm not even sure why I am entering into it. :)

This is a BAD problem to solve using polynomials. Just flat out BAAAAD. Polynomial models are nasty things, that have no good properties like monotonicity. At least none that are easily enforced. And high order polynomials are things to be avoided like the plague.

Ok, having said that, SOME aspects of the problem are easily solvable using polynomials, as long as you use mathematics.

Given your "constraint" is that for x==0 or y==0, then z MUST be 1? The solution is simple.

- The polynomial model must have a constant term of 1.

- There can be no terms in the model of the form c*x^n or c*y^n.

- All non-constant terms in the model must have at least one factor of x and one factor of y.

The second and third requirements are the same of course, just stated differently. So a valid model might be

z(x,y) = 1 + c1*x*y + c2*x^2*y + c3*x*y^2 + ...

As you can see, when either x or y is zero, then those other terms all drop out.

Conversely, suppose the model had a term of the form cn*x^n? Thus a pure power of x but without y in it? Suppose that x was nonzero, but y was zero? Then it must be true that z(x,0) will not be identically 1. (Think about it.) It may be equal to 1 for SOME values of x, but not for ALL values of x.

Therefore you could achieve your boundary condition goal by a simple choice of polynomial model. Estimating that model is as easy. Just use a model with no constant term at all, then fit z-1 as a function of x and y, but with NO constant term, and ONLY multinomial terms in the form I describe.

Ok, the next issue is how to solve the monotonicity problem. Sadly, for polynomials, you simply cannot do so. At least not easily.

The constraint that

diff(z(x,y),x) >= 0

is true, even if only over the support of your data, is a very difficult one to enforce for polynomials, even in one dimension. It is a problem that is highly nonlinear in the coefficients. Worse, the constraint is not itself one that tools like fmincon would use, since this is a constraint that must be valid everywhere.

Note that there is a solution for low order polynomials (up to cubic) in ONE dimension for this problem, but even for a cubic polynomial, the constraint is a nonlinear one. So, for example, I use that constraint form for spline models that must remain monotone. But even that gets nasty for two dimensional problems.

Ok, so given all that, is there a solution? I'd argue that a simple one is that posed by a tool like gridfit, although gridfit itself will NOT be amenable to this, since gridfit does not accept constraints. (You can find gridfit on the file exchange, as one of the tools I supply there.)

The idea is that gridfit, since it builds what is essentially a piecewise linear spline in two dimensions, WOULD be easily enough hacked to force the values of the function along the edges to be 1. You simply set those values to 1, and do not estimate the values at those nodes in the grid.

Next, monotonicity in x is nothing more than enforcing the constraint that

z(x+1,y) <= z(x,y)

for a function that is decreasing as a function of x. This can be done using a solver like lsqlin for the problem. Or you could transform the problem to be one that lsqnonneg can handle, IF you lack the optimization toolbox.

The end result of such a model would NOT be one where any functional form can be written down. It would simply be a smooth set of function values over a lattice, that emulates the data as well as possible. Could one write the code for all this easily? Sigh. I could, in perhaps a few hours. But for another to write it, that would require an understanding of how to build and solve large scale sparse linear algebra problems. It would require an understanding of the ideas behind regularization to enforce smoothness as a goal. It would require an understanding of how to use tools like lsqlin or lsqnonneg to enforce the required monotonicity constraints.

So the SIMPLE answer is to use a polynomial model with only low order terms, and ONLY the proper terms.

The complex answer is to use a tool like gridfit, but that might take some work to code for the less than expert.

An intermediate answer might be to use a nonlinear model, properly constructed to implicitly enforce the desired constraints. Sadly, that might be difficult to build without knowing the shape of the surface, because one would need to understand how to construct such a nonlinear model. One would need to understand the source of the data to choose a viable model.

Edit: I just noticed that the data was provided as an attachment.

x = MeshD.labx;
y = MeshD.laby;
z = MeshD.labz;

plot3(reshape(x,[],4),reshape(y,[],4),reshape(z,[],4),'o')
box on
grid on
view(27,35)

As you can see, the data is not truly monotonic itself in x.

reshape(z,[],4)
ans =
            1            1            1            1
      0.97741      0.98517      0.98386      0.99293
      0.96506      0.97776      0.97852      0.98762
      0.95648      0.97121        0.974      0.98295
      0.94889      0.96538      0.97039      0.97978
      0.95025      0.96783      0.96257      0.97726
      0.93876      0.95972      0.96145      0.97542
      0.93762      0.95852      0.95999      0.97366
      0.92935      0.95127      0.95278       0.9685
>> diff(reshape(z,[],4),[],2)
ans =
            0            0            0
    0.0077573   -0.0013096    0.0090745
     0.012699   0.00075837    0.0091021
     0.014732    0.0027827    0.0089524
      0.01649    0.0050132    0.0093895
     0.017583   -0.0052639     0.014695
     0.020965    0.0017266     0.013968
     0.020896    0.0014703     0.013671
     0.021919    0.0015129     0.015724

Were the data truly monotonic for z as a function of x, then there would be no negative results from that call to diff. So the data is somewhat noisy, as is all experimental data.

Worse, the data is not even terribly well behaved for your goals. LOOK AT YOUR DATA!

So here I've plotted the data, as well as a simple polynomial model on top as a surface. I estimated that model as

z(x,y) = 1 - 0.00018698*x*y - 2.5261e-07*x*y^2 + 3.7255e-10*x*y^3 + 3.7816e-05*x^2*y + 1.1786e-07*x^2*y^2

There is a serious problem here. See that this relatively low order model does nasty things in the corner at (x,y)=(3,500). It is not in fact monotonic at all. It is especially poor since we are using a polynomial to extrapolate. That is ALWAYS a BAD thing.

Almost as bad is what happens near x==0. See that IF you will force z(0,y)==1, then see what happens at z(0.5,y). Your function dips down rather far, then it goes up again to z(1,y). Then it drops off again for larger values of x.

So your data is rather strongly non-monotonic. In fact, it appears that I would not trust the data at all, if you require this relationship to be monotone decreasing for z as a function of x.

[Andreas Harter]

Well that was one hell of an answer, thank you very much for your time and effort. The problem seems to be more complex that I hoped for...

And damn, again I did a mistake in discribing my problem: The fit has to be monotonic towards the y-dimension NOT towards the x-dimension. z(1) < z(2) for y(1) > y(2) (if y or is increasing, z has to decrease monotonically). Sorry for the mess. With this constraint the data is still a bit noisy but ok. Of course, a monotonic behaviour of the x-dimension is not at all close to the lab data. l that was one hell of an answer! Thank you for your time and effort. I somehow hoped that there would be an easy answer for this problem but I slowly realize that it is way more complex that I thought it would be.

Edit: I tried the Curve Fitting Toolbox again and entered this term: 1+a*x*y+b*x^2*y+c*x*y^2+d*x^3*y+e*x*y^3+f*x^2*y^2

The edges are now okay and z is monotonically decreasing with regard to the y-dimension (although I did not declare any monotony constraints, this could just be luck). However R^2 is 0.7171 which is way to low. :/