# Approximate GCD

I have several pairs of bivariate polynomials that I want to find if they have common factors. The polynomials, however, have numerical errors because the coefficients are some algebraic numbers that I cannot obtain very easily so I represent them by floating point numbers. Increasing the precision of the coefficient will not work. So one way I do in Mathematica is to just plot the two polynomials on the real plane and see if they have common components or intersect in a tangent. This visual inspection is not ideal since I have several of these polynomials. So I was wondering if one can do such approximation by some form of "epsilon"-GCD in Mathematica. Any ideas? Here is an example of a code I write to check if two such polynomials have common factors:

{f, g} = {x0^3 - 8.05690292794*x0^2*x1 + 66.8206602296*x0^2 -
2.42037371481*x0*x1^2 + 17.8341472094*x0*x1 - 55.3041827489*x0 +
8.52288118623*x1^3 - 66.7864893146*x1^2 - 35.9578468971*x1 +
71.4328040275,
-1.7128257326*x0^2*x1 - 0.309381807141*x0^2 +
13.9547615636*x0*x1^2 - 111.58197742*x0*x1 - 20.6540588108*x0 +
1.0984988335*x1^3 - 23.4487086797*x1^2 + 116.275997087*x1 +
18.7247500823};

ContourPlot[{f == 0, g == 0}, {x0, -100, 100}, {x1, -10, 10}, ContourStyle -> {Blue, Dashed, Thick, Green}]

PolynomialGCD[f, g]


You will see that the output of PolynomialGCD will be 1. I want to do this for a lot of pairs of polynomial without visually inspecting them. The plot will give me an idea though that the two polynomials represent curves that share a line

Any hints/suggestion is appreciated.

• You will get several results by googling "approximate GCD polynomial in several variables". I doubt that there is any implementation in Mathematica (maybe @DanielLichtblau can clarify?). Maybe you can find implementations in another language, e.g. in Maple, and translate it. May 30, 2018 at 10:13

I have a few methods, all different (I believe) from Zhonggang's. I'll show one that works well for the example at hand. It has some restrictions.

(1) Really not well suited for more than 3 variables.

(2) It wants to know the expected degree of the GCD. This can be finessed of course, e.g. by trying several and deciding which result you like.

(3) You need to provide an "operating precision" and this may need to be greater than that of the input. In this case it is not necessary though.

This is from work that will appear. The main ideas are explained here:

https://www.researchgate.net/publication/313159586_Approximate_polynomial_GCD_by_approximate_syzygies_redux

Here is the code. I have four sub-methods for a key part, and give the two that seem to work best.

drlMatrix[n_] :=
Prepend[Table[-KroneckerDelta[j + k - (n + 1)], {j, n - 1}, {k, n}],
ConstantArray[1, n]]
totalDegree[poly_, vars_] :=
Module[{t}, Max[0, Exponent[poly /. Thread[vars -> t*vars], t]]]

powerprods[_, n_ /; n <= 0] := {1};
powerprods[vars_, deg_] :=
Expand[(1 + Total[vars])^deg] /. Plus -> List /. a_Integer*b_ :> b

normAndSign[poly_, vars_, prec_] := Module[
{dtl, coeffs},
dtl = First[
GroebnerBasisDistributedTermsList[poly, vars,
CoefficientDomain -> InexactNumbers[Precision[poly]],
MonomialOrder -> DegreeReverseLexicographic]];
coeffs = Flatten[dtl[[All, 2]]];
If[coeffs === {}, coeffs = {0}];
{Sqrt[N[coeffs, prec].coeffs], Sign[First[coeffs]]}
]
polyNorm[poly_, vars_, prec_] := First[normAndSign[poly, vars, prec]]

normalizePoly[poly_, vars_, prec_] :=
With[{pnorm = normAndSign[poly, vars, prec]},
Expand[pnorm[[2]]*poly/First[pnorm]]]
normalizeSyz[syz_, vars_, prec_] := Module[
{cnorm1, cnorm2},
cnorm1 = polyNorm[syz[[2]], vars, prec];
cnorm2 = polyNorm[syz[[3]], vars, prec];
Expand[syz/Sqrt[cnorm1^2 + cnorm2^2]]
]

