ENH: add numerically robust polygcd function for numerical calculation of polynomial greatest common divisor

GCD for polynomials is an important operation for computing sums of polynomial ratios A(x)/B(x) + C(x)/D(x). The naive approach to computing GCD(B(x),D(x)) is numerically ill-conditioned:

import numpy as np

def polygcd(a,b):
    '''return monic GCD of polynomials a and b'''
    pa = a
    pb = b
    M = lambda x: x/x[0]
    # find gcd of a and b
    while len(pb) > 1 or pb[0] != 0:
        # Danger Will Robinson! requires numerical equality
        q,r = np.polydiv(pa,M(pb))
        pa = pb
        pb = r
    return M(pa)

def polylcm(a,b):
    '''return (Ka,Kb,c) such that c = LCM(a,b) = Ka*a = Kb*b'''        
    gcd = polygcd(a,b)
    Ka,_ = np.polydiv(b,gcd)
    Kb,_ = np.polydiv(a,gcd)
    return (Ka,Kb,np.polymul(Ka,a))

I figured there was an easy workaround that would be more robust, but apparently there isn't. There is some promising literature on the subject (see below), but I'm woefully underskilled for porting these algorithms to numpy, much less understanding how they work:

Here's something that uses material from Jan Elias's thesis as a starting point. It's certainly not optimal but it works for me and maybe some of you can improve:

import numpy as np 
from scipy.linalg import toeplitz

def cauchy(p,k):
    return toeplitz(np.hstack([p,np.zeros(k-1)]),np.zeros(k))

def sylvester(pa,pb,k=1):
    na = len(pa)-k
    nb = len(pb)-k
    Sa = cauchy(pa,nb)
    Sb = cauchy(pb,na)
    return np.hstack([Sa,Sb])

def polygcd(pa,pb,thresh=1e-10):
    '''
    return polynomials (g,u,v) such that 
       g = approximate gcd(pa,pb)
       pa = g*u
       pb = g*v

    this utilizes a simplified (suboptimal) version of
    Jan Elias's thesis
    https://who.rocq.inria.fr/Jan.Elias/pdf/je_masterthesis.pdf'''

    # trivial case: at least one polynomial is degree 0:
    if len(pa) == 1 or len(pb) == 1:
        return (np.array([1]), pa, pb)   

    # compute rank of S:
    # (there are more efficient ways to do this, per the thesis)
    S = sylvester(pa,pb)
    _,s,_ = np.linalg.svd(S)
    k = (s<thresh).sum()

    # another trivial case: pa and pb are relatively prime
    if (k == 0):
        return (np.array([1]), pa, pb)

    # find nullspace of Sk, partition into [-v,u]
    # (there are probably more efficient ways to do this... not sure)
    Sk = sylvester(pa,pb,k)
    _,_,v=np.linalg.svd(Sk)
    uv = v[-1,:]

    # solve [Ck(u); Ck(v)]x = [pa; pb] for x
    v = -uv[:len(pb)-k]
    u = uv[-(len(pa)-k):]
    Cuv = np.vstack([cauchy(u,k+1), cauchy(v,k+1)])
    pab = np.hstack([pa,pb])
    x,_,_,_ = np.linalg.lstsq(Cuv,pab)    
    return x/x[0],u*x[0],v*x[0]

This looks interesting but I need some time to digest it. What is the numerical difficulty, leading zeros in the division/remainder? If so it may be useful to experiment with Chebyshev polynomials, although in that case the location of arbitrary zeros in the complex plane may not be quite as accurate. Hmm... I probably need to read the thesis.

Thinking about it some more, it appears to me that when using floating point there needs to be a scale parameter of some sort that will be problem dependent. For example, 1e-6*x**2 + x + 1 is very much linear near the origin, but for large x the quadratic term dominates. What sort of problem are you interesting in solving?

Great questions -- I'd like to respond but I got diverted from this topic at work for a few days so it will have to wait before I can give a good background.

The short answer is that I'm doing some rational transfer function manipulation, and was surprised to find out that even for small-order polynomials (e.g. 2nd-order), I'm running into numerical issues.

I read the thesis and scaling is covered in the pre-processing step, but I think the general approach is not suitable for most practical problems. For one thing, the distance between polynomials is measured as the norm of the difference in coefficients, which assumes that the coefficients are exact modulo some small error. For polynomials derived from data, that is unlikely to be the case, and indeed the degree needed is uncertain. There are other problems along that line. So an important question to answer is how do the polynomials arise in a particular problem. The scaling is also chosen to optimize the Sylvester matrix, and that may not be the proper criterion. That is to say, it is important to specify exactly what you are attempting to do, where the data comes from, and what use will be made of the result.

I have some working code for LCM and GCD of polynomials from some time ago and they work sufficiently well for most tasks (excluding the ill-conditioned corner cases). The code was written some time ago and I don't think it fits into my toolbox' overall theme so I was thinking about cleaning it up and putting it into NumPy.

Especially LCM code I have been using a lot for MIMO systems in transfer function representations to handle minimality. It seems to agree with other commercial software and articles in the literature (as much as I could tested). Would it be a viable addition to NumPy after some cleanup? Or does it have to be some C/Fortran implementation?

You find the GCD and LCM code here

I'd argue that this belongs in np.polynomial, and we should not be adding more features to np.poly1d

Navigation Menu

Repository navigation

ENH: add numerically robust polygcd function for numerical calculation of polynomial greatest common divisor #4829

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

ENH: add numerically robust polygcd function for numerical calculation of polynomial greatest common divisor #4829

Description

Activity

jason-s commented on Jun 27, 2014

charris commented on Jun 27, 2014

charris commented on Jun 27, 2014

jason-s commented on Jun 27, 2014

charris commented on Jun 28, 2014

ilayn commented on Dec 18, 2016

eric-wieser commented on Mar 11, 2017

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions