Constrained linear least squares not fitting data

Question

I am trying to fit a 3D surface polynomial of n-degrees to some data points in 3D space. My system requires the surface to be monotonically increasing in the area of interest, that is the partial derivatives must be non-negative. This can be achieved using Matlab's built in lsqlin function.

So here's what I've done to try and achieve this:

I have a function that takes in four parameters; x1 and x2 are my explanatory variables and y is my dependent variable. Finally, I can specify order of polynomial fit. First I build the design matrix A using data from x1 and x2 and the degree of fit I want. Next I build the matrix D that is my container for the partial derivatives of my datapoints. NOTE: the matrix D is double the length of matrix A since all datapoints must be differentiated with respect to both x1 and x2. I specify that Dx >= 0 by setting b to be zeroes.

Finally, I call lsqlin. I use "-D" since Matlab defines the function as Dx <= b.

function w_mono = monotone_surface_fit(x1, x2, y, order_fit)

% Initialize design matrix
A = zeros(length(x1), 2*order_fit+2);

% Adjusting for bias term
A(:,1) = ones(length(x1),1); 


% Building design matrix
for i = 2:order_fit+1
    A(:,(i-1)*2:(i-1)*2+1) = [x1.^(i-1), x2.^(i-1)];
end

% Initialize matrix containing derivative constraint.
% NOTE: Partial derivatives must be non-negative
D = zeros(2*length(y), 2*order_fit+1);


% Filling matrix that holds constraints for partial derivatives
% NOTE: Matrix D will be double length of A since all data points will have a partial derivative constraint in both x1 and x2 directions. 
for i = 2:order_fit+1
     D(:,(i-1)*2:(i-1)*2+1) = [(i-1)*x1.^(i-2), zeros(length(x2),1); ...
                                 zeros(length(x1),1), (i-1)*x2.^(i-2)];
end

% Limit of derivatives
b = zeros(2*length(y), 1);

% Constrained LSQ fit
options = optimoptions('lsqlin','Algorithm','interior-point');

% Final weights of polynomial
w_mono = lsqlin(A,y,-D,b,[],[],[],[],[], options);

end

So now I get some weights out, but unfortunately they do not at all capture the structure of the data. I've attached an image so you can just how bad it looks. .

I'll give you my plotting script with some dummy data, so you can try it.

%% Plot different order polynomials to data with constraints

x1 = [-5;12;4;9;18;-1;-8;13;0;7;-5;-8;-6;14;-1;1;9;14;12;1;-5;9;-10;-2;9;7;-1;19;-7;12;-6;3;14;0;-8;6;-2;-7;10;4;-5;-7;-4;-6;-1;18;5;-3;3;10];
x2 = [81.25;61;73;61.75;54.5;72.25;80;56.75;78;64.25;85.25;86;80.5;61.5;79.25;76.75;60.75;54.5;62;75.75;80.25;67.75;86.5;81.5;62.75;66.25;78.25;49.25;82.75;56;84.5;71.25;58.5;77;82;70.5;81.5;80.75;64.5;68;78.25;79.75;81;82.5;79.25;49.5;64.75;77.75;70.25;64.5];
y = [-6.52857142857143;-1.04736842105263;-5.18750000000000;-3.33157894736842;-0.117894736842105;-3.58571428571429;-5.61428571428572;0;-4.47142857142857;-1.75438596491228;-7.30555555555556;-8.82222222222222;-5.50000000000000;-2.95438596491228;-5.78571428571429;-5.15714285714286;-1.22631578947368;-0.340350877192983;-0.142105263157895;-2.98571428571429;-4.35714285714286;-0.963157894736842;-9.06666666666667;-4.27142857142857;-3.43684210526316;-3.97894736842105;-6.61428571428572;0;-4.98571428571429;-0.573684210526316;-8.22500000000000;-3.01428571428571;-0.691228070175439;-6.30000000000000;-6.95714285714286;-2.57232142857143;-5.27142857142857;-7.64285714285714;-2.54035087719298;-3.45438596491228;-5.01428571428571;-7.47142857142857;-5.38571428571429;-4.84285714285714;-6.78571428571429;0;-0.973684210526316;-4.72857142857143;-2.84285714285714;-2.54035087719298];

