Error when using lsqnonlin

I use this function to approximate measured data with a given function. I have earlier with succes performed the same procedure with simulated data.

Intensity represents the measured data.

              A_guess_01=0.8;
              kappa_guess_01=0.5;
              sigma_guess_01=0.8;
              A_guess_11=0.8;
              phasediff_guess_11=7*pi/4;
              kappa_guess_11=0.8;
              sigma_guess_11=1.4;
              Theta_turn=-pi/4;
              A_guess_11_sin=0;

              start_point(1,1) = A_guess_01;
              start_point(1,2) = kappa_guess_01;
              start_point(1,3) = sigma_guess_01;
              start_point(1,4)=A_guess_11;
              start_point(1,5)=phasediff_guess_11;
              start_point(1,6)=kappa_guess_11;
              start_point(1,7)=sigma_guess_11;
              start_point(1,8)=Theta_turn;
              start_point(1,9)=A_guess_11_sin;

              st=start_point;

              Hej=  @(p) abs((p(1)*besselj(0,p(2)*r)+besselj(1,p(6)*r).*(p(4)*cos(theta+p(8))+p(9)*sin(theta+p(8)))*exp(1i*p(5))).*(r/a<=1)...
                              + (p(1)*besselj(0, p(2)*a)/besselk(0, p(3)*a)*besselk(0,p(3)*r)+...                               (p(4)*besselj(1, p(6)*a)/besselk(1, p(7)*a)*cos(theta+p(8))...
+p(9)*besselj(1, p(6)*a)/besselk(1,p(7)*a)*sin(theta+p(8))).*besselk(1, p(7)*r).*exp(1i*p(5)))...
.*(r/a>1)).^2- Intensity;

   opts = optimoptions(@lsqnonlin,'DiffMaxChange', 0.05,'FinDiffType', 'central', 'Display','off','MaxFunEvals',2E7,'TolFun',1E-25, 'TolX',1E-25,'MaxIter',4E4);
          x0 = st; % arbitrary initial guess
          lb = 0.0*ones(size(st));
          ub = 5*ones(size(st));
          [p_estimated,resnorm,residuals,exitflag,output] = lsqnonlin(Hej,x0, lb,[], opts);

When I run this I get the following error.

      Error using snls (line 47)
      Objective function is returning undefined values at initial point. lsqnonlin cannot continue.

    Error in lsqncommon (line 149)
        [xC,FVAL,LAMBDA,JACOB,EXITFLAG,OUTPUT,msgData]=...

    Error in lsqnonlin (line 236)
    [xCurrent,Resnorm,FVAL,EXITFLAG,OUTPUT,LAMBDA,JACOB] = ...

    Error in Mode_decomposition_multiple_modes_intensity_trial (line 189)
            [p_estimated,resnorm,residuals,exitflag,output] = lsqnonlin(Hej,x0, lb,[],
            opts);%optimset('Display','off','MaxFunEvals',2E6,'TolFun',1E-9,'TolX',1E-9,'MaxIter',2E3));%,opts);

Both the input Intensity and the initial guess yields finite values.

All inputs to this error will be highly appreciated.

Answer 1

[Matt J]

Both the input Intensity and the initial guess yields finite values.

Not when I run your code. You have not defined r, theta, a, and perhaps other quantities used by Hej. Therefore calling Hej(x0) in isolation immediately yields an error.

[Stine Larsen]

Hi Matt,

You cannot run the code as I had not provided the input data. It is as follows

data=load('trial.mat');

X=data.X;
Y=data.Y;
% Here we consider the intensity in form of what is actually measured, thus
% we take the square root (to use the same format as the simulated modes).
Intensity=data.selected_data_1(:, :, 1);
%M=data.M;

% Making a square matrix

[size_min, index]=min(size(X));

move=floor((max(size(X))-size_min)/2);

if index==1
    X=X(:, (move+1):(size_min+move));
    Y=Y(:, (move+1):(size_min+move));
    Intensity=Intensity(:, (move+1):(size_min+move));
else
    X=X((move+1):(size_min+move), :);
    Y=Y((move+1):(size_min+move), :);
    Intensity=Intensity((move+1):(size_min+move), :);
end

%%Normalizing the input

Intensity=Intensity/max(max(Intensity));

%%Subtracting noise

figure
 contourf(X,Y, Intensity, 200, 'LineStyle', 'none')
            colormap hot
            axis equal
            colorbar

noise=max(reshape(Intensity(1:10,:), 1, 10*size(Intensity(:,1:10), 1)));

Intensity=Intensity-noise*ones(size(Intensity));

% Setting all negative values to 0

ll=find(Intensity<0);

Intensity(ll)=0;
[theta,r] = cart2pol(X,Y);

a=3;

The .mat file is attached.

The rest of the code has been submitted in previous posts.

Stine

[Matt J]

Nope, I still get an error

    Reference to non-existent field 'selected_data_1'.

    Error in test1 (line 9)
                    Intensity=data.selected_data_1(:, :, 1);

The data structure in your trial.mat file contains only the fields X,Y, and Intensity.

[Stine Larsen]

A new and correct (!) version of trial.mat has been uploaded.

[Matt J]

When I run with this revised trial.mat, and with the options revised as below,I get no errors.

 opts = optimoptions(@lsqnonlin,'DiffMaxChange', 0.05,'FinDiffType', 'central', 'Display','iter','MaxFunEvals',2E7,'TolFun',1E-10, 'TolX',1E-10,'MaxIter',4E4);

