# NSolve gives additional solutions that don't satisfy the equations!

I am trying to solve the following polynomial equations in Mathematica:

e= {-1 + c[1]^2 + s[1]^2, -1 + c[2]^2 + s[2]^2, -1 + c[3]^2 + s[3]^2, -1 +
(-0.70873 c[1] + c[2] - 0.70548 s[1]) y1[ 1], -1 + (-0.916596 c[2] + c[3] +
0.399814 s[2]) y1[ 2], -1 + (c[1] + 0.808085 c[3] + 0.589066 s[3]) y1[ 3],
-1 + (c[1] + 0.808085 c[3] + 0.589066 s[3]) y2[ 1], (-0.70548 c[1] +
0.70873 s[1] + s[2]) y1[ 1] - (-0.589066 c[3] + s[1] + 0.808085 s[3])
y2[1], -1 + (-0.70873 c[1] + c[2] - 0.70548 s[1]) y2[ 2], (0.399814 c[2] +
0.916596 s[2] + s[3]) y1[ 2] - (-0.70548 c[1] + 0.70873 s[1] + s[2]) y2[ 2],
-1 + (-0.916596 c[2] + c[3] + 0.399814 s[2]) y2[ 3], (0.589066 c[3] - s[1] -
0.808085 s[3]) y1[ 3] - (0.399814 c[2] + 0.916596 s[2] + s[3]) y2[3]}


using NSolve command. i.e.,

sol = NSolve[e==0, Variables[e]]


I am only interested in real solutions, so in principle I could also use

NSolve[e==0, Variables[e], Reals].


The command indeed solves this system and gives me 8 solutions, all real. However, when I insert them back to the equations to check if the solutions satisfy the equations, they don't!

Do[ Print[eatsol[i] = e /. Thread[var -> sol[[i]]]], {i, lrealsol}]


gives

{-1.,-1.,-1.,-1.,-1.,-1.,-1.,0.,-1.,0.,-1.,0.}
{-1.,-1.,-1.,-1.,-1.,-1.,-1.,0.,-1.,0.,-1.,0.}
{-1.,-1.,-1.,-1.,-1.,-1.,-1.,0.,-1.,0.,-1.,0.}
{-1.,-1.,-1.,-1.,-1.,-1.,-1.,0.,-1.,0.,-1.,0.}
{-1.,-1.,-1.,-1.,-1.,-1.,-1.,0.,-1.,0.,-1.,0.}
{-1.,-1.,-1.,-1.,-1.,-1.,-1.,0.,-1.,0.,-1.,0.}
{0.,6.93889*10^-17,7.97973*10^-17,0.,-1.11022*10^-16,-1.11022*10^-16,-1.11022*10^-16,-4.47201*10^-17,0.,1.4905*10^-17,-1.11022*10^-16,-1.82678*10^-18}
{0.,6.93889*10^-17,7.97973*10^-17,0.,-1.11022*10^-16,-1.11022*10^-16,-1.11022*10^-16,-4.47201*10^-17,0.,1.4905*10^-17,-1.11022*10^-16,-1.82678*10^-18}


Thus, the first 6 solutions don't even remotely satisfy the equations! I would be fine with some numerical errors etc. but it seems that the infinite solutions also appear here.

Do you know what is going on here? Is there any problem with NSolve?

Update: Please run the below code in Mathematica:

Clear["Global*"];
e = {-1 + c[1]^2 + s[1]^2, -1 + c[2]^2 + s[2]^2, -1 + c[3]^2 +
s[3]^2, -1 + (-0.5646917067485537 c[1] + c[2] -
0.8253019304045068 s[1]) y1[1], -1 + (0.5033628291371058 c[2] + c[3] -
0.8640751484929357 s[2]) y1[2], -1 + (c[1] - 0.47880555467239877 c[3] +
0.877920976406679 s[3]) y1[3], (-0.8253019304045068 c[1] +
0.5646917067485537 s[1] + s[2]) y1[1] - (-0.877920976406679 c[3] + s[1] -
0.47880555467239877 s[3]) y2[1], -1 + (c[1] - 0.47880555467239877 c[3] +
0.877920976406679 s[3]) y2[1], -1 + (-0.5646917067485537 c[1] + c[2] -
0.8253019304045068 s[1]) y2[2], (-0.8640751484929357 c[2] -
0.5033628291371058 s[2] + s[3]) y1[2] - (-0.8253019304045068 c[1] +
0.5646917067485537 s[1] + s[2]) y2[2], -1 + (0.5033628291371058 c[2] + c[3] -
0.8640751484929357 s[2]) y2[3], (0.877920976406679 c[3] - s[1] +
0.47880555467239877 s[3]) y1[3] - (-0.8640751484929357 c[2] -
0.5033628291371058 s[2] + s[3]) y2[3]};
var = Variables[e];
sol = NSolve[e == 0, var];
lsol = Length[sol];
Do[
Print[Chop[e /. sol[[i]]]],
{i, lsol}]

