Slim code for smoothing spline fit with predefined knots, concavity and monotonicity

I want a (spline based) function in which I can input a number of knots (or a set of knots) on noisy xy data, which may have 1 corner point of varying curvature.

• This function should also allow setting concavity/convexity constraints, and monotonically increasing / decreasing profile data.
• Code should be fairly short and suitable for compilation with MATLAB Coder (i.e. some Matlab commands are not applicable).

Available routines on the forum describes similar polynomial based solutions, which I tried and works well on the data but may need 7-8 order polynomials. Available solutions like SLMENGINE works well, but is a large code base and I need something short and "easy to work with for Matlab Coder", and hence avoid CSAPS, fnder, fnval ,etc. On the other hand, spline, ppval, mkpp,etc are ok to use with Coder.

Here is what I have managed so far. If I remove the constraints CSAPS works as expected, activating the C-constraints produce individually constrained segments as defined but they are not connected to form a continuous spline. Maybe the C constraints should be replaced by the A/B constraints but it is not clear to me how to do that.

function pp = fitSplineWithConstraints2(x, y, num_knots)
    % Fit a spline to data with monotonicity and concave shape constraints
    % Inputs:
    %   x - independent variable data
    %   y - dependent variable data
    %   num_knots - number of knots for the spline
    % Output:
    %   pp - piecewise polynomial structure

    % Define the knots for the spline
    knots = round(linspace(1, length(x), num_knots));
    xknots = x(knots);
    yknots = y(knots);
        
    % Initial guess for spline coefficients
    pp_init = spline(xknots, yknots);  
    p0 = pp_init.coefs(:);

    % Define the objective function (least squares error)
    % p here is pp.coefs:
    objective = @(p) sum((evaluateSpline(xknots, p, x) - y).^2);

    A = [];
    b = [];
    Aeq = [];
    beq = [];
    lb = [];
    ub = [];
    nonlcon = @(p) constraints(xknots, p, x);

    % Perform constrained optimization
    options = optimoptions('fmincon', 'Display', 'off','Algorithm','sqp');
    [p_opt,fval,exitflag,output] = fmincon(objective, p0, A, b, Aeq, beq, lb, ub, nonlcon, options);
    
   
    % Construct the piecewise polynomial structure
    pp = mkpp(xknots, reshape(p_opt, [], num_knots-1), 1);

    figure
    plot(xknots, yknots)
     hold on
     plot(xknots,ppval(pp_init, xknots),'--')
     plot(xknots,ppval(pp, xknots),'-.')

end

function y = evaluateSpline(knots, p, x)
    % Evaluate the spline at given points
    pp = mkpp(knots, reshape(p, [], length(knots)-1), 1);
    y = ppval(pp, x);
end

function [c, ceq] = constraints(knots, p, x)
    % Nonlinear constraints for monotonicity and concavity
    % Inputs:
    %   knots - spline knots
    %   p - spline coefficients
    %   x - independent variable data
    % Outputs:
    %   c - inequality constraints (c <= 0)
    %   ceq - equality constraints (ceq == 0)

    % Construct the piecewise polynomial structure
    pp = mkpp(knots, reshape(p, [], length(knots)-1),1);

    % Evaluate the first and second derivatives of the spline
    dp = fnder(pp, 1);
    ddp = fnder(pp, 2);

    % Evaluate the derivatives at the data points
    dp_vals = ppval(dp, x);
    ddp_vals = ppval(ddp, x);

    % Monotonicity constraint (1st derivative): 
    % set "-" for increasing
    c1 = -dp_vals;

    % Concavity constraint (2nd derivative): 
    % set "+" concave down
    c2 = ddp_vals;

    % Combine the constraints
    c = [c1;c2];
    ceq = [];
end

[Daniel]

Output example showing the discontinuities at the knots.

Answer 1

[Matt J]

One path would be to obtain the matrix form for the spline interpolator, which I do in the example below using func2mat from the FEX download,

https://www.mathworks.com/matlabcentral/fileexchange/44669-convert-linear-functions-to-matrix-form

Once you have that, the whole thing can be done as a linear least squares problem using lsqlin.

xknots = linspace(0, 5, 8); % Knot locations
xdata = linspace(0, 5, 200); %Finer sample x-locations

C=func2mat(  @(z) spline(xknots,z,xdata),  xknots); %spline interpolator as matrix


% Generate data for a convex monotonic function
f = @(x) x.^2 + log(1 + x);
ydata=f(xdata); 
  ydata=ydata+randn(size(ydata))*0.8; %samples+noise

yfit=doFit(C,ydata); %perform the fit

plot(xdata,ydata,'bx',xdata,yfit,'-r'); legend 'Data Points' 'Convex Monotonic Fit'

function [yfit,coeffs]=doFit(C,d)
%Fit with a convex monotonically increasing  spline

    ddx=diff(C,1,1);        %First order derivative operator, for monotonicity constraints
    d2dx2=diff(ddx,1,1);    %2nd order derivative operator, for convexity constraints

    Aineq=[-ddx;-d2dx2]; bineq=zeros(height(Aineq),1);

    opts=optimoptions('lsqlin','Display','none');
    coeffs=lsqlin(C,d(:),Aineq,bineq,[],[],[],[],[],opts);
    yfit=C*coeffs;

end

[Matt J]

which I do in the example below using func2mat from the FEX download

You could also improve performance somewhat using localfunc2mat from the same FEX submission, and/or by reusing the C matrix where possible.

[Daniel]

Thanks Matt, this is indeed a clean approach.