Computing Approximate Greatest Common Right Divisors of Differential Polynomials

TL;DR

This paper introduces an algorithm for computing approximate greatest common right divisors of differential polynomials, extending GCD concepts to non-commutative differential operators with proven local well-posedness and quadratic convergence.

Contribution

It presents the first algorithm for approximate GCRD of differential polynomials, including a SVD-based initial approximation and convergence analysis.

Findings

Algorithm effectively finds nearby differential polynomials with non-trivial GCRD.

Proven local well-posedness of the approximate GCRD problem.

Quadratic convergence of Newton iteration with a good initial guess.

Abstract

Differential (Ore) type polynomials with "approximate" polynomial coefficients are introduced. These provide an effective notion of approximate differential operators, with a strong algebraic structure. We introduce the approximate Greatest Common Right Divisor Problem (GCRD) of differential polynomials, as a non-commutative generalization of the well-studied approximate GCD problem. Given two differential polynomials, we present an algorithm to find nearby differential polynomials with a non-trivial GCRD, where nearby is defined with respect to a suitable coefficient norm. Intuitively, given two linear differential polynomials as input, the (approximate) GCRD problem corresponds to finding the (approximate) differential polynomial whose solution space is the intersection of the solution spaces of the two inputs. The approximate GCRD problem is proven to be locally well-posed. A method based on the singular value decomposition of a differential Sylvester matrix is developed to produce an initial approximation of the GCRD. With a sufficiently good initial approximation, Newton iteration is shown to converge quadratically to an optimal solution. Finally, sufficient conditions for existence of a solution to the global problem are presented along with examples demonstrating that no solution exists when these conditions are not satisfied.

Equations

$$\partial y(t)=y(t)\partial +y^{\prime }(t)\text{ for }y(t)\in \mathbb{R}(t).$$

$$f(\partial )=\sum _{0\le i\le M}{f}_i(t){\partial }^i\text{ acts on y(t) as }\sum _{0\le i\le M}{f}_i(t)\frac{d^i}{d{t}^i}y(t).$$

$$f=\frac{1}{f_{-1}}\sum _{0\le i\le M}{f}_i{\partial }^i,$$

$$\Vert p\Vert =\Vert p{\Vert }_2={\left(\sum _{0\le i\le d}{p}_i^2\right)}^{1/2}.$$

$$\Vert f\Vert =\Vert f{\Vert }_2={\left(\sum _{0\le i\le M}\Vert f_i{\Vert }_2^2\right)}^{1/2}.$$

$$\Vert S{\Vert }_F^2=\sum _{ij}\Vert S_{ij}{\Vert }^2.$$

$$\overrightarrow{\mathrm{deg}}(f)=({\mathrm{deg}}_t{f}_0,{\mathrm{deg}}_t{f}_1,\dots ,{\mathrm{deg}}_t{f}_M),\text{ for }{f}_0,\dots ,f_M\in \mathbb{R}[t].$$

$$\overrightarrow{\mathrm{deg}}(f)<\overrightarrow{\mathrm{deg}}(g)\text{ if }{\mathrm{deg}}_t{f}_i\le {\mathrm{deg}}_t{g}_i\text{ for }0\le i\le M.$$

$${\Psi }_K(f)=\frac{1}{f_{-1}}(f_0,f_1,\dots ,f_M,0,\dots ,0)\in \mathbb{R}(t{)}^{1\times K}.$$

$$\mathbf{p}=(p_0,p_1,\dots ,p_d)\in {\mathbb{R}}^{1\times (d+1)}$$

$$\mathbf{f}=({\mathbf{f}}_0,\dots ,{\mathbf{f}}_{\mathbf{M}})\in {\mathbb{R}}^{1\times L},$$

$$\mathbf{f}\in {\mathbb{R}}^{1\times (M+1)(d+1)}.$$

$${\mathcal{C}}_K^{\mathbf{R}}(h)=\left(\begin{array}{c}{\Psi }_{K+D+1}(h)\\ {\Psi }_{K+D+1}(\partial h)\\ ⋮\\ {\Psi }_{K+D+1}({\partial }^{K}h)\end{array}\right)\in \mathbb{R}[t{]}^{(K+1)\times (K+D+1)}$$

$${\mathcal{C}}_K^{\mathbf{L}}(h)=\left(\begin{array}{c}{\Psi }_{K+D+1}(h)\\ {\Psi }_{K+D+1}(h\partial )\\ ⋮\\ {\Psi }_{K+D+1}(h{\partial }^K)\end{array}\right)\in \mathbb{R}[t{]}^{(K+1)\times (K+D+1)},$$

$$f=\sum _{0\le i\le M}{f}_i{\partial }^i,f^{\ast }=\sum _{0\le i\le M-D}{f}_i^{\ast }{\partial }^i\text{ and }h=\sum _{0\le i\le D}{h}_i{\partial }^i,$$

$$(f_0,f_1,\dots ,f_M)=(f_0^{\ast },\dots ,f_{M-D}^{\ast }){\mathcal{C}}_{M-D}^{\mathbf{R}}(h).$$

$$(f_0,f_1,\dots ,f_M{)}^T={\mathcal{C}}_D^{\mathbf{L}}(f^{\ast })(h_0,h_1,\dots ,h_D{)}^T.$$

$$S=S(f,g)=\left(\begin{array}{c}{\mathcal{C}}_{N-1}^{\mathbf{R}}(f)\\ {\mathcal{C}}_{M-1}^{\mathbf{R}}(g)\end{array}\right)\in \mathbb{R}(t{)}^{(M+N)\times (M+N)}$$

$$w=(u_0,u_1,\dots ,u_{N-1},v_0,v_1,\dots ,v_{M-1})\in \mathbb{R}(t{)}^{1\times (M+N)}.$$

$$C_k(b)=\left(\begin{array}{ccc}{b}_0& & \\ b_1& \ddots & \\ ⋮& \ddots & b_0\\ b_m& & b_1\\ & \ddots & ⋮\\ & & b_m\end{array}\right)\in {\mathbb{R}}^{(m+k+1)\times (k+1)}.$$

$$\hat{S}\in {\mathbb{R}}^{(M+N)(\mu +1)\times (M+N)(\mu +d+1)}.$$

$${\mathrm{deg}}_{\partial }\text{mbox}gcrd(f,g)=\frac{\dim \text{mbox}{\mathrm{null}}_{\mathrm{\ell }}(\hat{S})}{\mu +d+1}.$$

$$(f_0,f_1,\dots ,f_M)=(f_0^{\ast },\dots ,f_{M-D}^{\ast }){\mathcal{C}}_{M-D}^{\mathbf{R}}(h).$$

$$v_{i}=\sum_{j=0}^{(M-D+1)d}v_{ij}t^{j},~{}~{}~{}\mbox{$i=-1\ldots{M-D}$,}$$

$$v_{-1}\cdot (f_0,\dots ,f_M)=(v_0,\dots ,v_{M-D}){\mathcal{C}}_{M-D}^{\mathbf{R}}(h)$$

$$\Phi (h,f^{\ast },g^{\ast })=\Vert f-f^{\ast }h{\Vert }_2^2+\Vert g-g^{\ast }h{\Vert }_2^2.$$

$$\Vert f-\tilde{f}{\Vert }_2^2+\Vert g-\tilde{g}{\Vert }_2^2$$

$$\Phi (p)\le \Phi (x)\le \Phi (q),$$

$$\varphi :\mathbb{R}[t][\partial {;}^{\prime }{]}^2\to {\mathbb{R}}^{(M+N+2)(d+1)}$$

$${\cal F}_{\Omega}=\left\{\begin{array}[]{cc}({\widetilde{f}},{\widetilde{g}})\in{\mathbb{R}}[t][\partial;^{\prime}]^{2}\text{ such that }&{\deg_{\partial}}{\widetilde{f}}=M,\\ &{\deg_{\partial}}{\widetilde{g}}=N,\\ &{\deg_{\partial}}{\widetilde{h}}\geq D,\\ &{\widetilde{h}}={\mbox{gcrd}}({\widetilde{f}},{\widetilde{g}}),\\ &\|{\widetilde{h}}\|\leq\Omega,\\ &\operatorname{{\mbox{lcoeff}}}_{t}(\operatorname{{\mbox{lcoeff}}}_{\partial}{\widetilde{h}})=1,\\ &\text{ and }\Xi(\phi({\widetilde{f}},{\widetilde{g}}))={\mbox{true}}\end{array}\right\}.$$

Full text

Given two differential polynomials, we present an algorithm to find nearby differential polynomials with a non-trivial GCRD, where nearby is defined with respect to a suitable coefficient norm. Intuitively, given two linear differential polynomials as input, the (approximate) GCRD problem corresponds to finding the (approximate) differential polynomial whose solution space is the intersection of the solution spaces of the two inputs.

The approximate GCRD problem is proven to be locally well-posed. A method based on the singular value decomposition of a differential Sylvester matrix is developed to produce an initial approximation of the GCRD. With a sufficiently good initial approximation, Newton iteration is shown to converge quadratically to an optimal solution. Finally, sufficient conditions for existence of a solution to the global problem are presented along with examples demonstrating that no solution exists when these conditions are not satisfied.

1 Introduction

The problem of computing the GCRD in a symbolic and exact setting dates back to Ore (1933), who presents a Euclidean-like algorithm. See (Bronstein and Petkovšek, 1994) for an elaboration of this approach. Li and Nemes (1997) introduces a differential-resultant-based algorithm which makes computation of the GCRD very efficient using modular arithmetic. The technique of Li and Nemes (1997) is an extension of ideas presented by Grigor’ev (1990) for computing GCRDs of differential operators.

The analogous approximate GCD problem for usual (commutative) polynomials has been a key topic of research in symbolic-numeric computing since its inception. A full survey is not possible here, but we note the deep connection between our current work and that of (Corless et al., 1995); see also (Karmarkar and Lakshman, 1996), (Sasaki and Sasaki, 1997), and (Zeng and Dayton, 2004). Also important to this current work is the use of so-called structured (numerical) matrix methods for approximate GCD, such as structured total least squares (STLS) and structured total least norm (STLN); see (Botting et al., 2005) and (Kaltofen et al., 2005). More directly employed later in this paper is the multiple polynomial approximate GCD method of Kaltofen et al. (2006). This latter paper also provides a nice survey of the current state of the art in approximate GCDs. Finally, we modify the proof of Kaltofen et al. (2007b), an optimization approach to computing the GCD of multiple, multivariate commutative polynomials, to prove the existence of a globally nearest GCRD.

The goal of this paper is to devise an efficient, numerically robust algorithm to compute the GCRD when the coefficients in ${\mathbb{R}}$ are given approximately. Given $f,g\in{\mathbb{R}}(t)[\partial;^{\prime}]$, we wish to find ${\widetilde{f}},{\widetilde{g}}\in{\mathbb{R}}(t)[\partial;^{\prime}]$, where ${\widetilde{f}}$ is * near * $f$ and ${\widetilde{g}}$ is

• near * $g$, such that ${\deg_{\partial}}{\mbox{gcrd}}({\widetilde{f}},{\widetilde{g}})\geq 1$, where near is taken with respect to a distributed Euclidean norm. That is, ${\widetilde{f}}$ and ${\widetilde{g}}$ have an exact, non-trivial GCRD.

