Minimal polynomial for the Wasserstein distance degree

This notebook contains the code explained in Seciton 4.2 of the paper “The algebraic degree of the Wasserstein distance” by Chiara Meroni, Bernhard Reinke, and Kexin Wang.

using Oscar
using HomotopyContinuation

Welcome to Nemo version 0.43.3

Nemo comes with absolutely no warranty whatsoever
  ___   ____   ____    _    ____
 / _ \ / ___| / ___|  / \  |  _ \   |  Combining ANTIC, GAP, Polymake, Singular
| | | |\___ \| |     / _ \ | |_) |  |  Type "?Oscar" for more information
| |_| | ___) | |___ / ___ \|  _ <   |  Manual: https://docs.oscar-system.org
 \___/ |____/ \____/_/   \_\_| \_\  |  Version 1.0.2

We implement here the algorithm for cubic polynomials. An analogous procedure can be used for polynomials of other degrees. For linear or quadratic polynomials, finding the minimal polynomial is quite straightforward. For degree \(4\) or higher, the elimination step might not terminate in reasonable time.

Cubic polynomials

d = 3;

The Wasserstein distance degree of two univariate polynomials \(p, q\) is the algebraic degree of the Wasserstein distance between \(p\) and \(q\). In degree three, the Wasserstein distance degree depends only on the number of their real and complex roots.

We define the functions that provide the minimal polynomial of WDdeg(\(p,q\)). Each of them constructs an ideal encoding the coefficients of the polynomials, and it eliminates an auxiliary variable. The output is a multiple of the minimal polynomial of the Wasserstein distance of \(p, q\). For sufficiently generic \(p, q\), the output is irreducible, hence it is the actual minimal polynomial.

function minimal_poly_33(coeff_p, coeff_q)
    ## Input:
        ## coeff_p = list of coefficients of p, polynomial with 3 real roots
        ## coeff_q = list of coefficients of q, polynomial with 3 real roots
    ## Output:
        ## a list of primary ideals among which there is the minimal polynomial of the Wasserstein distance of p,q
    R, a, b, z = polynomial_ring(QQ, "a" => 1:d, "b" => 1:d, "z" =>1:1)
    # a[1], a[2], a[3] denote the real roots of p
    # b[1], b[2], b[3] denote the real roots of q
    S, t = polynomial_ring(R, "t")
    pp = prod(t-a[i] for i=1:d)
    coeff_pp = collect(Oscar.coefficients(pp))
    qq = prod(t-b[i] for i=1:d)
    coeff_qq = collect(Oscar.coefficients(qq))
    I = ideal(R, [[coeff_p[i]-coeff_pp[i] for i=1:d+1];
                  [coeff_q[i]-coeff_qq[i] for i=1:d+1];
                   z[1] - sum([(a[i]-b[i])^2 for i=1:d])])
    J = eliminate(I, [a;b])
    return primary_decomposition(J)
end

minimal_poly_33 (generic function with 1 method)

function minimal_poly_13(coeff_p, coeff_q)
    ## Input:
        ## coeff_p = list of coefficients of p, polynomial with 1 real root
        ## coeff_q = list of coefficients of q, polynomial with 3 real roots
    ## Output:
        ## a list of primary ideals among which there is the minimal polynomial of the Wasserstein distance of p,q
    R, a, b, z = polynomial_ring(QQ, "a" => 1:d, "b" => 1:d, "z" =>1:1)
    # a[1], a[2]+ia[3], a[2]-ia[3] denote the roots of p
    # b[1], b[2], b[3] denote the real roots of q
    I = ideal(R, [a[1]+2*a[2]+coeff_p[3]//coeff_p[4], 2*a[1]*a[2]+a[2]^2+a[3]^2-coeff_p[2]//coeff_p[4], a[1]*(a[2]^2+a[3]^2)+coeff_p[1]//coeff_p[4],
                  b[1]+b[2]+b[3]+coeff_q[3]//coeff_q[4], b[1]*b[2]+b[2]*b[3]+b[3]*b[1]-coeff_q[2]//coeff_q[4], b[1]*b[2]*b[3]+coeff_q[1]//coeff_q[4],
                  z[1] - ((a[1]-b[1])^2 + (a[2]-b[2])^2 + (a[2]-b[3])^2 + 2*a[3]^2)])
    J = eliminate(I, [a;b])
    return primary_decomposition(J)
end

minimal_poly_13 (generic function with 1 method)

function minimal_poly_11(coeff_p, coeff_q)
    ## Input:
        ## coeff_p = list of coefficients of p, polynomial with 1 real root
        ## coeff_q = list of coefficients of q, polynomial with 1 real root
    ## Output:
        ## a list of primary ideals among which there is the minimal polynomial of the Wasserstein distance of p,q
    R, a, b, z = polynomial_ring(QQ, "a" => 1:d, "b" => 1:d, "z" =>1:1)
    # a[1], a[2]+ia[3], a[2]-ia[3] denote the roots of p
    # b[1], b[2]+ib[3], b[2]-ib[3] denote the roots of q
    # here the associated graphs Gamma have different cycle structures, hence we have multiple options for the minimal polynomial
    I1 = ideal(R, [a[1]+2*a[2]+coeff_p[3]//coeff_p[4], 2*a[1]*a[2]+a[2]^2+a[3]^2-coeff_p[2]//coeff_p[4], a[1]*(a[2]^2+a[3]^2)+coeff_p[1]//coeff_p[4],
                  b[1]+2*b[2]+coeff_q[3]//coeff_q[4], 2*b[1]*b[2]+b[2]^2+b[3]^2-coeff_q[2]//coeff_q[4], b[1]*(b[2]^2+b[3]^2)+coeff_q[1]//coeff_q[4],
                  z[1] - ((a[1]-b[1])^2 + 2*(a[2]-b[2])^2 + 2*(a[3]-b[3])^2)])
    J1 = eliminate(I1, [a;b])
    I2 = ideal(R, [a[1]+2*a[2]+coeff_p[3]//coeff_p[4], 2*a[1]*a[2]+a[2]^2+a[3]^2-coeff_p[2]//coeff_p[4], a[1]*(a[2]^2+a[3]^2)+coeff_p[1]//coeff_p[4],
                  b[1]+2*b[2]+coeff_q[3]//coeff_q[4], 2*b[1]*b[2]+b[2]^2+b[3]^2-coeff_q[2]//coeff_q[4], b[1]*(b[2]^2+b[3]^2)+coeff_q[1]//coeff_q[4],
                  3*z[1] - ((a[1]-b[2])^2 + a[3]^2 + (a[2]-b[2])^2 + (a[3]+b[3])^2 + (b[1]-a[2])^2 + b[3]^2)])
    J2 = eliminate(I2, [a;b])
    return [primary_decomposition(J1) , primary_decomposition(J2)]
