Approximate polynomial gcds

Problem 5.1, computation of approximate polynomial greatest common divisors (gcds). Given two polynomials, $u(x)=\sum_{i=0}^m u_ix^i$ and $v(x)=\sum_{i=0}^nv_ix^i$, a fixed polynomial norm $\Vert \cdot \Vert$, and a positive $\epsilon$, compute an $\epsilon$-gcd of u(x) and v(x) , that is, a nonunique polynomial $d^*(x)$ of the maximum degree $d_{\epsilon}$ that divides two polynomials $u^*(x)$ and $v^*(x)$ satisfying

 $$\ {\rm deg} \ u^* (x)\leq \ m, \ \ {\rm deg} \ v^* (x)\leq \ n ,$$

 $$\parallel u^* (x)-u(x)\parallel \leq \epsilon \parallel u(x) ... ...lel v^* (x)-v(x)\parallel \leq \epsilon \parallel v(x)\parallel .$$

The study of this problem is motivated by its applications to control linear systems, network theory and computer aided design (see some bibliography in [CGTW95], [EGL96], [CCC98], [P98], [P,a]).

Effective algorithms are known that estimate the weighted least-squares distances from u(x) and v(x) to two polynomials $u^*(x)$ and $v^*(x)$ satisfying (5.1) and divisible by a given polynomial $d^*(x)$ [CGTW95], [P,a]. The computation of such a candidate polynomial $d^*(x)$ is a harder problem, however. Only exponential time algorithms are known for its solution [Hou77], [KL96], [CCC98]; they are based on using the (weighted) least-squares norm and quadratic programming approach. There are several alternative techniques, which define the polynomial neighborhoods (5.2) based on the coefficient vector norms. These techniques give us heuristic solutions of Problem 5.1, that is, produce a common $\epsilon$-divisor $d_\epsilon (x)$ of u(x) and v(x) and, furthermore, for some input pairs u(x) and v(x) verify that the degree of such a divisor is maximum but for other input pairs leave the latter question open. These techniques include Euclidean algorithm [Sc85], [NS91], [HS96], [EGL96], [EGL97], [SaSa97], computations with subresultant matrices [CGTW95], [EGL96], [EGL97] and Lazard's algorithm [L81], having an equivalent matrix pencil version [KM94].

We propose two alternative approaches, based on Padé approximation and polynomial rootfinding (cf. [P98], [P,a]). The former one is similar to the subresultant approach. It leads to a little simpler computations (involving symmetric Hankel or Bezout matrices of smaller sizes, versus unsymmetric subresultant matrices) but to a little less efficient techniques for the insurance of the $\epsilon$-divisiblity. This approach relies on the observation that the polynomial $g(x)={\rm gcd} (u(x),v(x))$ satisfies the equations

 $$g(x)=\frac{u(x)}{w(x)}=\frac{v(x)}{z(x)}$$

where $(z(x),~ w(x))$ is the (m,n) Padé approximation of the power series $h(x)=\sum_{i=0}^{\infty} h_ix^i=\frac{u(x)}{v(x)}$. Here, with no loss of generality, we assume that $v(0)\ne 0$.

The computation of the gcd by this approach consists of the three successive stages of the computation of the following polynomials:

1.

$h_N(x)= h(x)$ mod $x^{N+1}, N=m+n$,

2.

z(x) and w(x) , being the (m,n) -Padé approximation of $h_N(x)$,

3.

g(x) of (5.3).

For the computation of an $\epsilon$-gcd, stages 2 and 3 are modified as follows [P98], [P,a]. At first, the pairs $(z_d(x), w_d(x)$) of polynomials being the (m-d, n-d) Padé approximations of $h_N(x)$ are computed for d recursively decreasing from a given upper bound $d^+$ on the degree of $\epsilon$-gcds (say, $d^+={\rm min}\{m,n\}$) down to an unknown nonnegative integer $d^-$; each computed polynomial $w_d(x)$ is tested for being an $\epsilon$-divisor of u(x) ; the process continues recursively for $d=d^+, ~d^+-1,\ldots ,$ until we either arrive at d=0 (in which case the value $d^-=0$ and the constant polynomial $g^-(x)=1$ are output) or satisfy the inequalities $\vert\vert u(x)-w_d(x)g^-(x)\vert\vert \leq \epsilon \vert\vert u(x)\vert\vert$, $\vert\vert v(x)-z(x)g^-(x)\vert\vert \leq \epsilon \vert\vert v(x)\vert\vert$for some polynomials $g^-(x)$ of degree at most d and z(x) of degree at most n-d (in which case $d^-=d$ and $g^-(x)$ are output).

