NSolve erroneously gives no solution to a polynomial system

Question

I have a polynomial system with three equations in three unknowns, the maximum degree is 26. Two equations are symmetric, i.e. eq1(x,y,z)=eq2(y,x,z). If I search for a symmetric solution with y=x, NSolve returns a number of solutions. However, if I apply NSolve to the original system, {} is returned after a while, without any errors. What could be going on here?

Edit: below are the polynomials in question.

eq1 = -400 z^2 (12500 y^3 z^5 (90 (-1+z) z+y (-130-1967 z+1701 z^2))+x^6 (25+420 z+164 z^2) (100 z^2+20 y z (10+189 z)+y^2 (100+3780 z+3321 z^2))^2+6250 x y^2 z^4 (180 (-1+z) z^2+4 y z (-350-3655 z+3213 z^2)+y^2 (-1300-35800 z-148177 z^2+120285 z^3))+5000 x^2 y^2 z^3 (60 z^2 (-25-200 z+126 z^2)+2 y z (-2375-54375 z-242440 z^2+176112 z^3)+y^2 (-3250-129775 z-1201160 z^2-2076368 z^3+1507086 z^4))+2 x^5 (25+420 z+164 z^2) (-10000 z^4-8000 y z^3 (5+92 z)+300 y^2 z^2 (-275-8545 z-45294 z^2+492 z^3)+20 y^3 z (-4250-188475 z-2030715 z^2-57492 z^3+278964 z^4)+y^4 (-32500-1942000 z-27773325 z^2-13103595 z^3+14755203 z^4+4901796 z^5))+500 x^3 y z^2 (-400 z^3 (10+89 z)+60 y z^2 (-525-10795 z-41594 z^2+11076 z^3)+4 y^2 z (-15000-535375 z-5230875 z^2-10618835 z^3+6404076 z^4)+y^3 (-32500-1699500 z-24450075 z^2-98053445 z^3-5080482 z^4+31443228 z^5))+25 x^4 z (-2000 z^4 (10+89 z)-800 y z^3 (200+4460 z+17477 z^2)+120 y^2 z^2 (-4500-153900 z-1409110 z^2-2323793 z^3+82656 z^4)+8 y^3 z (-90625-4401875 z-66277500 z^2-300577200 z^3-20833965 z^4+47108016 z^5)+y^4 (-325000-20987500 z-408812500 z^2-2528313000 z^3-1459453050 z^4+653500107 z^5+366000768 z^6)))
eq2 = -400 z^2 (10000 y^4 z^4 (-5 z (10+89 z)-2 y (25+420 z+164 z^2)+y^2 (25+420 z+164 z^2))+4000 x y^3 z^3 (-50 z^2 (10+89 z)-5 y z (200+4460 z+17477 z^2)-4 y^2 (125+4400 z+39460 z^2+15088 z^3)+y^3 (250+8925 z+81020 z^2+30996 z^3))+600 x^2 y z^2 (1875 (-1+z) z^4+500 y z^3 (-25-200 z+126 z^2)+50 y^2 z^2 (-525-10795 z-41594 z^2+11076 z^3)+5 y^3 z (-4500-153900 z-1409110 z^2-2323793 z^3+82656 z^4)+y^5 (2500+136500 z+2227025 z^2+11086740 z^3+4087044 z^4)+y^4 (-6875-329125 z-4766350 z^2-20412560 z^3-7221576 z^4+80688 z^5))+40 x^3 z (28125 (-1+z) z^5+625 y z^4 (-350-3655 z+3213 z^2)+250 y^2 z^3 (-2375-54375 z-242440 z^2+176112 z^3)+50 y^3 z^2 (-15000-535375 z-5230875 z^2-10618835 z^3+6404076 z^4)+5 y^4 z (-90625-4401875 z-66277500 z^2-300577200 z^3-20833965 z^4+47108016 z^5)+y^6 (25000+1837500 z+42668750 z^2+338995125 z^3+386232300 z^4+102937716 z^5)+y^5 (-106250-6496875 z-130624375 z^2-885247500 z^3-350209800 z^4+107736192 z^5+45750096 z^6))+x^4 (12500 z^5 (-130-1967 z+1701 z^2)+y^6 (25+420 z+164 z^2) (100+3780 z+3321 z^2)^2+6250 y z^4 (-1300-35800 z-148177 z^2+120285 z^3)+5000 y^2 z^3 (-3250-129775 z-1201160 z^2-2076368 z^3+1507086 z^4)+500 y^3 z^2 (-32500-1699500 z-24450075 z^2-98053445 z^3-5080482 z^4+31443228 z^5)+25 y^4 z (-325000-20987500 z-408812500 z^2-2528313000 z^3-1459453050 z^4+653500107 z^5+366000768 z^6)+2 y^5 (-812500-62200000 z-1515303125 z^2-12310874375 z^3-9689455125 z^4+4170740580 z^5+4478607612 z^6+803894544 z^7)))
eq3 = -400 (12500000 (-1+y) y^6 z^8 (25 z+2 y^2 (105+82 z)+5 y (5+84 z))+5000000 x (-1+y) y^5 z^7 (250 z^2+75 y z (10+147 z)+25 y^2 (20+693 z+4652 z^2)+y^3 (4200+60610 z+44772 z^2))+4 x^9 (105+82 z) (2500 z^4+500 y z^3 (20+273 z)+25 y^2 z^2 (600+16380 z+67481 z^2)+10 y^3 z (1000+40950 z+337405 z^2+201474 z^3)+y^4 (2500+136500 z+1687025 z^2+2014740 z^3+544644 z^4))^2+62500 x^2 y^4 z^5 (625 z^4 (-74+9 z)+500 y z^3 (-470-5607 z+378 z^2)+50 y^2 z^2 (-7550-193325 z-1082282 z^2+33228 z^3)+8 y^5 z (147000+4127900 z+25870635 z^2+17756198 z^3)+20 y^3 z (-10250-434275 z-5137625 z^2-16185017 z^3+61992 z^4)+y^4 (-16250-1576375 z-30283800 z^2-106933140 z^3+234123840 z^4+242064 z^5))+12500 x^3 y^3 z^4 (6250 z^5 (-50+9 z)+625 y z^4 (-3610-40630 z+6237 z^2)+500 y^2 z^3 (-11800-278675 z-1440017 z^2+168174 z^3)+50 y^3 z^2 (-135750-5012550 z-54425905 z^2-161386282 z^3+12110364 z^4)+8 y^6 z (1470000+61344500 z+729108450 z^2+2509338345 z^3+1543671402 z^4)+10 y^4 z (-323750-16777875 z-284394650 z^2-1725057270 z^3-2682939630 z^4+44504352 z^5)+y^5 (-406250-30198750 z-718189875 z^2-5601198800 z^3-5235243980 z^4+21387326112 z^5+86416848 z^6))+3125 x^4 y^2 z^3 (12500 z^6 (-74+9 z)+2500 y z^5 (-3610-40630 z+6237 z^2)+125 y^2 z^4 (-269600-6271800 z-31535234 z^2+4879809 z^3)+100 y^3 z^3 (-610500-22016600 z-231267240 z^2-676924767 z^3+95566338 z^4)+10 y^4 z^2 (-5582500-277080750 z-4514172950 z^2-27429959045 z^3-47636823957 z^4+5572844748 z^5)+4 y^7 z (14700000+814100000 z+14268765000 z^2+89719815000 z^3+175946555925 z^4+89184527498 z^5)+10 y^5 z (-2362500-151852500 z-3404834250 z^2-31477458950 z^3-107192270485 z^4-70258903867 z^5+4014973872 z^6)+y^6 (-3250000-266250000 z-7734616250 z^2-91274467875 z^3-353314512600 z^4+81822358510 z^5+682796135744 z^6+7739512272 z^7))+500 x^5 y z^2 (-2500000 z^7+62500 y z^6 (-470-5607 z+378 z^2)+12500 y^2 z^5 (-11800-278675 z-1440017 z^2+168174 z^3)+625 y^3 z^4 (-610500-22016600 z-231267240 z^2-676924767 z^3+95566338 z^4)+500 y^4 z^3 (-1077500-53309750 z-859760550 z^2-5217504375 z^3-9155373012 z^4+1373964516 z^5)+50 y^5 z^2 (-8225000-525798750 z-11613295500 z^2-108006433425 z^3-394956201355 z^4-373424805542 z^5+57710474136 z^6)+4 y^8 z (73500000+5073775000 z+117600525000 z^2+1077107192500 z^3+3626828197125 z^4+3688382662340 z^5+1139229494676 z^6)+10 y^6 z (-15375000-1227393750 z-34936781250 z^2-440165818750 z^3-2389380496125 z^4-4259671880765 z^5-543757517062 z^6+199369046880 z^7)+y^7 (-20312500-2011500000 z-71374225000 z^2-1119790209375 z^3-7315791538750 z^4-12681549935250 z^5+12305201074940 z^6+11922339291856 z^7+378251627040 z^8))+50 x^6 z (-6250000 z^8-2500000 y z^7 (30+431 z)+62500 y^2 z^6 (-7550-193325 z-1082282 z^2+33228 z^3)+12500 y^3 z^5 (-135750-5012550 z-54425905 z^2-161386282 z^3+12110364 z^4)+625 y^4 z^4 (-5582500-277080750 z-4514172950 z^2-27429959045 z^3-47636823957 z^4+5572844748 z^5)+500 y^5 z^3 (-8225000-525798750 z-11613295500 z^2-108006433425 z^3-394956201355 z^4-373424805542 z^5+57710474136 z^6)+50 y^6 z^2 (-53875000-4317337500 z-122350656250 z^2-1545553044375 z^3-8658001464750 z^4-17447002625280 z^5-4281527765016 z^6+1116768103056 z^7)+4 y^9 z (367500000+30385250000 z+876127612500 z^2+10549119450000 z^3+50543244451875 z^4+77207282058450 z^5+46560234887460 z^6+9670806703512 z^7)+50 y^7 z (-17562500-1748150000 z-61641381250 z^2-992616036250 z^3-7418598629625 z^4-21656294896665 z^5-10779809973576 z^6+837951801292 z^7+633852330048 z^8)+y^8 (-101562500-12916093750 z-561037093750 z^2-11032005028125 z^3-98895162378125 z^4-320995927116250 z^5+11541934933200 z^6+321942723499160 z^7+136511538445872 z^8+5502013053120 z^9))+x^8 (12500000 z^8 (-185+256 z)+5000000 y z^7 (-3700-43285 z+71528 z^2)+4 y^9 (105+82 z) (2500+136500 z+1687025 z^2+2014740 z^3+544644 z^4)^2+62500 y^2 z^5 (-16250-1576375 z-30283800 z^2-106933140 z^3+234123840 z^4+242064 z^5)+12500 y^3 z^4 (-406250-30198750 z-718189875 z^2-5601198800 z^3-5235243980 z^4+21387326112 z^5+86416848 z^6)+y^8 (25+420 z+164 z^2)^2 (-29225000-1286281250 z-16042309375 z^2+8122707000 z^3+54598605300 z^4+30273492096 z^5+803894544 z^6)+3125 y^4 z^3 (-3250000-266250000 z-7734616250 z^2-91274467875 z^3-353314512600 z^4+81822358510 z^5+682796135744 z^6+7739512272 z^7)+500 y^5 z^2 (-20312500-2011500000 z-71374225000 z^2-1119790209375 z^3-7315791538750 z^4-12681549935250 z^5+12305201074940 z^6+11922339291856 z^7+378251627040 z^8)+50 y^6 z (-101562500-12916093750 z-561037093750 z^2-11032005028125 z^3-98895162378125 z^4-320995927116250 z^5+11541934933200 z^6+321942723499160 z^7+136511538445872 z^8+5502013053120 z^9)+10 y^7 (-101562500-19520312500 z-1065143281250 z^2-25668581968750 z^3-289577442890625 z^4-1286535678011875 z^5-471989907683000 z^6+1653988149325800 z^7+1479996563635440 z^8+361065214964880 z^9+13535440402176 z^10))+10 x^7 (-6250000 z^8 (5+79 z)-12500000 y z^7 (20+663 z+4211 z^2)+125000 y^2 z^6 (-10250-434275 z-5137625 z^2-16185017 z^3+61992 z^4)+12500 y^3 z^5 (-323750-16777875 z-284394650 z^2-1725057270 z^3-2682939630 z^4+44504352 z^5)+3125 y^4 z^4 (-2362500-151852500 z-3404834250 z^2-31477458950 z^3-107192270485 z^4-70258903867 z^5+4014973872 z^6)+500 y^5 z^3 (-15375000-1227393750 z-34936781250 z^2-440165818750 z^3-2389380496125 z^4-4259671880765 z^5-543757517062 z^6+199369046880 z^7)+250 y^6 z^2 (-17562500-1748150000 z-61641381250 z^2-992616036250 z^3-7418598629625 z^4-21656294896665 z^5-10779809973576 z^6+837951801292 z^7+633852330048 z^8)+8 y^9 z (262500000+25286875000 z+872210062500 z^2+13019934956250 z^3+81022758665625 z^4+165148040303000 z^5+147358962336600 z^6+59875852069440 z^7+8997991630992 z^8)+20 y^7 z (-60156250-7708359375 z-334196875000 z^2-6638447356250 z^3-62930398093750 z^4-245935395612250 z^5-188063456267100 z^6+5474243584200 z^7+28342274761872 z^8+4103157149568 z^9)+y^8 (-101562500-19520312500 z-1065143281250 z^2-25668581968750 z^3-289577442890625 z^4-1286535678011875 z^5-471989907683000 z^6+1653988149325800 z^7+1479996563635440 z^8+361065214964880 z^9+13535440402176 z^10)))

