fitting a funciton with minimax error

I am trying to solve a constraint problem regarding minimax error. Basically, I want to fit data to a specific type of function where I minimize the maximum error as the polynomial fit oscillates. This might involve weighting the data differently, fixing some coeffs, etc. The excel file is attached.

I have this code:

clear; clc; close all;
mat =xlsread('C:\example2.xlsx','Sheet1','A2:C32');

density = mat(:,1);
eta = mat(:,2);
Z_MD = mat(:,3);
eta_c = 1/1.55;

I want to fit the x data (eta) vs. y data (Z_MD) to the following functional form:

So I need to solve for my Ak values. How I can I minimize the maximums of the relative error? Obviously, since it's a polynomial, the error will fluctuate. Can I use MATLAB to minimize the maximums?

Currently, when I use cftool to fit the data, when I plot the error of the fit with respect to Z_MD, the maximum error is not minimized, meaning as the polynomial fluctuates through the data points, the error is not bound a constant max error.

Edit: Note that it could be fitting eta and Z_MD. density and eta are the same thing, just multiplied by a constant basically.This is why I changed it to eta, so I want to minimize the maximum error of the polynomial fit to my equation that is fit to x = eta and y = Z_MD.

Answer 1

[John D'Errico]

mat =xlsread('example2.xlsx','Sheet1','A2:C32');
density = mat(:,1);
Z_MD = mat(:,3);
eta = mat(:,2);
eta_c = 1/1.55;

Your model is:

Z(eta) = 1 + 3*eta/(eta_c - eta) + sum(k*A_k*(eta/eta_c)^k)

where the sum operates over k = 1:4. The goal being a minimax fit. Since the first two terms in that model are not variable at all, we will just subtract them for the fit. But at first, plot EVERYTHING. Always plot your data. Look at it. Think about it.

plot(eta,Z_MD,'ro')

The result being a simple power type function. It might seem things will work at least reasonably well, as polynomials like that sort of stuff. So plot more. Think more. Look at what we see. First, remember that your function has a singularity in it. The singularity lives at eta==eta_c.

eta_c
eta_c =
      0.64516
      
plot(eta,Z_MD - (1 + 3*eta./(eta_c - eta)),'ro')

So this is the function you wish to fit with a polynomial.

The first problem is, it does not seem to pass throughzero at eta==0. It is less than zero down there.

Z(eta) - 1 - 3*eta/(eta_c - eta) = sum(k*A_k*(eta/eta_c)^k)

You should observe that the right hand side is identically zero when eta == 0. The left hand side (what I just plotted)is not so. So this will create some amount of lack of fit down there. But worse, what happens when eta approaches 0.64516? The left hand side has a singularity.

Does the right hand side? NO. There can be NO singularities in a polynomial. And the right hand side is purely a polynomial, and only a 4th degree one at that.

Do you really want to go through with this? Sigh.

First, ets see what happens when we just use a simple linear least squares. Backslash can handle that well enough.

A2 = [(eta./eta_c).^(1:4).*(1:4)]\(Z_MD - (1 + 3*eta./(eta_c - eta)))
A2 =
        5.302
      -20.991
       33.224
      -16.441

So A2 is the 2-norm solution. Given those coefficients, we can look at the fit.

Zpred = [(eta./eta_c).^(1:4).*(1:4)]*A2;
plot(eta,Z_MD - (1 + 3*eta./(eta_c - eta)),'ro')
hold on
grid on
plot(eta,Zpred,'b-')

Sigh. Double sigh. Well, it does look reasonably close to a minimax fit as it is. A 2-norm fit like this will often not be unreasonable. Lets look at the residuals.

plot(eta,Z_MD - (1 + 3*eta./(eta_c - eta)) - Zpred,'ro')
grid on
yline(0);

So indeed pretty close to a minimax fit, as I said. It only deviates by a bit. Is it really worth the effort, for a fit that does not fit that well in the first place?

I'd suggest that looking for a true minimax fit here is a bit extreme, given the significant problems in the fit. Do you really need it? Regardless, fminimax is much easier to use than you seem to be making it appear. Two lines suffice, although it would always be a good idea to use the 2-norm solution as a start point. So three lines.

mmfun = @(A) abs(Z_MD - 1 - 3*eta./(eta_c - eta) - [(eta./eta_c).^(1:4).*(1:4)]*A);
[Amm,~,maxres,exitflag] = fminimax(mmfun,A2)

Local minimum possible. Constraints satisfied.

fminimax stopped because the size of the current search direction is less than
twice the default value of the step size tolerance and constraints are 
satisfied to within the default value of the constraint tolerance.

<stopping criteria details>

Amm =
       5.7492
      -23.245
       36.892
      -18.375
maxres =
      0.15832
exitflag =
     4

So a little different from A2. But honestly, hardly worth the effort. Note that I was carful to use abs in mmfun, since fminimax minimizes the maxima, but not the absolute max. From the help:

fminimax finds a minimax solution of a function of several variables. 
    fminimax attempts to solve the following problem:
    min (max {FUN(X} )  where FUN and X can be vectors or matrices.

So you do need an abs in there.

Zpredmm = [(eta./eta_c).^(1:4).*(1:4)]*Amm;
plot(eta,Z_MD - (1 + 3*eta./(eta_c - eta)) - Zpredmm,'ro')
yline(0);
grid on

From a different perspective, can you honestly see the difference below, between the green and blue curves? Is that difference significant?

plot(eta,Z_MD,'ro',eta,(1 + 3*eta./(eta_c - eta)) + Zpredmm,'b-',eta,(1 + 3*eta./(eta_c - eta)) + Zpred,'g-')
grid on

Whatever floats your boat. Personally, I think focusing on a minimax approximation here is somewhat misguided, because there are essential flaws in that approximation that far exceed the tiny differences.

[Matt J]

Note that I was carful to use abs in mmfun, since fminimax minimizes the maxima, but not the absolute max.

It will when the AbsoluteMaxObjectiveCount input option is used, and this is probably necessary, to avoid differentiability issues in the abs() operator.

[Benjamin]

Wow, thanks for this extensive answer! I'm still digging through what you have, but have a couple of questions about this!

1) If I want to add a 5th term, do I just change the 1:4 to 1:5?

2) How can I set A1 = -0.4194, A2 = 0.5812, and A3 = 0.6439?