{-1.,-1.,-1.,-1.,-1.,-1.,0,-1.,-1.,0,-1.,0}
{-1.,-1.,-1.,-1.,-1.,-1.,0,-1.,-1.,0,-1.,0}
{-1.,-1.,-1.,-1.,-1.,-1.,0,-1.,-1.,0,-1.,0}
{-1.,-1.,-1.,-1.,-1.,-1.,0,-1.,-1.,0,-1.,0}
{-1.,-1.,-1.,-1.,-1.,-1.,0,-1.,-1.,0,-1.,0}
{-1.,-1.,-1.,-1.,-1.,-1.,0,-1.,-1.,0,-1.,0}
{0,0,0,0,0,0,0,0,0,0,0,0}
{0,0,0,0,0,0,0,0,0,0,0,0}


Basically, the problem is that the solutions given by the NSolve command are not actually the solutions of the equations!

-
Try e /. sol // Chop. –  b.gatessucks Oct 31 '12 at 20:58
Welcome to Mathematica.SE! I suggest the following: 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2)Read the FAQs! 3) When you see good Q&A, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. ALSO, remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign –  chris Oct 31 '12 at 20:59
or if you want to be super sure, 'ee=e// Rationalize[#, 0.00000001] &; sol = Solve[Thread[ee == 0], Variables[ee]]; ee /. sol // N[#, 20] & // Chop' –  chris Oct 31 '12 at 21:05
b.gatessucks and chris, thanks for your comments. My problem is different actually. The problem is that some of the equations at the first 6 solutions is -1, and not zero (or some small number). So Chop etc. is not going to help.d –  dbm Oct 31 '12 at 21:12
Folks, please stop downvoting this. it has actual merit and at this point is more or less clear. –  Daniel Lichtblau Oct 31 '12 at 23:36

The issue is numeric instability. If you rationalize the equations and solve at higher precision you will get small residuals.

re = Rationalize[e, 0];
var = Variables[e];
sol50 = NSolve[re == 0, var, WorkingPrecision -> 50];
lsol = Length[sol];

Max[Abs[re /. sol50]]

(* Out[208]= 4.83580900740416867634390903981.0824091485306464*^-49 *)


If you check sizes of variables you will notice many orders of magnitude separate them in some solutions. This is an indication that high precision may be needed in order to get significant digits for the smaller ones.

--- edit ---

Okay, I missed this subtlety of "infinite" solutions. They are indeed infinite, in an important sense (and this gives the numeric instability). In a certain way the system is not "generic". Specifically, coefficients satisfy one or more (unknown to me) relevant algebraic conditions that put the system on what is called the "discriminant variety" for the family of systems having the same Newton polytope. Remark: Since the system is generic for your purposes, that means your entire family of systems lives on this variety.

In such a situation the number of solutions can drop from the generic (in this case) 8. But that does not quite happen for the example above, because we have coefficients given only approximately, via machine doubles. So we actually have a system NEAR but not ON the discriminant variety.

So what happens to the 8 solutions when we are on this variety? In this case 6 have, loosely speaking, wandered off to infinity. Continuity of solutions as a function of coefficients implies that systems near the variety will have 6 "large" solutions. That's what NSolve is showing. When I round and solve to high precision I get a perfectly good solution set, albeit residuals of the "bad" solutions will be orders of magnitude larger than those of the good ones.

Homotopy continuation methods will not show behavior this for either of two reasons.

(1) They detect "diverging" solutions via heuristics in the tracking phase. In this case 6 would be dropped once they got too large in the tracking.

