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[
GroebnerBasis`DistributedTermsList[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[GroebnerBasis`DistributedTermsList[{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 =
GroebnerBasis`DistributedTermsList[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 =
Flatten[{Total[abscoeffs], Thread[objcoeffs <= abscoeffs],
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 = GroebnerBasis`DistributedTermsList[syz[[All, 1]], vars][[1]];
posns = Union[Flatten[dtl[[All, All, 1]], 1]];
len = Length[posns];
cols = Dispatch[Thread[posns -> Range[len]]];
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}}]
--- edit ---
We now have ApproximatePolynomialGCD in the Wolfram Function Repository.
ResourceFunction["ApproximatePolynomialGCD"][
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]
(* Out[219]= 1. + 0.0147812 x0 - 0.121578 x1 *)
I used this as one of the examples in the WFR function.
--- end edit ---