Linear differential polynomials and GCRD’s are key tools in finding closed form symbolic solution of systems of linear differential equations in modern computer algebra systems like Maple and Mathematica (see, e.g., (Salvy and Zimmermann, 1994) and (Abramov et al., 2005)). Equations with real (floating point) coefficients or parameters are regularly encountered and it is important to understand the stability of this fundamental tool in this case. Moreover, floating arithmetic is potentially much faster than managing large rational coefficients. We regard this paper as a positive and important initial exploration of this topic.

We commence with necessary preliminaries and well-known results that we expand upon in the remainder of this introductory section. In Section 2 we describe a linear algebra formulation of the approximate GCRD problem and that can be used in conjunction with truncated SVD (Giesbrecht and Haraldson, 2014; Haraldson, 2015; Corless et al., 1995) to compute nearby polynomials with an exact GCRD. Section 3 reformulates the approximate GCRD problem as a continuous unconstrained optimization problem. Sufficient conditions for existence of a solution are provided with an example showing that when this sufficient condition is not satisfied there is no solution. These results are complemented by showing that the Jacobian of the residuals has full rank and under ideal circumstances Newton iteration will converge quadratically. We generalize some results of Zeng and Dayton (2004) and Zeng (2011) to a non-commutative Euclidean domain showing that the problem is locally well-posed. In Section 4 we present our algorithms explicitly, discuss their complexity and evaluate the numerical robustness of our implementation on examples of interest.

A part of this work, presenting the SVD-based approach to approximate GCRD, but without the proof of existence a nearest solution or analysis of the corresponding optimization, is presented in the workshop paper (Giesbrecht and Haraldson, 2014). This is described in Section 4.1 of this current work.

1.1 Preliminaries

We review some well known results (Ore, 1933) and (Bronstein and Petkovšek, 1996) on differential polynomials.

The ring of differential (Ore) polynomials ${\mathbb{R}}(t)[\partial;^{\prime}]$ over the real numbers ${\mathbb{R}}$ provides a (non-commutative) polynomial ring structure to the linear ordinary differential operators. Differential polynomials have found great utility in symbolic computation, as they allow us to apply algebraic tools to the simplification and solution of linear differential equations; see (Bronstein and Petkovšek, 1994) for a nice introduction to the mathematical and computational aspects.

Let ${\mathbb{R}}(t)[\partial;^{\prime}]$ be the ring of differential polynomials over the function field ${\mathbb{R}}(t)$. ${\mathbb{R}}(t)[\partial;^{\prime}]$ is the ring of polynomials in $\partial$ with coefficients from the commutative field of rational functions, under the usual polynomial addition along with the non-commutative multiplication defined by

[TABLE]

Here $y^{\prime}(t)$ is the usual derivative of $y(t)$ with respect to $t$.

There is a natural action of ${\mathbb{R}}(t)[\partial;^{\prime}]$ on the space $\cal{C}^{\infty}[{\mathbb{R}}]$ of infinitely differentiable functions $y(t):{\mathbb{R}}\to{\mathbb{R}}$. In particular, for any $y(t)\in\cal{C}^{\infty}[{\mathbb{R}}]$,

[TABLE]

We maintain a right canonical form for all $f\in{\mathbb{R}}(t)[\partial;^{\prime}]$ by writing

[TABLE]

for polynomials $f_{-1},f_{0},\ldots,f_{M}\in{\mathbb{R}}[t].$ That is, with coefficients in ${\mathbb{R}}(t)$ always written to the left of powers of $\partial$. An analogous left canonical form exists as well.

A primary benefit of viewing differential operators in this way is that they have the structure of a left (and right) Euclidean domain. In particular, for any two polynomials $f,g\in{\mathbb{R}}(t)[\partial;^{\prime}]$, there is a unique polynomial $h\in{\mathbb{R}}(t)[\partial;^{\prime}]$ of maximal degree in $\partial$ such that $f=f^{*}h$ and $g=g^{*}h$ for $f^{*},g^{*}\in{\mathbb{R}}(t)[\partial;^{\prime}]$ (i.e., $h$ divides $f$ and $g$ exactly on the right). This polynomial $h$ is called the Greatest Common Right Divisor (GCRD) of $f$ and $g$ and it is unique up to multiplication by a unit (non-zero element) of ${\mathbb{R}}(t)$ (we could make this GCRD have leading coefficient 1, but this would introduce denominators from ${\mathbb{R}}[t]$, as well as potential numerical instability, as we shall see). An important geometric interpretation of GCRDs is that the GCRD $h$ of differential polynomials $f$ and $g$ is a differential polynomial whose solution space is the intersection of the solution spaces of $f$ and $g$.

Approximations require a norm, so we need a proper definition of the norm of a differential polynomial.

Definition 1.1**.**

We define the Euclidean norm for polynomials and a distributed coefficient norm for differential polynomials as follows:

For $p=\sum_{0\leq i\leq d}p_{i}t^{i}\in{\mathbb{R}}[t]$, define

[TABLE]

• 2.

For $f=\sum_{0\leq i\leq M}f_{i}\partial^{i}\in{\mathbb{R}}[t][\partial;^{\prime}]$, define

[TABLE]

We could extend the above definition of norm over ${\mathbb{R}}(t)$ and ${\mathbb{R}}(t)[\partial;^{\prime}]$. However it turns out that this is unnecessary and somewhat complicating. In practice, we perform most of our computations over ${\mathbb{R}}[t]$. In the cases where we are unable to avoid working over ${\mathbb{R}}(t)$, we simply solve an associate problem. This is done by clearing denominators and performing intermediate computations over ${\mathbb{R}}[t]$, then converting back to the representation over ${\mathbb{R}}(t)$. Note that the algebraic problem is always computing GCRDs and cofactors in ${\mathbb{R}}(t)[\partial;^{\prime}]$, and not the more intricate algebraic domain ${\mathbb{R}}[t][\partial;^{\prime}]$; see the discussion below.

Definition 1.2**.**

For any matrix $S\in{\mathbb{R}}[t]^{(M+N)\times(M+N)}$, we define the Frobenius norm $\|S\|_{F}$ by

[TABLE]

Main Problem: Approximate GCRD.  Given $f,g\in{\mathbb{R}}[t][\partial;^{\prime}]$ such that ${\mbox{gcrd}}(f,g)=1$ we wish to compute ${\widetilde{f}},{\widetilde{g}}\in{\mathbb{R}}[t][\partial;^{\prime}]$ with the same coefficient degree structure111The polynomial coefficients of $\partial^{i}$ have the same degree, i.e. $\deg{\widetilde{f}}_{i}\leq\deg f_{i}$ and $\deg{\widetilde{g}}_{i}\leq g_{i}$. as $f$ and $g$ such that $h={\mbox{gcrd}}({\widetilde{f}},{\widetilde{g}})$ with $D={\deg_{\partial}}h\geq 1$ and

• (i)

$\|f-{\widetilde{f}}\|^{2}_{2}+\|g-{\widetilde{g}}\|^{2}_{2}=\varepsilon$ is minimized, and

• (ii)

$D$ is the largest possible for the computed distance $\varepsilon$.

The differential polynomial $h$ is said to be an approximate GCRD of $f$ and $g$ if these conditions are satisfied. In general it is not easy to minimize $\varepsilon$, so instead we take a local optimization approach and compute an upper bound on this quantity. These upper-bounds will agree with the global minimum if $\varepsilon$ is sufficiently small. The algorithmic considerations will generally assume $D$ is fixed without loss of generality, since we can vary $D$ from 1 to $\min\{M,N\}$ to determine the (local) optimal value.

The approximate GCRD problem is a generalization of computing an ${\varepsilon{\text{-}}\mbox{GCD}}$ (Schönhage, 1985; Corless et al., 1995; Karmarkar and Lakshman, 1996; Emiris et al., 1997) in the commutative case. The requirement that the GCRD has maximal degree is difficult to certify outside the exact setting, however this usually is not a problem in our experiments. We prove that our formulation of the approximate GCRD problem has a solution with a minimal $\varepsilon$ (opposed to an infimum). Furthermore, if $D$ is fixed, then for a computed pair of nearby differential polynomials, we are able to certify that $\varepsilon$ is reasonably close to the optimal value through a condition number.

In our approach to the approximate GCRD problem we devise methods of performing division and computing an exact GCRD numerically. These tools are used in conjunction with our algorithm for computing a nearby pair of differential polynomials with an exact GCRD via the SVD. The nearby differential polynomials with an exact GCRD are used as an initial guess in a post-refinement Newton iteration.

It will also be necessary to define a partial ordering on differential polynomials. In later sections we will need to make use of this partial ordering to preserve structure.

Definition 1.3**.**

Let ${\overrightarrow{\deg}}:{\mathbb{R}}[t][\partial;^{\prime}]\to{\mathbb{Z}}^{M+1}$ be the degree vector function defined as

[TABLE]

For $f,g\in{\mathbb{R}}[t][\partial;^{\prime}]$ with ${\deg_{\partial}}f={\deg_{\partial}}g=M$ we write

[TABLE]

We define ${\overrightarrow{\deg}}(f)={\overrightarrow{\deg}}(g),{\overrightarrow{\deg}}(f)<{\overrightarrow{\deg}}(g),{\overrightarrow{\deg}}(f)\geq{\overrightarrow{\deg}}(g)$ and ${\overrightarrow{\deg}}(f)>{\overrightarrow{\deg}}(g)$ analogously.

We note that differential polynomials are written in a canonical ordering with highest degree coefficients appearing to the left in our examples. The degree vector function and most linearizations will appear in reverse order as a result. For convenience, we will assume that $\deg 0=-\infty$.

Definition 1.4**.**

Let $f\in{\mathbb{R}}[t][\partial;^{\prime}]$ where ${\deg_{\partial}}f=M$ is in standard form. The content of $f$ is given by $\operatorname{{\mbox{cont}}}(f)=\gcd(f_{0},f_{1},\ldots,f_{M})$. If $\operatorname{{\mbox{cont}}}(f)=1$, we say that the differential polynomial is primitive.

Proposition 1.5**.**

The ring ${\mathbb{R}}(t)[\partial;^{\prime}]$ is a non-commutative principal left (and right) ideal domain. For $f,g\in{\mathbb{R}}(t)[\partial;^{\prime}]$, with ${\deg_{\partial}}f=M$ and ${\deg_{\partial}}g=N$, we have the following properties (Ore, 1933).

1. (i)

${\deg_{\partial}}(fg)={\deg_{\partial}}f+{\deg_{\partial}}g$* , ${\deg_{\partial}}(f+g)\leq\max\{{\deg_{\partial}}f,{\deg_{\partial}}g\}$.* 2. (ii)

There exist unique $q,r\in{\mathbb{R}}(t)[\partial;^{\prime}]$ with ${\deg_{\partial}}r<{\deg_{\partial}}g$ such that $f=qg+r$ (right division with remainder). 3. (iii)