Basically by doing this, the correct behavior of the curve is captured as eta -> 0.

So I want these to be set to these values, and not let them vary.

Also note that at eta=0, then Z should =1. Perhaps these constraints I have mentioned above will help that be the case.

Then I can allow A4 to vary. This may not be a great fit, but then I can subsequently add more terms (A5, A6, etc) until I get the fit that I want.

Edit: Here is the code you have. I changed it so there are 5 terms. How can I easily set A1, A2, and A3 to values I want?

Edit: Also, I notice you plot the error in terms of Z_MD - Z_predmm - other terms. Is there a way to minimize the max such that the relative error is minimized?

clear; clc; close all;

mat =xlsread('C:\example2.xlsx','Sheet1','A2:C32');
density = mat(:,1);
Z_MD = mat(:,3);
eta = mat(:,2);
eta_c = 1/1.55;

figure(1)
plot(eta,Z_MD,'ro')
grid on;
xlabel('\eta');
ylabel('Z');

figure(2)
plot(eta,Z_MD - (1 + 3*eta./(eta_c - eta)),'ro')
grid on;


A2 = [(eta./eta_c).^(1:5).*(1:5)]\(Z_MD - (1 + 3*eta./(eta_c - eta)));

figure(3)
Zpred = [(eta./eta_c).^(1:5).*(1:5)]*A2;
plot(eta,Z_MD - (1 + 3*eta./(eta_c - eta)),'ro')
hold on
grid on
plot(eta,Zpred,'b-')

figure(4)
plot(eta,Z_MD - (1 + 3*eta./(eta_c - eta)) - Zpred,'ro')
grid on


