How do I hold some coefficients constant in a vector of initial guesses passed to lsqcurvefit?

I am working on a project for my year in industry as part of my university course. I have been tasked with automating a fairly user intensive process that involves many iterations of different initial guesses being passed to a levenberg marquardt solver (companies own software) to find the best fit. These initial guesses take the form K = [k1,o1;k2,o2,k3,o3;...;k22,o22] where k1:k22 are the actual initial guesses and o1:o22 is logical vector of the same length denoting whether or not the corresponding k value should be optimised or held constant during the optimisation.

I have my program working using lsqcurvefit for the situations where o1:o22 are all true but when I tried to add upper and lower bounds it tells me that the upper and lower bounds cannot be the same, which was how i intended to hold these values constant. Can anyone help me think of another way around this as I have invested a lot of time getting the program up to this point and have now been stuck on this particular problem for a week or so.

any help would be much appreciated, thank you.

Answer 1

[Matt J]

It's strange that it tells you that ub cannot equal lb, because I know that you can do this for other Optimization Toolbox functions, such as FMINCON and maybe LSQNONLIN, which you could try as alternatives.

However, since you say that this is a computationally intensive process that you are running many times, the most speed-optimal approach is to use o1:o22 to automatically rewrite the curve modelling function as a function of fewer variables. Lower-dimensional functions do crunch (and optimize) faster.

For example, if your curve was instead a simple 3-parameter quadratic

k(i)*x^2+k(2)*x+k(3)

you wouldn't want to include the polynomial terms where o(i)=0 because that would just waste computation. So you would do

function [F_k_xdata,k0]=CreateFitFun(o,k0_full)
%return handle to curve fit function and initial point in reduced space  

   A=bsxfun(@power,xdata(:),[2 1 0]);

    b=A(:,~o)*k0_full(~o);

    A=A(:,o); throw away fixed parameters  
    k0=k0_full(o);

   F_k_xdata=@(k,xdata) A*k(:)+b;

end

and then pass this F_k_data handle and initial point k0 to LSQCURVEFIT.

[Stevie]

thanks for answering.

firstly I should clarify that when i said user intensive i was referring to the old process and what i meant was that the user needed to manually copy and paste a lot of tedious information to and from excel for the program to find an answer. my project is to streamline it all within one project. the actual work the computer is doing isn't that intensive, just finding least squares fit and a bit of interpolating.

secondly when i say that the o1:o22 vector controls whether or not the corresponding k element should be optimised it is possible, indeed likely, that k is not just zero. it could be a fixed constant that would still need to be passed as a coefficient to the fitting function.

i did think about your solution in general though, splitting the coefficients for optimisation and the coefficients to be held constant into two separate vectors however this becomes very messy when i try to apply them to the correct locations inside the function.

if you say equal upper and lower bounds are possible in other functions then the error may also be to do with the fact that it claims Levenberg Marquardt solver does not accept upper or lower bounds at all, but it then does not even let me use the same upper and lower bounds with the trust region solver. and i also know of no reason why a levenberg marquardt solver could not have bounds applied to it.

[Matt J]

I didn't intend to assume k(~o)=0. I modified my original proposed code accordingly. As for messiness, I'd have to see your actual function to get a better idea of the difficulty. In any case, when you need speed, that's what you do...

Again, I cannot understand why you are limited from lb=ub, if the documentation doesn't forbid it, but you could try FMINCON if that's an easier path. It shouldn't have this restriction.

As for why Levenberg-Marquardt won't take bounds at all, you would have to solve a constrained quadratic sub-problem at every iteration and there is a likelihood of slow convergence. I can imagine TMW deeming it not worth the coding effort.