The GPGCD Algorithm with the Bézout Matrix
for Multiple Univariate Polynomials

With the progress of algebraic computation with polynomials and matrices, we are paying more attention to approximate algebraic algorithms. Algorithms for calculating approximate GCD, which are approximate algebraic algorithms, consider a pair of given polynomials $F$ and $G$ that are relatively prime in general, and find $\tilde{F}$ and $\tilde{G}$ which are close to $F$ and $G$, respectively, in the sense of polynomial norm, and have the greatest common divisor of a certain degree. These algorithms can be classified into two categories: 1) for a given tolerance (magnitude) of $\|F-\tilde{F}\|$ and $\|G-\tilde{G}\|$, maximize the degree of approximate GCD, and 2) for a given degree $d$, minimize the magnitude of $\|F-\tilde{F}\|$ and $\|G-\tilde{G}\|$.

In both categories, algorithms based on various methods have been proposed including the Euclidean algorithm ([1], [19], [20]), low-rank approximation of the Sylvester matrix or subresultant matrices ([5], [6], [10], [11], [12], [21], [24], [26], [27], [29]), Padé approximation ([17]), and optimizations ([3], [13], [15], [22], [28]). Among them, the second author of the present paper has proposed the GPGCD algorithm based on low-rank approximation of subresultant matrices by optimization ([24], [25], [26]), which belongs to the second category above. Based on the researches mentioned above, the authors of the present paper have proposed the GPGCD algorithm using the Bézout matrix ([2]), which is the previous research of this paper.

In this paper, we propose the GPGCD algorithm with the Bézout matrix for multiple polynomials, while subresultant matrices have been used in the original GPGCD algorithm for multiple polynomials ([25]). We show that the proposed algorithm is more efficient than the original one with maintaining almost the same accuracy for most of the cases.

The paper is organized as follows. In Section 2, we give a definition of the approximate GCD problem. In Section 3, we give a formulation of the transformation of the approximate GCD problem to the optimization problem using the Bézout matrix. In Section 4, we review the modified Newton method used for optimization. In Section 5, we illustrate the proposed algorithm and give a running time analysis. In Section 6, the results of experiments are shown.

Let $F_{1}(x),\dots,F_{n}(x)$ be univariate polynomials with real coefficients of degree at most $m$:

In this paper, for a polynomial $F(x)=f_{m}x^{m}+\dots+f_{0}x^{0}$, the norm $\|F\|$ denotes the 2-norm defined as $\|F\|_{2}:=(f_{m}^{2}+f_{m-1}^{2}+\dots+f_{0}^{2})^{\frac{1}{2}}$. For a vector $(a_{1},\dots,a_{n})\in\mathbb{R}^{n}$, the norm $\|(a_{1},\dots,a_{n})\|$ denotes the 2-norm defined as $\|(a_{1},\dots,a_{n})\|_{2}:=(a_{1}^{2}+\dots+a_{n}^{2})^{\frac{1}{2}}$.

Here, we give a definition of an approximate GCD.

Algorithms for calculating approximate GCD can be classified into two categories: 1) for a given tolerance $\epsilon$, make the degree of approximate GCD as large as possible, and 2) for a given degree $d$, minimize the magnitude $\big{\|}(\|F_{1}-\tilde{F}_{1}\|,\dots,\|F_{n}-\tilde{F}_{n}\|)\big{\|}$. In this paper, we focus on the second category of the approximate GCD algorithms, solving the following problem.

For solving the approximate GCD problem, we transfer the approximate GCD problem to a constrained minimization problem with the Bézout matrix.

For the Bézout matrix, we have a following theorem.

For polynomials $\tilde{F}_{1}(x),\dots,\tilde{F}_{n}(x)$ in Problem 2, let $\tilde{B}=\operatorname{Bez}(\tilde{F}_{1}(x),\dots,\tilde{F}_{n}(x))=(\bm{b}_{1},\dots,\bm{b}_{m})$. In the case that polynomials $\tilde{F}_{1}(x),\dots,\tilde{F}_{n}(x)$ have an exact GCD of degree $d$, Theorem 4 shows that $\bm{b}_{1},\dots,\bm{b}_{m-d}$ are linearly independent, and $\bm{b}_{m-d+1}$ can be represented as a linear combination of $\bm{b}_{1},\dots,\bm{b}_{m-d}$. So there exists a vector $\bm{y}\in\mathbb{R}^{m-d}$, such that $(\bm{b}_{1},\dots,\bm{b}_{m-d})\bm{y}=\bm{b}_{m-d+1}$.

Let $\bm{s}$ be the vector of the coefficients of the input polynomials: $\bm{s}:=(f_{10},\dots,f_{1m},\dots,f_{n0},\dots,f_{nm})$. Let the Bézout matrix of the input polynomials represented by $\bm{s}$ be $B(\bm{s})=(\bm{b}(\bm{s})_{1},\dots,\bm{b}(\bm{s})_{n})$. In the same way, let the Bézout matrix of the polynomials we find ($\tilde{F}_{1}(x),\dots,\tilde{F}_{n}(x)$) be $B(\bm{s}+\Delta\bm{s})$, where $\Delta\bm{s}=(\tilde{f}_{10}-f_{10},\dots,\tilde{f}_{1m}-f_{1m},\dots,\tilde{f}_{n0}-f_{n0},\dots,\tilde{f}_{nm}-f_{nm})$.