LSQNONLIN then shows the following iterative activity,

                                         Norm of      First-order 
 Iteration  Func-count     f(x)          step          optimality   CG-iterations
     0         19          323484                      5.74e+05
     1         38         4404.92        2.05096       5.26e+03            0
     2         57         2975.75       0.242177            599            0
     3         76         2773.48       0.237737            126            0
     4         95         2711.83        0.56732           37.8            0
     5        114         2700.84       0.273689           12.9            0
     6        133         2700.11       0.186991           12.8            0
     7        152         2700.09      0.0573804          0.947            0
     8        171         2700.09     0.00536302          0.013            0
     9        190         2700.09    5.15508e-05       0.000124            0

[Stine Larsen]

Thank you for trying. I still get the same mistake. I use Matlab 2014. Should that make a difference?

[Matt J]

When you evaluate Hej(x0) in isolation from the rest of the code, what value do you get?

[Stine Larsen]

I get only real numbers. In the figure Hej(x0) is plotted.

[Matt J]

Time for

 >> dbstop if error

I guess.

[Stine Larsen]

It stops at the following

% Check to see if the function vector or the user-supplied Jacobian at the initial point are well-defined.
% If not, then terminate immediately.
if any(~isfinite(fval)) 
    error(message('optimlib:snls:UsrObjUndefAtX0', caller));    
end

which makes no sense to me as the input is definitely finite.

[Stine Larsen]

Running

any(~isfinite(Hej(x0)))

yields a matrix purely with zeros.

[Matt J]

Insert

which -all Hej(x0)

right before the call to lsqnonlin. It will help make sure that you are calling the version of Hej that you think you are.

Also, implement Hej as an mfile function instead of an anonymous function. Then use DBSTOP again. When and if it stops, switch to the workspace of Hej and examine the values of x0 and Hej(x0).

[Stine Larsen]

I have now implemented the Hej-function as follows:

function          output=Hej(p)
global r
global theta
global a
global Intensity
output=abs((p(1)*besselj(0,p(2)*r)+besselj(1,p(6)*r).*(p(4)*cos(theta+p(8))+p(9)*sin(theta+p(8)))*exp(1i*p(5))).*(r/a<=1)...
                              + (p(1)*besselj(0, p(2)*a)/besselk(0, p(3)*a)*besselk(0,p(3)*r)+...                              
                              (p(4)*besselj(1, p(6)*a)/besselk(1, p(7)*a)*cos(theta+p(8))...
                                +p(9)*besselj(1, p(6)*a)/besselk(1,p(7)*a)*sin(theta+p(8))).*besselk(1, p(7)*r).*exp(1i*p(5)))...
                                .*(r/a>1)).^2- Intensity;

Implementing the lsqnonlin as

which -all Hej(x0)
          [p_estimated,resnorm,residuals,exitflag,output] = lsqnonlin(@Hej,x0, lb,[], opts)

shows that is the correct Hej-function. Examining x0 and Hej(x0) provides the same results as when it was an anomynous function.

I keep getting the same error as previosly.

[Stine Larsen]

In the debugging mode I get the following result for x0 and Hej(x0)

K>> Hej(x0)
Undefined function or variable 'x0'.

K>> x0
Undefined function or variable 'x0'.

[Matt J]

Because in the workspace of Hej, you have "p" and "output", not x0 and Hej(x0).

[Stine Larsen]

Okay. I am not quite sure what you are suggesting.

Any suggestions to recurring mistake, when not implenting the Hej-function as a separate fundtion.

[Matt J]

Put the following lines at the end of Hej.m, after the computation of the output

check = isfinite(output) & isreal(output),
keyboard

This will stop you in debug mode at the end of every call to Hej.m and will also display for you whether "output" has a legitimate value. Step through sucessive calls to Hej.m by using the "return" command in debug mode until the instant that check=0.

[Stine Larsen]

I have included the 2 new lines in Hej. In the first call of Hej, I get a matrix with a single zero. Running the following in debugging mode

k=find(check==0)

k =

        403602

[row, col]=find(check==0)

row =

     408

col =

     448

That is at origo.

I will try moving the vectors slightly so that nothing is to be evaluated at (0,0).

[Stine Larsen]

Moving away from origo there is no error.

Answer 2

[Alan Weiss]

Your settings of TolFun and TolX don't make sense. They should not be smaller than 1e-14, and should probably be much larger. See the documentation on tolerances.

It seems that you have complex numbers in your objective function. Maybe not, maybe the abs makes everything real. And I am also not sure, but it seems that you are summing the squares in your function. You shouldn't do that, you should pass the vector of values to lsqnonlin as documented (see the second sentence in Description). If you have a complex objective function, then you should use the Levenberg-Marquardt algorithm as explained here.

Now, finally, I come to your question. I am not sure, but it seems to me that your initial point is right on the edge of the region where your objective function is finite. Finite differences can step outside this region. So, if I am correct, all you need to do is set start_point(1,9) to a value above 0, maybe 1e-2.

Alan Weiss

MATLAB mathematical toolbox documentation

[Stine Larsen]

Hi Alan,

Thank you for the input!

I have tried with numerous settings for for TolX and TolFun including numbers larger than 1e-14 - all settings results in the same error.

I do not have complex numbers in my objective function, as you point out the abs makes them real. The Intensity matrix consists of purely real and positive numbers. My function will result in real-valued matrix (which for the given inputs may be plotted).

My initialpoint returns finite values (approximately between -1 and 1) for the object funtion. In principle the objective function should yield finite values as long as a>0.

I have tried changing start_point(1,9) - again with the same result.

New suggestions are very welcome.

Stine