mmfun = @(A) abs(Z_MD - 1 - 3*eta./(eta_c - eta) - [(eta./eta_c).^(1:5).*(1:5)]*A);
[Amm,~,maxres,exitflag] = fminimax(mmfun,A2);

figure(5)
Zpredmm = [(eta./eta_c).^(1:5).*(1:5)]*Amm;
plot(eta,Z_MD - (1 + 3*eta./(eta_c - eta)) - Zpredmm,'ro')
grid on

figure(6)
plot(eta,Z_MD,'ro',eta,(1 + 3*eta./(eta_c - eta)) + Zpredmm,'b-',eta,(1 + 3*eta./(eta_c - eta)) + Zpred,'g-')
grid on

[John D'Errico]

1) If I want to add a 5th term, do I just change the 1:4 to 1:5?

Yes. Nothing in what I did was specific to 4 terms.

Beware a little though. Higher order terms tend to do bad things in polynomials, if you get carried away. Order 5 is not too bad though, as eta is not large.

2) How can I set A1 = -0.4194, A2 = 0.5812, and A3 = 0.6439?

You do exactly the same as I did with the first two terms. I subtracted them from Z_MD. Now your model will only require terms for the 4th and 5th powers, so there are only two coefficients to estimate.

So a least squares estimate for the other terms looks like:

A123 = [-0.4194, 0.5812, 0.6439]';
Zfixed = 1 + 3*eta./(eta_c - eta) + (eta/eta_c).^[1 2 3].*[1 2 3]*A123;

A45 = [(eta./eta_c).^(4:5).*(4:5)]\(Z_MD - Zfixed)

mmfun = @(A) abs(Z_MD - Zfixed - [(eta./eta_c).^(4:5).*(4:5)]*A);
[Amm,~,maxres] = fminimax(mmfun,A45)

Amm =
       4.3708
      -4.3919
maxres =
      0.24247
      

So residuals that are close to double those you saw in my answer. But here we estimated only TWO coefficients, not 4.

Basically by doing this, the correct behavior of the curve is captured as eta -> 0.

So I want these to be set to these values, and not let them vary.

Also note that at eta=0, then Z should =1. Perhaps these constraints I have mentioned above will help that be the case.

No, they won't help, nor will they hurt. They simply won't matter.

In the function

Z = 1 + 3*eta/(eta_c - eta) + A1*(eta/eta_c) + ...

when eta is zero, the model is identically 1. You can pick any finite set of values for the coefficients A1,A2,A3,A4,A5,... and it will ALWAYS be 1.

However, it is important to note the data may not be consistent with that model.

Then I can allow A4 to vary. This may not be a great fit, but then I can subsequently add more terms (A5, A6, etc) until I get the fit that I want.

Do not go too far overboard. High order polynomials tend to do nasty things.

[Benjamin]

Attached is the residuals I get for the fminimax, does this look like what you got?

I also had another question. This seems to minimize the max error in terms of (Z_MD - fit). Is there a way to minimize the relative error, i.e. (fit - Z_MD)/Z_MD

Basically, I think it is currensly set up to minimize the max error in terms of the absolute error. How do I get it to do the same thing but for relative error?

[Matt J]

This seems to minimize the max error in terms of (Z_MD - fit). Is there a way to minimize the relative error, i.e. (fit - Z_MD)/Z_MD

Isn't it obvious that you would just change mmfun, accordingly? Although, I still think it is inadvisable to include the abs() in mmfun.

[Benjamin]

Are you suggesting I remove abs from mmfun?

Without removing abs,

I changed

mmfun = @(A) abs(Z_MD - Zfixed - [(eta./eta_c).^(4:5).*(4:5)]*A);

to

mmfun = @(A) abs(((Zfixed + [(eta./eta_c).^(4:5).*(4:5)]*A)-Z_MD)/Z_MD);

but after:

Zpredmm = [(eta./eta_c).^(4:5).*(4:5)]*Amm;

and I plot this:

plot(eta, ((Zfixed + Zpredmm)-Z_MD)./Z_MD,'ro')

I get the figure attached (rel_error).