end

minimal_poly_11 (generic function with 1 method)

We put the three functions together.

function minimal_poly(coeff_p, coeff_q)
    ## Input:
        ## coeff_p = list of coefficients of p
        ## coeff_q = list of coefficients of q
    ## Output:
        ## [(n,m), L] where
        ## n is the number of real roots of p
        ## m is the number of real roots of q
        ## L is a list of primary ideals among which there is the minimal polynomial of the Wasserstein distance of p,q
    @var t
    # we compute the number of real and complex roots of p and q
    p = t^3*coeff_p[1]+t^2*coeff_p[2]+t*coeff_p[3]+coeff_p[4]
    real_p = nreal(HomotopyContinuation.solve(System([p])))
    q = t^3*coeff_q[1]+t^2*coeff_q[2]+t*coeff_q[3]+coeff_q[4]
    real_q = nreal(HomotopyContinuation.solve(System([q])))
    # based on the number of real roots, we apply the appropriate function
    if real_p == 3
        if real_q == 3
            return [(real_p, real_q), minimal_poly_33(coeff_p, coeff_q)]
        elseif real_q == 1
            return [(real_p, real_q), minimal_poly_13(coeff_q, coeff_p)]
        end
    elseif real_p == 1
        if real_q == 3
            return [(real_p, real_q), minimal_poly_13(coeff_p, coeff_q)]
        elseif real_q == 1
            return [(real_p, real_q), minimal_poly_11(coeff_p, coeff_q)]
        end
    end
end

minimal_poly (generic function with 1 method)

Tests

We test our algorithm on various coefficient lists. We start with two polynomials that have all real roots.

coeff_p = [-32 47 -15 1] # this represents the polynomial p = -32 + 47*x - 15*x^2 + x^3
coeff_q = [-951 311 -32 1] # this represents the polynomial q = -951 + 311*x - 32*x^2 + x^3

1×4 Matrix{Int64}:
 -951  311  -32  1

mp = minimal_poly(coeff_p, coeff_q)

Tracking 3 paths... 100%|███████████████████████████████| Time: 0:00:14
  # paths tracked:                  3
  # non-singular solutions (real):  3 (3)
  # singular endpoints (real):      0 (0)
  # total solutions (real):         3 (3)

2-element Vector{Any}:
 (3, 3)
 Tuple{MPolyIdeal{QQMPolyRingElem}, MPolyIdeal{QQMPolyRingElem}}[(Ideal with 1 generator, Ideal with 1 generator)]

The output confirms that both \(p\) and \(q\) have \(3\) real roots.

mp[2][1][1]

Ideal generated by
  z[1]^6 - 1278*z[1]^5 + 660151*z[1]^4 - 175286252*z[1]^3 + 25035024823*z[1]^2 - 1809005947030*z[1] + 51596679643601

As predicted in Table 2 of [arXiv:2401.12735], the degree of the minimal polynomial is \(6\).

For a polynomial with one real root, and one with three real roots we have:

coeff_p = [-56 87 -94 1]
coeff_q = [3 4 -5 1]

1×4 Matrix{Int64}:
 3  4  -5  1

mp = minimal_poly(coeff_p, coeff_q)

2-element Vector{Any}:
 (1, 3)
 Tuple{MPolyIdeal{QQMPolyRingElem}, MPolyIdeal{QQMPolyRingElem}}[(Ideal with 1 generator, Ideal with 1 generator)]

Oscar.degree(mp[2][1][1])

9

Whereas for two polynomials with one real root each, we have:

coeff_p = [-7 -3 1 2]
coeff_q = [4 -5 -1 3]

1×4 Matrix{Int64}:
 4  -5  -1  3

mp = minimal_poly(coeff_p, coeff_q)

2-element Vector{Any}:
 (1, 1)
 Vector{Tuple{MPolyIdeal{QQMPolyRingElem}, MPolyIdeal{QQMPolyRingElem}}}[[(Ideal with 1 generator, Ideal with 1 generator)], [(Ideal with 1 generator, Ideal with 1 generator)]]

In this case we get \(2\) candidates for the minimal polynomial.

mp[2]

2-element Vector{Vector{Tuple{MPolyIdeal{QQMPolyRingElem}, MPolyIdeal{QQMPolyRingElem}}}}:
 [(Ideal with 1 generator, Ideal with 1 generator)]
 [(Ideal with 1 generator, Ideal with 1 generator)]

However, the two candidates have the same degree.

(Oscar.degree(mp[2][1][1][1]), Oscar.degree(mp[2][2][1][1]))

(18, 18)