There exists $h\in{\mathbb{R}}(t)[\partial,^{\prime}]$ of maximal degree in $\partial$ with $f=f^{*}h$ and $g=g^{*}h$. $h$ is called the GCRD (Greatest Common Right Divisor) of $f$ and $g$, written ${\mbox{gcrd}}(f,g)=h$. $f^{*}$ and $g^{*}$ are called the left co-factors of $f$ and $g$. The GCRD is unique up to multiplication from a unit belonging to ${\mathbb{R}}(t)$. 4. (iv)

There exist $\sigma,\tau\in{\mathbb{R}}(t)[\partial;^{\prime}]$ such that $\sigma f=\tau g=\ell$ for $\ell$ of minimal degree. $\ell$ is called the LCLM (Least Common Left Multiple) of $f$ and $g$, written ${\mbox{lclm}}(f,g)=\ell$. The LCLM is unique up to multiplication from a unit belonging to ${\mathbb{R}}(t)$. 5. (v)

${\deg_{\partial}}{\mbox{lclm}}(f,g)={\deg_{\partial}}f+{\deg_{\partial}}g-{\deg_{\partial}}{\mbox{gcrd}}(f,g)$.

In an algebraic context we can clear denominators of our inputs and assume without loss of generality that our GCRD belongs to ${\mathbb{R}}[t][\partial;^{\prime}]$. We will also assume our inputs and output are primitive. Again, this is not algebraically necessary but will be important for the convergence of our subsequent optimization formulation (see Section 3.2). It is important to note that the co-factors of the GCRD need not belong to ${\mathbb{R}}[t][\partial;^{\prime}]$ even if we have $f,g,h\in{\mathbb{R}}[t][\partial;^{\prime}]$ such that ${\mbox{gcrd}}(f,g)=h$. This is not unexpected, as a similar situation occurs when computing GCD’s over ${\mathbb{Z}}[x]$, where cofactors in the GCD of primitive polynomials may well lie in ${\mathbb{Q}}[x]\setminus{\mathbb{Z}}[x]$. In essence, this is a computational technique to narrow the input domain, not a change to the problem being considered.

A related but considerably more difficult problem is computing ideal bases and factorizations completely within ${\mathbb{R}}[t][\partial;^{\prime}]$. This has been dealt with algebraically and in terms of exact computation by a number of authors, though not with respect to approximate coefficients; see for example (Heinle and Levandovskyy, 2016; Giesbrecht et al., 2016; Bell et al., 2017).

Most of our results involve transforming a representation of $f\in{\mathbb{R}}(t)[\partial;^{\prime}]$ into a representation over ${\mathbb{R}}(t)^{1\times K}$ for $K\geq{\deg_{\partial}}f$. We make extensive use of the following map.

Definition 1.6**.**

For $f\in{\mathbb{R}}(t)[\partial;^{\prime}]$ of degree $M$ in $\partial$ as in (1.1), and $K>M$, we define

[TABLE]

That is, $\Psi_{K}$ maps polynomials in ${\mathbb{R}}(t)[\partial;^{\prime}]$ of degree (in $\partial$) less than $K$ into ${\mathbb{R}}(t)^{1\times K}$.

It will be useful to linearize (differential) polynomials, that is, express them as an element of Euclidean space. For $p\in{\mathbb{R}}[t]$ with $\deg_{t}p=d$ we write

[TABLE]

For $f=f_{0}+f_{1}\partial+\cdots f_{M}\partial^{M}\in{\mathbb{R}}[t][\partial;^{\prime}]$ with ${\deg_{\partial}}f=M$ and $\deg_{t}f_{i}=d_{i}$ we write

[TABLE]

where $L=(d_{0}+1)+\cdots+(d_{M}+1)$. If $d\geq\max\{d_{i}\}$ we will sometimes pad each ${\bf{f_{i}}}$ with zeros to have precisely $d+1$ coefficients, and by a slight abuse of notation regard

[TABLE]

We will not do this unless specifically stated.

2 Computing the GCRD via Linear Algebra

In this section we demonstrate how to reduce the computation of the GCRD to that of linear algebra over ${\mathbb{R}}(t)$, and then over ${\mathbb{R}}$ itself. This approach has been used in the exact computation of GCRDs (Li and Nemes, 1997) and differential Hermite forms (Giesbrecht and Kim, 2013), and has the benefit of reducing differential, and more general Ore problems, to a system of equations over a commutative field. Here we will show that it makes our approximate version of the GCRD problem amenable to numerical techniques. We note that for computing approximate GCRDs of differential polynomials, much as for computing approximate GCDs of standard commutative polynomials, the Euclidean algorithm is numerically unstable, and thus we employ resultant-based techniques, as described below.

Since ${\mathbb{R}}(t)[\partial;^{\prime}]$ is a right (and left) Euclidean domain (Ore, 1933), a GCRD may be computed by solving a Diophantine equation corresponding to the Bézout coefficients. Using the subresultant techniques of Li (1998), we are able to transform the non-commutative problem over ${\mathbb{R}}(t)[\partial;^{\prime}]$ into a commutative linear algebra problem over ${\mathbb{R}}(t)$. This is done through a Sylvester-like resultant matrix. By using resultant-like matrices we are able to express the Bézout coefficients as a linear system over ${\mathbb{R}}(t)$ and compute a GCRD via nullspace basis computation.

Lemma 2.1**.**

Suppose $f,g\in{\mathbb{R}}(t)[\partial;^{\prime}]$ with ${\deg_{\partial}}f=N$ and ${\deg_{\partial}}g=M$. Then ${\deg_{\partial}}{\mbox{gcrd}}(f,g)\geq 1$ if and only if there exist $u,v\in{\mathbb{R}}(t)[\partial;^{\prime}]$ such that ${\deg_{\partial}}u<M$, ${\deg_{\partial}}v<N$, and $uf+vg=0$.

Proof.

This follows immediately from Proposition 1.5. ∎∎

Using Lemma 2.1 we can solve a Bézout-like system to compute a GCRD of two differential polynomials. This is characterized by the differential Sylvester matrix, based on the subresultant method of Li and Nemes (1997).

Definition 2.2**.**

Suppose $h\in{\mathbb{R}}[t][\partial;^{\prime}]$ has ${\deg_{\partial}}=D$. For any $K\in{\mathbb{N}}$, the matrix

[TABLE]

is the $K^{th}$ right differential convolution matrix of $h$. We note that the entries of ${\mathcal{C}}^{{\bf R}}_{K}(h)$ are written in their right canonical form, where the $\partial$’s appear to the right (polynomials in ${\mathbb{R}}[t]$ appear to the left). We note that $\deg_{t}\partial^{i}h=\deg_{t}h$, so the degree in $t$ of all entries of ${\mathcal{C}}^{{\bf R}}_{K}(h)$ is at most $\deg_{t}h$.

We analogously define the $K^{th}$ left differential convolution matrix of $h$ as ${\mathcal{C}}^{{\bf L}}_{K}(h)$ as

[TABLE]

where elements are written in their left canonical form, where the $\partial$’s appear to the left (polynomials in ${\mathbb{R}}[t]$ always appear to the right).

Both right and left differential convolution matrices can be used to perform multiplication. Suppose $f^{*}\in{\mathbb{R}}(t)[\partial;^{\prime}]$, $h\in{\mathbb{R}}(t)[\partial;^{\prime}]$ and $f=f^{*}h\in{\mathbb{R}}[t][\partial;^{\prime}]$, with

[TABLE]

with $f_{i},h_{i}\in{\mathbb{R}}[t]$ and $f^{*}_{i}\in{\mathbb{R}}(t)$. We can express the product of $f^{*}$ and $h$ as

[TABLE]

Similarly, we may write

[TABLE]

In keeping with our canonical ordering, we express our results in terms of right differential convolution matrices. We carefully observe that both the right and left differential convolution matrices described correspond to right multiplication. Left multiplication can be formulated in a similar manner.

Let $f,g\in{\mathbb{R}}(t)[\partial;^{\prime}]$ with ${\deg_{\partial}}f=M$ and ${\deg_{\partial}}g=N$. Then by Lemma 2.1 we have that ${\deg_{\partial}}{\mbox{gcrd}}(f,g)\geq 1$ if and only if there exist $u,v\in{\mathbb{R}}(t)[\partial;^{\prime}]$ such that ${\deg_{\partial}}u<N,{\deg_{\partial}}v<M$ and $uf+vg=0$. We can encode the existence of $u,v$ as an $(M+N)\times(M+N)$ matrix over ${\mathbb{R}}(t)$ in what we will call the differential Sylvester matrix.

Definition 2.3**.**

The matrix

[TABLE]

is the differential Sylvester matrix of $f$ and $g$.

This matrix (Li and Nemes, 1997) is analogous to the Sylvester matrix of real polynomials; see (von zur Gathen and Gerhard, 2013, Chapter 6). As expected, many useful properties of the Sylvester matrix over real polynomials still hold with the differential Sylvester matrix. These similarities become evident when we consider

[TABLE]

Then $uf+vg=0$ implies that $wS=0$, hence $w$ is a non-trivial vector in the (left) nullspace of $S$. In particular, this solution is equivalent to saying that $S$ is singular. Clearing denominators of $f$ and $g$, we may assume that $u,v\in{\mathbb{R}}[t][\partial;^{\prime}]$, i.e., they have polynomial coefficients, which implies that $S\in{\mathbb{R}}[t]^{(M+N)\times(M+N)}$. Moreover, for $f,g\in{\mathbb{R}}[t][\partial;^{\prime}]$ with $\deg_{t}f\leq d$ and $\deg_{t}g\leq d$ then $\deg_{t}S_{ij}\leq d$.

We summarize these results in the following lemma.

Lemma 2.4**.**

Suppose $f,g\in{\mathbb{R}}[t][\partial;^{\prime}]$, where ${\deg_{\partial}}f=M$, ${\deg_{\partial}}g=N$, $\deg_{t}f\leq d$ and $\deg_{t}g\leq d$.

1. (i)

$S=S(f,g)$* is singular if and only ${\deg_{\partial}}{\mbox{gcrd}}(f,g)\geq 1$.* 2. (ii)

${\deg_{\partial}}{\mbox{gcrd}}(f,g)=\dim\operatorname{{\mbox{null}_{\ell}}}(S)$, where $\operatorname{{\mbox{null}_{\ell}}}(S)$ is the left nullspace of $S$. 3. (iii)

For any $w=(u_{0},\ldots,u_{N-1},v_{0},\ldots,v_{M-1})\in{\mathbb{R}}(t)^{1\times(M+N)}$ such that $wS=0$, we have $uf+vg=0$, where $u=\sum_{0\leq i<N}u_{i}\partial^{i}$ and $v=\sum_{0\leq i<M}v_{i}\partial^{i}$. 4. (iv)

Suppose that ${\deg_{\partial}}{\mbox{gcrd}}(f,g)\geq 1$. Then there exists $w\in{\mathbb{R}}[t]^{1\times(M+N)}$ such that $wS=0$ and $\deg_{t}w\leq\mu=2(M+N)d$.

Proof.

Part (i) – (iii) follow from Lemma 2.1 and the discussion above. Part (iv) follows from an application of Cramer’s rule and a bound on the degree of the determinants of a polynomial matrix. ∎∎

2.1 Linear Algebra over ${\mathbb{R}}$