Edit 2: As is noted in the comments, the system above is huge, and it takes too long to try different solution methods. Interestingly, while working on the initial problem, I came across a much simpler system, but with the same NSolve issue. Here are the equations:

foc = {-x^4 y^2+2 x^3 y^3+x^2 y^4-2 x^4 y z+4 x^3 y^2 z-x^4 y^2 z+10 x^2 y^3 z-2 x^3 y^3 z+4 x y^4 z-x^2 y^4 z-2 x^4 z^2+8 x^2 y^2 z^2+8 x y^3 z^2+2 y^4 z^2,x^4 y^4+2 x^3 y^5-x^2 y^6+8 x^4 y^3 z+18 x^3 y^4 z-2 x^4 y^4 z+4 x^2 y^5 z-4 x^3 y^5 z-2 x y^6 z-2 x^2 y^6 z+16 x^4 y^2 z^2+44 x^3 y^3 z^2-4 x^4 y^3 z^2+28 x^2 y^4 z^2-10 x^3 y^4 z^2-8 x^2 y^5 z^2-2 y^6 z^2-2 x y^6 z^2+12 x^4 y z^3+40 x^3 y^2 z^3-2 x^4 y^2 z^3+38 x^2 y^3 z^3-6 x^3 y^3 z^3+8 x y^4 z^3-7 x^2 y^4 z^3-2 y^5 z^3-4 x y^5 z^3-y^6 z^3+3 x^4 z^4+12 x^3 y z^4+14 x^2 y^2 z^4+4 x y^3 z^4-y^4 z^4,y (2 x^2 y^3+8 x^2 y^2 z+4 x y^3 z+10 x^2 y z^2+10 x y^2 z^2-2 x^2 y^2 z^2+y^3 z^2-x y^3 z^2+4 x^2 z^3+6 x y z^3-4 x^2 y z^3+2 y^2 z^3-4 x y^2 z^3-2 x^2 z^4-3 x y z^4-y^2 z^4)}

