Problem with fit Model with odd order polynomial

Hi, I have the data as below.

Its fits OK'ish to a 3rd order (odd polynomial), but then becomes terrible when trying a 5th order (odd) polynomial

3rd order fit below

and here is the 5th order

Here is my code:

        ax=app.UIAxes3;
        [x,y] = getDataFromGraph(app,ax,1);  % My function which grabs the data from the plot

        % 3rd order (odd)
        ft=fittype('a*x^3+b*x');            %    ft=fittype('a*x^5+b*x^3+c*x');
        [fitobj,~,~,~]=fit(x,y,ft);         %[fitobj, goodness, output, convmsg]=fit(x,N(:,i),ft)
        coeff1=fitobj.a;
        coeff2=fitobj.b;        

        yfit=coeff1*x.^3+coeff2*x;             
        r=y-yfit; 
        hold(ax,'on');
        plot(ax,x,yfit,'r--');

        %5th order (odd)
        ft=fittype('a*x^5+b*x^3+c*x');
        [fitobj,~,~,~]=fit(x,y,ft);    %[fitobj, goodness, output, convmsg]=fit(x,N(:,i),ft)
        coeff1=fitobj.a;
        coeff2=fitobj.b;   
        coeff3=fitobj.c;   
        yfit=coeff1*x.^5+coeff2*x.^3+coeff3*x;     %Remember dot notation
        plot(ax,x,yfit,'g--');

Answer 1

[Matt J]

Don't use a custom fit type. Use the builtin poly5 fit type, with bounds on the even degree coefficients.

[x,y] = readvars('DistortionTableData.csv');  % My function which grabs the data from the plot

lb=[-inf,0,-inf,0,-inf,0];

% 3rd order (odd)
ft=fittype('poly3');                        
[fitobj,~,~,~]=fit(x,y,ft,'Lower',lb(1:4),'Upper',-lb(1:4),'Normalize','on')      

fitobj = 
     Linear model Poly3:
     fitobj(x) = p1*x^3 + p2*x^2 + p3*x + p4
       where x is normalized by mean 1.285e-14 and std 2024
     Coefficients (with 95% confidence bounds):
       p1 =       0.286  (0.2638, 0.3082)
       p2 =           0  (fixed at bound)
       p3 =     -0.1792  (-0.2227, -0.1357)
       p4 =           0  (fixed at bound)

plot(fitobj,x,y)

%5th order (odd)
ft=fittype('poly5');
[fitobj,~,~,~]=fit(x,y,ft,'Lower',lb,'Upper',-lb,'Normalize','on')

fitobj = 
     Linear model Poly5:
     fitobj(x) = p1*x^5 + p2*x^4 + p3*x^3 + p4*x^2 + p5*x + p6
       where x is normalized by mean 1.285e-14 and std 2024
     Coefficients (with 95% confidence bounds):
       p1 =     0.04478  (0.01573, 0.07382)
       p2 =           0  (fixed at bound)
       p3 =      0.1371  (0.03795, 0.2362)
       p4 =           0  (fixed at bound)
       p5 =    -0.08362  (-0.1592, -0.008065)
       p6 =           0  (fixed at bound)

plot(fitobj,x,y)

[Jason]

Thanks Matt. The reason im not using poly5 is i need even powers of x to be zero

And why use data normalisation?

Thanks

[Matt J]

But as I've demonstrated, we can constrain the even powers to zero with bounds.

I suggested data normalization because your x -data is of very different order of magnitude than your y-data. However, my tests don't show much difference in this case.

[Jason]

Thankyou

Is there any reason why you use -inf in lb and then used -lb for the upper bound.

Could you equally have done

lb=[inf,0,inf,0,inf,0];

With

fitobj,~,~,~]=fit(x,y,ft,'Lower',-lb,'Upper',lb)

[Matt J]

No, there's no difference.

[Jason]

Thanks

[Matt J]

I changed my mind in light of the conversation in Torsten's answer thread. You should definitely use normalization. I've now edited my original answer accordingly.

Answer 2

[Torsten]

In the case of the polynomial of degree 5, the design matrix is rank-deficient:

        
        T = readmatrix("DistortionTableData.csv");
        x = T(:,1);
        y = T(:,2);
        
        % 3rd order (odd)
        A = [x.^3,x];
        rank(A)

ans = 2

        b = y;
        sol = A\b;
        a = sol(1);
        b = sol(2);
        figure(2)
        hold on
        plot(x,a*x.^3+b*x)
        plot(x,y)
        hold off

●

        %5th order (odd)
        A = [x.^5,x.^3,x];
        rank(A)

ans = 2

        b = y;
        sol = A\b;

Warning: Rank deficient, rank = 2, tol = 4.322070e+05.

        a = sol(1);
        b = sol(2);
        c = sol(3);
        figure(4)
        hold on
        plot(x,a*x.^5+b*x.^3+c*x)
        plot(x,y)
        hold off

[Matt J]

But you wouldn't use such a direct inversion as A\b. You would use a QR decomposition:

[x,y] = readvars('DistortionTableData.csv');  % My function which grabs the data from the plot

 A = [x.^5,x.^3,x];
 
 [Q,R]=qr(A,0);
 
 sol = R\(Q.'*y)

sol = 3×1

1.0e-04 * 0.0000 0.0000 -0.4130

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

 rank(R)

ans = 3

[Jason]

Thanks Torsten. Whats this method of finding the coefficients called so i can read further (i.e sol=A/b)

[Torsten]

@Matt J

Did you look at cond(R) ?

According to the flow chart, "mldivide" uses the QR solver in case of a non-square matrix A:

https://uk.mathworks.com/help/matlab/ref/double.mldivide.html

@Jason

Whats this method of finding the coefficients called so i can read further (i.e sol=A/b)

It's just solving the overdetermined linear system of equations for the unknown parameter vector in the least-squares sense:

a*x.^5 + b*x.^3 + c*x = y

[Matt J]

Did you look at cond(R) ?

We know without looking at cond( R ) that it will be the same as cond(A), but since R is triangular, the error amplification will not be as bad:

[x,~] = readvars('DistortionTableData.csv');  % My function which grabs the data from the plot

c=[1;2;3]; %ground truth coefficients

A = [x.^5,x.^3,x];

[errMLD, errQR]=compareSolvers(A,c)

Warning: Rank deficient, rank = 2, tol = 4.322070e+05.

errMLD = 3.0000

errQR = 0.0017

In any case, the best thing to do here would probably be to normalize the x-data, which improves cond(A) greatly:

x=x/4000;
A = [x.^5,x.^3,x];

cond(A)

ans = 38.6801

[errMLD, errQR]=compareSolvers(A,c)

errMLD = 9.9301e-16

errQR = 1.2608e-14

function [errMLD, errQR]=compareSolvers(A,c)

 y=A*c;

 [Q,R]=qr(A,0);
 
 errMLD=norm(A\y-c);      %error using direct mldivide
 errQR=norm(R\(Q'*y)-c);  %error using QR
 
end