Such an algorithm outputs a common $\epsilon$-divisor of u(x) , v(x) having degree $d^-$. If other methods show that $d^-$ is actually an upper (and consequently the tight) bound on the degree of $\epsilon$-gcd, then the computed $\epsilon$-divisor $g^-$ is a certified $\epsilon$-gcd of u(x) and v(x) .

The arithmetic cost of performing the above computations is bounded by $O_A(( d^+ - d^- + 1) n \log n)$, for a generic input pair (u(x) , v(x)), and by $O_A((d^+-d^- +1 )n\log^2 n)$, for any pair (u(x) , v(x)) , provided that the Padé approximations are computed by means of the extended Euclidean algorithm [BGY80], [BP94]. An alternative way is to obtain them via the computations with the associated structured (Hankel or Bezout) matrices, that is, via the computation of the rank r of such an $n\times n$ matrix $H_n$ or $B_n$ followed by solving a nonsingular linear system of r equations with a structured coefficient matrix [BGY80], [BP94]. In the context of $\epsilon$-gcd computation, one computes the rank numerically and obtains an upper bound $d^+$ on the degree of $\epsilon$-gcds equal to the number of the eigenvalues of $H_n$ and $B_n$ lying in the interval $(-\sigma(\epsilon), \sigma(\epsilon))$ for a positive function $\sigma(\epsilon)$. This way of computation is slower but allows numerical stabilization and defines an upper bound $d^+$ on the degree of an $\epsilon$-gcd.

In our second approach, we rely on the observation that the behavior of the gcd of the perturbed pairs of polynomials is directly governed by the perturbation of their zeros and only indirectly by the perturbation of their coefficients. Thus, in this approach, we first recall Theorem 3.1 to compute some higher precision approximations, $y^*_i$ and $z^*_j$, to all the zeros $y_i$ of u(x) and $z_j$ of v(x) , $i=1,\ldots,m$; $j=1,\ldots,n$. Then we define a bipartite graph with the two sets of vertices, $Y=\{y^*_i\}$ and $Z=(z^*_j)$, and with the edges connecting $y^*_{i(k)}$ and $z^*_{j(k)}$ if and only if $\vert y^*_{i(k)}-z^*_{j(k)}\vert<\delta(\epsilon)$ for a specified positive $\delta(\epsilon)$. We compute the maximum cardinality matching $M=\{(y^*_{i(k)},z^*_{j(k)}\}$ in the graph and the candidate $\epsilon$-gcd, $g(x)=\prod_k(x-(y^*_{i(k)}+z^*_{j(k)})/2)$, where the product is over all k such that $(y^*_{i(k)},z^*_{j(k)})\in M$. Finally, we test g(x) for being a common $\epsilon$-divisor of u(x) and v(x) and compare its degree with the known upper bound $d^+$, set equal to min$\{ m, n \}$ or obtained by different approaches.

The computational cost of our second approach is quite low; it is defined by Theorem 3.1 and by the well known bound of $O(N^{2.5})$comparisons, N=m+n , on the cost of the matching computation. We also show that the latter bound can be decreased to O(mn) comparisons (by replacing the matching stage by the computation of the connected components in a bipartite graph) if elastic choice of $\epsilon$is allowed (cf. [P98], [P,a]).

This approach leads to a natural extension of the definition of $\epsilon$-gcd where the norm (of the perturbation) of a polynomial is defined in terms of the absolutely maximum (perturbation of) its zeros rather than its coefficients. In this case our approach always outputs correctly an $\epsilon$-gcd. Furthermore, the resulting solution of such an $\epsilon$-gcd problem can be extended to a heuristic algorithm for the problem of computing a Hankel matrix of a lower rank $r^-$approximating a given Hankel matrix of rank r , with $r^- < r$, this problem has important applications to signal and image processing [Cad88], [CaW90], [DeM94], [vdV96].