Please note that these are not symmetric anymore. NSolve[foc,Reals] gives {} without any errors. So does NSolve[Join[foc, {x >= 1, y >= 1, z >= 1}], Reals] (the constraints give the domain I am interested in). Running GroebnerBasis before NSolve still produces no results but the following error message appears from NSolve:

Infinite solution set has dimension at least 1.
Returning intersection of solutions with
(171802 x)/178835-(113492 y)/178835-(121484 z)/178835 == 1.

The best part of it, the solution exists as can be easily seen from the following contour plot (and by playing with the upper limits it seems the solution is unique):

ContourPlot3D[Evaluate@Thread@(foc == 0), {x, 1, 10}, {y, 1, 10}, {z, 1, 10}]

Solve, given the constraints on the domain, can find that unique solution. Somehow, NSolve fails in this case...

Edit 3: belisarius has pointed out that foc/.{x->0,y->0} = {0,0,0}. It is also the case that foc/.{x->0,z->0} = {0,0,0} and foc/.{y->0,z->0} = {0,0,0}. I cannot trivially delete these roots from the original equations, because the original equations do not factor so. However, these roots do seem to present the problem when using NSolve after GroebnerBasis as is evident from the error message. So, I've tried dropping all the small common factors from the Groebner basis.

gb = GroebnerBasis[foc,{x,y,z}];
NSolve[FactorList[#][[-1,1]]&/@gb,Reals]

gives the set of the remaining solutions, with the one I was looking for being among them.

So, it seems that NSolve gets caught in the positive dimensionality part of the solution, even if explicitly told the region (x>1,y>1,z>1), where that positive dimensionality is not present. When run after GroebnerBasis, NSolve at least admits it, but when run on the original equations, it seems to hide the fact.

(I accept belisarius' answer as the comment below helped me to at least have an idea of what might be going wrong.)

Edit 4: Just for fun:

addConstraints[foc_] := {Sequence @@ foc, x >= 1, y >= 1, z >= 1}
Reduce[addConstraints@Thread[foc == 0], {x, y, z}, Backsubstitution -> True] // N // Timing
{6.87911, x == 6.18083 && y == 4.58414 && z == 3.14352}

gb := GroebnerBasis[foc, {x, y, z}];
NSolve[addConstraints[FactorList[#][[-1, 1]] & /@ gb], Reals] // Timing
{2.97643, {{x -> 6.18083, y -> 4.58414, z -> 3.14352}}}

Solve[addConstraints@Thread[foc == 0]] // N // Timing
{6.84882, {{x -> 6.18083, y -> 4.58414, z -> 3.14352}}}

Answer 1

One way to deal with the singularities (that lead to dimensional components of the solution space) when x, y and/or z are zero is through a change of variables that lets excessive zeros be factored out.

changevar = {y -> x s, z -> x s t};
Short /@ Factor[foc /. changevar]

We can use FactorList to get at the irreducible factors and delete constants and simple variables. We find six real solutions, only one of which has all variables positive. (It was not clear to me whether the lower bound of 1 on the solutions was imposed to try to get rid of the zeros or was part of the system.)

nfoc = Flatten@ DeleteCases[
    FactorList[#][[All, 1]] & /@ (foc /. changevar),
    _?NumericQ | s | x | t,
    {2}];
nsol = NSolve[nfoc == {0, 0, 0}, {x, s, t}, Reals, WorkingPrecision -> 25];
sol = Thread[{x, y, z} -> #] & /@ ({x, y, z} /. changevar /. nsol)
(*
  {{x ->  20.6512,   y ->  1.14137,  z -> -0.452436},
   {x -> -18.1434,   y -> 28.1209,   z -> -6.50644},
   {x ->   6.18083,  y ->  4.58414,  z ->  3.14352},
   {x ->   1.82148,  y -> -4.7767,   z -> -1.25265},
   {x ->   0.885016, y ->  1.62512,  z -> -0.298749},
   {x ->  -1.20071,  y -> -0.747601, z ->  1.57008}}
*)

Verify:

foc /. sol
(* 
  {{0.*10^-18, 0.*10^-18, 0.*10^-20}, {0.*10^-14, 0.*10^-11, 0.*10^-15},
   {0.*10^-18, 0.*10^-16, 0.*10^-18}, {0.*10^-20, 0.*10^-18, 0.*10^-20},
   {0.*10^-23, 0.*10^-22, 0.*10^-23}, {-0.*10^-22, 0.*10^-21, 0.*10^-22}}
*)

Answer 2

Too much fuss :)

rf1 = Reduce[And @@ Thread[foc == 0] && x > 0 && y > 0 && z > 0, {x, y, z}, 
             Quartics -> True, Backsubstitution -> True];

Show[ContourPlot3D[Evaluate@Thread[foc == 0], {x, 1, 7}, {y, 1, 6}, {z, 1, 6},
                  Mesh -> None], 
     Graphics3D[{Green, Sphere[{x, y, z}, .3] /. (ToRules@rf1 // N)}]]

rf1 // N
(* x == 6.18083 && y == 4.58414 && z == 3.14352 *)