Search for numerical solutions of underdetermined systems of quadratic equations I'm looking for an algorithm that can quickly generate an approximate real number solutions of an underdetermined quadratic system. My task is to explore its algebraic variety. I'm interested in a special cases of quadratic systems that arise from mechanics:

set of variables $X = \{ x_i, y_i, z_i | 1 \leqslant i \leqslant N\}$ (3*N variables for $N \in \mathbb{N}$),

a system consists of two kinds of equations:

1) $(x_i-x_j)^2+(y_i-y_j)^2+(z_i-z_j)^2 + p_{ij} = 0$ for $i, j \in [1,N]$ and $p_{ij} \in \mathbb{R}$.

2) a more general case of 1: $(x_i-x_k)(x_j-x_k)+(y_i-y_k)(y_j-y_k)+(z_i-z_k)(z_j-z_k) + a_{ijk} = 0 $ for $i, j, k \in [1,N]$ and $a_{ijk} \in \mathbb{R}$.

It's guaranteed that variety of solutions is stable under small changes of coefficients.
Here is an example with six-dimensional variety of solutions:

(x1-x2)^2 + (y1-y2)^2 + (z1-z2)^2 - 1 = 0
(x3-x2)^2 + (y3-y2)^2 + (z3-z2)^2 - 1 = 0
(x3-x4)^2 + (y3-y4)^2 + (z3-z4)^2 - 1 = 0
(x1-x4)^2 + (y1-y4)^2 + (z1-z4)^2 - 1 = 0
(x2-x1)*(x4-x1) + (y2-y1)*(y4-y1) + (z2-z1)*(z4-z1) - 1/3 = 0
(x1-x2)*(x3-x2) + (y1-y2)*(y3-y2) + (z1-z2)*(z3-z2) - 1/3 = 0
(x2-x3)*(x4-x3) + (y2-y3)*(y4-y3) + (z2-z3)*(z4-z3) - 1/3 = 0
(x1-x4)*(x3-x4) + (y1-y4)*(y3-y4) + (z1-z4)*(z3-z4) - 1/3 = 0


First, the algorithm is allowed to process this underdetermined system for some time. After that, it should get as input a set of constraints, like $x_1 = 0$ or $z_N = 1$, that turn the system into determined and quickly generate real solutions of it. For instance, one can find 8 real solutions of the system above after constraining all of $\{x_4,y_4,z_4,y_3,z_3,z_2\}$ to zero. 
I need the algorithm to be scalable so it can be applied for N=25-30, however high precision of solutions isn't required, 3-4 orders will do.

I've come up with a few approaches:


I tried to calculate Gröbner bases for such systems in Magma:


    P&lt; x1,y1,z1,... &gt; := PolynomialRing(RationalField(),3*N);
    I := ideal&lt;P | my equations &gt;;
    Groebner(I);


It works for zero-dimensional ideals, but coefficients were big even for N=7. In cases with dimensions >0, calculation took too long (with enough RAM) and I have been able to compute Gröbner bases only for small systems.
So I tried compute approximate Gröbner bases by running OpenF4 algorithm with datatype of coefficients changed on double, but the error grew significantly on each step of the computation.


Also I consider a calculation of Gröbner basis in grevlex order with subsequent transformation to lexicographic using FGLM-like algorithm.
I found this article https://arxiv.org/abs/1703.07810 and this looks useful for my purpose. I'm not familiar with such methods, but probably this is exactly what I need.


Here are my questions:

Is there a way to calculate an approximate Gröbner bases in reasonable time using my idea or any other approach, maybe for other order with subsequent transformation to lexicographic? 

Can this problem be solved without Gröbner bases using methods like №3 or any other?

Any comments on my ideas and new suggestions would be highly appreciated.
This is a bit late for you, perhaps, but it might be worth looking at cylindrical algebraic decomposition. CAD takes a general statement in the real closed fields (i.e. an arbitrary boolean statement involving equalities and inequalities), and returns a finite set of "solutions", each of which is a "procedure" to find those solutions by assigning the unknown variables.

I haven't tried it on your example (because I don't have any software handy to do it!), but as an example, here is the CAD of $x^2 + y^2 = 1$:


$-1 &lt; x &lt; 1\,\wedge\,y = \sqrt{1-x^2}$
$-1 &lt; x &lt; 1\,\wedge\,y = -\sqrt{1-x^2}$
$x = -1\,\wedge\,y = 0$
$x = 1\,\wedge\,y = 0$


There are four "solutions" in this decomposition, and each one gives you a way to find it, by first determining an assignment to $x$, and using that to get an assignment to $y$.

This would work extremely well if you knew which variables were going to be constrained, because you could specify that they appear first in each "solution". In this case, if you added the constraint that $x = 0$, you would just go through each solution and find that $y=1$ and $y=-1$ are the only solutions, since the last two don't apply.