🎲 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 common zeros are still (at most) isolated points, but you cannot hand such a system to a total-degree start system, which needs a square target.
The fix 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 has extra solutions too. The workflow is always: 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, then 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 list(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 sorted(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 squaring (every row pushed to the
maximum degree 3) would track \(3^2 = 9\). You may also pass your own exact matrix,
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 split into separate variable groups, randomization works group by group and the multihomogeneous start system tracks far fewer paths.
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.
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¶
The whole tutorial as one runnable script – assemble nothing, just run it:
"""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 list(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()