Let $S\in{\mathbb{R}}[t]^{(M+N)\times(M+N)}$ be the differential Sylvester matrix of $f,g\in{\mathbb{R}}[t][\partial;^{\prime}]$ of degrees $M$ and $N$ respectively in $\partial$, and degrees at most $d$ in $t$. From Lemma 2.4 we know that if a GCRD of $f$ and $g$ exists, then there is a $w\in{\mathbb{R}}[t]^{1\times(M+N)}$ such that $wS=0$, with $\deg_{t}w\leq\mu=2(M+N)d$.

Definition 2.5**.**

The $k^{th}$ convolution matrix of $b\in{\mathbb{R}}[t]$ with $\deg b=m$ is defined as

[TABLE]

Let $a\in{\mathbb{R}}[t]$ with $\deg a=\mu$ and define the mapping $\Gamma:{\mathbb{R}}[t]\to{\mathbb{R}}^{(\mu+1)\times(\mu+d+1)}$ by $\Gamma(a)=C_{d}(a)^{T}$. $\Gamma(a)$ is the left multiplier matrix of $a$ with respect to the basis $\langle 1,t,\ldots,t^{\mu+d}\rangle$.

A differential convolution matrix generalizes the convolution matrix in the role of linearizing multiplication between differential polynomials.

Definition 2.6**.**

Given the $(M+N)\times(M+N)$ differential Sylvester matrix $S$, we apply $\Gamma$ entry-wise to $S$ to obtain ${\widehat{S}}\in{\mathbb{R}}^{(M+N)(\mu+1)\times(M+N)(\mu+d+1)}$; each entry of $S$ in ${\mathbb{R}}[t]$ is mapped to a block entry in ${\mathbb{R}}^{(\mu+1)\times(\mu+d+1)}$ of ${\widehat{S}}$. We refer to ${\widehat{S}}$ as the inflated differential Sylvester matrix of $f$ and $g$.

Lemma 2.7**.**

Let $f,g\in{\mathbb{R}}[t][\partial;^{\prime}]$ have differential Sylvester matrix $S\in{\mathbb{R}}[t]^{(M+N)\times(M+N)}$ and inflated differential Sylvester matrix

[TABLE]

There exists a $w\in{\mathbb{R}}[t]^{1\times(M+N)}$ such that $wS=0$, if and only if there exists a ${\widehat{w}}\in{\mathbb{R}}^{(\mu+d+1)\times(M+N)(\mu+1)}$ such that ${\widehat{w}}{\widehat{S}}=0$. More generally,

[TABLE]

Proof.

This follows directly from the definition of $\Gamma$ and Lemma 2.4. ∎∎

We note that ${\widehat{S}}$ is no longer a square matrix. This will not pose too many problems as we will see in the following sections.

2.2 Division Without Remainder

While multiplication of differential polynomials with approximate numerical coefficients is straightforward, division is somewhat more difficult. We will generally require a division without remainder, for the computation of which we use a least squares approach. Given $f,h\in{\mathbb{R}}[t][\partial;^{\prime}]$ as in (2.1), we wish to find an $f^{*}\in{\mathbb{R}}(t)[\partial;^{\prime}]$ such that $\|f-f^{*}h\|$ is minimized. We will assume as usual that ${\deg_{\partial}}f=M$, ${\deg_{\partial}}h=D$ and $\deg_{t}f,\deg_{t}h\leq d$.

Much as in the (approximate polynomial) commutative case, we do this by setting the problem up as a linear system and then finding a least squares solution. Let us assume for now that $f=f^{*}h$ is exact, so this can be expressed as a linear system over ${\mathbb{R}}(t)$ by writing

[TABLE]

This system of equations is over-constrained (over ${\mathbb{R}}(t)$), but we note that the sub-matrix formed from the last $M-D+1$ columns of ${\mathcal{C}}^{{\bf R}}_{M-D}(h)$ is lower triangular, with diagonal entry $h_{D}\in{\mathbb{R}}[t]$. Thus, any exact quotient $h\in{\mathbb{R}}[t][\partial;^{\prime}]$ such that $f=f^{*}h$, in lowest terms, must have denominators dividing $h_{D}^{M-D+1}$, and in particular have denominators of degree at most $(M-D+1)\deg_{t}h_{D}\leq(M-D+1)d$. Equivalently, $h_{D}^{M-D+1}f^{*}\in{\mathbb{R}}[t]^{M-D+1}$. By applying Cramer’s rule on the last $M-D+1$ columns of ${\mathcal{C}}^{{\bf R}}_{M-D}(h)$, the degrees of the numerators in $f^{*}$ must be at most $(M-D+1)d$. Using this information we can formulate an associated problem with coefficients from ${\mathbb{R}}[t]$ and avoid performing linear algebra over ${\mathbb{R}}(t)$.

Now let $v_{-1},v_{0},\ldots,v_{M-D}$ be generic polynomials in $t$, with indeterminate coefficients of degree at most $(M-D+1)d$. I.e.,

[TABLE]

for indeterminates $v_{ij}$ with $v_{-1}\neq 0$. Then we are seeking to solve the linear system of equations

[TABLE]

for the $v_{ij}$. For each entry $f_{i}$ we have $(M-D+1)d+d+1$ equations; this is the degree $(v_{0},\ldots,v_{M-D}){\mathcal{C}}^{{\bf R}}_{M-D}(h)$ plus one, and we get one equation per coefficient. Hence there are $(M+1)((M-D+1)d+d+1)$ equations in $(M-D+2)(M-D+1)d$ unknowns. We then use a standard linear least squares solution to find the $v_{i}$ which minimizes the residual, and thus minimizes $\|f-f^{*}h\|$.

It may be desirable to find the lowest degree $v_{-1}$ which meets this criteria, for which we can use a simple binary search for a lower degree with reasonable residual (or alternatively use an SVD-based identification procedure).

Finally, a more straightforward approach to solving (2.2) is to simply use the solution from the last $M-D+1$ columns of ${\mathcal{C}}^{{\bf R}}_{M-D}(h)$. The last $M-D+1$ columns of ${\mathcal{C}}^{{\bf R}}_{M-D}(h)$ are lower triangular, with diagonal entries consisting of $h_{D}\in{\mathbb{R}}[t]$. While this does not yield a solution to the least squares normal equations, it is usually sufficiently good in practice, and considerably easier to formulate.

3 Optimization-based Formulation of Approximate GCRD

First we standardize some notation and assumptions. We assume that $f,g,{\widetilde{f}},{\widetilde{g}},h\in{\mathbb{R}}[t][\partial;^{\prime}]$ and $f^{*},g^{*}\in{\mathbb{R}}(t)[\partial;^{\prime}]$. Moreover, we assume that ${\widetilde{f}}=f^{*}h$ and ${\widetilde{g}}=g^{*}h$ and $h={\mbox{gcrd}}({\widetilde{f}},{\widetilde{g}})$. Intuitively, $f,g$ are our “input polynomials” and we will be identifying “nearby” ${\widetilde{f}},{\widetilde{g}}$ with a non-trivial GCRD $h$. Note that $f^{*},g^{*}$ have rational function coefficients. Later we will find it useful to clear fractions and work with a primitive associate.

We also assume degree bounds as follows: ${\deg_{\partial}}f={\deg_{\partial}}{\widetilde{f}}=M$, $\deg_{t}f,\deg_{t}{\widetilde{f}}\leq d$, ${\deg_{\partial}}g,{\deg_{\partial}}{\widetilde{g}}=N$, $\deg_{t}g,\deg_{t}{\widetilde{g}}\leq d$, ${\deg_{\partial}}h=D$, ${\deg_{\partial}}f^{*}=M-D$ and ${\deg_{\partial}}g^{*}=N-D$.

Using the method of Giesbrecht and Haraldson (2014), essentially the generalization of the SVD-based method of Corless et al. (1995) to differential polynomials, we will make an initial guess for ${\widetilde{f}},{\widetilde{g}}$; details are described in Section 4.1 of this paper. We then use optimization techniques to hone in on polynomials with minimal distance. While the techniques in that paper are not particularly effective at providing a nearest solution, they do provide a suitable initial guess, which we employ here.

We next describe how to formulate an objective function $\Phi$ that, when minimized, corresponds to a solution to the approximate GCRD problem. Define the objective function $\Phi:{\mathbb{R}}[t][\partial;^{\prime}]\times{\mathbb{R}}(t)[\partial;^{\prime}]^{2}\to{\mathbb{R}}$ as

[TABLE]

In keeping up with our notation from earlier, we observe that ${\widetilde{f}}=f^{*}h$ and ${\widetilde{g}}=g^{*}h$ in the context of the objective function $\Phi$, as $f$ and $g$ will typically be relatively prime. To compute guesses for the co-factors given $h$, we will perform an approximate division without remainder using the method of Section 2.2. We only require an initial guess for $f^{*}$ and $g^{*}$ to minimize $\Phi$, so this factorization doesn’t need to be exact, in the event that ${\mbox{gcrd}}(f,g)=h$.

We show that $\Phi$ has an attainable global minimum under appropriate assumptions. More precisely, there exist non trivial ${\widetilde{f}}$ and ${\widetilde{g}}$ such that

[TABLE]

is minimized. Furthermore, we will show that the approximate GCRD problem is locally well-posed.

3.1 Existence of Solutions

Lemma 3.1**.**

Let $f,h\in{\mathbb{R}}[t][\partial;^{\prime}]$, with monic leading coefficients, be not necessarily primitive, such that $f=f^{*}h$ for $f^{*}\in{\mathbb{R}}[t][\partial;^{\prime}]$ with ${\deg_{\partial}}f=M$ and ${\deg_{\partial}}h=D$. Then $\|f^{*}\|$ is bounded above.

Proof.

It follows that $f^{*}$ is bounded by the computing the Cramer solution to (2.2) using the last $M-D+1$ columns of ${\mathcal{C}}^{{\bf R}}_{M-D}(h)$. ∎

As an observation, we relax the assumption that $f$ is primitive (we work with an associate instead) in order to guarantee that $f^{*}\in{\mathbb{R}}[t][\partial;^{\prime}]$. This can be taken without loss of generality as the quantity $\|\operatorname{{\mbox{cont}}}(f)\|_{2}^{2}$ is bounded above and away from zero (as its leading coefficient is monic). Thus we may divide by it without affecting the quality of the results, as $\|f-f^{*}h\|$ is still well defined.

We will make use of the following well known fact from (Rudin, 1976, Theorem 4.16).

Fact 3.2**.**

Suppose that $\Phi$ is a continuous real function on a compact metric space $X$. Then there exist points $p$ and $q$ in $X$ such that

[TABLE]

for all $x\in X$. Precisely, $\Phi$ attains its minimum and maximum values at $p$ and $q$ respectively.

We first state a general version of the theorem where a logical predicate $\Xi:{\mathbb{R}}^{k}\to\{\text{true},\text{false}\}$ (for some $k$) can be chosen to impose additional constraints on the problem. For the rest of this section let

[TABLE]

