Surface fit with (partial) differential equation constraint

For a 3rd-order surface fit, let's say an sfit object z=f(x,y), I would like for the change of z with respect to x to never be negative (partial dz/dx >= 0). Is it possible to set up such a constraint when fitting or maybe afterwards, and if so how would one proceed?

I need this for a pressure ratio (x2) vs mixing ratio relation. When increasing the pressure upstream, I need to assume the mixing ratio can only increase (or at least stay constant).

Edit: I guess what I'm looking for is a tool for or means of monotonic surface fitting, based on scattered non-monotonic x-y-z data.

Answer 1

[John D'Errico]

There are several ways one could do this, and YES, it is possible. There are some issues one must resolve of course.

The simplest solution is to use a tool like my gridfit. That tool is simply a piecewise linear regression spline, defined in terms of the function values at each node in the grid. So were you to build a 100x100 grid, there are 10000 unknowns to estimate using a regression scheme. The constraint that dz/dx>=0 everywhere is simple to apply here of course, since one simply forms one inequality constraint equation for each pair of neighboring x-connected nodes. So there would be 99*100 such linear inequality constraints. The resulting system might be rather large for LSQLIN to solve I fear. However, you could probably solve a 20x20 problem easily enough, so only 400 unknowns, and 380 linear inequality constraints.

As I said, the above is the simple solution, although it might take some CPU time depending on how large of a grid you chose to populate.

Next, we can consider a 3rd order model in x and y. I'll assume:

z(x,y) = a00 + a10*x + a01*y + a20*x^2 + a11*x*y + a02*y^2 + a30*x^3 + a21*x^2*y + a12*x*y^2 + a03*y^3

There are 10 coefficients here. The problem, if we wish to constrain the slope, dz/dx to be positive everywhere inside the bounds, is the constraint becomes a nonlinear one. Worse, since it must be true everywhere in that domain, it gets a bit nasty. A simple solution is to define a set of NECESSARY constraints, for positivity. That is, at each of a large set of points in the (x,y) domain of interest, we will set a linear constraint such that

dz/dx >= 0

Once a point (x_k,y_k) is chosen, that derivative becomes a linear inequality constraint in the parameters.

If we choose a large enough set of points for that set of linear inequality constraints, then you can see that there MAY be some small region in the domain of interest where the constraint may SLIGHTLY fail, but only by a small amount if our test set is chosen to be large enough. The result is a linear least squares problem, with 10 unknowns and as many linear inequality constraints as you choose to generate. The more constraints you choose, then better the result.

Again, the above scenario is a NECESSARY, but NOT SUFFICIENT constraint on slope positivity. It should get you pretty close however.

Another approach can be found from the work of Fritsch and Carlson. (This is the idea that originally spawned PCHIP by the way.) Here we can form a SUFFICIENT constraint on the linear system for the coefficients of the cubic model. The problem is that the constraint system would be nonlinear but still convex, so the trick is to convert that nonlinear constraint into a set of LINEAR inequalities that ensure the solution has the desired slope.

The difference between the last two options (necessary versus sufficient) is a subtle one. he former may result in a solution that may have some small areas of slightly negative slope. The latter may possibly result in a slightly too conservative fit, but it always ensures slope positivity.

Of course, there are no codes written that I know of to directly solve any of these variations, but it is not terribly difficult. In fact, option2 is almost trivial to cobble together, and option 1 not too much harder. The downside of the polynomial model options is they are simply polynomial models, and not terribly high order ones at that. And since low order polynomials tend to be somewhat poor in their abilities to fit, there will be limits here too.

My personal suggestion is to throw together option 2. See how well it does, if there are regions of negative slope encountered. (Note that this will not generate a GLOBALLY positive slope, only a slope that is everywhere positive in the region of interest.) If the third order model was deemed insufficient after seeing the fit quality and lack of fit, then I'd suggest option 1 for a reasonably small grid, perhaps 20x20, maybe a bit larger. See how fine you could go before causing LSQLIN to have a hissy fit.

[Star Strider]

That function is linear in the parameters. Will it work with a nonlinear function?

[John D'Errico]

Well, if the function is nonlinear in the parameters, then no, it won't work. However all the examples I've seen in this thread show models that are linear in the parameters, so lsqlin would work. Of course a third order polynomial model is nonlinear in the independent variables, but linear in the parameters. This is not a problem for LSQLIN.

Otherwise, one would be forced to use a tool like FMINCON to involve the linear inequality constraints.

[Image Analyst]

Does he have cross terms (x*y, etc.) in the model? If not, then why not use John's fine "polyfitn"?

[Daan]

In my case I only need to constrain the slope in the domain of my x-y data, so no need to consider the extrapolated bit for now. I guess this is what is intended with trying Option 2 first? Of course, someone else might require broader constraints and/or more accuracy. I will try to put something together, but evidently I'm no expert on the subject :)

Thank you btw for elaborating, I will also have a look at your gridfit tool.

Image Analyst, I did use polyfitn, a great tool indeed! I can obtain for example good third-order fits with it. The result is of course a polynomial with cross terms. I tried to avoid those by fitting a 'monotonic' surface such as a (quarter of an) elliptic paraboloid, x^2/a^2 + y^2/b^2 = z/c, even higher order, using anonymous functions, but have failed so far to get decent fits.

[John D'Errico]

The problem is polyfitn will not handle the constraint, although it has no problem with cross terms. And gridfit is also unconstrained, so only the basic idea would work there too.

Option 2 is pretty basic, but remember, it has limited abilities in terms of fitting. I'll see what I can hack together today.

[Daan]

Thanks again! Learning a great deal here.

I recently found that I should be looking for a relation of the form z=f(a,b,c), i.e. more than two independent variables. This is not a problem for polyfitn, its specialty actually, but does this change things? Or is the problem just getting more extensive? Sorry if I'm asking something trivial here.

[John D'Errico]

Polyfitn has no problem with three independent variables. It can handle as many as you choose. But again it has no ability to constrain the slope.

[Daan]

Right, but the strategies for fitting with one or more slope constraints you set forth above still hold for higher dimensions?

Answer 2

[Paul Shepherd]

I'm not entirely sure I understand the question, but let's say the equation you are fitting is z=k1*x+k2*y. If you are doing the fitting with the fit() function, you can generate a fit options structure using the fitoptions() function.

Once you have the fit options structure, you can assign a minimum value for your coefficients of 0: fit_opts_struct.Lower =[0 0];

If you are trying to fit something more like z=k11*x^2 + k12*x + k13 + k21*y^3 ...

then you are somewhat limited by the fact that the Lower array does not allow terms to interact. To do that, you would need to consider building an algorithm that fits k12 recursively based on the current value of k11, and changing the limit on k12 each time you try a new value of k11.

Hope this helps, Paul

[Daan]

Such an algorithm would be nice. Hypothetically I could try to build some sort of iterative loop which checks for certain derivatives (for example dz/dx of my z=f(x,y) surface fit) to be positive, but I wouldn't know how to proceed with the subsequent control of the coefficients.