🎲 Randomize an overdetermined system¶
Most tutorials solve a square system – as many equations as unknowns. But polynomial models are often overdetermined: more equations than unknowns. The zeros we compute in numerical algebraic geometry are still isolated points, but you cannot hand an overdetermined system directly to a path tracker, which needs a square target (to do linear algebra, for example).
The numerical algebraic geometric fix to overdetermination is randomization: replace the \(N\) functions with \(n\) generic combinations (\(n\) = the number of unknowns). Every common zero of the original system is a zero of each combination, so the genuine solutions survive – but the square randomized system may have extra solutions. The workflow is: randomize, solve the square system, then keep only the points that also satisfy the original system.
This tutorial does that twice, to show two different things. First, an ordinary system in one block of variables, and second a system whose variables split into groups – where the grouping cuts the path count.
Part 1 – one block of variables¶
Three curves in the plane: the unit circle, the line \(x = y\), and the vertical lines \(2x^2 = 1\). Along \(x = y\) the circle gives \(2x^2 = 1\), which is the third equation, so all three meet at \((\pm\tfrac{1}{\sqrt2}, \pm\tfrac{1}{\sqrt2})\) – two real points you can read off by hand.
x, y = bertini.Variable('x'), bertini.Variable('y')
original = bertini.System()
original.add_variable_group(bertini.VariableGroup([x, y]))
original.add_function(x*x + y*y - 1) # the unit circle (degree 2)
original.add_function(x - y) # the line x = y (degree 1)
original.add_function(2*x*x - 1) # the lines x = ±1/√2 (degree 2)
assert original.degrees() == [2, 1, 2] # three equations, two unknowns: overdetermined
randomize() returns a new square system; the original is left untouched.
For a single variable group it sorts the functions by descending degree and forms \(R = [\,I
\mid C\,]\) with a random block \(C\), so the randomized degrees are the \(n\) largest
original degrees and the total-degree path count is their product – here \(2 \times 2 = 4\),
the minimum a randomization can achieve.
randomized = original.randomize()
assert randomized.num_functions() == 2 # squared up
assert original.num_functions() == 3 # original unchanged
assert randomized.degrees() == [2, 2]
Note
The descending sort is what keeps the count minimal: the degree-1 line is folded into the degree-2 rows, rather than becoming a degree-1 target that a folded degree-2 function would inflate. The payoff grows with the degree spread. For degrees \((3, 2, 2)\) randomization tracks \(3 \times 2 = 6\) paths, where a degree-blind or dense squaring (every row pushed to the maximum degree 3) would track \(3^2 = 9\).
Note
You may also pass your own exact matrix R, calling
original.randomize(R), in which case the functions are kept in their given order.
Solve the square system, then filter: evaluate the original system at each computed point and keep the ones with a tiny residual.
bertini.random.set_random_seed(1) # reproducible generic coefficients + gamma
zd = bertini.ZeroDimSolver(randomized, endgame='cauchy', mptype='adaptive', startsystem='binomial')
zd.solve()
solutions = zd.all_solutions()
assert len(solutions) == 4 # the two we want, plus two extraneous
def satisfies_original(point):
coords = np.array([complex(c) for c in point])
return max(abs(complex(v)) for v in original.eval(coords)) < 1e-7
true_solutions = [np.array([complex(c) for c in s]) for s in solutions if satisfies_original(s)]
r = 1.0 / np.sqrt(2.0)
found = sorted((round(p[0].real, 4), round(p[1].real, 4)) for p in true_solutions)
assert found == sorted([(round(r, 4), round(r, 4)), (round(-r, 4), round(-r, 4))])
Two of the four points survive the filter: \((\tfrac{1}{\sqrt2}, \tfrac{1}{\sqrt2})\) and \((-\tfrac{1}{\sqrt2}, -\tfrac{1}{\sqrt2})\) – exactly where all three curves meet. The other two solve the random combination but not the original equations, so the filter drops them.
Part 2 – variables in separate groups¶
Now something different: when the unknowns are split into separate variable groups, randomization works group by group and the multihomogeneous start system can track far fewer paths than a total-degree start system.
The vehicle is a bilinear system: each function is degree 1 in \(x\) and degree 1 in \(y\) – with \(x\) and \(y\) in their own groups. Three equations whose only common solution is \((1, 1)\): from \(x = y\) and \(x + y = 2\) we get \(x = y = 1\), and indeed \(xy = 1\).
x, y = bertini.Variable('x'), bertini.Variable('y')
original = bertini.System()
original.add_variable_group(bertini.VariableGroup([x])) # x is its own group
original.add_variable_group(bertini.VariableGroup([y])) # y is its own group
original.add_function(x*y - 1) # xy = 1
original.add_function(x + y - 2) # x + y = 2
original.add_function(x - y) # x = y
randomized = original.randomize()
assert randomized.num_functions() == 2
Solve with the multihomogeneous start system, which exploits the grouping. Two bilinear equations on \(\mathbb{P}^1 \times \mathbb{P}^1\) meet in the multihomogeneous Bézout number of paths – here \(2\) (the coefficient of \(z_x z_y\) in \((z_x + z_y)^2\)) – whereas a total-degree start would track \(2^2 = 4\). Same answer, half the work. Simple example; the difference can be much more dramatic.
bertini.random.set_random_seed(3)
zd = bertini.ZeroDimSolver(randomized, endgame='cauchy', mptype='adaptive', startsystem='mhom')
zd.solve()
def satisfies(point):
coords = np.array([complex(c) for c in point])
return max(abs(complex(v)) for v in original.eval(coords)) < 1e-7
true_solutions = [np.array([complex(c) for c in s]) for s in zd.all_solutions() if satisfies(s)]
assert len(true_solutions) == 1
p = true_solutions[0]
assert abs(p[0] - 1) < 1e-7 and abs(p[1] - 1) < 1e-7
One solution survives: \((1, 1)\), exactly as read off by hand. The difference from Part 1 is
the grouping – declaring x and y as separate variable groups let the multihomogeneous
start system see the bilinear structure and track fewer paths to the same answer.
Complete example¶
"""Randomize an overdetermined system -- Bertini 2 tutorial.
Squares up overdetermined systems by generic randomization, then filters.
Run: python randomize.py
"""
import numpy as np
import bertini
def part1_one_block():
"""Part 1 -- one block of variables: three curves in the plane."""
x, y = bertini.Variable('x'), bertini.Variable('y')
original = bertini.System()
original.add_variable_group(bertini.VariableGroup([x, y]))
original.add_function(x*x + y*y - 1) # the unit circle (degree 2)
original.add_function(x - y) # the line x = y (degree 1)
original.add_function(2*x*x - 1) # the lines x = ±1/√2 (degree 2)
assert original.degrees() == [2, 1, 2] # three equations, two unknowns: overdetermined
# randomize returns a new square system; the original is untouched
randomized = original.randomize()
assert randomized.num_functions() == 2 # squared up
assert original.num_functions() == 3 # original unchanged
assert sorted(randomized.degrees()) == [2, 2]
# solve the square system, then filter against the original
bertini.random.set_random_seed(1) # reproducible generic coefficients + gamma
zd = bertini.ZeroDimSolver(randomized, endgame='cauchy', mptype='adaptive', startsystem='binomial')
zd.solve()
solutions = zd.all_solutions()
assert len(solutions) == 4 # the two we want, plus two extraneous
def satisfies_original(point):
coords = np.array([complex(c) for c in point])
return max(abs(complex(v)) for v in original.eval(coords)) < 1e-7
true_solutions = [np.array([complex(c) for c in s]) for s in solutions if satisfies_original(s)]
r = 1.0 / np.sqrt(2.0)
found = sorted((round(p[0].real, 4), round(p[1].real, 4)) for p in true_solutions)
assert found == sorted([(round(r, 4), round(r, 4)), (round(-r, 4), round(-r, 4))])
def part2_variable_groups():
"""Part 2 -- variables in separate groups: a bilinear system."""
x, y = bertini.Variable('x'), bertini.Variable('y')
original = bertini.System()
original.add_variable_group(bertini.VariableGroup([x])) # x is its own group
original.add_variable_group(bertini.VariableGroup([y])) # y is its own group
original.add_function(x*y - 1) # xy = 1
original.add_function(x + y - 2) # x + y = 2
original.add_function(x - y) # x = y
randomized = original.randomize()
assert randomized.num_functions() == 2
# solve with the multihomogeneous start system, which exploits the grouping
bertini.random.set_random_seed(3)
zd = bertini.ZeroDimSolver(randomized, endgame='cauchy', mptype='adaptive', startsystem='mhom')
zd.solve()
def satisfies(point):
coords = np.array([complex(c) for c in point])
return max(abs(complex(v)) for v in original.eval(coords)) < 1e-7
true_solutions = [np.array([complex(c) for c in s]) for s in zd.all_solutions() if satisfies(s)]
assert len(true_solutions) == 1
p = true_solutions[0]
assert abs(p[0] - 1) < 1e-7 and abs(p[1] - 1) < 1e-7
def main():
part1_one_block()
part2_variable_groups()
if __name__ == '__main__':
main()