Search for numerical solutions of underdetermined systems of quadratic equations

Question

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:

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

    P< x1,y1,z1,... > := PolynomialRing(RationalField(),3*N);
    I := ideal<P | my equations >;
    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.

2. Also I consider a calculation of Gröbner basis in grevlex order with subsequent transformation to lexicographic using FGLM-like algorithm.

3. 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.

Answer 1

When you have only type-1 equations, this is the graph layout problem. There are a variety of algorithms for finding approximate solutions to that problem; you might try seeing which ones can be adapted to your setting.

In particular, gradient descent looks like a good approach. For any particular assignment to the variables, let $e_i$ denote the error in the $i$th equation (e.g., the squared difference between the left-hand side and the right-hand side). Let $\Psi = e_1 + \dots + e_n$ be the sum of errors for all $n$ equations. Then you could try using gradient descent (or any other optimization method) to find an assignment to the variables that minimizes the objective function $\Psi$. Note that solutions might be sensitive to the initialization for the variables, so you might want to repeat many times with different random initializations. This might give you reasonable solutions.