How can I distinguish a genuine solution of polynomial equations from a numerical near miss? 
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.
 A: Interval/ball arithmetic may help, actually.
It can be used to prove existence of solutions to multivariate systems like this one. The main idea is: reformulate your system as a fixed-point system $x = g(x)$, $g: \mathbb{R}^n \to \mathbb{R}^n$. If you manage to find a certain product of intervals (hypercube) $X = [a_1,b_1] \times [a_2,b_2] \times \dots \times [a_n, b_n]$ such that $g(X) \subseteq X$ (which you can check using interval arithmetic: interval arithmetic computes an enclosure $g(X) \subseteq \tilde{g}(X)$, so if you have $\tilde{g}(X) \subseteq X$ you have won), then by the Brouwer fixed-point theorem your system has an exact solution inside $X$.
There are tricks to get a 'good' map $g$ such that the inclusion can be verified; for instance, centering on an approximate computed solution and multiplying by an approximated inverse of the Jacobian. Look for the interval Newton method, or the Krawczyk method. Rump's paper on Acta Numerica is a good introduction to these topics.
(EDIT: As the comments below mention, OP's problem is more general than this one though, since the system is overdetermined.)
Interval arithmetic does not mean "replace double x with interval<double> x and run your algorithm again"; this is a common misconception. Most of the times it works best as an a posteriori verification method, starting from an approximate solution that has been computed with a completely different (numerical, unverified) method.
There are tricky and ill-conditioned cases where many digits of precision will be needed in your arithmetic (think $(x-1)^2+\varepsilon=0$), but the main idea works well, especially for simple solutions. (Rump's review paper contains methods to find double solutions as well, essentially augmenting the system with its Jacobian.)
With an overdetermined system like the one OP asks about, one could convert multiple equations to a single one with $f(x)=g(x)=0 \iff f(x)^2+g(x)^2=0$ (since the question is about real solutions), but this gives a problem with multiple solutions, so it's not simply a matter of applying that fixed-point method.
A: There is a technique called Smale's $\alpha$ theory that is used to certify solutions to polynomial equations. This technique is implemented by Haunstein and Sottile, it is sometimes considered to be computationally heavy. There are some recent alternatives:
https://dl.acm.org/doi/10.1145/3326229.3326235
In general, developing certification techniques that are computationally tractable is an active research area. You might also take a look at this a bit older paper:
http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.756.9975&rep=rep1&type=pdf