be the combined coefficient vector function, i.e. for arbitrary $f,g\in{\mathbb{R}}[t][\partial;^{\prime}]$ we write $\phi(f,g)=({\bf{f}},{\bf{g}})$, where ${\bf{f}}$ and ${\bf{g}}$ are padded with zeros to have the desired dimensions.

The following lemma and its proof are analogous to (Kaltofen et al., 2007b, Theorem 2), which in turn generalizes the univariate argument of (Kaltofen et al., 2007a, Theorem 1).

Theorem 3.3** (Existence of Global Minima).**

Let $f,g\in{\mathbb{R}}[t][\partial;^{\prime}]\backslash\{0\}$, let $d=\max\{\deg_{t}f,\deg_{t}g\}$, ${\deg_{\partial}}f=M$, ${\deg_{\partial}}g=N$ and $D\leq\min\{M,N\}$. Furthermore, let $\Xi:{\mathbb{R}}^{(M+N+2)(d+1)}\to\{{\mbox{true}},{\mbox{false}}\}$ be a predicate on $\phi(f,g)$. We assume that the preimage $\Xi^{-1}({\mbox{true}})$ is a topologically closed set in ${\mathbb{R}}^{(M+N+2)(d+1)}$ with respect to the Euclidean norm. For a given $\Omega\in{\mathbb{R}}_{>0}$ we define the set of possible solutions by

[TABLE]

Suppose that ${\cal F}_{\Omega}\neq\emptyset$. Then the minimization problem

[TABLE]

has an attainable global minimum.

Proof.

Without loss of generality, we assume that $M\leq N$. Then we iterate the minimization over all $\ell\in{\mathbb{Z}}_{\geq 0}$ such that $D\leq\ell\leq M$ and coefficients $\rho\subset{\mathbb{R}}[t]^{\ell}$. Let ${\cal H}_{\ell,\rho}$ denote the set of all differential polynomials over ${\mathbb{R}}[t][\partial;^{\prime}]$ of degree $\ell$ with coefficients from $\rho$. We optimize over the continuous real objective function

[TABLE]

for $h\in{\cal H}_{\ell,\rho}$, ${\deg_{\partial}}f^{*}\leq M-D$ and ${\deg_{\partial}}g^{*}\leq N-D$. We fix the leading coefficient of $h$ with respect to $\partial$ to be monic, that is $\operatorname{{\mbox{lcoeff}}}_{t}\operatorname{{\mbox{lcoeff}}}_{\partial}h=1$.

Since the leading coefficient of $h$ is monic, we can write $G={\mbox{gcrd}}(f^{*}h,g^{*}h)$ with ${\deg_{\partial}}G\geq D$. Since $G$ is a multiple of $h$, we normalize $G$ so that $\operatorname{{\mbox{lcoeff}}}_{t}\operatorname{{\mbox{lcoeff}}}_{\partial}G=1$, i.e. the leading coefficient of $G$ is also monic. The restriction on $h$ that the leading coefficient of $h$ is monic enforces that ${\deg_{\partial}}G\geq D$. Furthermore, we restrict the domain of our function $\Phi$ to those $h,f^{*}$ and $g^{*}$ for which $(f^{*}h,g^{*}h)\in{\cal F}_{\Omega}$. If there is no such common factor $h$ and co-factors $f^{*}$ and $g^{*}$, then this pair of $\ell$ and $\rho$ does not occur in the minimization (3.2). By assumption we have that ${\cal F}_{\Omega}\neq\emptyset$, so there must be at least one possible case. We note that if $(0,0)\in{\cal F}_{\Omega}$, then $f^{*}=g^{*}=0$.

Now suppose that for the given $\ell$ and $\rho$, there are ${\widetilde{h}}\in{\cal H}_{\ell,\rho}$ and ${\widetilde{f}}^{*},{\widetilde{g}}^{*}$ satisfying ${\deg_{\partial}}{\widetilde{f}}^{*}\leq M-\ell$ and ${\deg_{\partial}}{\widetilde{g}}^{*}\leq N-\ell$ such that $({\widetilde{f}}^{*}{\widetilde{h}},{\widetilde{g}}^{*}{\widetilde{h}})\in{\cal F}_{\Omega}$. We shall prove that the function $\Phi$ has a value on a closed and bounded set (i.e., compact with respect to the Euclidean metric) that is smaller than elsewhere. Hence $\Phi$ attains a global minimum by Fact 3.2.

Clearly any solution ${\widetilde{h}}\in{\cal H}_{\ell,\rho}$ and ${\widetilde{f}}^{*},{\widetilde{g}}^{*}$ with $({\widetilde{f}}^{*}h,{\widetilde{g}}^{*}h)\in{\cal F}_{\Omega}$ but with $\Phi({\widetilde{h}},{\widetilde{f}}^{*},{\widetilde{g}}^{*})>\Phi(h,f^{*},g^{*})$ can be discarded. So the norm of the products $\|{\widetilde{f}}^{*}{\widetilde{h}}\|_{2}$ and $\|{\widetilde{g}}^{*}{\widetilde{h}}\|_{2}$ can be bounded from above. We have that $\|{\widetilde{h}}\|$ is bounded above by Lemma 3.1 because it is a right factor of ${\widetilde{G}}={\mbox{gcrd}}({\widetilde{f}}^{*}{\widetilde{h}},{\widetilde{g}}^{*}{\widetilde{h}})$ with $\|{\widetilde{G}}\|\leq\Omega$. We note that ${\widetilde{h}}$ has a monic leading coefficient, so $\|{\widetilde{h}}\|\geq 1$. We have that $\|{\widetilde{f}}^{*}\|$ and $\|{\widetilde{g}}^{*}\|$ (or the appropriate associate) are both bounded above by Lemma 3.1.

Thus we can restrict the domain of $\Phi$ to values that lie within a sufficiently large closed ball $B$. The function $\zeta$ that maps $(h,f^{*},g^{*})$ to the combined coefficient vector $\phi(f^{*}h,g^{*}h)$ of $f^{*}h$ and $g^{*}h$ is continuous. We minimize over $\zeta^{-1}(\Xi^{-1}({\mbox{true}})\cap\zeta(B))$, which is a compact set. ∎∎

For the less general version of the theorem, given arbitrary $f,g\in{\mathbb{R}}[t][\partial;^{\prime}]$, we define

[TABLE]

We observe that ${\cal S}$ is a closed subset of ${\mathbb{R}}^{(M+N+2)(d+1)}$, where ${\deg_{\partial}}f=M$, ${\deg_{\partial}}g=N$ and $d=\max\{\deg_{t}f,\deg_{t}g\}$. The set ${\cal S}$ corresponds to the combined coefficient vectors of ${\widetilde{f}}$ and ${\widetilde{g}}$ that have the same degree structure as $f$ and $g$.

Corollary 3.4**.**

Let $f,g\in{\mathbb{R}}[t][\partial;^{\prime}]\backslash\{0\}$, let $d=\max\{\deg_{t}f,\deg_{t}g\}$, ${\deg_{\partial}}f=M$, ${\deg_{\partial}}g=N$ and $D\leq\min\{M,N\}$. For a given $\Omega\in{\mathbb{R}}_{>0}$ we define the set of possible solutions by

[TABLE]

Suppose that ${\cal F}_{\Omega}\neq\emptyset$. Then the minimization problem

[TABLE]

has an attainable global minimum.

We note that Theorem 3.3 does not guarantee a unique minimum of $\Phi$, merely that $\Phi$ has an attainable minimum (as opposed to an infimum). The choice of $h,f^{*}$ and $g^{*}$ that we optimize over is important. If $\operatorname{{\mbox{lcoeff}}}_{t}\operatorname{{\mbox{lcoeff}}}_{\partial}h$ vanishes or $\|h_{0}\|,\ldots,\|h_{D-1}\|$ are quite large, then $f^{*}$ and $g^{*}$ can be ill-conditioned in the approximate GCRD problem. Furthermore, choosing overly large, small or poor degree structure in $t$ for $h$ can result in a $\Phi$ that cannot be minimized for the specified structure, but would otherwise have a minimum for a different choice of $h$.

Example 3.5**.**

Consider $f=\partial^{2}-2\partial+1$ and $g=\partial^{2}+2\partial+2$ (see (Kaltofen et al., 2007b) for an example with complex perturbations). Then $f$ and $g$ do not have a degree 1 approximate GCRD. That is, we show that there does not exist ${\widetilde{f}},{\widetilde{g}}\in{\mathbb{R}}[t][\partial;^{\prime}]$ where ${\deg_{\partial}}{\mbox{gcrd}}({\widetilde{f}},{\widetilde{g}})=1$ and $\|f-{\widetilde{f}}\|_{2}^{2}+\|g-{\widetilde{g}}\|_{2}^{2}$ is minimized.

The real monic Karmakar-Lakshman distance (Karmarkar and Lakshman, 1996, 1998) of

[TABLE]

occurs when the rational function

[TABLE]

is minimized for $h_{0}\in{\mathbb{R}}$. The minimum value (if it exists) of this function corresponds to the approximate GCRD $h=\partial-h_{0}$. The infimum is $2$, which is unattainable. There is no attainable global minimum.

The non-monic real Karmakar-Lakshman distance is $2$, which is achieved if and only if the leading coefficient vanishes. The minimum occurs when the rational function

[TABLE]

is minimized. The minimum value of this function corresponds to the approximate GCRD $h=h_{1}\partial+1$.

In particular, if we consider ${\widetilde{f}}=(-2\partial+1)(\varepsilon\partial+1)$ and ${\widetilde{g}}=(2\partial+2)(\varepsilon\partial+1)$, then $\|f-{\widetilde{f}}\|_{2}^{2}+\|g-{\widetilde{g}}\|_{2}^{2}$ becomes arbitrarily near $2$ as $\varepsilon\to 0$.

There is no real degree $1$ approximate GCRD, as

[TABLE]

is not defined in the monic case. In the non-monic case, if a minimum exists then it occurs when $\operatorname{{\mbox{lcoeff}}}_{\partial}h$ vanishes, so the minimum value is not defined either.

This example illustrates that not all $f,g\in{\mathbb{R}}[t][\partial;^{\prime}]$ have an approximate GCRD. Furthermore, we see that the requirement that $\operatorname{{\mbox{lcoeff}}}_{t}\operatorname{{\mbox{lcoeff}}}_{\partial}h=1$ and $\|h\|$ is bounded, from Theorem 3.3 are required, even if there are no additional constraints imposed.

Now it remains to show that it is possible to obtain a (locally) unique solution to $\Phi$. One of many equivalent conditions for uniqueness of an exact GCRD, is to require it to be primitive and have a monic leading coefficient. Numerically, to obtain a unique solution of the approximate GCRD problem, we impose the same constraints, making solutions locally unique.

3.2 Convergence of Newton Iteration and Conditioning

From Theorem 3.3 and Corollary 3.4 we know a solution to the approximate GCRD problem exists. We now show that a standard Newton iteration will converge quadratically when starting with an estimate sufficiently close to an approximate GCRD. We first describe the Jacobian of the residuals and show that the Jacobian has full rank. This leads to a first-order approximation of the Hessian matrix showing that it is locally positive definite around a global minimum when the residual is sufficiently small. The implication is that Newton’s method will converge quadratically (Boyd and Vandenberghe, 2004). If we consider structured perturbations, then we are able to obtain results similar to that of Zeng and Dayton (2004) to the overall conditioning of the system.