setCoefficientPrecision[a_?NumberQ, prec_] :=
If[Abs[a] < 10^(-prec), 0, SetPrecision[a, prec]]
setCoefficientPrecision[a_?NumberQ*b_?(! NumberQ[#] &), prec_] :=
setCoefficientPrecision[a, prec]*b
setCoefficientPrecision[(a_Plus | a_Times | a_List), prec_] :=
Map[setCoefficientPrecision[#, prec] &, a]
setCoefficientPrecision[a_, _] := a

getRelations[opoly1_, opoly2_, prec_, degmin_] := Catch[Module[
{poly1, poly2, vars, p2deg, mat, c, coords, newpolys, allvars, mgb,
degbnd, syz, degs, monoms},
vars = Variables[{opoly1, opoly2}];
poly1 = normalizePoly[opoly1, vars, prec];
poly2 = normalizePoly[opoly2, vars, prec];
mat = {{poly1, 1, 0}, {poly2, 0, 1}};
coords = Array[c, 3];
newpolys = mat.coords;
vars = Reverse[
Last[GroebnerBasisDistributedTermsList[{poly1, poly2}, vars,
MonomialOrder -> DegreeReverseLexicographic, Sort -> True]]];
allvars = Join[vars, RotateLeft[coords]];
degbnd = Join[ConstantArray[0, Length[vars]], {1, 1, 1}];
mgb = GroebnerBasis[newpolys, allvars, DegreeBound -> {degbnd, 1},
CoefficientDomain -> InexactNumbers[prec],
MonomialOrder -> drlMatrix[Length[allvars]]];
If[Head[mgb] === GroebnerBasis, Throw[$Failed]]; syz = mgb /. c[j_] :> UnitVector[3, j]; degs = Map[totalDegree[#, vars] &, syz[[All, 2]]]; p2deg = totalDegree[poly2, vars]; syz = Pick[syz, Map[p2deg - # >= degmin &, degs]]; degs = Map[totalDegree[#, vars] &, syz[[All, 2]]]; syz = Map[normalizeSyz[#, vars, prec] &, syz]; If[Length[vars] <= 3, syz = Flatten[Table[ monoms = powerprods[vars, (p2deg - degs[[j]] - degmin)]; Map[syz[[j]]*# &, monoms] , {j, Length[syz]}], 1] ]; {vars, Expand[syz]}]] getObjective[syz_, coords_, vars_] := Module[ {residuals, objcoeffs}, residuals = setCoefficientPrecision[syz[[All, 1]], MachinePrecision]; objcoeffs = GroebnerBasisDistributedTermsList[Expand[coords.residuals], vars, MonomialOrder -> DegreeReverseLexicographic][[1, All, 2]]; objcoeffs = Apply[Times, Tally[objcoeffs], {1}]; objcoeffs = Map[{# /. rr_?NumberQ*c[aa_] :> Sign[rr]*c[aa], #} &, objcoeffs]; objcoeffs = GatherBy[objcoeffs, First]; objcoeffs = Map[Total, objcoeffs[[All, All, 2]]]; objcoeffs ] recoverSyzygy[result_, coords_, syz_] := Catch[Module[{min, vals}, If[! ListQ[result] || Length[result] != 2, Throw[$Failed]];
{min, vals} = result;
Expand[(coords /. vals).syz]
]]

approxPolynomialSyzygyByLP[opoly1_, opoly2_, prec_, degmin_] := Module[
{vars, syz, objcons, c, coords, objcoeffs, a2, abscoeffs, result},
{vars, syz} = getRelations[opoly1, opoly2, prec, degmin];
coords = Array[c, Length[syz]];
objcoeffs = getObjective[syz, coords, vars];
abscoeffs = Array[a2, Length[objcoeffs]];
objcons =
Thread[-objcoeffs <= abscoeffs], Total[coords] == 1}];
SetOptions[LinearProgramming, Method -> "CLP",
Tolerance -> Automatic];
result =
FindMinimum[Evaluate[objcons], Join[coords, abscoeffs],
MaxIterations -> 10000];
recoverSyzygy[result, coords, syz]
]

approxPolynomialSyzygyByUnconstrainedOpt[opoly1_, opoly2_, prec_,
degmin_] := Module[
{vars, syz, dtl, posns, len, mat, cols, vec, result},
{vars, syz} = getRelations[opoly1, opoly2, prec, degmin];
dtl = GroebnerBasisDistributedTermsList[syz[[All, 1]], vars][[1]];
posns = Union[Flatten[dtl[[All, All, 1]], 1]];
len = Length[posns];
mat = Flatten[
MapIndexed[{#1[[1]] /. cols, #2[[1]]} -> N[#1[[2]]] &,
dtl, {2}]];
mat = SparseArray[mat, {len, Length[dtl[[All, All, 1]]]}];
vec = -mat[[All, 1]];
mat = mat[[All, 2 ;; -1]];
result = Prepend[LeastSquares[mat, vec], 1];
normalizeSyz[Expand[result.syz], vars, prec]
]

approxPolynomialGCD[opoly1_, opoly2_, prec_, degmin_,
syzFunction_: approxPolynomialSyzygyByLP] := Catch[Module[
{poly1, poly2, vars, vlen, syz, t, newdeg, newpoly, nextsyz,
result, residual},
vars = Variables[{opoly1, opoly2}];
poly1 = normalizePoly[opoly1, vars, prec + 2];
poly2 = normalizePoly[opoly2, vars, prec + 2];
vlen = Length[vars];
syz = syzFunction[poly1, poly2, prec, degmin];
If[syz === $Failed, Throw[$Failed]];
newpoly = setCoefficientPrecision[syz[[2]], prec + 2];
newdeg =
Exponent[newpoly /. Thread[vars -> RandomReal[1, vlen]*t], t];
nextsyz = syzFunction[newpoly, poly2, prec, newdeg];
result = setCoefficientPrecision[nextsyz[[2]], prec];
residual = (polyNorm[
normalizePoly[
Expand[setCoefficientPrecision[syz[[3]], prec]*result, prec],
vars, prec] - poly1, vars, prec] +
polyNorm[
normalizePoly[Expand[newpoly*result], vars, prec] - poly2,
vars, prec]);
{residual, result}
]]


The example:

{f, g} = {x0^3 - 8.05690292794*x0^2*x1 + 66.8206602296*x0^2 -
2.42037371481*x0*x1^2 + 17.8341472094*x0*x1 - 55.3041827489*x0 +
8.52288118623*x1^3 - 66.7864893146*x1^2 - 35.9578468971*x1 +
71.4328040275, -1.7128257326*x0^2*x1 - 0.309381807141*x0^2 +
13.9547615636*x0*x1^2 - 111.58197742*x0*x1 - 20.6540588108*x0 +
1.0984988335*x1^3 - 23.4487086797*x1^2 + 116.275997087*x1 +
18.7247500823};


Here we use the default method. We guess the degree is 1 since the plot looks like the two curves share a line. The result is of the form {residual, gcd} where "residual" is some measure of discrepancy from being a perfect gcd.

{res, agcd} = approxPolynomialGCD[f, g, MachinePrecision, 1]

(* Out[203]= {4.39798952826*10^-13,
0.709809760293 + 0.0104918057935 x0 - 0.0862970808038 x1} *)


The other method is in good agreement.

{res2, agcd2} =
approxPolynomialGCD[f, g, MachinePrecision, 1,
approxPolynomialSyzygyByUnconstrainedOpt]

(* Out[204]= {7.16133018154*10^-12, -0.709809760289 -
0.0104918057905 x0 + 0.0862970808035 x1} *)


Check how close these are to one another, and to the line in the prior response, by normalizing so constant term in each is unity.

h = -216.787714755602 - 3.20437211328366*x0 + 26.3565647925793*x1;

In[232]:= Expand[{agcd2/First[agcd2], agcd/First[agcd], h/First[h]}]

(* Out[232]=
{1. + 0.0147811517642 x0 - 0.121577760171 x1,
1. + 0.0147811517682 x0 - 0.121577760171 x1,
1. + 0.0147811517682 x0 - 0.121577760171 x1} *)


I plot these with some smallish offsetting in an effort to make the overlapping parts all appear.

ContourPlot[{agcd == -.02, f == -.04, g == 0}, {x0, -100,
100}, {x1, -10, 10},
ContourStyle -> {{Thickness[.005], Orange}, {Blue,
Dashed}, {Thickness[.003], Green}}]


I'll post this as a separate method (since it really is one). Comes from work found here (including the code shown below). The idea is to use a fairly standard method of computing the LCM of a pair of polynomials, using a Groebner basis. Then we allow some tolerancing up to some precision, and hope the GB computation handles that in a reasonable way. [Remark: While the method I already gave also uses a numeric GB, it is in no way subject to this tolerancing issue. An optimization step handles the combining of polynomials needed to get to the GCD.]

The code this time is quite short.

setCoefficientPrecision[a_?NumberQ, prec_] :=
If[Abs[a] < 10^(-prec), 0, SetPrecision[a, prec]]
setCoefficientPrecision[a_?NumberQ*b_?(! NumberQ[#] &), prec_] :=
setCoefficientPrecision[a, prec]*b
setCoefficientPrecision[(a_Plus | a_Times | a_List), prec_] :=
Map[setCoefficientPrecision[#, prec] &, a]
setCoefficientPrecision[a_, _] := a

floatPolynomialLCM[poly1_, poly2_, tol_] :=
Module[{vars, mat, v, cvars, newpolys, rels, gb, rul},
vars = Variables[{poly1, poly2}];
mat = {{1, 1, 1}, {poly1, 0, 0}, {0, poly2, 0}};
cvars = Array[v, 3];
newpolys = mat.cvars;
rels = Flatten[Union[Outer[Times, cvars, cvars]]];
newpolys = Join[newpolys, rels];
gb = GroebnerBasis[newpolys, Prepend[vars, Last[cvars]],
Most[cvars], MonomialOrder -> EliminationOrder, Tolerance -> tol,
CoefficientDomain -> InexactNumbers[Precision[newpolys]],
Sort -> True];
rul = Map[(# -> {}) &, rels];
gb = Flatten[gb /. rul];
First[gb] /. Last[cvars] -> 1]

floatPolynomialGCD[p1_, p2_, tol_] :=
Expand[PolynomialReduce[p1*p2, floatPolynomialLCM[p1, p2, tol],
CoefficientDomain -> InexactNumbers][[1, 1]]]


The example:

{f, g} = {x0^3 - 8.05690292794*x0^2*x1 + 66.8206602296*x0^2 -
2.42037371481*x0*x1^2 + 17.8341472094*x0*x1 - 55.3041827489*x0 +
8.52288118623*x1^3 - 66.7864893146*x1^2 - 35.9578468971*x1 +
71.4328040275, -1.7128257326*x0^2*x1 - 0.309381807141*x0^2 +
13.9547615636*x0*x1^2 - 111.58197742*x0*x1 - 20.6540588108*x0 +
1.0984988335*x1^3 - 23.4487086797*x1^2 + 116.275997087*x1 +
18.7247500823};


Find the approx GCD.

agcd3 = floatPolynomialGCD[f, g, 10^(-6)]

(* Out[77]= -115.879043771 - 1.7128257326 x0 + 14.088314591 x1 *)


Normalizing shows it is quite close to what was found in the other responses.

Expand[agcd3/First[agcd3]]

(* Out[78]= 1. + 0.014781151767 x0 - 0.121577760158 x1 *)


Here is an approach based (tweaked somewhat loosely) on AMVGCD in Gao et al. (2004). The tweaks are using Tolerance option to MatrixRank to determine the degree of the gcd and using Series to perform the approximate division.

The function gensys generates a linear system to represent $0 = u\,f + v\,g$ and returns the polynomials $u$ and $v$, their coefficients (as symbolic unknowns), and the coefficient matrix of the linear system. (The system $0 = u\,f + v\,g$ is the syzygy of Daniel's answer.) The argument deficiency determines the total degree of $u$ and $v$; it is subtracted from the degrees of $g$ and $f$ [sic], respectively. In the construction of ux and vx, there is a tacit assumption that f and g are bivariate polynomials. There is no need for this assumption, but it was simpler to code.

ClearAll[gensys];
gensys[f_, g_, deficiency_: 1, vars_: {u, v}] :=
Module[{df, dg, args, ux, vx, unknowns, mat, coeffs},
df = Max[Total@*First /@ CoefficientRules[f]] - deficiency; (* deg. of f *)
dg = Max[Total@*First /@ CoefficientRules[g]] - deficiency; (* deg. of g *)
args = Variables[{f, g}];
ux = Sum[vars[[1]][i, j] args[[1]]^i args[[2]]^j, {i, 0, dg}, {j, 0, dg - i}];
vx = Sum[vars[[2]][i, j] args[[1]]^i args[[2]]^j, {i, 0, df}, {j, 0, df - i}];
unknowns = Flatten@{
Table[vars[[1]][i, j], {i, 0, dg}, {j, 0, dg - i}],
Table[vars[[2]][i, j], {i, 0, df}, {j, 0, df - i}]};
coeffs = CoefficientArrays[ux f + vx g, args];  (* leads to the coeff. mat. *)
coeffs = coeffs /. sa_SparseArray :> sa["NonzeroValues"];
mat = Last@CoefficientArrays[Flatten@coeffs, unknowns];
{ux, vx, unknowns, mat}
];


There are two passes. The first approximates the degree of the gcd kk.

$tolerance = 10^-10; {ux, vx, unknowns, mat} = gensys[f, g]; kk = k /. First@Solve[ {k (k + 1)/2 == Length@unknowns - MatrixRank[mat, Tolerance ->$tolerance],
k >= 0}, k]
(*  1  *)


The degree kk should be checked. It may be possible that the rank of mat might not match a value of k (k + 1)/2 for an integer k. (Note this formula also assumes there are only two variables.) If kk == 0, then stop and the gcd is 1; otherwise, proceed to the next pass with kk as the deficiency for gensys. Now we should have exactly one singular value approximately equal to zero. The corresponding right singular vector gives an approximate solution to the linear system $0 = u\,f + v\,g$.

{ux, vx, unknowns, smat} = gensys[f, g, kk];
{left, ss, right} = SingularValueDecomposition[smat];
tmp = Reverse@{f, g}/{ux, vx} /. Thread[unknowns -> Last@Transpose@right];
gcds = With[{polys = Series[tmp, {x0, 0, 1}, {x1, 0, 1}] // Normal // Expand},
polys/polys[[All, 1]] // Expand // Chop
]

(*  {1. + 0.014781151768271154 x0 - 0.12157776017117022 x1,
1. + 0.014781151768222383 x0 - 0.12157776017103819 x1}   *)


The gcd is given by both f/vx and g/ux. You need only calculate one of them, but I showed both.

I just want to answer my own question which points out to an implementation that is not in Mathematica. Yes, I know there were approximate GCD algorithms and papers about them but I wanted an implementation that I could use since I don't think I will have the luxury of implementing myself for my purpose. I found one, unfortunately not in Mathematica but in Matlab. But I can parse the results to Mathematica to verify. The implementation is here: http://homepages.neiu.edu/~zzeng/naclab.html

Now from my last example I can do this in Matlab:

h = PolynomialGCD(f,g)


and this gives

h = -216.787714755602 - 3.20437211328366*x0 + 26.3565647925793*x1;


and this is exactly the coinciding the line in the figure in my original post (one can plot this in Mathematica to verify).

• If need a programmatic way to do that "within" Mathematica and if you have a MatLab installation (and a licence) at disposal, then you should try MatLink. May 30, 2018 at 15:17