Automatic Generator of Minimal Problem Solvers
Zuzana Kukelova
1
, Martin Bujnak
1,2
and Tomas Pajdla
1
1
Center for Machine Perception
Czech Technical University, Prague
2
Microsoft Corporation
kukelova,bujnam1,pajdla@cmp.felk.cvut.cz
Abstract. Finding solutions to minimal problems for estimating epipolar geom-
etry and camera motion leads to solving systems of algebraic equations. Often,
these systems are not trivial and therefore special algorithms have to be designed
to achieve numerical robustness and computational efficiency. The state of the
art approach for constructing such algorithms is the Gr
¨
obner basis method for
solving systems of polynomial equations. Previously, the Gr
¨
obner basis solvers
were designed ad hoc for concrete problems and they could not be easily applied
to new problems. In this paper we propose an automatic procedure for gener-
ating Gr
¨
obner basis solvers which could be used even by non-experts to solve
technical problems. The input to our solver generator is a system of polynomial
equations with a finite number of solutions. The output of our solver generator
is the Matlab or C code which computes solutions to this system for concrete
coefficients. Generating solvers automatically opens possibilities to solve more
complicated problems which could not be handled manually or solving existing
problems in a better and more efficient way. We demonstrate that our automatic
generator constructs efficient and numerically stable solvers which are compara-
ble or outperform known manually constructed solvers. The automatic generator
is available at http://cmp.felk.cvut.cz/minimal
3
.
1 Introduction
Many problems can be formulated using systems of algebraic equations. Examples are
the minimal problems in computer vision, i.e. problems solved from a minimal num-
ber of point correspondences, such as the five point relative pose problem [20], the
six point focal length problem [18], six point generalized camera problem [19], the
nine point problem for estimating para-catadioptric fundamental matrices [9], the eight
point problem for estimating fundamental matrix and single radial distortion parameter
for uncalibrated cameras [11], the six point problem for estimating essential matrix and
single radial distortion parameter for calibrated cameras [12, 4], the nine point prob-
lem for estimating fundamental matrix and two different radial distortion parameters
for uncalibrated cameras [12, 4]. These are important problems with a broad range of
applications [10].
3
This work has been supported by EC projects FP6-IST-027787 DIRAC and MRTN-CT-2004-
005439 VISIONTRAIN and grants MSM6840770038DMCMIII, STINT Dur IG2003-2 062
and MSMT KONTAKT 9-06-17.
Often, polynomial systems which arise are not trivial. They consist of many polyno-
mial equations in many unknowns and of higher degree. Therefore, special algorithms
have to be designed to achieve numerical robustness and computational efficiency. The
state of the art method for constructing such algorithms, the solvers, is the Gr
¨
obner basis
method for solving systems of polynomial equations. It was used to solve all previously
mentioned computer vision problems.
The Gr
¨
obner basis solvers in computer vision are mostly designed for concrete prob-
lems and in general consist of three key steps. In the first step, the problem is solved
using a computer algebra system, e.g. Macaulay 2 or Maple, for several (many) random
coefficients from a finite field. This is done using a general technique for solving poly-
nomial equations by finding a Gr
¨
obner basis for original equations [6]. In this phase,
the basic parameters of the problem are identified, such as whether there exists a finite
number of solutions for general coefficients and how “hard” is to obtain the Gr
¨
obner
basis, and which monomials it contains. This procedure is automatic and relies on gen-
eral algorithms of algebraic geometry such as Buchberger [1] or F4 [8] algorithm. The
computations are carried out in a finite field in order to avoid expensive growth of co-
efficients and to avoid numerical instability.
In the second step, a special elimination procedure, often called “elimination tem-
plate” is generated. This elimination template says which polynomials from the ideal
should be added to the initial polynomial equations to obtain a Gr
¨
obner basis or at least
all polynomials needed for constructing a special matrix, called “action matrix”, and
thus solving the initial equations. The goal of this step is to obtain a computationally
efficient and numerically robust procedure. Until now, this step has been mainly man-
ual requiring to trace the path of elimination for random coefficient form the finite field,
checking redundancies and possible numerical pitfalls, and writing down a program in
a procedural language such as Matlab or C.
In the last step, the action matrix is constructed from the resulting equations and
the solutions to the original problem are obtained numerically as the eigenvalues or
eigenvectors of the action matrix [7].
The first and the third step are standard and well understood. It is the second step
which still involves considerable amount of craft and which makes the process of solver
generating rather complex and virtually impenetrable for a non-specialist. Moreover, for
some problems it need not be clear how their elimination templates were generated and
therefore non-specialists often use them as black boxes, are not able to reimplement
them, improve them or create similar solvers for their own new problems.
In this paper we present an automatic generator of the elimination templates for
Gr
¨
obner basis solvers to an interesting class of systems of polynomial equations which
appear in computer vision problems. We have to accept that there is no hope to obtain an
efficient and robust solver for completely general systems of polynomial equations since
the general problem is known to be NP-complete and EXPSPAPCE-complete [15]. On
the other hand, all known computer vision problems share the property that the elimi-
nation path associated with their solution is the same for all interesting configurations
of coefficients of a given problem.
Consider, for instance, the classical 5 point problem for estimating the essential
matrix from 5 point correspondences [16, 20]. In general, the elimination path used to
solve the 5 point problem depends on actual image coordinates measured. Fortunately,
for non-degenerate image correspondences, i.e. for those which lead to a finite number
of essential matrices, the elimination path is always the same. Therefore, it is enough
to find the path for one particular non-degenerate configuration of coefficients and then
use the same path for all non-degenerate configuration of coefficients. The paths for
degenerate configurations of coefficients may be very different and there may be very
many of them but we need not consider them.
We propose an automatic generator that finds one elimination path. It can find any
path. The choice of the path is controlled by the particular coefficients we choose to
generate the elimination template. We demonstrate that our automatic generator con-
structs efficient and numerically stable solvers which are comparable or outperform
known manually constructed solvers.
The input into our automatic generator is the system of polynomial equations which
we want to solve with a particular choice of coefficients from Z
p
that choose the par-
ticular elimination path. For many problems, the interesting “regular” solutions can be
obtained with almost any random choice of coefficients. Therefore, we use random co-
efficients from Z
p
. The output from the generator is the Matlab or C code solver that
returns solutions to this system of polynomial equations for concrete coefficients from
R. In online computations, only this generated solver is called.
In the next two sections we review the Gr
¨
obner basis method for solving systems
of polynomial equations and the solvers based on this method. Section 5 is dedicated
to our automatic procedure for generating Gr
¨
obner basis solvers. Then, we demonstrate
the results of our automatic generator on some minimal problems.
2 Solving systems of polynomial equations
Our goal is to solve a system of algebraic equations
f
1
(x) = ... = f
m
(x) = 0 (1)
which are given by a set of m polynomials F = {f
1
, ..., f
m
| f
i
∈ C [x
1
, ..., x
n
]} in n
variables x = (x
1
, ..., x
n
) over the field of complex numbers.
Such system of algebraic equations can be written in a matrix form
M X = 0, (2)
where X is a vector of all monomials which appear in these equations and M is a coef-
ficient matrix. In the next we mostly consider this matrix representation of systems of
equations and, for example, by Gauss-Jordan (G-J) elimination of equations we under-
stand G-J elimination of the corresponding coefficient matrix M.
Solving systems of algebraic polynomial equations is a very challenging problem.
There doesn’t exist one robust, numerically stable and efficient method for solving such
systems in general case. Therefore, special algorithms have to be designed for specific
problems.
The general Gr
¨
obner basis method for solving systems of polynomial equations can
be quite inefficient in some cases but it was recently used successfully as the basis for
efficient solvers of computer vision minimal problems. We now describe this method
and the solvers based on this method.
3 Gr
¨
obner basis method
The polynomials F = {f
1
, ..., f
m
| f
i
∈ C [x
1
, ..., x
n
]} define ideal I, which is the set of
all polynomials that can be generated as polynomial combinations of initial polynomials
F
I = {Σ
m
i=1
f
i
p
i
| p
i
∈ C [x
1
, ..., x
n
]} , (3)
where p
i
are arbitrary polynomials from C [x
1
, ..., x
n
].
We can define division by an ideal I in C [x
1
, ..., x
n
] as the division by the set F of
generators of I. There are special sets of generators, Gr
¨
obner bases, of the ideal I, for
which this division by the ideal I is well defined in the sense that the remainder on the
division doesn’t depend on the ordering of the polynomials in the Gr
¨
obner basis G. This
means that the remainder of an arbitrary polynomial f ∈ C [x
1
, ..., x
n
] on the division
by G in a given monomial ordering is uniquely determined. Furthermore, f ∈ I if and
only if the reminder of f on the division by G is zero (f
G
= 0). This implies that
f
G
+ g
G
= f + g
G
and f
G
g
G
G
= f g
G
Gr
¨
oebner bases generate the same ideal as the initial polynomial equations and
therefore have the same solutions. However, it is important that they are often easier
to solve (e.g. the reduced Gr
¨
obner basis w.r.t. the lexicographic ordering contains poly-
nomial in one variable only). Computing such basis and “reading off” the solutions
from it is one standard method for solving systems of polynomial equations. Although
this sounds really nice and easy the reality is much harder. The problem is that the com-
putation of Gr
¨
obner bases is in general an EXPSPACE-complete problem, i.e. that large
space is necessary for storing intermediate results. Fortunately in many specific cases,
the solution to a system of polynomial equations computed using Gr
¨
obner bases can be
obtained much faster.
For solving systems of polynomial equations, the most suitable ordering is the lex-
icographic one which results in a system of equations in a “triangular form” with one
equation in one variable only. Unfortunately, computation of such Gr
¨
obner basis w.r.t.
the lexicographic ordering is very time consuming and for most of the problems can
not be used. Therefore in many cases a Gr
¨
obner basis G under another ordering, e.g.
the graded reverse lexicographical ordering (grevlex), which is often easier to compute,
is constructed. Then, other properties of this basis are used to obtain solutions to the
initial system of polynomial equations.
Thanks to the property that the division by the ideal I is well defined for the Gr
¨
obner
basis G, we can consider the space of all possible remainders on the division by I. This
space is know as a quotient ring and we will denote it as A = C [x
1
, ..., x
n
] /I. It is
known that if I is a radical ideal [6] and the set of equations F has a finite number of
solutions N , then A is a finite dimensional space with dim(A) = N. Now we can use
nice properties of a special action matrix defined in this space, to find solutions to our
system of equations (1).
Consider the multiplication by some polynomial f ∈ C [x
1
, ..., x
n
] in the quotient
ring A. This multiplication defines a linear mapping T
f
from A to itself. Since A is
a finite-dimensional vector space over C, we can represent this mapping by its matrix
M
f
with respect to some monomial basis B =
n
x
α
|x
α
G
= x
α
o
of A, where x
α
is
a monomial x
α
= x
α
1
1
x
α
2
2
...x
α
n
n
and x
α
G
is the reminder of x
α
on the division by
G. The action matrix M
f
can be viewed as a generalization of the companion matrix
used in solving one polynomial equation in one unknown. It is because solutions to our
system of polynomial equations (1) can be easily obtained from the eigenvalues and
eigenvectors of this action matrix [7].
4 Gr
¨
obner basis solver
The Gr
¨
obner basis method for solving systems of polynomial equations based on ac-
tion matrices was recently used to solve many minimal problems in computer vision.
The solvers to all these minimal problems are very similar and are based on the same
concepts. Many minimal problems in computer vision, including all mentioned above,
have the convenient property that the monomials, appearing in the set of initial gen-
erators F are always same irrespectively from the concrete coefficients arising from
non-degenerate image measurements. Therefore, the leading monomials of the corre-
sponding Gr
¨
obner basis, and thus the monomials in the basis B of the quotient ring A,
are generally the same and can be found once in advance. This is an important observa-
tion which helps to solve all these problems efficiently.
The first step of these solvers is to analyze the particular problem, i.e. whether it
is solvable and how many solutions there are, in a randomly chosen finite prime field
Z
p
(Z/ hpi) with p >> 7. Coefficients in Z
p
can be can be represented accurately and
efficiently. It speeds up computations, minimizes memory requirements and especially
avoids numerical instability. Computing with floating point approximations of the co-
efficients may lead to numerical instability since it may be difficult to determine when
the coefficients become zero.
Next, the Gr
¨
obner basis G, and the basis B are found for many random coefficients
from Z
p
. Thanks to algebraic geometry theorem [21], we know that if the bases G
and B remain stable for many different random coefficients, i.e. if B consists of the
same monomials, they are generically equivalent to the bases of the original system of
polynomial equations with rational coefficients.
With this information in hand the solver can be created. The solver typically con-
sists of hand made elimination templates [18, 19, 9, 11] (or one template [2–5]) that
determine which polynomials from the ideal I should be added to the initial equations
to obtain the Gr
¨
obner basis G or at least all polynomials needed for constructing the
action matrix M
f
. These elimination templates are the crucial part of the solver.
An important observation has been made in [11]. It was observed that the action
matrix can be constructed without computing a complete Gr
¨
obner basis G. All we need
for creating the action matrix M
f
is to construct polynomials from the ideal I with
leading monomials from the set (f · B)\ B and the remaining monomials from B. This
fact was in some way implicitly used in previous solvers [18, 19, 9] but hasn’t been
fully articulated.
Consider that the basis B =
©
x
α(1)
, . . . , x
α(N)
ª
of A, which has been found once
in advance by computations in Z
p
. Then, the polynomials we need for constructing the
action M
f
matrix are of the form
q
i
= f x
α(i)
+ h
i
, (4)
![Fig. 3. The log10 relative errors of the radial distortion parameters λ1 (Left) and λ2 (Right) for 1000 runs of three different solvers. Original solver with the basis selection [4] (Red, darker), original solver without the basis selection (Blue, dotted) and our generated solver (Green, lighter) on the synthetic dataset without noise .](/figures/fig-3-the-log10-relative-errors-of-the-radial-distortion-24lksd4h.webp)
![Fig. 2. The log10 relative errors of the focal length for 10000 runs of two solvers. Original solver [18] (Red, darker) and our generated solver (Green, lighter) on the synthetic dataset without noise.](/figures/fig-2-the-log10-relative-errors-of-the-focal-length-for-2vhagw2s.webp)
