why do I get an error with fsolve in R2012b if in R2007a it works fine ?

I'm facing with different results returned by the algorithm fsolve between the releases R2007a and R2012b. In my case, when I run the code below in Matlab R2007a:

opts = optimset('Display','iter','FunValCheck','on','MaxIter',200,'MaxFunEvals',400,'TolX',1e-6,'TolFun',1e-6);
%
[filename,pathname] = uigetfile('*.dat','Select the data file to test fsolve(.)');
fstr = [pathname filename];
a = dlmread(fstr,' '); % load envelope data file with 1000 i.i.d. samples
%
Z = a;
ER = mean(Z);
ER2 = mean(Z.^2);
ER4 = mean(Z.^4);
x0 = [2.25; 1.00];
%
fl = @(x)ekm6(x,ER,ER2,ER4); % function handle
[x,fval,exitflag,output] = fsolve(fl,x0,opts)

which calls the function ekm6.m:

function f = ekm6(x,ER,ER2,ER4)
%
f = [(1 + 2*x(1))/(x(2)*(1 + x(1))^2) - (ER4/ER2^2 - 1);
    (gamma(x(2) + 1/2)*(hypergeom(-1/2,x(2),-x(1)*x(2))/gamma(x(2))))/sqrt((1 + x(1))*x(2)) - ER/sqrt(ER2)];

I get the result below (the expected result):

                                         Norm of      First-order   Trust-region
 Iteration  Func-count     f(x)          step         optimality    radius
     0          3      0.00312567                        0.0294               1
     1          6    2.70862e-005       0.103684        0.00224               1
     2          9    2.37394e-009      0.0123891      2.06e-005               1
     3         12    1.43461e-017    0.000183517       1.6e-009               1
Optimization terminated: first-order optimality is less than options.TolFun.

x =

       2.2729
       1.1127

...

output = 

       iterations: 3
        funcCount: 12
        algorithm: 'trust-region dogleg'
    firstorderopt: 1.5979e-009
          message: 'Optimization terminated: first-order optimality is less than options.TolFun.'

NOTE: the file loaded by command dlmread() is a sequence of 1000 i.i.d. random variables (attached)

On the other hand, when the same code above is applied in Matlab R2102b, I get the error:

Error using lsqfcnchk/checkfun (line 136)
User function '@(x)ekm6(x,ER,ER2,ER4)' returned a complex value when evaluated;
 FSOLVE cannot continue.

Error in C:\Program Files\MATLAB\R2012b\toolbox\shared\optimlib\finDiffEvalAndChkErr.p>finDiffEvalAndChkErr (line 26)

Error in C:\Program Files\MATLAB\R2012b\toolbox\shared\optimlib\finitedifferences.p>finitedifferences (line 128)

Error in trustnleqn (line 97)
    [JACfindiff,~,~,numFDfevals] = finitedifferences(x,funfcn{3},[],[],[],Fvec,[],[], ...

Error in fsolve (line 403)
    [x,FVAL,JACOB,EXITFLAG,OUTPUT,msgData]=...

Error in test_fsolve (line 17)
[x,fval,exitflag,output] = fsolve(fl,x0,opts)

The optimization toolbox is properly installed in R2012b since I've successfully tested the Example 1 in the fsolve() help page.

Is there any change in the fsolve() code of R2012b that leads to this error ? How can I run this code without error with R2012b ? I would appreciate any help, cause I've not found similar question in the Mathworks support page. Best regards.

Answer 1

[Matt J]

I doubt the key question is the difference between the MATLAB versions. The key question is probably whether ekm6 is indeed generating complex values, and if so, why. Modify as follows (or set a conditional breakpoint) to trap the occurrence of a complex f.

function f = ekm6(x,ER,ER2,ER4)
%
f = [(1 + 2*x(1))/(x(2)*(1 + x(1))^2) - (ER4/ER2^2 - 1);
    (gamma(x(2) + 1/2)*(hypergeom(-1/2,x(2),-x(1)*x(2))/gamma(x(2))))/sqrt((1 + x(1))*x(2)) - ER/sqrt(ER2)];

is ~isreal(f), keyboard; end

[Antonio Ribeiro]

Thanks Matt, for suggesting me to debug ekm6.m. In fact, the function f returns a complex number when executed on R2012b, which does not happen in R2007a. I edited ekm6.m with your suggestion:

function f = ekm6(x,ER,ER2,ER4)
%
f = [(1 + 2*x(1))/(x(2)*(1 + x(1))^2) - (ER4/ER2^2 - 1);
    (gamma(x(2) + 1/2)*(hypergeom(-1/2,x(2),-x(1)*x(2))/gamma(x(2))))/sqrt((1 + x(1))*x(2)) - ER/sqrt(ER2)];
%
TF = isreal(f);
if TF == false
    keyboard;
end
end

I clicked the run button in test_solve.m, after more than 10 min e many system memory used, the keyboard gave me the values for ekm6.m:

K>> x(1)

ans =

      2.2500

K>> x(2)

ans =

      1.0000

K>> f

f =

     0.0554          
    -0.0074 + 0.0000i

K>> hypergeom(-1/2,x(2),-x(1)*x(2))

ans =

     1.8957 +  1.691e-72i

The problem is with hypergeom that returns a complex number when calculated for the iterative variables x(1) e x(2) in fsolve. Interestingly, this does not happen with R2007a.

If I calculate the hypergeom in the work space with the numerical starting values there is no problem

>> hypergeom(-1/2,1,-2.25*1)

ans =

      1.8957

Also I made x symbolic as below but the error remained

opts = optimset('Display','iter','FunValCheck','on','MaxIter',200,'MaxFunEvals',400,'TolX',1e-6,'TolFun',1e-6);
%
[filename,pathname] = uigetfile('*.dat','Select the data file to test fsolve(.)');
fstr = [pathname filename];
a = dlmread(fstr,' '); % load envelope data file with s samples
%
Z = a;
ER = mean(Z);
ER2 = mean(Z.^2);
ER4 = mean(Z.^4);
x0 = [2.25; 1.00];
%
syms x
fl = @(x)ekm6(x,ER,ER2,ER4);
[x,fval,exitflag,output] = fsolve(fl,x0,opts)

I will post a new question specifically about the hypergeom function. Regards.