Starting from a geometric problem, I came up with a system of multivariate (many lines and points) polynomial equations where some equations are redundant (because they correspond to redundant geometric constraints. However, the constraints are not sufficient to define a geometry so the system is underdetermined.
In the case of no redundand geometric constraints, if I have $n$ unknowns and $m$ equations ($m$ is always $<n$ if there are no redundant equations), using Python SymPy I'm able to express $m$ variables in terms of the remaining $n-m$ free variables. I actually have to try a few times choosing random samples of $m$ variables to solve for, but the computation is fast.
However, most often, I have redundant constraints, and in such case I don't know how many variables I have to solve for (since I don't know the number of independent equations), so I should try samples of different length but the computation unsurprisingly takes too much time.
How can I find an equivalent system with no redundant equations? May Gröbner basis help?
Practical example
As a practical example please find below the Python with SymPy code to build such underdetermined system from a geometric problem:
import random
import sympy as sy
sy.init_printing(use_unicode=True, wrap_line=False)
from func_timeout import func_timeout, FunctionTimedOut # pip install func-timeout /// conda install -c conda-forge func_timeout
# cross product for 2D entities in homogeous coordinates
# see https://staff.fnwi.uva.nl/r.vandenboomgaard/IPCV20162017/LectureNotes/MATH/homogenous.html
def cross(a,b):
return a[1]*b[2]-a[2]*b[1],-a[0]*b[2]+a[2]*b[0],a[0]*b[1]-a[1]*b[0]
ne = 8 # number of edges
E = sy.Matrix(sy.MatrixSymbol('E', ne, 3)) # matrix of unknown parameters of edges, in homogenous coordinates
for idx in range(E.shape[0]):
E[idx,2] = 1 # set initial value of 3rd homogeneus coordinate
ort_E = sy.Matrix(tuple((E[e,1],-E[e,0],0) for e in range(E.shape[0]))) # edges orthogonal to originals passing through origin (0,0)
prl = [(1,3),(2,4),(3,5),(5,7),(6,8)] # parallel edges (human writable indices of edges starting at 1)
# orthogonal edges (human writable indices of edges starting at 1)
ort = [(1,2),(5,6)] # strictly necessary conditions
# ort = [(1,2),(2,3),(3,4),(4,1),(5,6),(6,7),(7,8),(8,1)] # 6 redundant conditions
prl_eqs = [cross(E[e1-1,:],E[e2-1,:])[2] for e1,e2 in prl] # 2 edges in homogeneous coordinates are parallel when the 3rd coordinate of their intersection point is zero
ort_eqs = [cross(ort_E[e1-1,:],E[e2-1,:])[2] for e1,e2 in ort] # 2 edges in homogeneous coordinates are orthogonal when the 3rd coordinate of the intersection point of <edge-orthogonal-to-the-1st> with <2nd-edge> is zero
eqs = [sy.simplify(eq) for eq in prl_eqs + ort_eqs] # list of equations derived from all conditions
print('\nNumber of edges: ',ne)
print('Number of edges DOF (degrees of freedom): ',2*ne)
print('\nNumber of equations (derived from conditions on edges): ',len(eqs))
print('Number of unknowns in equations: ',len(sum(eqs).atoms(sy.matrices.expressions.matexpr.MatrixElement)))
G = sy.groebner(eqs)
print('\nLength of Groebner basis: ',len(G))
for g in G: print(g)
print('\nTry direct solving on 7 random unknowns.\n')
eqs_unknowns = list(sum(eqs).atoms(sy.matrices.expressions.matexpr.MatrixElement))
step = 0
for _ in range(200): # with no redundant equations usually a solution is found within 100 iterations
step += 1
solving_unk = random.sample(eqs_unknowns,7)
print('Step ',step,', trying solving on ', solving_unk,end=' \r')
try:
sol = func_timeout(20, sy.solve, args=(eqs,solving_unk)) # fun_timeout runs `sy.solve(eqs,solving_unk)`
if isinstance(sol,dict) and all(x!=0 for x in sol.values()):
print('\nSolution found: ',sol)
break
except FunctionTimedOut:
continue
except Exception as e:
# print('Solving terminated by error: ',e)
continue
else:
print('\nNo solution found.')
To switch between no-redundant/redundant cases you just have to comment/uncomment lines 20-21 accordingly. With no redundant constraint/equations we have:
Number of edges: 9
Number of edges DOF (degrees of freedom): 18
Number of equations (derived from conditions on edges): 7
Number of unknowns in equations: 16
Length of Groebner basis: 19
and usually a solution (just using Sympy solve) is found within 100 iterations of 7 random chosen unknowns to solve the system for. With 6 redundant constraint/equations we have:
Number of edges: 9
Number of edges DOF (degrees of freedom): 18
Number of equations (derived from conditions on edges): 13
Number of unknowns in equations: 16
Length of Groebner basis: 32
and no solutions is found (at least after thousands iterations) for any 7 random chosen unknowns to solve system for. That's why I'm asking if Gröbner basis can help.
Just to know, the constraints define the following geometry:
If you want to have a square system and then discard the problem of redundant equations, consider that you have $n$ equations $$f_i=f_i(x_1,x_2,\cdots,x_n)=0 \qquad \text{for}\qquad i=1,2,\cdots m ~>~n$$
Consider now the norm $$\Phi=\Phi(x_1,x_2,\cdots,x_n)=\frac 12\sum_{k=1}^m \Big[f_i(x_1,x_2,\cdots,x_n)\Big]^2$$ and differentiate with respect to $x_i$ (for $i=1,2,\cdots,n$) to face $n$ equations $$g_i(x_1,x_2,\cdots,x_n)=\frac {\partial \Phi}{\partial x_i}=\sum_{k=1}^m f_k ~\frac {\partial f_k}{\partial x_i}=0$$