In this section we assume without loss of generality that $f^{*},g^{*}\in{\mathbb{R}}[t][\partial;^{\prime}]$ are primitive, and that $f$ and $g$ may no longer be primitive to simplify computations. We need to clear fractions of rational functions to apply our coefficient norms, and to linearize $h,f^{*}$ and $g^{*}$ as vectors of real numbers.

The residual of the approximate GCRD is

[TABLE]

where

[TABLE]

Intuitively, $\eta$ represents the number of components of $({\bf{f}},{\bf{g}})\in{\mathbb{R}}^{1\times\eta}$. Let $\nu$ be the number of variables needed to represent the coefficients of $h,f^{*}$ and $g^{*}$, i.e. $({\bf{h}},{\bf{f^{*}}},{\bf{g^{*}}})\in{\mathbb{R}}^{1\times\nu}$.

Recall that when $f=f^{*}h$, we can linearize this relationship with differential convolution matrices, by writing

[TABLE]

If $f_{i}$ is a coefficient of $f$ with $\deg_{t}f_{i}=d$, then we may write

[TABLE]

This relationship may be linearized over ${\mathbb{R}}$ through the use of convolution matrices. Writing

[TABLE]

we now have a direct method of computing ${\bf{f_{i}}}$ in terms of the coefficients of $f^{*}$ and $h$.

If we differentiate ${\bf{f^{*}h}}$ with respect to an entry from ${\bf{f^{*}}}$, then we will obtain the corresponding (linearized) row of ${\mathcal{C}}^{{\bf R}}_{M-D}(h)$. Similarly, differentiating ${\bf{f^{*}h}}$ with respect to an entry of ${\bf{h}}$ will give us a (linearized) column of ${\mathcal{C}}^{{\bf L}}_{D}(f^{*})$. This relationship becomes clear when we observe that

[TABLE]

Differentiating ${\bf{g^{*}h}}$ with respect to variables from ${\bf{g^{*}}}$ and ${\bf{h}}$ will produce similar results.

The Jacobian of $r(h,f^{*},g^{*})$ for arbitrary $h,f^{*}$ and $g^{*}$ may be expressed (up to column permutation) in block matrix form as

[TABLE]

where the block matrices are linearized accordingly. In our formulation of the approximate GCRD problem we normalize $\operatorname{{\mbox{lcoeff}}}_{t}\operatorname{{\mbox{lcoeff}}}_{\partial}h$ so that it is a predetermined constant, which results in essentially the same Jacobian as described above.

The only difference in the Jacobians, is that the $\nu^{th}$ column would become the zero column if differentiated with respect to $\operatorname{{\mbox{lcoeff}}}_{t}\operatorname{{\mbox{lcoeff}}}_{\partial}h$, since $\operatorname{{\mbox{lcoeff}}}_{t}\operatorname{{\mbox{lcoeff}}}_{\partial}h$ is constant. When normalized for computational purposes, the Jacobian belongs to ${\mathbb{R}}^{\eta\times\nu-1}$ instead (the last column is deleted). In the general case when ${\mbox{gcrd}}(f^{*},g^{*})=1$, $J$ is rank deficient by 1 and the $\nu^{th}$ column is a linear combination of the other columns. The following lemma, similar to (Zeng, 2011, Lemma 4.1), formalizes this statement.

Lemma 3.6**.**

Let $r$ be the residual described earlier with Jacobian $J$. Suppose that $\operatorname{{\mbox{lcoeff}}}_{t}\operatorname{{\mbox{lcoeff}}}_{\partial}h$ is a fixed non-zero constant. If ${\mbox{gcrd}}(f^{*},g^{*})=1$, then all non-zero columns of $J$ are linearly independent.

Proof.

Let $\vec{e}_{\nu}\in{\mathbb{R}}^{1\times\nu}$ be a unit vector whose last component is $1$. We write

[TABLE]

We shall prove the equivalent statement that the matrix

[TABLE]

has full rank.

Suppose the converse holds, then there exists $q_{1},q_{2},p\in{\mathbb{R}}[t][\partial;^{\prime}]$ with ${\deg_{\partial}}q_{1}\leq M-D,{\deg_{\partial}}q_{2}\leq N-D$ and ${\deg_{\partial}}p\leq D$ such that their combined coefficient vector satisfies

[TABLE]

Expressing this as multiplication over ${\mathbb{R}}[t][\partial;^{\prime}]$, we have that

[TABLE]

We conclude that ${\mbox{gcrd}}(f^{*}p,g^{*}p)=p$, as ${\mbox{gcrd}}(f^{*},g^{*})=1$. If $p=0$ or $q_{1}=0$ or $q_{2}=0$, then we are done (as ${\mathbb{R}}[t][\partial;^{\prime}]$ is a domain). Suppose that $p\neq 0$ and $q_{1}\neq 0$ and $q_{2}\neq 0$. Accordingly we must also have that ${\mbox{gcrd}}(q_{1}h,q_{2}h)={\mbox{gcrd}}(q_{1},q_{2})\alpha h=p$ for some $\alpha\neq 0$. Since ${\deg_{\partial}}p\leq{\deg_{\partial}}h$ it follows that ${\mbox{gcrd}}(q_{1},q_{2})=1$ so $p=\alpha h$.

Since $p=\alpha h$ we must have that $\alpha f^{*}=q_{1}$ and $\alpha g^{*}=q_{2}$. Now,

[TABLE]

On the other hand,

[TABLE]

This occurs if and only if $\alpha=0$. But in this case $p=0$ as well, so

[TABLE]

It follows that the only vector in the null space is the zero vector, hence $\begin{pmatrix}J\\ \vec{e}_{\nu}\end{pmatrix}$ has full rank. Since any subset of linearly independent vectors is also linearly independent, we have that when $\operatorname{{\mbox{lcoeff}}}_{t}\operatorname{{\mbox{lcoeff}}}_{\partial}h$ is a fixed non-zero constant that $J$ has rank $\nu-1$. ∎∎

Note that from the proof we see that if $\operatorname{{\mbox{lcoeff}}}_{t}\operatorname{{\mbox{lcoeff}}}_{\partial}h$ were not fixed, then the vector $({\bf{f^{*}}},{\bf{g^{*}}},{\bf{h}})^{T}$ forms a basis for the nullspace of $J$. Intuitively, if we did not fix $\operatorname{{\mbox{lcoeff}}}_{t}\operatorname{{\mbox{lcoeff}}}_{\partial}$ in advance, then there would be infinitely many tuples of $(h,f^{*},g^{*})$ with the same degree structure over ${\mathbb{R}}[t]$ that minimized $\Phi$, since for any $\alpha\neq 0$ we have

[TABLE]

In other words, we need to normalize $h$ in advance to obtain a unique solution.

Corollary 3.7**.**

Let $r$ be the residual defined earlier in this section with $\operatorname{{\mbox{lcoeff}}}_{t}\operatorname{{\mbox{lcoeff}}}_{\partial}h$ a non-zero constant. If $r=0$, then the Hessian matrix $\nabla^{2}\Phi(h,f^{*},g^{*})$ is positive definite.

Proof.

Let $J$ be the Jacobian of $r$. $J$ has full rank, so $J^{T}J$ has full rank and is positive semidefinite. If $r=0$, at the global minimum we have that $2J^{T}J=\nabla^{2}\Phi$, and $\nabla^{2}\Phi(h,f^{*},g^{*})$ is positive definite. ∎∎

When there is no residual, the Hessian $\nabla^{2}\Phi(h,f^{*},g^{*})$ is positive definite. It follows that if $f$ and $g$ are perturbed by a sufficiently small amount, then $\nabla^{2}\Phi$ remains locally positive definite, and Newton iteration will converge to the (local) global minimum with an initial guess that is sufficiently close.

We are able to obtain a condition number for a structured perturbation through the Jacobian of the residuals. Since $J$ has full rank, the smallest singular value $\sigma_{\nu-1}$ of $J(r(h,f^{*},g^{*}))$ is strictly positive. If we consider structured perturbations, then we are able to show that the approximate GCRD problem is (locally) well-posed.

In the next lemma, we make use of the fact that for any $f\in{\mathbb{R}}[t][\partial;^{\prime}]$, we have that $\|f\|_{2}=\|{\bf{f}}\|_{2}$.

Lemma 3.8**.**

Let $f,g,h,f^{*},g^{*}\in{\mathbb{R}}[t][\partial;^{\prime}]$ be such that $\Phi(h,f^{*},g^{*})<\varepsilon$ for some $\varepsilon>0$, with $\operatorname{{\mbox{lcoeff}}}_{t}\operatorname{{\mbox{lcoeff}}}_{\partial}h$ a fixed non-zero constant. Suppose ${\widehat{f}},{\widehat{g}},\widehat{h},\widehat{f^{*}},\widehat{g^{*}}\in{\mathbb{R}}[t][\partial;^{\prime}]$ possess the same degree structures as $f,g,h,f^{*}$ and $g^{*}$ and that

[TABLE]

Then,

[TABLE]

Proof.

Let $J=J(r(h,f^{*},g^{*}))$ be the Jacobian of the residuals from earlier in this section. We have that

[TABLE]

Ignoring high order terms and using the well known fact that for a (left) pseudo inverse $J^{+}$ of $J$, that $\|J^{+}\|_{2}=\frac{1}{\sigma_{\nu-1}}$ gives us

[TABLE]

A straightforward application of the triangle inequality gives

[TABLE]

∎

Corollary 3.9**.**

Suppose that $h_{opt},f^{*}_{opt},g^{*}_{opt}\in{\mathbb{R}}[t][\partial;^{\prime}]$ are a locally unique global minimum of $\Phi$ in some neighborhood around $h,f^{*}$ and $g^{*}$. If

[TABLE]

for $\varepsilon>0$, then

[TABLE]

If we compute different approximate GCRD pairs of $f$ and $g$ (using different optimization techniques or initial guesses), then we are able to bound the size of the perturbations of $f^{*},g^{*}$ and $h$ based on how near they are. Furthermore, this corollary allows us to certify an upper bound on the distance between our computed approximate GCRD tuple and the actual global minimum.

4 Implementation of Approximate GCRD

This section discusses the particulars and implementation of the algorithms. The algorithms are described in a Maple-like pseudo code, with Matlab style matrix indexing. All of the algorithms have been implemented in the Maple programming language. For convenience, the notation and assumptions introduced at the start of Section 3 will hold, unless otherwise stated. Additionally, we will assume that content from differential polynomials can be removed numerically, as computed quantities are typically not primitive due to round-off errors.

The matrices $S=S(f,g)\in{\mathbb{R}}[t]^{(M+N)\times(M+N)}$ will be the differential Sylvester matrix of $f$ and $g$, and ${\widehat{S}}={\widehat{S}}(f,g)\in{\mathbb{R}}^{(M+N)(\mu+1)\times(M+N)(\mu+d+1)}$ will be the inflated differential Sylvester matrix of $f$ and $g$, where $\mu=2(M+N)d$.