From the above, we can transfer Problem 2 to a constrained minimization problem with the objective function:

and the constraints:

with variables:

We consider a constrained minimization problem of minimizing an objective function $f(\bm{x}):\mathbb{R}^{s}\rightarrow\mathbb{R}$ which is twice continuously differentiable, subject to the constraints $\bm{g}(\bm{x})=(g_{1}(\bm{x}),\dots,g_{t}(\bm{x}))^{T}=\bm{0}$, where $g_{i}(\bm{x})$ is a function of $\mathbb{R}^{s}\rightarrow\mathbb{R}$ and is also twice continuously differentiable. For solving the constrained minimization problem, we use the modified Newton method by Tanabe ([23]), which is a generalization of the Gradient Projection method ([18]), as in the original GPGCD algorithm ([26]). For $\bm{x}_{k}$ which satisfies $\bm{g}(\bm{x}_{k})=0$, we calculate the search direction $\bm{d}_{k}$ and the Lagrange multipliers $\bm{\lambda}_{k}$ by solving the following linear system

where $J_{\bm{g}}(\bm{x})$ is the Jacobian matrix represented as

We give an algorithm for calculating approximate GCD of multiple polynomials using the Bézout matrix in this section.

We have implemented our GPGCD algorithm on a computer algebra system Maple 2021. For $i=1,\dots,n$, the test polynomials $F_{i}(x)$ are generated as

Here, $\tilde{F}_{i}(x)$ and $H(x)$ are polynomials of degrees $\deg(F_{i})-d$ and $d$, respectively. Note that $F_{i}(x)$ are relatively prime polynomials. They are generated as polynomials that their coefficients are floating-point numbers and their absolute values are not greater than 10.¹⁾¹⁾1)Coefficients are generated with the Mersenne Twister algorithm [14] by built-in function Generate with RandomTools:-MersenneTwister in Maple, which approximates a uniform distribution on $[-10,10]$. The noise polynomials $\hat{F}_{i}(x)$ are polynomials of degree $\deg(F_{i})-1$, which are randomly generated with coefficients given as the same as for $\tilde{F}_{i}(x)$ and $H(x)$.

We have generated 5 test groups of test polynomials, each group comprising 100 tests. We set the degree of the input polynomials for all groups as 10, and the degree of the approximate GCD as $3,4,5,6,7$, respectively (shown as in Table 1). We set the number of polynomials for each tests as $n=10$, and set the norm of the noise $e$ as $e=0.01$ in our tests. The stop criterion $\epsilon$ and the stop width $\alpha$ in Algorithm 1 are set as $0.1$ and $1$, respectively.

We have carried out the tests on CPU Intel(R) Xeon(R) Silver 4210 at 2.20GHz with RAM 256GB under Linux 5.4.0 with Maple 2021.

We have proposed an algorithm using the Bézout matrix based on the GPGCD algorithm for multiple polynomials.

The experimental results show that, in most of the cases, the minimizer of the proposed algorithm is closed to the minimizer of the GPGCD-Sylvester-multiple algorithm. For the cases in which the minimizers of both algorithms are closed to each other, total computing time, the number of iterations, average computing time per iteration of the GPGCD-Bézout-multiple algorithm are smaller than those of the GPGCD-Sylvester-multiple algorithm, respectively. Thus, we see that the proposed algorithm is more efficient than the GPGCD-Sylvester-multiple algorithm. By comparing the computing time per iteration, the experimental results show that when the degree of the approximate GCD is higher, the computing time per iteration of the proposed algorithm is smaller. This is consistent with the arithmetic running time analysis shown in Section 5.4.

The previous research has shown that the GPGCD-Bézout algorithm was not so accurate in the sense of the norm of perturbations compared with the GPGCD-Sylvester algorithm. However, new experimental results show that, with the improvement for calculating the approximate polynomials as shown in (16), the GPGCD-Bézout-multiple algorithm calculate approximate GCDs with almost the same accuracy as the GPGCD-Sylvester-multiple algorithm in the sense of the norm of perturbations, for the tests in which the minimizers of both algorithms are closed to each other.

On the other hand, the experimental results show that the number of tests in which the perturbation of the GPGCD-Sylvester-multiple algorithm is smaller than that of the proposed algorithm is larger than that of the opposite. Thus, the GPGCD-Sylvester-multiple algorithm shows higher stability than the proposed algorithm. Improving stability of the proposed algorithm will be one of the topics of further research.

This research was partially supported by JSPS KAKENHI Grant Number JP20K11845. The author (Chi) was supported by the research assistant fellowship from Graduate School of Pure and Applied Sciences, University of Tsukuba.