% Used to plot the surface in all points in the grid
X1 = meshgrid(-10:1:20);
X2 = flipud(meshgrid(30:2:90).');

figure;
for i = 1:4

    w_mono = monotone_surface_fit(x1, x2, y, i);

    y_nr = w_mono(1)*ones(size(X1)) + w_mono(2)*ones(size(X2));

    for j = 1:i
        y_nr = w_mono(j*2)*X1.^j + w_mono(j*2+1)*X2.^j;
    end

    subplot(2,2,i);
    scatter3(x1, x2, y); hold on;

    axis tight;

    mesh(X1, X2, y_nr);
    set(gca, 'ZDir','reverse');

    xlabel('x1'); ylabel('x2');
    zlabel('y');
%     zlim([-10 0])
end

I think it may have something to do with the fact that I haven't specified anything about the region of interest, but really I don't know. Thanks in advance for any help.

Answer 1

Alright I figured it out.

The main problem was simply an error in the plotting script. The value of y_nr should be updated and not overwritten in the loop.

Also I figured out that the second derivative should be monotonically decreasing. Here's the updated code if anybody is interested.

%% Plot different order polynomials to data with constraints

x1 = [-5;12;4;9;18;-1;-8;13;0;7;-5;-8;-6;14;-1;1;9;14;12;1;-5;9;-10;-2;9;7;-1;19;-7;12;-6;3;14;0;-8;6;-2;-7;10;4;-5;-7;-4;-6;-1;18;5;-3;3;10];
x2 = [81.25;61;73;61.75;54.5;72.25;80;56.75;78;64.25;85.25;86;80.5;61.5;79.25;76.75;60.75;54.5;62;75.75;80.25;67.75;86.5;81.5;62.75;66.25;78.25;49.25;82.75;56;84.5;71.25;58.5;77;82;70.5;81.5;80.75;64.5;68;78.25;79.75;81;82.5;79.25;49.5;64.75;77.75;70.25;64.5];
y = [-6.52857142857143;-1.04736842105263;-5.18750000000000;-3.33157894736842;-0.117894736842105;-3.58571428571429;-5.61428571428572;0;-4.47142857142857;-1.75438596491228;-7.30555555555556;-8.82222222222222;-5.50000000000000;-2.95438596491228;-5.78571428571429;-5.15714285714286;-1.22631578947368;-0.340350877192983;-0.142105263157895;-2.98571428571429;-4.35714285714286;-0.963157894736842;-9.06666666666667;-4.27142857142857;-3.43684210526316;-3.97894736842105;-6.61428571428572;0;-4.98571428571429;-0.573684210526316;-8.22500000000000;-3.01428571428571;-0.691228070175439;-6.30000000000000;-6.95714285714286;-2.57232142857143;-5.27142857142857;-7.64285714285714;-2.54035087719298;-3.45438596491228;-5.01428571428571;-7.47142857142857;-5.38571428571429;-4.84285714285714;-6.78571428571429;0;-0.973684210526316;-4.72857142857143;-2.84285714285714;-2.54035087719298];

% Used to plot the surface in all points in the grid
X1 = meshgrid(-10:1:20);
X2 = flipud(meshgrid(30:2:90).');

figure;
for i = 1:4

    w_mono = monotone_surface_fit(x1, x2, y, i);

    % NOTE: Should only have 1 bias term
    y_nr = w_mono(1)*ones(size(X1));

    for j = 1:i
        y_nr = y_nr + w_mono(j*2)*X1.^j + w_mono(j*2+1)*X2.^j;
    end

    subplot(2,2,i);
    scatter3(x1, x2, y); hold on;

    axis tight;

    mesh(X1, X2, y_nr);
    set(gca, 'ZDir','reverse');

    xlabel('x1'); ylabel('x2');
    zlabel('y');
%     zlim([-10 0])
end

And here's the updated function

function [w_mono, w] = monotone_surface_fit(x1, x2, y, order_fit)

% Initialize design matrix
A = zeros(length(x1), 2*order_fit+1);

% Adjusting for bias term
A(:,1) = ones(length(x1),1); 

% Building design matrix
for i = 2:order_fit+1
    A(:,(i-1)*2:(i-1)*2+1) = [x1.^(i-1), x2.^(i-1)];
end

% Initialize matrix containing derivative constraint.
% NOTE: Partial derivatives must be non-negative
D = zeros(2*length(y), 2*order_fit+1);

for i = 2:order_fit+1
    D(:,(i-1)*2:(i-1)*2+1) = [(i-1)*x1.^(i-2), zeros(length(x2),1); ...
        zeros(length(x1),1), -(i-1)*x2.^(i-2)];
end

% Limit of derivatives
b = zeros(2*length(y), 1);

% Constrained LSQ fit
options = optimoptions('lsqlin','Algorithm','active-set');
w_mono = lsqlin(A,y,-D,b,[],[],[],[],[], options);
w = lsqlin(A,y);

end

Finally a plot of the fitting (Have used a new simulation, but fit also works on given dummy data).