The presentation and theoretical analysis of the algorithms is presented in a bottom-up manner, reflecting their dependencies. Asymptotic upper bounds on the number of floating point operations required are provided. Furthermore, we discuss whether the output of the algorithm can be certified in some manner, when applicable.

We demonstrate the robustness of our algorithms in practice. Specific examples are provided to thoroughly demonstrate the steps of the algorithms. We investigate interesting families of input. In particular, we investigate exact inputs with an exact GCRD, and perturbed differential polynomials with varying errors and noise introduced. The test cases of differential polynomials of interest to us have

• •

low degree in $t$ and high degree in $\partial$ (unbalanced in $\partial$),

• •

high degree in $t$ and low degree in $\partial$ (unbalanced in $t$), and

• •

proportional degrees in $t$ and $\partial$ (balanced degrees).

4.1 Algorithms for Approximate GCRD

We adapt techniques from the exact setting to a numerical setting to compute an exact GCRD numerically. These algorithms compute the rank of the differential Sylvester matrix and a least squares solution to a polynomial linear system, corresponding to the Bézout coefficients. We describe an algorithm for finding nearby differential polynomials introduced in (Giesbrecht and Haraldson, 2014), whose (inflated) differential Sylvester matrix is nearly singular. Using the least squares numeric GCRD algorithm, we can compute an approximate GCRD candidate from the nearly singular differential Sylvester matrix. From this candidate, we extract a guess for the co-factors numerically and proceed with post-refinement Newton iteration.

4.1.1 Numerical Computation of a GCRD

Before we can compute a GCRD numerically, the rank of the differential Sylvester matrix needs to be determined. Our numeric rank algorithm is an adaptation of the rank algorithm used by Corless et al. (1995). There are

[TABLE]

trivial singular values222The inflated differential Sylvester matrix has more columns than rows, however the nullspace of the columns contains the information pertaining to the GCRD. The trivial singular values are the zero singular values occuring from there being more columns than rows., and $\mu/2<\mu+d+1$, the column block size. These trivial singular values need to be accounted for when annihilating small singular values. In the full rank case, we should not underestimate the rank of $S$ by inferring from ${\widehat{S}}$, as there are strictly fewer trivial singular values than the column block size.

Algorithm 1 computes a reasonable guess for the degree in $\partial$ of an approximate GCRD, although it is not generally certifiable. When ${\mbox{gcrd}}(f,g)$ is non-trivial (no errors present in the input coefficients), we compute (generically) the degree of the GCRD of $f$ and $g$. In the exact setting, we can now formulate a linear algebra problem over ${\mathbb{R}}[t]$ to compute a GCRD. We present two solutions to this problem. Algorithm 2 solves this problem using linear algebra over ${\mathbb{R}}(t)$. Algorithm 3 linearizes the problem over ${\mathbb{R}}$ and computes a least squares solution.

In the implementation of Algorithm 2, we take special care to ensure that $\operatorname{{\mbox{lcoeff}}}_{t}\operatorname{{\mbox{lcoeff}}}_{\partial}h$ does not vanish when $h$ is normalized. If $\operatorname{{\mbox{lcoeff}}}_{t}\operatorname{{\mbox{lcoeff}}}_{\partial}h$ vanishes, then this could be an indication that the input is ill-conditioned or content removal of $h$ failed. In either case, it is possible that this instance of the approximate GCRD problem will not have an attainable global minimum in accordance with Theorem 3.3.

4.1.2 Nearby Differential Polynomials with GCRD Algorithm

The matrix ${\widehat{S}}$ is highly structured, as it is composed of block Toeplitz matrices. When we consider the matrix ${\widehat{S}}+\Delta{\widehat{S}},$ the nearest (unstructured) matrix of prescribed rank deficiency, we have considerable flexibility in how we recover the coefficients of ${\widetilde{f}}$ and ${\widetilde{g}}$, nearby differential polynomials with an exact, non-trivial GCRD as in (3.1). However, the matrix ${\widehat{S}}+\Delta{\widehat{S}}$ is not generally an inflated differential Sylvester matrix, but it is probably reasonably close to one (see Giesbrecht and Haraldson (2014); Haraldson (2015)). We recall that the mapping $\Gamma:{\mathbb{R}}[t]\to{\mathbb{R}}^{(\mu+1)\times(\mu+d+1)}$ generates the (rectangular) Toeplitz blocks of ${\widehat{S}}$. To recover the coefficients of ${\widetilde{f}}$ and ${\widetilde{g}}$ one must make a suitable definition for the mapping $\Gamma^{-1}:{\mathbb{R}}^{(\mu+1)\times(\mu+d+1)}\to{\mathbb{R}}[t]$. We use $\Gamma^{-1}$ to find ${\widetilde{f}},{\widetilde{g}}\in{\mathbb{R}}[t][\partial;^{\prime}]$ such that ${\widehat{S}}({\widetilde{f}},{\widetilde{g}})\approx{\widehat{S}}(f,g)+\Delta{\widehat{S}}(f,g)$.

Regardless of our choice of $\Gamma^{-1}$, this method of recovering ${\widetilde{f}}$ and ${\widetilde{g}}$ can lead to a differential Sylvester matrix that does not have the desired numeric rank, as determined by Algorithm 1. The perturbation $\Delta{\widehat{S}}$ is unstructured while $\Gamma({\widetilde{f}}_{i})$ and $\Gamma({\widetilde{g}}_{j})$ are (highly structured) Toeplitz matrices. Consequently, some non-zero terms of $\Delta{\widehat{S}}$ are ignored.

4.1.3 Numeric Right Division

Numeric right division without remainder between two differential polynomials is a rational function linear algebra problem. The (approximate) quotient is a solution to a linear system, in a least squares sense. We present a naive algorithm that works well in practice and a more rigorous linear least squares variant.

The solution to this system may not be in (approximate) lowest terms. In our implementation we use approximate GCD and real linear least squares to resolve this. We note that total least squares can also be employed to prevent the need of an approximate GCD computation to put the rational function coefficients in lowest terms.

4.1.4 Improved GCRD via Optimization: Newton’s Method

Using Algorithm 5, we can compute an initial guess for an approximate GCRD, $h_{init}$. We can perform right division without remainder numerically to compute initial guesses for the co-factors, $f^{*}_{init}$ and $g^{*}_{init}$. We now have enough information to set up a post-refinement Newton iteration, to hopefully compute an approximate GCRD. When the co-factors have polynomial coefficients, the products $f^{*}h$ and $g^{*}h$ are always polynomial. This makes Newton iteration a very straightforward procedure, as the objective function

[TABLE]

is easily computed. However, when the co-factors have rational function coefficients, the quantities $f^{*}h$ and $g^{*}h$ usually have rational function coefficients due to round-off error. We can clear fractions and compute the least squares solution of an equivalent associate problem.

The normalization we impose, that $\operatorname{{\mbox{lcoeff}}}_{t}\operatorname{{\mbox{lcoeff}}}_{\partial}h$ is fixed, ensures the solution is (locally) unique, by Corollary 3.7. We note that this normalization can be changed. However one must ensure that the normalization vector is not orthogonal to $({\bf{f^{*}}},{\bf{g^{*}}},{\bf{h}})$. We now generalize the Newton iteration for the instance when the co-factors have rational function coefficients.

4.2 Analysis of Algorithms

In this section we assess the computational cost in terms of the number of floating point operations or flops. Where applicable, we discuss the numerical stability of the algorithms and whether or not their output can be certified. The algorithms are analyzed in the order they were presented. The assumption that content can be removed numerically is not without loss of generality; content removal can be unstable if implemented poorly.

In our implementation we remove content by (re)formulating our solutions as a solution to a (total) least squares problem. This can be done by performing the SVD on a generalized Sylvester matrix of several univariate polynomials (Kaltofen et al., 2006) to infer the degree of the content. Computing the degree of the content this way generalizes the method of Corless et al. (1995) to several polynomials. In our implementation the only important information is the degree of an approximate GCD, so we assume that the run-time of approximate GCD is cubic in the number of variables. One could compute an approximate GCD of several polynomials and perform a least squares division, however post-refinement would likely be needed. We generally assume that unstructured linear algebra techniques are used on the problems, however structured methods could lead to a modest asymptotic improvement.

4.2.1 Analysis of Algorithm 1 – DeflatedRank

The number of flops Algorithm 1 requires is dominated by the cost of performing the SVD on ${\widehat{S}}$. The SVD requires $O((M+N)^{3}(\mu+d+1)^{3})=O((M+N)^{6}d^{3})$ flops, using standard arithmetic. As mentioned earlier, this algorithm is generally not certified to produce the degree of an approximate GCRD.

4.2.2 Analysis of Algorithm 2 – NumericGCRD

The number of flops Algorithm 2 requires is ultimately bounded by the cost of computing the rank of $S$ using Algorithm 1. The cost of Algorithm 1 is $O((M+N)^{6}d^{3})$ flops. The cost of computing a GCRD given the degree in $\partial$ is $O((M+N)^{3})$ operations over ${\mathbb{R}}(t)$ which corresponds to $O((M+N)^{3}d^{2})$ flops. The cost of the approximate GCD and division to remove content depends on the specific method used, but is usually negligible when compared to the rank computation.

This algorithm is not numerically stable for large degree inputs in $t$ and $\partial$. Performing linear algebra over ${\mathbb{R}}(t)$ leads to considerable degree growth in $t$, and removing (approximate) content with a division further perturbs the coefficients of the GCRD. The output of this algorithm is not certified to be correct in most instances.

4.2.3 Analysis of Algorithm 3 – NumericGCRDviaLS

There are $(M+N)(\mu+d+1)$ equations and $(M+N)(\mu+1)+(D+1)(\mu+d+1)=O((M+N)^{2}d)$ unknowns. The cost of computing the least squares solution is $O((M+N)^{6}d^{3})$ flops. The cost of inferring the content by looking at singular values of the (generalized) Sylvester matrix is bounded by $O((D)^{3}(\mu+d+1)^{3})=O((M+N)^{6}d^{3})$ flops. The total number of flops required for this algorithm is $O((M+N)^{6}d^{3})$.

This algorithm relies on solving a real linear least squares problem. As such, this algorithm is numerically stable, provided that the underlying least squares problem is reasonably conditioned, and solved in a reasonable way. One such method of solving the least squares problem is the SVD and arising pseudo-inverse. We are able to certify the correctness of the answer obtained via least squares, provided that the underlying approximate GCD algorithm computes the degree of the content correctly.

4.2.4 Analysis of Algorithm 4 – DeflatedPerturbation

The number of flops Algorithm 4 requires is $O((M+N)^{2}d^{2})$, assuming that $\Gamma^{-1}$ uses the weighted block average. We use this in our implementation. This algorithm is not certified to provide meaningful output.

4.2.5 Analysis of Algorithm 5 – NearbyWithGCRD

The number of flops Algorithm 5 requires is dominated by the cost of computing the singular values of ${\widehat{S}}$, which is $O((M+N)^{6}d^{3})$ flops.

This algorithm is exactly the same as Algorithm 2 when ${\widehat{S}}$ has the desired rank deficiency. In the event that the input is approximate, the quality of our answer depends on the largest singular value of ${\widehat{S}}$ that we annihilate. This algorithm is not certified to provide meaningful output, but if used in conjunction with Algorithm 3, the output can be certified as a least squares approximation to the solution of the Bézout coefficients.

