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

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 < x < 1\,\wedge\,y = \sqrt{1-x^2}$
• $-1 < x < 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.

Answer 2

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.