When I plot the old version:

plot(eta,Z_MD - Zfixed - Zpredmm,'ro')

but still using this version of mmfun,

I get the figure attached (abs_error).

It seems like it is still minimizing the max absolute error right? At least that is what the graphs seems to say.

The graph of the relative error varies from -0.02 to 0.06, whereas the one of absolute error varies -0.25 to 0.25. So while I think that I minimized the relative error, why does it still seem to minimizing the absolute error instead?

[Benjamin]

Any chance you could comment on my last questions?

Should I remove abs?

And why does my new mmfun not minimize the max relative error? It still seems to minimize the absolute error.

Also, just curious, how does:

A45 = [(eta./eta_c).^(4:7).*(4:7)]\(Z_MD - Zfixed);

solve for the coeffs? I mean clearly it does, but you don't have A_k anywhere in the actual equation. How does it know to create these coefficients?

[Catalytic]

You should remove abs and instead use AbsoluteMaxObjectiveCount (that's what it's meant for). I don't know why John's version works without it. Freakish luck, I suppose.

Also, you shouldn't set AbsoluteMaxObjectiveCount blindly to some value like 1. Read about what it does and how to use it!

Answer 2

[Catalytic]

Consider using fminimax.

[Benjamin]

Can you should me the general way to set this up with my equation? Here is the equation in MATLAB. These A values are set here, but I want them to be able to change in order to minimize the max error.

Why does the following not work? I'm not sure how to get my polynomial fit in the form I showed in my question with the proper A values. This just gives me very small numbers in fval.

clear; clc; close all;
mat =xlsread('C:\example2.xlsx','Sheet1','A2:C32');

eta = mat(:,2);
Z_MD = mat(:,3);


x0 = [-0.4194 0.5812 0.6439 0.4730];
[x,fval] = fminimax(@Z_fit,x0);

function Z = Z_fit(eta)
    x0 = [-0.4194 0.5812 0.6439 0.4730];
    eta_c = 0.6452;
    rpf = eta./eta_c;
    Z = (1+3.0*eta/(eta_c-eta)+x0(1)+rpf.*(2.0*x0(2)+rpf.*(3.0*x0(3)+rpf.*4.0*x0(4)))) * (1.0/4.0 * 1/eta_c) ;
end

[Catalytic]

Can you should me the general way to set this up with my equation?

If I showed you how to set it up with your equation, it wouldn't be "general". Here is a relevant example, though, from the fminimax documentation https://www.mathworks.com/help/optim/ug/fminimax.html#mw_af58b8ea-61b8-4369-a11e-a34798d9535d

Do not simply copy/paste everything and expect it to apply out-of-the-box to your problem. In particular, once you have read the example, you should see that you will have to adapt the way that AbsoluteMaxObjectiveCount is used in your case.

[Benjamin]

The below code is what I put together. I'm not sure if what I'm doing is correct or not though. What should I be changing AbsoluteMaxObjectiveCount to? The fit to my data is not stellar, though its not absolutely miserable either. Any chance you could explain what I am missing here? The maximum error does not seem to be minimized here.

clear; clc; close all;
mat =xlsread('C:\example2.xlsx','Sheet1','A2:C32');

eta = mat(:,2);
Z_MD = mat(:,3);
eta_c = 1/1.55;
rpf = eta./eta_c;


%Initial Guess
x0 = [-0.4194 0.5812 0.6439 0.4730];

%Set options
options = optimoptions('fminimax','AbsoluteMaxObjectiveCount',1);

%Define function
fun = @(eta)[(1+3.0*eta/(eta_c-eta)+x0(1)+rpf.*(2.0*x0(2)+rpf.*(3.0*x0(3)+rpf.*4.0*x0(4)))) * (1.0/4.0 * 1/eta_c)];


A = []; % No constraints
b = [];
Aeq = [];
beq = [];
lb = [];
ub = [];
nonlcon = [];

%Solve for values of x1 to minimize the max error.
x1 = fminimax(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon,options);

[John D'Errico]

@Benjamin: read my answer.