4.2.6 Analysis of Algorithm 6 – NaiveNumericRightDivision

The number of flops Algorithm 6 requires depends on the method used to solve the linear system. The particular system is highly structured so we can solve it by backwards substitution directly, which costs $O((M-D)^{2})$ operations over ${\mathbb{R}}(t)$. This corresponds to $O((M-D)^{2}d^{2})$ flops. An upper bound on the degree required for approximate GCD computations is $(M-D+1)d$. The total cost of each approximate GCD computation is at most $O(((M-D)d)^{3})$ flops. There are at most $M-D+1$ approximate GCD computations performed, so the total cost of the algorithm is $O((M-D)^{4}d^{3})$ flops.

The output of this answer is generally only certifiable if the residual of a least squares division is zero, i.e. the coefficients are exact. If we assume that $\operatorname{{\mbox{lcoeff}}}_{t}\operatorname{{\mbox{lcoeff}}}_{\partial}h=1$ and $\|h\|$ is not arbitrarily large, then the backwards substitution is well conditioned. The approximate GCD computations and following divisions can perturb the coefficients, so the algorithm can be unstable for poorly conditioned inputs. This is especially problematic when $\operatorname{{\mbox{lcoeff}}}_{\partial}h$ is poorly conditioned.

4.2.7 Analysis of Algorithm 7 –

NumericRightDivisionViaLS

If $f^{*}$ has polynomial coefficients, then $\deg_{t}f^{*}\leq\deg_{t}f\leq d$ as $f$ and $h$ have polynomial coefficients as well. If $f^{*}$ has rational function coefficients, we recall from Section 2.2 that there are $O(M(M-D)d)$ equations and $O((M-D)^{2}d)$ unknowns. The cost of solving this linear least squares problem is $O((M(M-D)d)^{3})\subseteq O(M^{6}d^{3})$ flops.

The output of this algorithm is certified as a linear least squares solution. Like Algorithm 6, the conditioning of this algorithm is strongly related to the conditioning of $h_{D}$.

4.2.8 Analysis of Algorithms

8–9 – NewtonIteration

and ModifiedNewtonIteration

We transform the problem of computing a GCRD to that of optimizing $\Phi:{\mathbb{R}}[t][\partial;^{\prime}]\times{\mathbb{R}}(t)[\partial;^{\prime}]^{2}\to{\mathbb{R}}$. We can assume without loss of generality that $f^{*}$ and $g^{*}$ have polynomial coefficients, as we can solve an equivalent associate problem instead. The dominating cost of the Newton iteration is solving a linear system to get the next value which requires $O(\nu^{3})$ operations, where $\nu$ is the number of variables needed to represent the coefficients of $h,f^{*}$ and $g^{*}$.

Newton iteration can fail for many reasons, (it is, afterall, a locally convergent method) however our Newton iteration usually fails because:

$\nabla^{2}\Phi$ is positive semidefinite at a point in the iteration, the stationary point is a saddle point;

• 2.

The initial guess is poorly chosen and $\nabla^{2}\Phi$ is indefinite at a point.

In the event that Newton iteration fails we can perform a Gauss-Newton iteration instead. Despite Gauss-Newton iteration having at least linear convergence, $J^{T}J$ is positive definite, so saddle points are no longer a problem if the optimal residual is sufficiently small. According to Corollary 3.7, if the residual is sufficiently small then Newton iteration will converge to a global minimum.

4.3 Examples and Experimental Results

This section contains some examples of our implementation.333A proof-of-concept implementation of the algorithms is available at https://www.scg.uwaterloo.ca/software/ApproxOreFoCM-2019.tgz. The (inflated) differential Sylvester matrix is ill-conditioned for large degree inputs in $t$ and $\partial$. This ill-conditioning occurs because the columns (rows) become unbalanced due to the falling factorials, where some columns have a Frobenius norm factorially larger than others. We restrict ourselves to modest examples with minimal coefficient growth. Computations are done using the default precision in Maple, which is approximately 10 decimal points of accuracy.

Example 4.1** (No Noise, many factors).**

[TABLE]

We compute initial guesses (removing content numerically where appropriate):

[TABLE]

The quality of this initial guess is

[TABLE]

The condition number for the Hessian matrix evaluated at our initial guess is $18354.38336$ and our smallest eigenvalue is $.00314$. Since $\nabla^{2}\Phi$ is locally positive definite, we know that we will converge to a unique (local) minimum. The minimum we converge to is $2.33030\times 10^{-20}$.

The exact GCRD in this example is $h=(\partial+4t-1)(\partial-1)(\partial-1)$.

Example 4.2** (Noise).**

In this example we introduced normalized noise of size $10^{-5}$ to $f$ and $g$.

[TABLE]

We compute initial guesses (removing content numerically where appropriate):

[TABLE]

The quality of this initial guess is

[TABLE]

The condition number for the Hessian matrix evaluated at our initial guess is $21971.20356$ and our smallest eigenvalue is $.00818$. Since $\nabla^{2}\Phi$ is locally positive definite, we know that we will converge to a unique (local) minimum. The minimum we converge to is $1.06759\times 10^{-10}$, which is roughly the amount of noise we added.

Example 4.3** (GCRD via LS).**

In this example we added a noise factor of $10^{-4}$ to $f$ and $g$. Performing Linear Algebra over ${\mathbb{R}}(t)$ produced completely unacceptable answers, so we used a Least Squares algorithm to compute an approximate GCRD.

[TABLE]

Using a Least Squares variant of our Numeric GCRD algorithm, we are able to compute (without removing content):

[TABLE]

The quality of this initial guess is

[TABLE]

The condition number for the Hessian matrix evaluated at our initial guess is $148.62547$ and our smallest eigenvalue is $.04615$. Since $\nabla^{2}\Phi$ is locally positive definite, we know that we will converge to a unique (local) minimum. The minimum we converge to is $9.53931\times 10^{-9}$.

4.4 General Examples

We provide results that demonstrate the robustness of our algorithms. We consider differential polynomials whose degrees in $t$ and $\partial$ are balanced and unbalanced. The coefficients of the inputs were generated using the Maple routine randpoly(). The inputs $f$ and $g$ were normalized so that $\|f\|=\|g\|=1$. We introduced normalized noise to the coefficients of $f$ and $g$, so that the relative error is size of the perturbation. Precisely, if $f+\Delta f$ and $g+\Delta g$ are perturbed from $f$ and $g$ by the quantities $\Delta f$ and $\Delta g$, then the relative error in the coefficients of $f$ and $g$ is given by $\|\Delta f\|_{2}=\|\Delta g\|_{2}$.

We recall that the Newton iteration optimizes $\|f-{\widetilde{f}}\|_{2}^{2}+\|g-{\widetilde{g}}\|_{2}^{2}$, which is the sum of the squares of the errors. The initial error and error from post-refinement are expressed as the sum of square errors accordingly. In our experiments, the Initial Error is the quantity $\|f-f_{init}\|_{2}^{2}+\|g-g_{init}\|_{2}^{2}$ and the error after post-refinement, Newton Error is the quantity $\|f-f_{opt}\|_{2}^{2}+\|g-g_{opt}\|_{2}^{2}$. In all of the examples when there were no perturbations in the coefficients of $f$ and $g$, our numeric GCRD algorithm and post-refinement procedures were able to compute an exact GCRD to machine precision. When perturbations imposing a relative error of $10^{-8}$ in the coefficients of $f$ and $g$ were introduced, we were able to compute a solution to the approximate GCRD problem in every example.

Introducing perturbations imposing a relative error of order $10^{-4}$ and $10^{-2}$ into the coefficients of $f$ and $g$ prevented computation of an approximate GCRD in some examples. Instead, we provide examples of the largest perturbation in the coefficients of $f$ and $g$ that we were able to compute an approximate GCRD. Instances that are denoted as “FAIL” occur when the post-refinement did not converge. The implementation of Newton’s method is not globalized to converge to a stationary point, hence the iterates may diverge. In our examples iterates diverge because the Hessian matrix is indefinite at an initial guess.

4.4.1 Balanced Degrees in $t$ and $\partial$

The following results of experiments were conducted on differential polynomials whose degrees in $t$ and $\partial$ were proportional, or balanced.

[TABLE]

4.4.2 Unbalanced Degrees in $\partial$

The following results of experiments were conducted on differential polynomials whose degrees in $\partial$ were relatively larger than their degree in $t$.

[TABLE]

4.4.3 Unbalanced Degrees in $t$

The following results of experiments were conducted on differential polynomials whose degrees in $t$ were relatively larger than their degree in $\partial$.

[TABLE]

5 Conclusion

In this paper we have formally defined an approximate GCRD problem for differential polynomials, and given an approach to a robust numerical solution. We have seen that, under reasonable assumptions the approximate GCRD problem is well posed. In particular, we show that Newton iteration will converge to an optimal solution if the residual is sufficiently small. We employ the earlier results in (Giesbrecht and Haraldson, 2014), analogous to SVD-based approximate GCD methods like Corless et al. (1995), to compute a reasonable initial estimate for the Newton iteration. The results were presented for real differential polynomials, however the results generalize in a very straight forward way to the instance of complex differential polynomials.

We believe that some aspects of our problems could also be approached from a structured low-rank approximation viewpoint Kaltofen et al. (2005); Schost and Spaenlehauer (2016). In particular, the work of Schost and Spaenlehauer (2016) can be used to obtain an initial low-rank differential Sylvester matrix in which co-factors and a GCRD can be extracted for post-refinement. This holds more generally than differential polynomials, and a particular example to consider is the shift operator, commonly associated with linear difference equations.

Another area of future work is in the certification of the degree of an approximate GCRD. We can obtain a reasonable guess by enumerating over the degrees of all possible approximate GCRDs, similar to the Structured Total Least Norm approach adopted for multivariate polynomial approximate GCD Kaltofen et al. (2006). A possible direction would be to look at the differential subresultant sequence and the singular values of their inflated block matrices Emiris et al. (1997).

The differential polynomials defined in this paper are special case of more general Ore polynomials, which have broader application in the solution of differential and difference equations. In particular, we could potentially apply our methods in the context of $q$-differentiation (Jackson differentiation) or derivations on exponential polynomials. Ultimately, any Ore structure will have a well-defined Sylvester-like matrix (see, e.g., Giesbrecht and Kim (2013)). However, the numerical properties of different derivations may well be quite difficult or even problematic, and may well introduce poles or other significant sources of numerical instability.

We also hope, the results of this paper are a foundation for extending the approximate polynomial toolbox to other problems with differential polynomials and more general linear differential operators. Much like approximate GCD, the approximate GCRD is both a stepping stone and a key tool towards operations like approximate factorization and (functional) solution of differential polynomials. More immediately, computation of an approximate GCRD enables computation of a corresponding approximate LCLM, and multiple GCRD’s, and to multiple differential variables (i.e., iterated Ore polynomials), which provide an effective method for dealing with linear PDEs.

Acknowledgements

The authors would like to thank George Labahn for his comments. The authors would also like to thank the two anonymous referees for their careful reading and comments.