Cross-posted from MSE, where this question was asked over a year ago with no answers.
Suppose I have a large system of polynomial equations in a large number of real-valued variables. \begin{align} f_1(x_1, x_2, &\dots, x_m) = 0 \\ f_2(x_1, x_2, &\dots, x_m) = 0 \\ &\vdots \\ f_n(x_1, x_2, &\dots, x_m) = 0 \\ \end{align} (In my particular case, I have about $n \approx 1000$ equations of degree $10$ in about $m \approx 200$ variables.) By numerical means, I've found an approximate solution vector $(\tilde{x}_1, \tilde{x}_2, \dots, \tilde{x}_m)$ at which the value of each $f_j$ is very small: $$\lvert f_j(\tilde{x}_1, \tilde{x}_2, \dots, \tilde{x}_m) \rvert < 10^{-16} \quad \forall j = 1, \dots, n.$$ This leads me to believe that a genuine solution of my system exists somewhere in a small neighborhood of $(\tilde{x}_1, \tilde{x}_2, \dots, \tilde{x}_m)$, and that the small residuals I see are due to round-off error in finite (IEEE double) precision arithmetic. However, it could conceivably be the case that the zero loci of my polynomials $f_j$ come very close to each other (within $10^{-16}$) but do not mutually intersect. How can I rigorously tell which is the case?
I could, of course, further refine my solution using quadruple- or extended-precision arithmetic to push the residuals even closer to zero, but this would only provide supporting empirical evidence.
If it helps, all of my polynomials have integer coefficients and can be evaluated exactly on integer and rational inputs. However, my approximate solution $(\tilde{x}_1, \tilde{x}_2, \dots, \tilde{x}_m)$ is probably irrational.
In principle, there are methods of computational algebraic geometry (Groebner bases, cylindrical decomposition, etc.) that can algorithmically decide the existence of a true mathematical solution to a polynomial system, but my system is completely out of reach of all such algorithms I know. Buchberger's algorithm, for example, has doubly-exponential time complexity in the number of input variables.
Note that interval/ball arithmetic won't help, because even if I can show that each $f_j$ exactly assumes the value $0$ in a small neighborhood of $(\tilde{x}_1, \tilde{x}_2, \dots, \tilde{x}_m)$, it could be the case that a different point zeroes out each $f_j$, and no single point simultaneously zeroes out all of them.