Using monodromy to recover symmetries of polynomial systems

Galois/monodromy groups attached to parametric systems of polynomial equations provide a method for detecting the existence of symmetries in solution sets. Beyond the question of existence, one would like to compute formulas for these symmetries, towards the eventual goal of solving the systems more efficiently. We describe and implement one possible approach to this task using numerical homotopy continuation and multivariate rational function interpolation. We describe additional methods that detect and exploit a priori unknown quasi-homogeneous structure in symmetries. These methods extend the range of interpolation to larger examples, including applications with nonlinear symmetries drawn from vision and robotics.


Introduction
Structured systems of nonlinear equations appear frequently in applications like computer vision and robotics.Although the word "structure" can be interpreted in many ways, one of its aspects that is strongly connected to the complexity of solving is the algebraic degree of the problem to be solved.In many contexts, this may simply refer to the number of solutions of a system (usually counted over the complex numbers).However, if we adopt this definition without scrutiny, we may fail in certain special cases to detect additional structure such as symmetry.
To answer more refined questions involving structure, one can often consider a Galois/monodromy group naturally associated to the problem of interest.In this case, "problem" refers to a parametric family of problem instances which must be solved for different sets of parameter values.In our work, we are primarily interested in geometric Galois groups arising from algebraic extensions of functions fields of varieties over the complex numbers. 1urrently, a number of heuristic methods for computing Galois/monodromy groups using numerical homotopy continuation methods have been proposed and implemented, eg.[2,3].It is also fairly well-understood how Galois/monodromy groups encode important structural properties such as decomposability, or the existence of problem symmetries which may be expressed as rational functions known as deck transformations.Thus, Galois/monodromy group computation provides us a useful toolkit for detecting the existence of special structure.However, one key challenge remains: once we know that our problem does have such special structure, can we use this information to solve systems more efficiently?
Our work focuses on a natural first step towards addressing this challenge: given the data of a numerical Galois/monodromy group computation, can we recover formulas for the deck transformations?
In this paper, we describe and implement a novel framework that solves this problem, namely recovering symmetries.Our framework combines techniques for numerically computing Galois/monodromy groups with a new scheme for floating-point interpolation of multivariate rational functions which represent the underlying symmetries.This leads to an improvement upon the prior interpolation techniques described in [4].While the novelty of [4] was the combination of monodromy techniques with interpolation, this paper provides a substantial new contribution in the form of expanding the interpolation methods to detect a quasi-homogeneous structure and exploit this when interpolating (see Section 5).As a result, previously intractable problems are now solvable.For example, improvements are made on the five-point relative pose problem from [4] (see Section 6.1).New problems, such as the Perspective-3-Point (P3P) and radial camera relative pose, can also now be solved, as described in Section 6.2 and Section 6.3, respectively.
In Section 2, we provide some context for our approach by considering related previous works.In Section 3, we establish terminology and useful background facts.In Section 4, we describe our a basic approach to interpolating deck transformations and illustrate it on simple examples.In Section 5, we describe a more sophisticated variant of this approach, which exploits the quasihomogeneous structure of certain symmetries to reduce the size of the associ-also why we cannot pick inputs for the interpolation problem arbitrarily.With that said, we point out that assuming exact inputs could also be relevant if, say, certified homotopy continuation (see [16,17,18,19]) is used, augmented by some additional postprocessing.
Additional methodology introduced in this version of the paper relies on detecting structure that was a priori unknown quasi-homogeneous structure.These techniques rely on discovering scaling symmetries by computing the Smith Normal Form of an integer matrix, a well-studied computational problem.Detecting such scaling symmetries with the integer linear algebra techniques has been done before in multiple different contexts [20,21,22,23,24].The advances in this paper highlight the use of the Smith Normal Form for both the variables and parameters, instead of simply the variables.

