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

Question

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!

Answer 1

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.8358090074041686763439090398`1.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 ---

Answer 2

ss = NSolve[Thread[e == 0], Variables@e, WorkingPrecision -> 50];
Select[Abs /@ Join @@ (e /. ss ), # > 10^-9 &]
(*
 {}
*)

Answer 3

I wanted to comment a "Thank You" to @dbm for asking and @Daniel Lichtblau for answering, but don't have the necessary reputation yet.

There's also a way to bypass this problem in certain cases that wasn't explicitly stated here, so I figured I'd mention it.

I ran into this same issue when randomly generating momenta (also for a particle physics problem), and so had $k_1 + k_2 + k_3 + k_4 = 0$, where $k_i$ are four-component vectors also satisfying $k_i^2 = 0$. These were the explicit algebraic constraints that were putting me on or off (not sure which) the discriminant variety, per @Daniel Lichtblau's answer:

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.

The tiny numerical error I got from solving the constraint equations $$ k_1 + k_2 + k_3 + k_4 = 0 \quad \textrm{and} \quad k_i^2 = 0 $$ numerically was enough to result in the infinite solutions described by @dbm. The solution in my specific case (and likely in many cases where the algebraic constraints can be handled explicitly), was to impose the constraints analytically instead of numerically. That reduced the numerical error enough that NSolve returned only the non-infinite solutions, all of which satisfied the original equations.