FSOLVE seems to find complex roots correctly (although it shouldn't)

I am using an oldish version of MatLab, r2009a, where the manual entry on FSOLVE states that, "when x has complex variables, the variables must be split into real and imaginary parts".

However, FSOLVE seems to handle complex solutions quite well. Here's a simple example:

 opt=optimset('Maxiter',200,'TolX',1e-7,'Tolfun',1e-7);
 fsolve(@(x)1i*x^2-2i*x+5i,0.5i,opt)

 Optimization terminated: first-order optimality is less than options.TolFun.

 ans =

   1.0000 + 2.0000i

I've also checked the actual problem I'm working with -- an eigenvalue problem for a second-order linear differential operator -- and FSOLVE works for it perfectly well too (I've solved the problem twice: by applying FSOLVE to the complex problem "as is" and by splitting it into real and imaginary parts).

So, my question is: in which cases FSOLVE canNOT handle complex roots?

Answer 1

[Alan Weiss]

The description of when fsolve works with complex functions is in Complex Numbers in Optimization Toolbox Solvers. In short, the fsolve algorithms have not changed, but MathWorks has only relatively recently qualified its algorithms to work with complex numbers under the condition that the objective functions are analytic, and the algorithms often work better using a complex initial point.

Alan Weiss

MATLAB mathematical toolbox documentation

[Eugene Benilov]

Many thanks, Alan.

[John D'Errico]

A complex start point is not always sufficient, as I show in my answer. The idea comes down to a basin of attraction. Good starting values are often the answer to getting the right result from an optimization.

Answer 2

[John D'Errico]

Fsolve CAN work, under some circumstances.

Consider the function x^3 + 2*x^2 + 3*x + 4.

roots([1 2 3 4])
ans =
    -1.6506 +          0i
   -0.17469 +     1.5469i
   -0.17469 -     1.5469i

A real root, and a complex conjugate pair.

roots([1 1 1 1])
ans =
           -1 +          0i
   1.5266e-16 +          1i
   1.5266e-16 -          1i
>> roots([1 2 3 4])
ans =
      -1.6506 +          0i
     -0.17469 +     1.5469i
     -0.17469 -     1.5469i

Create a function, throw it into fsolve.

fun = @(x) x.^3 + 2*x.^2 + 3*x + 4
fun =
  function_handle with value:
    @(x)x.^3+2*x.^2+3*x+4

fsolve(fun,0)

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the default value of the function tolerance, and
the problem appears regular as measured by the gradient.

<stopping criteria details>

ans =
      -1.6506

In fact, it will matter not where on the real line I start fsolve out, it will converge to the unique real root.

ezplot(fun,[-2,3])
grid on

So in a sense, fsolve will ALWAYS fail to find any complex roots on this problem. One reason for that is fsolve will never be pushed off the real line.

Sometimes, if the function returns a complex result, fsolve will start looking in that direction.

fun2 = @(x) x.^3 + 2*x.^2 + 3*x + 4 + eps*i;
fun2(.1)
ans =
        4.321 + 2.2204e-16i

fsolve(fun,0)

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the default value of the function tolerance, and
the problem appears regular as measured by the gradient.

<stopping criteria details>

ans =
      -1.6506

Not here though, even though every function evaluation will always return a result in the complex plane.

What if we start out by looking in the complex plane?

fsolve(fun,.1 + .1*i)

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the default value of the function tolerance, and
the problem appears regular as measured by the gradient.

<stopping criteria details>

ans =
      -1.6506 + 1.6812e-14i

Nope. So even if the initial value lies in the complex plane, fsolve need not always find the complex root. So what am I doing "wrong" here? Can I get fsolve to find a complex root here?

Well, yes. Really, this can come down to the basins of attraction of this function.

fsolve(fun,.1 + i)

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the default value of the function tolerance, and
the problem appears regular as measured by the gradient.

<stopping criteria details>

ans =
     -0.17469 +     1.5469i

Apparently, if I start fsolve out sufficiently near a complex root, it will converge to the desired complex solution.

So what is a basin of attraction?

It is the locus of all points that if you start fsolve out from any of those points, it will converge to a given root of the function, to within some slop based on the convergence tolerance. A basin of attraction need not be a convex set, in fact, they generally are not terribly simple regions. Basins of attraction need not be closed sets, since they need not contain the boundary of the set. Some basins of attraction may converge to a point of divergence, so a point at infinity. And some basins may converge to a stable point that is not a root. Or there might be nasty divergences that are simply cycles, so no convergence ever occurs. Finally, a basin of attraction is algorithm dependent. Different algorithms may have very different basins.

I could spend a fair amount of time on this, carefully mapping the basins of attraction for the three roots of fun. Sigh. Well, it might be "fun" ...

fun = @(x) x.^3 + 2*x.^2 + 3*x + 4;
R = roots([1 2 3 4]);

[r,c] = meshgrid(-10:.1:10);

basinmap = NaN(size(r));

opts = optimset('fsolve');
opts.Display = 'none';

H = waitbar(0);
for ii = 1:numel(r)
  waitbar(ii/numel(r),H)
  res = fsolve(fun,r(ii) + i*c(ii),opts); 
  [~,basinmap(ii)] = min(abs(R - res));
end
delete(H)

H = pcolor(basinmap);
set(H,'edgecolor','none')

You should see that the three basins are not simple regions, with hairy, fuzzy edges delineating them. Worse, that are not even contiguous sets, with isolated spots in one region that lie in the basin of another root. Do you see those yellow and cyan spots inside the dark blue region?

The fundamental issue is that of starting values here. Good starting values will often make or break an optimizer.

[Eugene Benilov]

Thanks, John, for a detailed response. I agree that, since FSOLVE finds only one root at a time, a lot depends on the initial guess.

I was more interested, however, in whether or not FSOLVE is supposed to find complex roots in principle.

[John D'Errico]

Yes, it is "supposed" to find them. Perhaps I would say it better as it is ABLE to find them.

All of the computations in fsolve support complex arithmetic. And linear algebra works on complex numbers. And since fsolve uses linear algebra, the conclusion is yes. Fsolve is ABLE to give complex results.

The entire point of my response is if you want it to find a complex solution, then you need to start the search in the basin of attraction of the desired solution.

Said differently, on many functions, the real line is almost ALWAYS in the basin of attraction of purely real solutions. So if you start with a real start point, don't ever expect it to find a complex root. It might happen if your function returns complex results because complex numbers tend to beget more complex numbers, but that is not the expectation. I showed the example of fun2 in my answer, which always returns a slightly complex result, but the solution did not change.

Or, look at the plot I generated! The x and y axes are the real and imaginary parts of the start point. ANYTHING in the dark blue region converged to the real solution. So even starting with a complex start point yielded the real root, IF you started from that location in the complex plane.