Background
In this work, we are interested in solving polynomial systems whose solutions correspond to points in a generic fiber of a branched cover of complex algebraic varieties.Here we collect some definitions and theoretical facts that we need to work within this framework.The section concludes with Proposition 3.10 and Corollary 3.11, which justify the correctness of our general interpolation setup used in Sections 4 and 5.
Definition 3.1.Let X and Z be irreducible algebraic varieties of dimension m over the complex numbers.A branched cover is a dominant, rational map f : X Z.The varieties X and Z are called the total space and the base space of the cover, respectively.The number of (reduced) points in the preimage over a generic z ∈ Z is called the degree of f, denoted deg(f ).
Essentially, the base space Z in Definition 3.1 can be thought of as a space of parameters or observations.The fiber f −1 (z) over some particular z ∈ Z should usually be understood as the solutions of a particular problem instance specified by z.Oftentimes, Z may be assumed to be an affine space C m , and in this case we write p ∈ C m for parameter values.The assumptions that f is dominant and dim X = m imply that there is a finite, nonzero number of solutions for almost all parameters.Counting solutions over C, that number is deg(f ).Additionally, the total space X is often either 1. an irreducible variety consisting of problem-solution pairs, for some system of polynomials f 1 , . . ., f k ∈ C[x, p], with projection f : X → C m given by f (x, p) = p, or 2. an affine space of unknowns X = C m , and f : C m C m .Cases (1) and (2) for the total space X given above are closely related.Indeed, (2) reduces to (1) if we take X to be the graph of f.Conversely, it can often be the case that the variety X has a unirational parametrization p : C m X.In this case, ( for all x in a nonempty Zariski-open subset of X.The maps g and h give a decomposition of f. The projection f : X → C 23 given by f (a, . . ., z) → (a, . . ., w) is a branched cover of degree 32.If we let Y be the set of all (a, . . ., w, x, ŷ) ∈ C 25 such that then g : X → Y given by g(a, . . ., w, x, y, z) = (a, . . ., w, xz, yz) and h : Y → Z given by h(a, . . ., w, x, ŷ) = (a, . . ., w) show that f is decomposable in the sense of Definition 3.2.Here we have deg(h) = 8 and deg(g) = 4.
The Galois/monodromy group is an invariant that allows us to decide whether or not a branched cover is decomposable, without actually exhibiting a decomposition.We recall the basic definitions here.For a branched cover f : X Z, fix a dense Zariski-open subset U ⊂ Z such that f −1 (z) consists of d = deg(f ) points.Over a regular locus, the branched cover f restricts to a d-sheeted covering map in the usual sense given by f −1 (U ) → U.For any basepoint z ∈ U, we may construct via path-lifting a group homomorphism from the fundamental group π 1 (U ; z) to the symmetric group S d .
More precisely, if γ : [0, 1] → U is any map that is continuous with resepct to the Euclidean topology, then the unique lifting property [25,Prop. 1.34] implies that there are precisely d continuous lifts γ1 , . . ., γd : [0, 1] → π −1 (U ) satisfying f • γi (t) = γ i (t) for all i = 1, . . ., d and t ∈ [0, 1].In particular, γi (0), γi (1) ∈ f −1 (z), and there is a permutation σ γ that sends each γi (1) to γi (0).One may check that this permutation is independent of the chosen representative γ of the homotopy class [γ] ∈ π 1 (U ; z).Thus, for our chosen U and z we may define the monodromy representation, This gives a group homomorphism, whose image is a subgroup of S d , which turns out to be independent of the choice of U and z.Definition 3.5.The Galois/monodromy group of a branched cover f is the subgroup of S d given by the image of the map (2).
The abstract structure of the Galois/monodromy group, although interesting, is not our main focus.Instead, we will be mainly interested in the action of this group given by (2).Since X is irreducible, this action is transitive (see eg. [26,Lemma 4.4,p87]).
The monodromy action also provides a clean characterization of decomposable branched covers.Recall that the action of a group G on a finite set B is said to be imprimitive if there exists a nontrivial partition that for any g ∈ G and B i there exists a B j with g • B i = B j .If B has d elements and G is a finite, transitive subgroup of S d , it follows that the subsets B i must all have the same size.The sets B 1 , . . ., B k are called blocks of the imprimitive action, and are said to form a block system.Proposition 3.6.(See eg.[11,Proposition 1].)A branched cover is decomposable if and only if its Galois/monodromy group is imprimitive.
In general, a transitive, imprimitive permutation group has a block system B 1 , . . ., B k whose blocks all have the same size l, and is thus permutationisomorphic to a subgroup of the wreath product S l ≀ S k .Unlike the previous example, there are a number of surprising cases of decomposable branched covers where the Galois/monodromy group is a proper subgroup of the associated wreath product: for instance, the five-point problem of Section 6.1.
We point out that Proposition 3.6 dates back, at least in some form, to work of Ritt on polynomial decompositions [27].This work is directly related to decomposition problems for polynomials and rational functions studied in computer algebra (see eg. [28,29]).
However, the main focus in this paper is not decomposability per se.Rather, we are interested in a property that is usually stronger: the existence of symmetries.A natural, and general, notion of symmetry can be obtained by studying the embedding of function fields f * : C(Z) C(X) induced by a branched cover.The field extension C(X)/C(Z), although not usually a Galois extension, may nevertheless a have a nontrivial group of automorphisms.These automorphisms correspond to rational maps Ψ : X X with f • Ψ = f.Topologically, these comprise the group of deck transformations of f.Proposition 3.8 below explains the relationship between deck transformations and decomposability, and provides an analogue of Proposition In the final results of this section, Proposition 3.10 and Corollary 3.11, we use the terminology generic path for a given branched cover f : X Z.This means a path α : [0, 1] → U where U is some suitably small set, either a regular locus in Z or its preimage in X.In the former case, we write α x for the unique lift of a path α through f starting at x ∈ f −1 (α(0)).Proposition 3.10.Let f : X Z be a branched cover with a fixed generic point x ∈ X.Then the value of a deck transformation Ψ ∈ Deck(X/Z) at a generic point x ′ ∈ X is completely determined via path-lifting by where it sends x.Explicitly, where α is a generic path in X from x to x ′ (see Figure 1).
Proof.We refer to the proof of [25,Prop. 1.33] and the general definition of a lift given on [25, p. 60].The deck transformation Ψ is a lift of f to X in the sense of this definition.This means the proof of Proposition 1.33 can be applied to construct a deck transformation Ψ ′ with Ψ ′ (x) = Ψ(x).This construction uses lifts of a generic path α to construct Ψ ′ , with the additional property that Ψ ′ (x ′ ) = (f • α) Ψ(x) (1).The unique path-lifting property then implies that Ψ(x ′ ) = Ψ ′ (x ′ ).
A consequence of Proposition 3.10 is that the correspondence between solutions for a fixed set of parameters under a fixed deck transformation Ψ is preserved under path-lifting.Corollary 3.11.Let f : X Z be a branched cover and Ψ ∈ Deck(X/Z).Let z ∈ Z be a generic point and β be a generic path in Z starting at z (see Figure 2).Then for x ∈ X z we have In other words, the points in two lifts of β-one starting at x, the other starting at Ψ(x)-are conjugate under Ψ (see Figure 2).Proof.By Proposition 3.10, Ψ( (1), which, in turn, is

Basic method -dense interpolation
Consider a branched cover encoding problem-solution pairs (x, p), As mentioned in the introduction, we may compute the Galois/monodromy group of f using numerical homotopy continuation.This is possible provided that we make the following assumptions about how our branched cover is given as input.
Assumption 4.1.For the branched cover defined in (5), assume that n rational functions f 1 , . . ., f n vanishing on X are known, and that we have access to a sampling oracle that produces generic (x * , p * ) ∈ X such that the n × n Jacobian ∂f ∂x (x * , p * ) has rank n.Assumption 4.1 is often satisfied in practice, including cases where even a set-theoretic description of X is not known.Additionally, we assume that homotopy continuation-specifically, coefficient parameter homotopy-can be used to track d known solutions for fixed, generic parameter values p * (corresponding to f −1 (p * )) to d solutions for some other parameter values p ∈ C m (corresponding to f −1 (p)).These parameter homotopies are the basis of the unspecified subroutines in lines 1 and 9 of Algorithm 1.
An important observation is that we can interpolate each of the coordinate functions ψ j (x, p) in ( 6) independently.We assume that the rational function ψ j contains only monomials up to total degree D. Since these monomials may or may not involve the parameters p, we distinguish the parameterdependent and parameter-independent settings, in which we take the number of monomials t to be either (for parameter-independent ψ j (x)) Our task is then to recover two vectors of unknown coefficients such that ψ j can be represented on X as In Equation ( 7), the vectors α k , β k ∈ Z n+m ≥0 range over a suitable set of multidegrees, depending on whether we are in the parameter-dependent or parameter-independent setting.If we know that ( this gives a homogeneous linear constraint on a and b, Suppose we have already computed permutations generating the monodromy group based at parameter values p 1 ∈ C m , and let for some element of the centralizer σ ∈ Cent S d (Mon(f, p 1 )) corresponding to Ψ. Now, Corollary 3.11 implies that we may obtain additional sample points satisfying (8) by tracking parameter homotopies using the system f 1 , . . ., f n .Specifically, we may track the solution curves with initial values x 1 , x ′ 1 from p 1 to generic p i ∈ C m for i = 1, . . ., 2t, which then allows us recover the coordinate functions of Ψ. Proposition 4.1 (Correctness of Algorithm 1).Suppose that ψ j in (6) can be represented as the quotient of polynomials with degree ≤ D and t monomials each.For a sufficiently generic sample ) for all i, suppose a ⊤ b ⊤ is a solution to the 2t linear equations given by ( 8) for i = 1, . . ., 2t, which lies outside the span of all solutions with a = 0 or b = 0. Then the rational function obtained by restricting ψ a,b (x, p) to X equals ψ j .
Proof.The assumption that a ⊤ b ⊤ is a nontrivial linear combination of solutions with a, b ̸ = 0 ensures that ψ a,b is a well-defined, nonzero rational function on X.Such a function of the form ( 7) is determined by its values on 2t generic points of X.Since ψ j , by assumption, is also such a function, the 2t linear constraints (8) force ψ j and ψ a,b to agree on X.
Thus, to interpolate ψ j , we may determine from the linear equations ( 8) a 2t × 2t Vandermonde-type coefficient matrix A. We represent the nullspace of A by the column-span of a matrix N with 2t rows.Although Proposition 4.1 can be viewed as a uniqueness statement, the matrix N will generally have more than one column, even for generic samples (x 1 , p 1 ), . . ., (x 2t , p 2t ) ∈ X.The "extra" columns of N appear for two reasons: 1.There may exist different representatives of ψ j on X of the form (7), whose coefficient vectors are linearly independent.2. The nullspace of A may contain spurious solutions not satisfying the hypothesis a = 0 or b = 0 in Proposition 4.1.For instance, fixing b = 0 we may interpolate polynomial functions of the form t k=1 a k • (x, p) α k vanishing on X.In the same way, fixing a = 0 we interpolate polynomial functions of the form t k=1 b k • (x, p) β k vanishing on X.For some applications it may be necessary to pick a sparse representative from the nullspace of A. In general, finding the sparsest vector in the nullspace of a matrix is NP-hard [30].Nevertheless, in many cases we may find a relatively good sparse representative by looking at the reduced row echelon form of N ⊤ for some particular ordering of its columns and picking one with the fewest zeros subject to the additional constraints a, b ̸ = 0. We illustrate some of the choices involved on two simple examples.
The Galois/monodromy group and deck transformation group are both S 2 .When interpolating a nontrivial deck transformations of degree D = 1, we obtain the reduced row echelon form for N ⊤ below.
We see that Ψ(x, p) has 2 different representatives 1 x and −x − p, which both agree on X.There is no clear choice of "best representative".In terms of sparsity, the representative 1 x is superior.However, one might instead prefer −x − p since it is a polynomial.
which has a unique non-identity deck transformation Ψ = (ψ 1 , ψ 2 ).If we interpolate parameter-dependent deck transformations, we may find matrices A 1 and A 2 representing Ψ which are 8 × 8.The reduced row echelon forms of the transposed nullspaces are and for ψ 2 (x, y, p) we have If we are not interested in the sparsest representative, then we may take In this example, it is possible to find the sparsest polynomial representative for ψ 1 by solving an auxiliary linear system.In other words, we compute a linear combination of rows a ⊤ b ⊤ = r ⊤ N 1 such that b ⊤ = 1 0 ⊤ and a ⊤ contains the minimum number of zeros.First, to obtain b ⊤ = 1 0 ⊤ , we solve a linear system obtained from the right 4 × 4 block of N 1 , The general solution of this system is given by Using r to form a linear combination of rows now from the left 4 × 4 block of N 1 , we obtain a ⊤ = −1 r r + 1 r + 1 .
To maximize the sparsity, we may set r = −1 to obtain a ⊤ b ⊤ = −1 −1 0 0 1 0 0 0 which encodes the function Our pseudocode in Algorithm 1 outlines a degree-by-degree procedure for interpolating the full set of deck transformations up to a given degree D * .To implement such a procedure, there are many design choices that could improve performance or meet the needs of a particular task.Among the design choices, we note that the monodromy, parameter homotopy, and get_representative subroutines on respective lines 1, 9, and 17 are left unspecified.Our implementation relies on HomotopyContinuation.jl for the first two of these subroutines.For get_representative, our implementation chooses the sparsest row in the rref matrix N j .For the final output of line 19, we heuristically truncate "small" entries of N j of size < 10 −5 .
Finally, we note the following improvements to the pseudocode in Algorithm 1, which we have used in our implementation.
1. Computing the monodromy group and centralizer in lines 1-2 is an offline task which only needs to be performed once for a given family of systems.
Input: F = (f 1 , . . ., f n ) and (x * , p * ) as in Assumption 4.1, representing f as in (5); an upper bound for the total degree D * of monomials in each interpolant Output: Partially-specified rational maps representing the group of deck transformations, {Ψ 1 , . . ., Ψ q } = Deck(f ), with all coordinate functions representable in degree ≤ D * specified (x (1) In practice, we might only need to recover generators of the deck transformation group.The needed modifications are trivial, since deck transformations are interpolated independently.3. To restart the computation at a higher degree limit D * , one can use previously-computed samples from X.In principle, one can also draw > 2t samples and compute the nullspace of the resulting rectangular matrices A j .4. To minimize the number of calls to the parameter homotopy subroutine, one can attempt to track samples in "batches": since every fiber consists of d points and each point gives 1 constraint on ψ j , then we need to obtain a complete set of d solutions for r different sets of parameters (including p) such that In our experience, this strategy can work well, but comes with the additional caveat that the samples need not satisfy the genericity conditions of Proposition 4.1, since multiple parameter values are duplicated.We encoutered one (ultimately benign) instance of this phenomenon in our study of Alt's problem Section 6.4.In this example, we had d = 8652 > 2t = 650, and this strategy resulted in many more spurious rows in N due to all samples using the same parameter values.

Continuous symmetries and multigraded interpolation
As we have already seen, the deck transformations of a branched cover provide us with one useful way to formalize the study of symmetries of a parametric family of polynomial systems.However, this is by no means the only useful notion of symmetry-deck transformations, although they may depend on both parameters and unknowns, only act nontrivially on the unknowns.In this section, we consider symmetries that act nontrivially on both parameters and unknowns.We assume some basic notions about algebraic groups acting rationally on algebraic varieties, eg. as treated in [31, §1].After a general discussion of branched covers which are equivariant with respect to a general rational action, we specialize to the case of scaling symmetries, also known as quasi-torus actions.The deck transformations of a branched cover equipped with these symmetries are quasi-homogeneous with respect to an associated multigrading.This leads to Algorithm 2, a multigraded refinement of Algorithm 1, which can produce Vandermonde matrices of considerably smaller size.

Equivariant Branched Covers
Definition 5.1.Let f : X Z be a branched cover and G be an algebraic group acting rationally on both X and Z.We say for all (x, g) in some Zariski-open subset of X × G where both sides of (10) are defined.If G is irreducible and dim G > 0, we say that the elements of G are continuous fiber-respecting symmetries of f .
Since G is a smooth variety, we note that its irreducibility is equivalent to its connectedness in the analytic topology.
The branched covers associated to systems of algebraic equations occurring in applications are typically equivariant with respect to some underlying symmetries of the problem.For example, the pose estimation problems considered in Sections 6.1 to 6.3 are naturally associated to branched covers which are invariant under certain actions of the special Euclidean group SE C (3), the complexified group of rotations and translations in 3-space, as well as the scaling symmetries that are the primary focus of this section.
The next result is likely known.We include a proof for completeness.
Proposition 5.2.Let f : X Z be a branched cover with a group of continuous fiber-respecting symmetries G, and let Ψ ∈ Deck(f ) be any deck transformation.Then Ψ is G-equivariant-that is, for all (g, x) in some dense Zariski-open subset of X × G.
Proof.Let U G ⊂ G, U X ⊂ X be nonempty Zariski-open subsets such that both the group action map and Ψ are defined for all pairs in U X × U G .Let θ g : [0, 1] → U G be a smooth path in G from id G to g.For any x ∈ U X , we define a path U X from x to g • x by the rule To show the two paths γ x,g and γ Ψ(x),g have the same projection by f , consider the chain of equalities = θ g (t)•f (Ψ(x)) where ( 1) and ( 3) follow from G-equivariance of f and (2) follows from the fact that f • Ψ = f .Now, using Corollary 3.11, simply note In what follows, we will restrict our attention to special cases where the group G is abelian which leads directly to algorithmic improvements.However, in view of the ubiquity of equivariant branched covers for more general groups, it would be interesting to study any further applications of this structure to the problem of deck transformation recovery.

Scaling symmetries from Smith Normal Forms
Consider again a branched cover of problem-solution pairs f : X → C m as in Assumption 4.1, where X ⊂ V(F ) ⊂ C n+m is an irreducible affine variety locally defined by the square polynomial system . . .
with all a ij nonzero.Each set of exponents {α i1 , . . ., α ik i } in ( 12) is called the monomial support of f i , and denoted supp(f i ).
For any λ ∈ C * and u ∈ Z n+m , and (x, p) ∈ C n+m , we define We shall be interested in detecting certain scaling symmetries on X of the form λ • (x, p) = λ u ⊙ (x, p).We consider both continuous symmetries C * ×X X, as well as discrete symmetries Z s ×X X, where Z s denotes the integers modulo s.In the latter case, this means λ s = 1.Such symmetries can be detected using well-known methods of integer linear algebra.From the system F in (12), we form a matrix of shifted exponent vectors and compute its Smith Normal Form, The blocking in ( 15) is chosen so that every row of U 1 lies in the left-nullspace of A. Similarly, each row of U 2 becomes a left null vector of A after reduction modulo some d i .We first consider continuous symmetries arising from U 1 .
Definition 5.3.Fix f, F, and U 1 ∈ Z r×(n+m) as above, and consider the rows of U 1 , denoted u T 1 , . . ., u T r ∈ Z 1×(n+m) .We define an associated group of continuous scaling symmetries G con to be the image of the homomorphism As an abstract group, we have G con ∼ = (C * ) r , and this group acts on C n+m and C m via the coordinate-wise product ⊙ as in (13).In fact, G con also acts on X.To see this, first observe for any point (x, p) ∈ X ⊂ V(F ) that orbits G con • (x, p) are all contained in V(F ) by construction.On the other hand, (1, . . ., 1) • (x, p) ⊂ X, so by connectivity we may conclude G con • (x, p) ⊂ X.Thus, the branched cover f is equivariant with respect to the action of the connected algebraic group G con .By Proposition 5.2, it follows that any deck transformation Ψ ∈ Deck(f ) must commute with the action of G con .This implies that the coordinate functions of Ψ are quasi-homogeneous.Definition 5.4.A rational function p(x) q(x) ∈ C(X) is said to be quasi-homogeneous with respect to a continuous scaling u if there exists d ∈ Z with We may verify that a quasi-homogeneous rational function can be represented as the quotient of two quasi-homogeneous polynomials.
Proposition 5.5.Consider a nonzero rational function on an irreducible affine variety X, represented as p(x) q(x) where p and q are both polynomials that do not vanish on X.If the function is quasi-homogeneous with respect to a free scaling u, then it can be represented as a(x) b(x) on X, where both a(x) and b(x) are both quasi-homogeneous polynomials with respect to u.Moreover, we may assume all monomials in a (resp.b) occur in p (resp.q), supp(a(x)) ⊆ supp(p(x)), supp(b(x)) ⊆ supp(q(x)).
Proof.The univariate Laurent polynomial ring C[x][λ −1 ] is equipped with a natural Z-grading with respect to degrees in λ.Let be expressions of p(λ u ⊙ x) and q(λ u ⊙ x) in homogeneous components with respect to this grading.Note that the integers r i (resp.s i ) are distinct, and without loss of generality we may assume that none of the p i or q i vanish on X.From the identities Clearly we have f (λ, x) = 0 for all x ∈ X and λ ∈ C * .Let . Suppose that there exists k ∈ [h 1 ] such that for all j ∈ [h 2 ] we have r k ̸ = t j .Then we can rewrite (16) as for suitable polynomials f i (x) with e i are all distinct.Since a univariate Laurent polynomial that evaluates to zero at infinitely many points is the zero polynomial, we obtain a contradiction Thus, for every i ∈ [h 1 ] there exists j ∈ [h 2 ] such that r i = t j .By symmetry, for every j ∈ [h 2 ] there exists i ∈ [h 1 ] such that r i = t j .In particular, Again, since r i are all distinct, we have Thus, for any (i, j) ∈ H, the rational function p(x)/q(x) has an equivalent representation p i (x)/q j (x).Taking (a, b) = (p i , q j ) gives the desired conclusion.
We now turn our attention to discrete scaling symmetries arising from the submatrix U 2 appearing in the Smith Normal Form (15).These will form a finite group G dis , defined analogously to G con .However, to do this we must consider some technicalities not present in the continuous case.Write where the blocking is chosen so that each U 2,i ∈ Z r i ×(n+m) corresponds to a distinct elementary divisor d i in D.
For each of the blocks U 2,i in (17), let the rows of this matrix after reduction modulo d i be denoted by . Write U i for the group of d i -th roots of unity in C * , and consider the homomorphism ((λ 1,1 , . . ., λ 1,r 1 ), . . ., Definition 5.6.The group of discrete scaling symmetries associated to f and F consists of all scalings in the image of ( 18) that map X to itself and commute with all deck transformations in Deck(f ).
Abstractly, G dis is a finite abelian group, As in the continuous case, f is G dis -equivariant.However, both requirements that elements of G dis map X to itself and commute with Deck(f ) are now nontrivial.
Example 5.7.Consider the system, with (n, m) = (4, 3), defined by The basic idea behind constructing the system F in this example is to take a Galois cover with Galois group S 3 , and apply a birational change of coordinates so that one of its deck transformations becomes a scaling.The following Macaulay2 [32] code was used: The variety V(F ) ⊂ C 7 has the irreducible decomposition Let X be the first of these irreducible components.We see that the scaling on C 7 that sends x 1 → −x 1 and fixes all other coordinates does not map X to itself.Moreover, by construction we have that the scaling on X that sends x 4 → −x 4 and fixes all other coordinates does not commute with the full deck transformation group Deck(f ) ∼ = S 3 .Since neither symmetry belongs to G dis , this shows why both requirements of Definition 5.3 are necessary.
Example 5.7 underscores the need for a procedure for computing G dis .We now summarize a "probability-one" procedure based on homotopy continuation that accomplishes this task.For each elementary divisor d i , consider all Z d i -linear combinations of the modulo-d i reduced rows of U 2,i .For each linear combination u ∈ Z n+m d i , we check if the associated scaling commutes with each deck transformation.This is done with a probability-one homotopy test-generate random intermediate parameter values p 1 ∈ C m and track along two linear segment parameter homotopies-first from p 0 to p 1 , then from p 1 to λ u ⊙ p 0 .2By Corollary 3.11, the following holds with probabilityone: the discrete scaling u commutes with Ψ ∈ Deck(f ) if and only if the endpoint obtained by tracking Ψ(x 0 ) along the homotopies is the same as the endpoint obtained by first tracking the start point x 0 and then applying Ψ.This probability-one test succeeds for all Ψ ∈ Deck(f ) if and only if u determines an element of G dis .

Quasi-homogeneous interpolation
Retaining the notation established earlier in the section, we now consider the multidegree map, ie. the group homomorphism where r ′ i ≤ r i and the Z d i -rowspan of each U ′ 2,i is contained in the Z d i -rowspan of U 2,i .We assume that the multidegree map is compatible with f in the sense that the deck transformations commute with the group of scaling symmetries in G con × G dis , obtained by applying Hom(•, C * ) to the image of (19).The procedure for computing G dis described above allows us to easily determine a set of maximally compatible U ′ 2,1 , . . ., U ′ 2,k from U 2,1 , . . ., U 2,k .Our approach to interpolation of quasi-homogeneous deck transformations is summarized by the pseudocode of Algorithm 2. The overall structure is much the same as Algorithm 1.The main difference, responsible for the improved performance, is that a potentially much smaller basis of monomials for interpolation is chosen, so as to incorporate the quasi-homogeneity of Input: F = (f 1 , . . ., f n ), (x * , p * ), D * as in Algorithm 1, and the f -compatible multidegree map specified by U = (U 1 , U ′ 2,1 , . . ., U ′ 2,k ) as in (19).Output: As in Algorithm 1 (x (1) , . . ., x (d) ), Mon(f, p * ) ← run_monodromy(F, x * , p * ) {σ 1 , . . ., σ q } ← Cent S d (Mon(f, p * )) // WLOG q > 1, σ 1 = id for i ← 2; i ≤ q; i ← i + 1 do for j ← 2; j ≤ q; j ← j + 1 do Deck(f ) encoded by the symmetry group G con × G dis .An even more refined strategy, not pursued here, would be to consider quasi-homogeneity at the level of individual deck transformations, or even their coordinate functions.Nevertheless, experiments in Section 6 show that working with just the symmetries in G con × G dis already produces considerable savings.
On line 8 of Algorithm 2, the set of all monomials of degree ≤ D * is divided into monomial classes.The sizes of these classes govern complexity of each interpolation step.Each class consists of key-value pairs, where the values are monomials whose corresponding key records the multidegree of the numerator of some candidate rational function representation.The multidegree of the denominator in such a representation is uniquely determined by the class.Pseudocode for computing the multidegree of the denominator from a given class is provided in Algorithm 3 (according to the G con × G dis -equivariance of Deck(f ) (11)).
, and j indexing the rational function ψ j in Ψ ∈ Deck(f ) Output: The corresponding denominator multidegree, ) Algorithm 3: Subroutine called on line 14 of Algorithm 2.

Examples and Experiments
Our implementation of Algorithms 1 and 2 is written in the Julia programming language.It depends on the following packages: 1. HomotopyContinuation.jl [34], 2. AbstractAlgebra.jl,part of the Nemo system [35], and 3. GAP.jl, part of the OSCAR system [36].
Basic instructions for using the code (which is currently under active development) can be found at the link: https://multivariatepolynomialsystems.github.io/DecomposingPolynomialSystems.jl/dev/.
All timings reported were obtained with a 2022 Mac M1 with 8GB RAM.

Five-point relative pose
One of the most well-known minimal problems in computer vision is the classical five-point problem as shown in Figure 3.While many solvers exist for this problem [37], and the symmetry is well-known, this section aims to show how the methods in this paper can recover this symmetry without any a priori knowledge.
We consider two slightly different formulations of the five-point problem.Section 6.1.1 presents an "inhomogeneous" formulation that appeared previously in the conference paper [4].For this formulation, we were able to recover only the coordinate functions of the deck transformation which are parameter-independent.Section 6.1.2shows how incorporating scaling symmetries in Algorithm 2 leads to a much more tractable problem, in which the full deck transformation can be recovered.A summary of these experiments is given in Table 1 below-we refer to the corresponding subsections for details.

Inhomogeneous formulation
We first consider the following setup.There are 5 correspondences between 2D image points x 1 ↔ y 1 , . . ., x 5 ↔ y 5 .These 2D data points are 2×1 vectors, and are assumed to be images of 5 world points under two calibrated cameras, where the two camera frames differ by a rotation R and a translation t.The task for this problem is to solve for the relative orientation [R | t] ∈ SE R (3) between the two cameras and each of the five points in 3D space, as measured by their depths with respect to the first and second camera frames.
Writing α 1 , . . ., α 5 for the depths with respect to the first camera and β 1 , . . ., β 5 for the depths with respect to the second camera, solutions to the five-point problem must satisfy a system of polynomial equations and inequations: The unknown depths and translation t are defined in projective space, meaning t, α 1 , . . ., α 5 , β 1 , . . ., β 5 can only be recovered up to a common scale factor.One option to remove this ambiguity is to treat these unknowns as homogeneous coordinates on a 12-dimensional projective space, then for generic data x 1 , . . ., x 5 , y 1 , . . ., y 5 , there are at most finitely many solutions in (R, t, α 1 , . . ., α 5 , β 1 , . . ., β 5 ) ∈ SO C (3)×P 12  C to the system (20).This finiteness is what creates the minimal problem structure.In practice, these solutions may be computed by working in a fixed affine patch of P = P 12 C such that the inequality t ̸ = 0 is satisfied (e.g., a ⊤ t = 1 for a random a ∈ C 3 ).
There are exactly 20 solutions over the complex numbers for generic data in Z = (C 2 ) 5 × (C 2 ) 5 .The solutions to (20) are naturally identified with the fibers of a branched cover f : X → Z, where With our chosen formulation, the branched cover f has a single deck transformation Ψ known as the twisted pair, whose coordinate functions are We see that Ψ consists of coordinate functions of total degree at most 3.The effect of this deck transformation on solutions to the five-point problem is illustrated in Figure 4.The coordinate functions Ψ R , Ψ t are parameterindependent, whereas the Ψ α i , Ψ β i are parameter-dependent.
We ran Algorithm 1 on the formulation (20) with the upper bound for the total degree D * = 3.However, when running Algorithm 1, we considered only the parameter-independent setting, for which t = 22+3 3 = 2300.In the parameter-dependent setup, we would have 2t = 2 22+20+3 3 = 28380 coefficients to interpolate.This exceeded the capacity of our machine.
The computation described above succeeded in recovering Ψ R and Ψ t in (21).For the coordinate functions Ψ α i , Ψ β i , no reasonable representative was found-all rows of N were such that a ≈ 0 or b ≈ 0. These coordinate functions remain "missing" in Algorithm 1.This is expected-the expressions for these functions in (21) are parameter-dependent.
Further details for this formulation are given in Table 1.
Figure 4: Twisted pair symmetry of the five-point problem.

Quasi-homogeneous formulation
In the quasi-homogeneous formulation we consider the image points to be the points in P 2 , i.e. we introduce 10 more parameters and consider the system of polynomial equations and inequations: where x 1 , . . ., y 5 are now 3 × 1 parameter vectors.The continuous scaling symmetries of the formulation ( 22) are given by where i = 1, . . ., 5, and hence G con ∼ = (C * ) 11 .The free scalings are represented by In this example, all discrete scalings detected by the Smith Normal Form commute with Ψ. Hence r 1 = r ′ 1 in the notation of Section 5.3.This commutativity can be detected with the probability-one homotopy test described in Section 5.2.This can also be understood a priori because these symmetries are instances of a continuous, non-scaling symmetry expressed as in (24), where R 1 , R 2 ∈ SO C (3) may be arbitrary rotations.
The last column of Table 1 shows that the quasi-homogeneous approach allows us to run parameter-dependent interpolation of the twisted pair (21).Using line 8 of Algorithm 2 we partition the monomials up to total degree 3 in both unknowns and parameters into classes w.r.t. the multidegree map given by (U 1 , U ′ 2,1 ).We obtain 14339 classes of monomials, the largest of which has size 36.This partitioning makes parameter-dependent interpolation feasible.
We close our discussion of the five-point problem with a remark that the depths α i , β i and world points X i can easily be recovered from known rotations and translations (R, t) (assuming generic data) using linear algebra.In the homogeneous formulation, this is based on the relation Thus, the five-point problem illustrates that interpolating deck transformations may be easier after eliminating certain variables.On the other hand, our success with interpolating the twisted pair on depths and world points showcases the utility of techniques that exploit quasi-homogeneity, like Algorithm 2, and selecting a formulation amenable to these techniques.

Perspective 3-Point
Another famous algebraic problem associated with computer vision is the so-called P3P problem.In fact, early studies of this problem date back centuries to Lagrange and Grunert-see [38] for a more complete history.Similar to our discussion of the five-point problem, we may consider an inhomogeneous formulation of the problem as in Section 6.1.1, and a quasi-homogeneous formulation, In the above, the parameters consist three image points x i (in C 2 or P 2 , respectively) and three world points X i (in C 3 or P 3 .)Unlike the five-point problem, in P3P the world points are known-this makes the former a relative pose problem, and the latter an absolute pose problem.In addition to the unknwon rotations, translations, depths, and world points, there is an additional unknown vector n ∈ C 3×1 which defines the normal to the plane spanned by X 1 , . . ., X 3 .We include the normal in our formulation because it reduces the total degree of the problem's single deck transformation, which is given by The formulation (27) allows us to decompose the monomials into smaller classes, since its group of scalings is larger that of (26).Its continuous part is isomorphic to (C * ) 7 and is given by The discrete part has generators analogous to (24): Thus, G con ∼ = (C * ) 7 and   26)and (27).Both algorithms were able to successfully interpolate (28) using the a priori knowledge that this deck transformation is parameter-independent.As expected, running Algorithm 2 on the quasi-homogeneous formulation is more efficient.Additionally, Algorithm 2, unlike Algorithm 1, succeeds when run with the parameter-dependent setting.Intriguingly, the largest monomial class is the same in both parameter dependent and independent cases.
Once again, we find the importance of a carefully-chosen formulation can be key to recovering deck transformations.One initial experiment in which n was eliminated from the quasi-homogeneous formulation Equation (27), in which case the deck transformation has coordinate functions of total degree D * = 5.Although the sizes of Vandermonde matrices returned by Algorithm 2 appeared to be reasonable, we were not able to robustly interpolate these degree-5 multivariate rational functions.We speculate this is due to the well-known fact that large Vandermonde matrices are ill-conditioned [39,40].

Radial camera relative pose
Another problem in computer vision that has been recently tackled [41] is the problem of 3D reconstruction or relative pose from 4 images made by a radial camera.A radial camera may be understood as a linear map P : P 3 P 1 , thus given by a 2 × 4 matrix.The radial camera P associates a world point in P 3 with the radial line passing through the center of distortion in an image and the projection of the world point under the usual pinhole model (as for the five-point relative pose problem).The center of distortion may be assumed to be [0 : 0 : 1] ∈ P 2 , so that the equation of the radial line is parametrized by the projected image point [u : v : 1] ∈ P 2 as a direction vector, thus giving a point l = [u : v] ∈ P 1 .With these assumptions, a pinhole camera P pin : P 3 P 2 can be associated with a radial camera P as follows: Assuming that P pin is calibrated results in the radial camera matrix is the Cayley parametrization of the first 2 rows of a 3 × 3 camera rotation matrix and t ∈ C 2 are the first 2 elements of the camera translation vector.As explained in [41], it is enough to have 4 cameras and 13 world points to achieve finite number of solutions.The problem is then formulated by We may choose a world coordinate system by fixing (according to [41, Section 3.2]) (Here e 2 = 0 1 ⊤ .)We may eliminate the first two unknowns for every world point X i using the relation we obtain A formulation with 3584 complex solutions: Here X i denotes the third coordinate of X i .The enormous amount of solutions indicates that the problem has to be checked for decomposability (symmetry existence) in order to design for it more efficient solvers [41].
The continuous scaling symmetries of the formulation (33) are given by: where i = 1, . . ., 13, j = 1, . . ., 4 and hence G con ∼ = (C * ) 52 .The group of discrete scalings is isomorphic to Z 2 2 .As explained in [41], exploiting these symmetries and further structure of the branched cover allows for a more practical solution than naively tracking 3584 homotopy paths.The group of deck transformations of this problem is isomorphic to By running Algorithm 2 we were able to recover all 16 elements of Deck(f ).The details of this experiment are summarized in Table 3.As we can see from Table 3, running parameter-dependent version of Algorithm 2 (even though every deck transformation is parameter-independent) results into larger Vandermonde matrices, which in turn forces the algorithm to track more paths.We also notice that running monodromy and tracking paths is the bottleneck for this problem, which may be attributed to a large number of variables and solutions.
In the end, we recover the following four generators of Deck(f ), In [41], these formulas for Ψ i had to be worked out carefully by hand.Algorithm 2 furnishes an automatic derivation.We now turn our attention to Alt's problem.This is a classic problem of kinematic synthesis which was first solved using homotopy continuation in work of Morgan, Sommese, and Wampler [42].Several more recent works have used monodromy to verify their result, eg.[43,44,45].Here we explain how this problem can be modeled using a branched cover, and show how its well-known symmetry group can be recovered in our approach.
The formulation we use follows [42], employing the standard convention of isotropic coordinates.A vector in the plane is represented by two variables x, x ∈ C. For the purpose of solving polynomial systems, x and x are treated as independent complex variables; for any physically meaningful solutions, these coordinates will be related by complex conjugation.With this convention, angles T = e iθ are modeled by points on the hyperbola T T = 1.
In Figure 5, the vectors x and y point from the coupler point p 0 to the upper joints of the four bar, and vectors a and b point from P 0 towards the ground pivots.The four-bar mechanism has four revolute joints: two connecting the left "crank" and right "rocker" bars to the ground pivots, and another two connecting these bars to the base of the coupler triangle.The motion of the mechanism is induced by rotating the crank bar about its ground pivot.Atop the coupler triangle sits the coupler point (⋆), which traces out a curve as the mechanism moves.Without loss of generality, we may assume (0, 0) is a point on this curve.
Alt's problem can be stated as follows: given nine task positions p 0 = 0, p 1 , . . ., p 8 ∈ C, determine the mechanism parameters x, y, a, b and angles Q j , T j , S j such that the coupler point moves from p 0 to p i for i = 1, . . ., 8.
Here T j = e iλ j , S j = e iµ j as in Figure 5, and Q j = e iθ j gives the rotation of the coupler point (⋆) as it moves from p 0 to p j .
Referring to Figure 5, we may write down for each j = 1, . . ., 8 four loop-closure equations, and their conjugates.Consequently, the orientation of the coupler point may be written as a rational function in the mechanism parameters and the other angles, T j (a, b, x, y, Q j , S j ) = (y − x) −1 (S j (y The rocker angle S j is an algebraic function of degree 2 in the quantities x = (x, x, y, ȳ, a, ā, b, b) and the crank angle Q j .That is, for some A, B, C ∈ Q[x, Q j ].We note that for generic, fixed values of mechanism parameters x, this equation defines an elliptic curve in the affine plane of (Q j , S j ) ∈ C 2 .Since the discriminant of the quadratic (36) is square-free, we may define an irreducible variety X ′ = {(x, Q 1 , . . ., S 8 ) ∈ C 24 | (36), A(x, Q j ) ̸ = 0 hold, j = 1, . . ., 8}.
Using (34), each coupler point can now be expressed in terms of rational functions on X ′ , say p j (x, Q j , S j ), and pj (x, Qj , Sj ) for the conjugate.We may then take as an irreducible variety of problem-solution pairs X ⊂ C 8 × C 16 be the closed image of X ′ under the map (x, D, Q) → (x, p(x, Q, S), p(x, Q, S)).This gives a branched cover f : X → C 16 .
Although not yet formally proved, there is strong evidence that deg(f ) = 8652.Following the elimination strategy described in [42], we obtain a system of 8 equations We note that extending Ψ Rob to the eliminated variables {Q j , S j , T j } yields parameter-dependent coordinate functions.
The first numerical evidence that deg(f ) = 8652 was given in [42].Later on, the lower bound deg(f ) ≥ 8652 was certified by Hauenstein and Sottile using Smale's α-theory [46].A rigorous proof that this bound is tight remains an open problem.More recently, Sottile and Yahl have posed the problem of determining the Galois/monodromy group of the branched cover f [1, §7.3].From equations (37), we heuristically computed permutations in Mon(f ) using monodromy loops.This produced 4 permutations using default settings.Using Proposition 3.8, we determine that the deck transformation group is isomorphic to S 3 , generated by a transposition and 3-cycle corresponding to (38) and (39), respectively.
We remark that, in this formulation of Alt's problem, the methods of Section 5 detect that the group of scaling symmetries G con × G dis is trivial.Thus, Algorithm 2 yields no improvement over Algorithm 1 in this case.It would, however, be interesting to investigate these symmetries for other formulations of Alt's problem.As we can see from Table 4, running both parameter-independent and parameter-dependent formulations using Algorithm 1 yield similar results, where both successfully interpolate the deck transformations.For this example, 38 Algorithm 2 produces the same results as Algorithm 1 and thus is omitted from Table 4.

Conclusion
In summary, we have proposed a novel method for recovering hidden symmetries of commonly-occuring parametric polynomial systems.Despite its heuristic nature, our experiments demonstrate that the method is capable of delivering results, even on examples with a relatively larger number of solutions like Sections 6.3 and 6.4.Additionally, we have introduced a novel quasi-homgeneous interpolation framework in Section 5, delivering muchimproved results compared to the baselines in our conference paper [4].
One obvious avenue for further research is to test more examples and develop better heuristics.It would be highly of interest to develop more robust numerical interpolation methods, eg. by exploiting bases other than the standard monomials, and test them on examples similar to those in Section 6.There is also potential for fruitful contact with more traditional methods of symbolic computation.In addition to our comments in Section 2, we point out that some hybrid symbolic-numerical methods may be useful in practice.For instance, it seems plausible that one could (1) run Algorithm 1 or 2 until recovering coordinate functions for the deck transformations on some subset of variables y ⊂ x, then (2) eliminate the remaining variables x \ y and use parametric Gröbner bases to solve for their coordinate functions using the interpolated expressions from step (1).Such hybrid methods might also be useful for recovering deck transformations when a decomposition as in Definition 3.2 is already known, or vice-versa.
As noted in Section 3, decomposability and the existence deck transformations are closely related, but not equivalent.Indeed, the radial camera relative pose problem of Section 6.3 has a deck transformation group of order 16, which implies the degree 3584 branched cover associated to this problem decomposes as a composition of covers of degree 16 and 224.The latter cover, despite not having any deck transformations, does in fact decompose further into covers of degree 4, 2 and 28.Developing numerical methods for determining the maps and intermediate varieties appearing in such a decomposition is a highly appealing next step.

2 . 3 . 4 .
and f : X → Z given by f (a, b, c, d, x) = (a, b, c, d).The projection f is a decomposable branched cover in the sense of Definition 3.2.To see this, take Y = V(a(y 3 − 3y) + b(y 2 − 2) + cy + d) ⊂ C 5 , and define g : X Y by g(a, b, c, d, x) = (a, b, c, d, x 2 +1 x ), and h : Y → Z by h(a, b, c, d, y) = (a, b, c, d).The degrees of maps satisfy 6 = deg(f ) = deg(h • g) = deg(h) • deg(g) = 3 • Example The following example is based on [11, §2.3.2], and belongs to a general class of examples where decomposability can be detected via equations' Newton polytopes.Let Z = C 23 , and X ⊂ C 26 be the vanishing locus of the three equations below:
Definition 3.2.A branched cover f : X Z is said to be decomposable if there exist two branched covers g : X Y and h : Y 3.6 for detecting the existence of deck transformations.Proofs are in [12, §2.1].For the branched cover f of Example 3.4, the centralizer of its Galois/monodromy group S 4 ≀ S 8 in S 32 is trivial.Thus, this decomposable branched cover has no nontrivial deck transformations.

Table 1 :
Timings/details for recovering deck transformation of the five-point problem.

Table 2 :
Timings/details for recovering deck transformation of the P3P problem.

Table 3 :
Timings/details for recovering deck transformations of the radial relative pose problem.

Table 4 :
Timings/details for recovering deck transformations of the Alt's problem.