How can I more efficiently define a square root with the branch cut along the positive imaginary axis?

Question

I have a system of equations I need to solve involving square roots. I don't always have a good guess what the solutions are, and so I am fortunate that NSolve can handle these equations. I want to be able to try moving the branch cut of the square root function to the positive imaginary axis since I want continuous functions as I cross the positive and negative real axis. I can easily accomplish this as follows:

Define my new square root function as

MySqrt[z_] := If[Re[z] < 0 && Im[z] >= 0, -Sqrt[z], +Sqrt[z]]

Now for some equations I pull out of a hat:

EQNS = {x + MySqrt[y] - 2/MySqrt[z] + z - 1, 
         y - x + MySqrt[x] + 1/MySqrt[z] - 2, 
           MySqrt[x] + 1/MySqrt[y] - 2 y + z};
NSolve[EQNS == 0, {x, y, z}] // Timing

Yields:

  Out[990]={0.92901, 
    {{x -> 0.698929 + 3.41798 I, y -> 1.03321 + 0.455899 I, z -> -0.301342 - 0.075144 I},
     {x -> 2.72138, y -> 1.42833, z -> 0.37027}}}

Now to check this against the standard definition:

EQNS = {x + Sqrt[y] - 2/Sqrt[z] + z - 1, 
       y - x + Sqrt[x] + 1/Sqrt[z] - 2, 
       Sqrt[x] + 1/Sqrt[y] - 2 y + z};
NSolve[EQNS == 0, {x, y, z}] // Timing

Yields:

Out[988]={0.279246, {
{x -> 0.698929 + 3.41798 I, y -> 1.03321 + 0.455899 I, z -> -0.301342 - 0.075144 I}, 
{x -> 0.698929 - 3.41798 I, y -> 1.03321 - 0.455899 I, z -> -0.301342 + 0.075144 I}, 
{x -> 0.698929 + 3.41798 I, y -> 1.03321 + 0.455899 I, z -> -0.301342 - 0.075144 I}}}

Which brings me to my question: Why does my definition for the square root make it so that it takes NSolve so much longer? It is about a factor of 4 here, which isn't a huge deal, but I am solving even more complicated equations with more variables and I am finding the increase to be a factor of 10 or greater. Can anyone suggest a better way to define a new (faster?) square root function? Or make any other suggestions to make this process faster?

Answer 1

Multiplying a number by $e^{\pi i/2} = i$ is equivalent to rotating the phase by $90^\circ$. So, if you pre-multiply your number by $i$ and take the square root, you will be using the same branch cut as the default. Then, to compensate, you will need to multiply by $e^{-\pi i/4}$ to rotate the phase backwards by $45^\circ$. Here is a function that does this:

isqrt[d_] := -Sqrt[I d](-1)^(3/4)

First, let's compare isqrt with MySqrt:

sample = RandomComplex[{-5 - 5 I, 5 + 5 I}, 100];
MinMax @ Abs[isqrt[sample] - MySqrt /@ sample]

{0., 6.28037*10^-16}

So, basically the same (I check the norm of the differences to account for differences in numerics). Next, let's check the timing:

EQNS={x+MySqrt[y]-2/MySqrt[z]+z-1,y-x+MySqrt[x]+1/MySqrt[z]-2,MySqrt[x]+1/MySqrt[y]-2 y+z};
NSolve[EQNS==0,{x,y,z}] //Timing

EQNS={x+isqrt[y]-2/isqrt[z]+z-1,y-x+isqrt[x]+1/isqrt[z]-2,isqrt[x]+1/isqrt[y]-2 y+z};
NSolve[EQNS==0,{x,y,z}] //Timing

{0.783682, {{x -> 0.698929 + 3.41798 I, y -> 1.03321 + 0.455899 I, z -> -0.301342 - 0.075144 I}, {x -> 2.72138, y -> 1.42833, z -> 0.37027}}}

{0.228611, {{x -> 0.698929 + 3.41798 I, y -> 1.03321 + 0.455899 I, z -> -0.301342 - 0.075144 I}, {x -> 2.72138, y -> 1.42833, z -> 0.37027}}}

The reason that MySqrt is so slow is because MySqrt[y] contains an "unevaluated" If object, while isqrt[y] is a simple arithmetic object:

MySqrt[y]
isqrt[y]

If[Re[y] < 0 && Im[y] >= 0, -Sqrt[y], +Sqrt[y]]

-(-1)^(3/4) Sqrt[I y]

NSolve probably has to use more complicated algebraic methods when presented with these If objects.