(2) In some cases one already has a correct result for a related system from the same part of the discriminant variety. In this case a common tactic is to use what is called the "Cheater's Homotopy", using these old solutions to go to solutions for the new system. Here we know in advance that we have the right number of starting points and the tracking will work just fine (barring other numeric issues).

NSolve has a system option that can be used, in a heuristic manner, to elicit better behavior. Or worse, depending on mood. The idea is implementation dependent: we compute a numerical Groebner basis, and a nondefault setting will tell that code to set to zero any coefficients that are below some threshold smaller than the average size in the rest of the polynomial. I show how to use this below, with 10^(-6) for that tolerance. Remark: if I go to 10^(-4) it will hang. That's one way the heuristic can mess up. Another is that it can give garbage results. It is very much a situation where one needs to make sure the results (when they appear) seem plausible. Even finding a viable tolerance is a trial-and-error process, informed perhaps by the nature of the problem. Here I know we have approximated coefficients (only) to several digits, so I expect that 10^(-8) or so will work for my purpose. And it did, as did 10^(-6).

SetSystemOptions["NSolveOptions" -> {"Tolerance" -> 10^(-6)}];
sol = NSolve[e == 0, var]

(* {{c[1] -> -0.840799, c[2] -> 0.0280179, c[3] -> 0.877839,
s[1] -> 0.541347, s[2] -> -0.999607, s[3] -> -0.478956,
y1[1] -> 17.8457, y1[2] -> 0.56958, y1[3] -> -0.594672,
y2[1] -> -0.594672, y2[2] -> 17.8457,
y2[3] -> 0.56958}, {c[1] -> 0.840799, c[2] -> -0.0280179,
c[3] -> -0.877839, s[1] -> -0.541347, s[2] -> 0.999607,
s[3] -> 0.478956, y1[1] -> -17.8457, y1[2] -> -0.56958,
y1[3] -> 0.594672, y2[1] -> 0.594672, y2[2] -> -17.8457,
y2[3] -> -0.56958}} *)


Here is the example where I first ran into this near-discriminant-variety computational mishap four years ago.

http://homepages.math.uic.edu/~jan/

It took some head-scratching between myself and the web page author to figure out what had gone wrong.

--- end edit ---

-
That worked! However, if I still think there is more fundamental problem here: NSolve gives solutions at infinity. If you see the actual solutions in sol, some variables take values of the order 10^16. These are solutions at infinity for a reasonable numerical precision say 10^13 or so. So there is a catch22 here. If we choose WorkingPrecision->50, the solutions at infinity are considered as affine solutions, if we choose the default WorkingPrecision then these solutions still appear but don't satisfy the equations! May be I am explaining it not properly?! –  dbm Oct 31 '12 at 23:04
I'll edit my response tomorrow to address that issue of "infinite" solutions. Meanwhile you might mention where the system came from or the nature of the problem. That might help when I add an explanation. –  Daniel Lichtblau Oct 31 '12 at 23:39
Thanks for your comment and your request to others to stop down-voting. The problem comes from particle theory where I have a nonlinear parametric system of equations and I want to find all real solutions of it. In particular, I want to see if the no. of real solutions is constant for the parametric system. The present system is at a particular parameter-point (a generic point). I started having doubt about the infinite solutions was from my cross-checking the solutions with homotopy continuation method which gave only two solutions always instead of 8 by NSolve. Thanks again. –  dbm Oct 31 '12 at 23:57
Edited to add explanation. –  Daniel Lichtblau Nov 1 '12 at 14:19
@Licthblau, Excellent! Thanks a lot for your detailed explanation. In short, we need to be extra careful while using NSolve for systems with approximate coefficients. I already knew discriminant variety and homotopy continuation method. But I didn't know about the near discriminant variety issue, and specially that it could affect solving systems so drastically. Thanks a lot again! dbm368 –  dbm Nov 1 '12 at 16:40
ss = NSolve[Thread[e == 0], Variables@e, WorkingPrecision -> 50];
Select[Abs /@ Join @@ (e /. ss ), # > 10^-9 &]
(*
{}
*)
`
-
Please see the update in my original massage. –  dbm Oct 31 '12 at 22:38
@dbm Please see my answer –  belisarius Oct 31 '12 at 22:55