Adjacent 2-by-2 Minors¶
This section documents the computation of the irreducible factors of a pure dimensional solution set, using the example of the 2-by-2 adjacent minors of a general 2-by-n matrix. The 2-by-2 minors of a 2-by-n matrix of indeterminates originates in algebraic statistics.
from phcpy.families import adjacent_minors
help(adjacent_minors)
shows
adjacent_minors(rows, cols)
Returns all adjacent 2-by-2 minors of a general
matrix of dimensions rows by cols.
This system originated in a paper on lattice walks and
primary decomposition, written by P. Diaconis, D. Eisenbud,
and B. Sturmfels, published by Birkhauser in 1998 in
Mathematical Essays in Honor of Gian-Carlo Rota,
edited by B. E. Sagan and R. P. Stanley,
volume 161 of Progress in Mathematics, pages 173--193.
See also the paper by S. Hosten and J. Shapiro on
Primary decomposition of lattice basis ideals, published in 2000
in the Journal of Symbolic Computation, volume 29, pages 625-639.
Let us look at the simplest nontrivial case: the adjacent 2-by-2 minors of a 2-by-3 matrix.
pols = adjacent_minors(2, 3)
for pol in pols:
print(pol)
shows
x_1_1*x_2_2-x_2_1*x_1_2;
x_1_2*x_2_3-x_2_2*x_1_3;
We have two polynomials in six variables. Therefore, we expect the solution set to be four dimensional. The two polynomials are quadrics. So, the degree of the solution set is expected to be four. The question is whether the four dimensional solution set of degree four is irreducible or not.
computing a witness set¶
A witness set of a positive dimensional solution set consists of
the original polynomial system, augmented with as many linear equations as the dimension of the solution set; and
generic points, as many as the degree of the solution set, computed as solutions of the augmented polynomial system.
The embedding adds extra slack variables, which are zero at the generic points.
The witness set data structure reduces the computation of a positive dimensional solution set to computing the isolated solutions of one embedded polynomial system.
from phcpy.sets import double_embed
epols = double_embed(6, 4, pols)
for pol in epols:
print(pol)
shows
+ x_1_1*x_2_2 - x_2_1*x_1_2 + (-9.99086911101846E-01-4.27240455127090E-02*i)*zz1 + (4.14818420957611E-01-9.09904213439104E-01*i)*zz2 + (-8.98599566975477E-01 + 4.38769664210604E-01*i)*zz3 + (-9.87637111749394E-01 + 1.56757569180295E-01*i)*zz4;
- x_2_2*x_1_3 + x_1_2*x_2_3 + (-9.50915060423189E-01 + 3.09452012209265E-01*i)*zz1 + (8.94934777859038E-01-4.46196978226426E-01*i)*zz2 + (4.28909177508030E-01-9.03347617171477E-01*i)*zz3 + (9.99997394966914E-01 + 2.28255545098161E-03*i)*zz4;
zz1;
zz2;
zz3;
zz4;
+ (2.52642772862563E-01-1.64562207779935E-01*i)*x_1_1 + (2.38172268997908E-01-1.84886617118381E-01*i)*x_2_2 + (-1.91482090565462E-01-2.32902769201594E-01*i)*x_2_1 + (7.30954525688645E-02 + 2.92516915276440E-01*i)*x_1_2 + (-3.84959957593650E-02-2.99043724594892E-01*i)*x_2_3 + (-2.71276842664000E-01-1.31597741406691E-01*i)*x_1_3 + (-2.64408104251273E-04-3.01511228642393E-01*i)*zz1 + (2.72300014364985E-01 + 1.29467343704581E-01*i)*zz2 + (-2.99939451765406E-01-3.07476207820803E-02*i)*zz3 + (-2.95299363009813E-01 + 6.08882346195858E-02*i)*zz4 - 3.01511344577764E-01;
+ (-4.58249807046168E-01-5.40906757392526E-02*i)*x_1_1 + (4.14862726733922E-02 + 4.11506360074024E-01*i)*x_2_2 + (-2.32018034572266E-01 + 9.48639573945220E-02*i)*x_2_1 + (-2.80184386911322E-01-4.28818448828079E-01*i)*x_1_2 + (-1.20422832345040E-01-1.42772536416847E-01*i)*x_2_3 + (-1.75000947987830E-01-2.52077630538672E-02*i)*x_1_3 + (1.00853655248805E-01-1.59726738098879E-01*i)*zz1 + (1.47157089788836E-01 + 9.31710575065957E-02*i)*zz2 + (-3.18937022210321E-02-2.67242905524245E-01*i)*zz3 + (-1.70655540108348E-01 + 3.05332146451031E-02*i)*zz4+(-9.36346290234469E-02 + 2.17662695080044E-01*i);
+ (-8.19176308116668E-02 + 9.31779891695395E-02*i)*x_1_1 + (-2.23021019174348E-01 + 9.08742650844521E-02*i)*x_2_2 + (-2.55417845632647E-01 + 1.92893995983262E-01*i)*x_2_1 + (2.85202344419800E-02 + 4.06697473728353E-01*i)*x_1_2 + (-3.26980478298163E-01 + 5.39551406506719E-02*i)*x_2_3 + (-4.72487633926490E-03 + 3.03185406062837E-01*i)*x_1_3 + (2.49469923033290E-01 + 7.22359362814556E-02*i)*zz1 + (3.57054498446798E-01 + 1.84302643405729E-01*i)*zz2 + (2.44627852456931E-01-2.19409644826369E-02*i)*zz3 + (-6.45211959942027E-02-2.81391298463502E-01*i)*zz4+(1.78511770624327E-02-2.88585493131170E-01*i);
+ (-2.47934564140415E-01 + 3.78474051361417E-01*i)*x_1_1 + (2.57696009847760E-01-1.97353959339253E-01*i)*x_2_2 + (-1.32713498879035E-01-6.68083554253052E-02*i)*x_2_1 + (1.71277270038441E-01-1.69671423207279E-01*i)*x_1_2 + (2.45498609773496E-01-3.00565440004147E-01*i)*x_2_3 + (-1.33200278933953E-01 + 9.57961946142634E-02*i)*x_1_3 + (2.53856038005847E-01 + 2.28882119288700E-01*i)*zz1 + (3.08364764282092E-03-3.24002619621010E-01*i)*zz2 + (2.00774472927742E-01 + 2.60472818035180E-01*i)*zz3 + (-2.51374429242122E-01 + 8.28490424271253E-03*i)*zz4+(-7.60100938321149E-02-1.82187404809906E-01*i);
Now we compute the second part of the witness set, using the blackbox solver.
from phcpy.solver import solve
esols = solve(epols)
for (idx, sol) in enumerate(esols):
print('Solution', idx+1, ':')
print(sol)
shows four generic points on the four dimensional solution set:
Solution 1 :
t : 1.00000000000000E+00 0.00000000000000E+00
m : 1
the solution for t :
x_1_1 : 1.26098995657825E+00 1.23960253907080E-01
x_2_2 : 6.08716792661783E-02 -3.72692628829057E-01
x_2_1 : -1.49366498111755E+00 -1.72487935488432E+00
x_1_2 : 1.17926528430272E-01 1.73403649424959E-01
zz1 : 2.01661798954576E-33 -1.36621089687380E-31
zz2 : 4.18934314423693E-32 1.19013821639201E-31
zz3 : 9.09330358138888E-33 -3.97471598577921E-32
zz4 : 9.41147475091108E-32 -5.25062461176711E-32
x_1_3 : -5.43007038294243E-01 4.77552159454532E-01
x_2_3 : 1.30126856409899E+00 4.91781165629461E-02
== err : 4.606E-15 = rco : 1.659E-02 = res : 1.221E-15 =
Solution 2 :
t : 1.00000000000000E+00 0.00000000000000E+00
m : 1
the solution for t :
x_1_1 : 7.84232277025150E-01 2.84185732210647E-02
x_2_2 : 1.09164662176987E-01 -6.94770373095709E-01
x_2_1 : 3.30315020994409E-02 -9.45083305517748E-01
x_1_2 : 5.76431528154133E-01 9.13299754350158E-02
zz1 : -1.29868775020467E-32 -1.23245957744023E-32
zz2 : -1.02927640426824E-32 -1.56443669914841E-32
zz3 : -1.96403692130211E-32 -7.02597925370710E-33
zz4 : 0.00000000000000E+00 0.00000000000000E+00
x_1_3 : -6.32506552290242E-01 1.49836349424459E-01
x_2_3 : 1.81539704759152E-01 7.61970172619267E-01
== err : 1.716E-15 = rco : 2.536E-02 = res : 1.499E-15 =
Solution 3 :
t : 1.00000000000000E+00 0.00000000000000E+00
m : 1
the solution for t :
x_1_1 : -1.07029652224429E+00 -2.84414365427269E-02
x_2_2 : 4.26584686826893E-01 4.61346027545957E-01
x_2_1 : -8.74524024788936E-02 -7.16513058512809E-01
x_1_2 : 7.70137841287209E-01 -5.24903704207380E-01
zz1 : -2.60172998762874E-31 -7.95558581219627E-32
zz2 : 1.94241870643223E-31 -7.60476764044638E-32
zz3 : -1.18161218057147E-31 6.47719980338975E-32
zz4 : 3.21269292411328E-31 9.43964867510126E-32
x_1_3 : 1.88771487997930E+00 2.52881775952061E+00
x_2_3 : -1.49855102312210E+00 1.51018382382141E+00
== err : 5.217E-15 = rco : 1.384E-02 = res : 1.749E-15 =
Solution 4 :
t : 1.00000000000000E+00 0.00000000000000E+00
m : 1
the solution for t :
x_1_1 : 1.03149817185821E+00 1.48840565152077E-01
x_2_2 : -1.97069785649029E-31 -4.19546322472859E-32
x_2_1 : -1.96423581053286E+00 -2.24826401526688E+00
x_1_2 : -6.05197132478501E-32 6.74819240986552E-32
zz1 : -4.35083051092473E-32 8.73329103580651E-32
zz2 : 6.61155389590362E-32 7.61859548953403E-32
zz3 : 2.10093447239302E-32 -1.28492740537369E-32
zz4 : 0.00000000000000E+00 0.00000000000000E+00
x_1_3 : -1.40341064895811E-01 1.23712580175709E+00
x_2_3 : 1.45874399778069E+00 6.42372628489355E-02\
== err : 4.557E-15 = rco : 1.108E-02 = res : 1.332E-15 =
As expected, we find four solutions, equal to the degree of the solution set.
monodromy breakup¶
The numerical irreducible decomposition of a pure dimensional solution set is a list of tuples, each tuple represents an irreducible component, with two elements
the list of labels to the generic points in the witness set; and
the certificate of the linear trace.
The monodromy breakup algorithm refines the witness set, partitioning the generic points in the witness set corresponding to the irreducible components. The stop test in the monodromy looping algorithm is provided by the linear trace, which serves as a certificate for the numerical computations.
The breakup algorithm runs in verbose mode:
from phcpy.factor import double_monodromy_breakup
deco = double_monodromy_breakup(epols, esols, 4, verbose=True)
In verbose mode, the progress of the algorithm is printed:
... running monodromy loops in double precision ...
... initializing the grid for the linear trace ...
The diagnostics of the trace grid :
largest error on the samples : 1.056863766729658e-14
smallest distance between the samples : 1.1622396157009902
... starting loop 1 ...
new permutation : [2, 1, 3, 4]
number of factors : 4 -> 3
the decomposition :
factor 1 : ([1, 2], 0.2775442892800747)
factor 2 : ([3], 0.2775442892800704)
factor 3 : ([4], 9.992007221626409e-16)
the permutation : 2 1 3 4 : 4 -> 3
calculated sum at samples : -6.39839733359069E-02 1.09505498482705E-01
value at the linear trace : 1.07936961323818E-02 -9.32611213290812E-02
calculated sum at samples : -1.22710696579410E-01 2.89111792077992E-02
value at the linear trace : -1.97488366047695E-01 2.31677799019584E-01
calculated sum at samples : 2.68414486121697E-01 -5.27451814780890E-01
value at the linear trace : 2.68414486121697E-01 -5.27451814780890E-01
Certifying with linear trace test...
calculated sum at samples : -6.39839733359069E-02 1.09505498482705E-01
value at the linear trace : 1.07936961323818E-02 -9.32611213290812E-02
The witness points 1 2 do not define a factor.
The factorization cannot be certified.
... starting loop 2 ...
new permutation : [1, 2, 3, 4]
number of factors : 3 -> 3
the decomposition :
factor 1 : ([1, 2], 0.2775442892800747)
factor 2 : ([3], 0.2775442892800704)
factor 3 : ([4], 9.992007221626409e-16)
... starting loop 3 ...
the permutation : 1 2 3 4 : 3 -> 3
calculated sum at samples : -6.39839733359069E-02 1.09505498482705E-01
value at the linear trace : 1.07936961323818E-02 -9.32611213290812E-02
calculated sum at samples : -1.22710696579410E-01 2.89111792077992E-02
value at the linear trace : -1.97488366047695E-01 2.31677799019584E-01
calculated sum at samples : 2.68414486121697E-01 -5.27451814780890E-01
value at the linear trace : 2.68414486121697E-01 -5.27451814780890E-01
Certifying with linear trace test...
calculated sum at samples : -6.39839733359069E-02 1.09505498482705E-01
value at the linear trace : 1.07936961323818E-02 -9.32611213290812E-02
The witness points 1 2 do not define a factor.
The factorization cannot be certified.
new permutation : [2, 3, 1, 4]
number of factors : 3 -> 2
the decomposition :
factor 1 : ([1, 2, 3], 4.163336342344337e-15)
factor 2 : ([4], 9.992007221626409e-16)
the permutation : 2 3 1 4 : 3 -> 2
calculated sum at samples : -1.86694669915317E-01 1.38416677690504E-01
value at the linear trace : -1.86694669915313E-01 1.38416677690503E-01
calculated sum at samples : 2.68414486121697E-01 -5.27451814780890E-01
value at the linear trace : 2.68414486121697E-01 -5.27451814780890E-01
Certifying with linear trace test...
calculated sum at samples : -1.86694669915317E-01 1.38416677690504E-01
value at the linear trace : -1.86694669915313E-01 1.38416677690503E-01
The witness points 1 2 3 define a factor.
calculated sum at samples : 2.68414486121697E-01 -5.27451814780890E-01
value at the linear trace : 2.68414486121697E-01 -5.27451814780890E-01
The witness points 4 define a factor.
The factorization is certified.
calculated sum at samples : -1.86694669915317E-01 1.38416677690504E-01
value at the linear trace : -1.86694669915313E-01 1.38416677690503E-01
calculated sum at samples : 2.68414486121697E-01 -5.27451814780890E-01
value at the linear trace : 2.68414486121697E-01 -5.27451814780890E-01
As a summary, the contents of deco is written as
from phcpy.factor import write_decomposition
write_decomposition(deco)
which produces
factor 1 : ([1, 2, 3], 4.163336342344337e-15)
factor 2 : ([4], 9.992007221626409e-16)
There are two irreducible factors, one of degree three, and another of degree one. The floating-point certificates are close to machine precision.