---
abstract: 'We describe how to solve simultaneous [Padé]{}approximations over a power series ring ${{\mathsf{K}}}[[x]]$ for a field ${{\mathsf{K}}}$ using ${O^{\scriptscriptstyle \sim}\!}(n^{\omega - 1} d)$ operations in ${{\mathsf{K}}}$, where $d$ is the sought precision and $n$ is the number of power series to approximate. We develop two algorithms using different approaches. Both algorithms return a reduced sub-bases that generates the complete set of solutions to the input approximations problem that satisfy the given degree constraints. Our results are made possible by recent breakthroughs in fast computations of minimal approximant bases and Hermite [Padé]{}approximations.'
author:
- |
    Johan Rosenkilde, né Nielsen\
    Technical University of Denmark\
    Denmark\
    jsrn@jsrn.dk
- |
    Arne Storjohann\
    University of Waterloo\
    Canada\
    astorjoh@uwaterloo.ca
- |
    Johan Rosenkilde, né Nielsen\
    \
    \
    Arne Storjohann\
    \
    \
bibliography:
- 'bibtex.bib'
title: 'Algorithms for Simultaneous [Padé]{}Approximations [^1] '
---

Introduction {#sec:intro}
============

The Simultaneous [Padé]{}approximation problem concerns approximating several power series $S_1,\ldots,S_n \in {{\mathsf{K}}}[[x]]$ with rational functions $\frac {\sigma_1}\lambda,\ldots,\frac {\sigma_n} \lambda$, all sharing the same denominator $\lambda$. In other words, for some $d \in {\mathbb Z\xspace}_{\geq 0}$, we seek $\lambda \in {{\mathsf{K}}}[x]$ of low degree such that each of $${{\textnormal{rem}}}(\lambda S_1,\  x^d) , {{\textnormal{rem}}}(\lambda S_2,\  x^d),\ \ldots,\ {{\textnormal{rem}}}(\lambda S_n,\  x^d)$$ has low degree. The study of Simultaneous [Padé]{}approximations traces back to Hermite’s proof of the transcendence of $e$ [@hermite_sur_1878]. Solving Simultaneous [Padé]{}approximations has numerous applications, such as in coding theory, e.g. [@feng_generalization_1991; @schmidt_collaborative_2009]; or in distributed, reliable computation [@clement_pernet_high_2014]. Many algorithms have been developed for this problem, see e.g. [@beckermann_uniform_1992; @olesh_vector_2006; @sidorenko_linear_2011; @nielsen_generalised_2013] as well as the references therein. Usually one cares about the regime where $d \gg n$. Obtaining $O(n d^2)$ is classical through successive cancellation, see [@beckermann_uniform_1994] or [@feng_generalization_1991] for a Berlekamp–Massey-type variant. Using fast arithmetic, the previous best was ${O^{\scriptscriptstyle \sim}\!}(n^\omega d)$, where $\omega$ is the exponent for matrix multiplication, see \[ssec:cost\]. That can be done by computing a minimal approximant basis with e.g. [@giorgi_complexity_2003; @GuptaSarkarStorjohannValeriote11]; this approach traces back to [@barel_general_1992; @beckermann_uniform_1992]. Another possibility which achieves the same complexity is fast algorithms for solving structured linear systems, e.g. [@bostan_solving_2008]; see [@chowdhury_faster_2015] for a discussion of this approach.

A common description is to require $\deg \lambda < N_0$ for some degree bound $N_0$, and similarly $\deg {{\textnormal{rem}}}(\lambda S_1,\, x^d) < N_i$ for $i = 1,\ldots,n$. The degree bounds could arise naturally from the application, or could be set such that a solution must exist. A natural generalisation is also to replace the $x^d$ moduli with arbitrary $g_1,\ldots,g_n \in {{\mathsf{K}}}[x]$. Formally, for any field ${{\mathsf{K}}}$:

\[prob:sim\_pade\] Given a tuple $({\bm{S}}, {\bm{g}}, {\bm{N}})$ where

-   ${\bm{S}} = (S_1,\ldots,S_n) \in {{\mathsf{K}}}[x]^n$ is a sequence of polynomials,

-   ${\bm{g}} = (g_1,\ldots,g_n) \in {{\mathsf{K}}}[x]^n$ is a sequence of moduli polynomials with $\deg S_i < \deg g_i$ for $i=1,\ldots,n$,

-   and ${\bm{N}} = (N_0,\ldots,N_n)
     \in {\mathbb Z\xspace}_{\geq 0}^{n+1}$ are degree bounds satisfying $1\leq N_0 \leq \max_i \deg g_i$ and $N_i \leq \deg g_i$ for $i=1,\ldots,n$,

find, if it exists, a non-zero vector $(\lambda, \phi_1, \ldots, \phi_n)$ such that

1.  $\lambda S_i \equiv \phi_i \mod g_i$ for $i = 1,\ldots, n$, and \[p1item1\]

2.  $\deg \lambda < N_0$ and $\deg \phi_i < N_i$ for $i=1,\ldots,n$.

We will call any vector $(\lambda, \phi_1, \ldots, \phi_n)$ as above *a solution* to a given Simultaneous [Padé]{}approximation problem. Note that if the $N_i$ are set too low, then it might be the case that no solution exists.

\[ex:simpade\] Consider over ${\mathbb F_{2}\xspace}[x]$ that $g_1 = g_2 = g_3 = x^5$, and ${\bm{S}} = (S_1,S_2,S_3) = 
      \left(x^{4} + x^{2} + 1,\,x^{4} + 1,\,x^{4} + x^{3} + 1\right)$, with degree bounds ${\bm{N}} = (5, 3, 4, 5)$. Then $\lambda_1 = x^4 + 1$ is a solution, since $\deg \lambda_1 < 5$ and $$\lambda_1 {\bm{S}} \equiv
      \left(x^{2} + 1,\ 1,\ x^{3} + 1\right)
      \mod x^5 \ .$$ $\lambda_2 = x^{3} + x$ is another solution, since $$\lambda_2 {\bm{S}} \equiv
      \left(x,\ x^{3} + x,\ x^4+x^3 + x\right)
      \mod x^5 \ .$$ These two solutions are linearly independent over ${\mathbb F_{2}\xspace}[x]$ and span all solutions.

Several previous algorithms for solving \[prob:sim\_pade\] are more ambitious and produce an entire *basis* of solutions that satisfy the first output condition $\lambda S_i \equiv \phi_i \mod g_i$ for $i=1,\ldots,n$, including solutions that do not satisfy the degree bounds stipulated by the second output condition. Our algorithms are slightly more restricted in that we only return the sub-basis that generates the set of solutions that satisfy both output requirements of \[prob:sim\_pade\]. Formally:

\[prob:sim\_pade\_basis\] Given an instance of \[prob:sim\_pade\], find a matrix $A \in {{\mathsf{K}}}[x]^{* \times (n+1)}$ such that:

-   Each row of $A$ is a solution to the instance.

-   All solutions are in the ${{\mathsf{K}}}[x]$-row space of $A$.

-   $A$ is $(-{\bm{N}})$-row reduced[^2].

The last condition ensures that $A$ is minimal, in a sense, according to the degree bounds ${\bm{N}}$, and that we can easily parametrise which linear combinations of the rows of $A$ are solutions. We recall the relevant definitions and lemmas in \[sec:preliminaries\].

We will call such a matrix $A$ a *solution basis*. In the complexities we report here, we cannot afford to compute $A$ explicitly. For example, if all $g_i = x^d$, the number of field elements required to explicitly write down all of the entries of $A$ could be $\Omega(n^2d)$. Instead, we remark that $A$ is completely given by the problem instance as well as the first column of $A$, containing the $\lambda$ polynomials.[^3] Our algorithms will therefore represent $A$ row-wise using the following compact representation.

For a given instance of \[prob:sim\_pade\_basis\], a *solution specification* is a tuple $({\bm{\lambda}},{\bm{\delta}}) \in
{{\mathsf{K}}}[x]^{k \times 1} \times {\mathbb Z\xspace}_{<0}^k$ such that the *completion* of ${\bm{\lambda}}$ is a solution basis, and where ${\bm{\delta}}$ are the $(-{\bm{N}})$-degrees of the rows of $A$.

The *completion* of ${\bm{\lambda}} = (\lambda_1,\ldots,\lambda_k){^\top}$ is the matrix $$\begin{bmatrix}
    \lambda_1 & {{\textnormal{rem}}}(\lambda_1 S_1,\ g_1) & \ldots & {{\textnormal{rem}}}(\lambda_1 S_n,\ g_n) \\
    \vdots & & \ddots & \vdots \\
    \lambda_k & {{\textnormal{rem}}}(\lambda_k S_1,\ g_1) & \ldots & {{\textnormal{rem}}}(\lambda_k S_n,\ g_n) \\
  \end{bmatrix}
  \ .$$

Note that ${\bm{\delta}}$ will consist of only negative numbers, since any solution ${\bm{v}}$ by definition has $\deg_{-{\bm{N}}} {\bm{v}} < 0$.

A solution specification for the problem in \[ex:simpade\] is $$({\bm{\lambda}}, {\bm{\delta}}) = \big( [x^4 + 1,\ x^3 + x]{^\top},\ (-1, -1) \big)  \ .$$ The completion of this is $$A = \begin{bmatrix}
           x^4 + 1 & x^{2} + 1 & 1 & x^{3} + 1 \\
           x^3 + x & x  & x^{3} + x & x^4+x^3 + x
         \end{bmatrix}$$ One can verify that $A$ is $(-{\bm{N}})$-row reduced.

We present two algorithms for solving \[prob:sim\_pade\_basis\], both with complexity $O\big(n^{\omega-1}\, {{\mathsf{M}}}(d)\,(\log d)\,(\log
d/n)^2\big)$, where $d = \max_i \deg g_i$ and ${{\mathsf{M}}}(d)$ is the cost of multiplying two polynomials of degree $d$, see \[ssec:cost\]. They both depend crucially on recent developments that allow computing minimal approximant bases of non-square matrices faster than for the square case [@zhou_efficient_2012; @jeannerod_computation_2016]. We remark that from the solution basis, one can also compute the expanded form of one or a few of the solutions in the same complexity, for instance if a single, expanded solution to the simultaneous [Padé]{}problem is needed.

Our first algorithm in \[sec:dual\] assumes $g_i = x^d$ for all $i$ and some $d \in {\mathbb Z\xspace}_{\geq 0}$. It utilises a well-known duality between Simultaneous [Padé]{}approximations and Hermite [Padé]{}approximations, see e.g. [@beckermann_uniform_1992]. The Hermite [Padé]{}problem is immediately solvable by fast minimal approximant basis computation. A remaining step is to efficiently compute a single row of the adjoint of a matrix in Popov form, and this is done by combining partial linearisation [@GuptaSarkarStorjohannValeriote11] and high-order lifting [@storjohann_high-order_2003].

Our second algorithm in \[sec:intersect\] supports arbitrary $g_i$. The algorithm first solves $n$ single-sequence [Padé]{}approximations, each of $S_1,\ldots,S_n$. The solution bases for two problem instances can be combined by computing the intersection of their row spaces; this is handled by a minimal approximant basis computation. A solution basis of the full Simultaneous [Padé]{}problem is then obtained by structuring intersections along a binary tree.

Before we describe our algorithms, we give some preliminary notation and definitions in \[sec:preliminaries\], and in \[sec:subroutines\] we describe some of the computational tools that we employ.

Both our algorithms have been implemented in Sage v. 7.0 [@stein_sagemath_????] (though asymptotically slower alternatives to the computational tools are used). The source code can be downloaded from <http://jsrn.dk/code-for-articles>.

Cost model {#ssec:cost}
----------

We count basic arithmetic operations in ${{\mathsf{K}}}$ on an algebraic RAM. We will state complexity results in terms of an exponent $\omega$ for matrix multiplication, and a function ${{\mathsf{M}}}(\cdot)$ that is a multiplication time for ${{\mathsf{K}}}[x]$ [@von_zur_gathen_modern_2012 Definition 8.26]. Then two $n\times n$ matrices over ${{\mathsf{K}}}$ can be multiplied in $O(n^{\omega})$ operations in ${{\mathsf{K}}}$, and two polynomials in ${{\mathsf{K}}}[x]$ of degree strictly less than $d$ can be multiplied in ${{\mathsf{M}}}(d)$ operations in ${{\mathsf{K}}}$. The best known algorithms allow $\omega < 2.38$ [@coppersmith_matrix_1990; @LeGall14], and we can always take ${{\mathsf{M}}}(d) \in O(n (\log n) (\operatorname{loglog}n))$ [@CantorKaltofen].

In this paper we assume that $\omega > 2$, and that ${{\mathsf{M}}}(d)$ is super-linear while ${{\mathsf{M}}}(d) \in O(d^{\omega-1})$. The assumption ${{\mathsf{M}}}(d) \in O(d^{\omega-1})$ simply stipulates that if fast matrix multiplication techiques are used then fast polynomial multiplication should be used also: for example, $n \, {{\mathsf{M}}}(nd) \in O(n^{\omega} \, {{\mathsf{M}}}(d))$.

Preliminaries {#sec:preliminaries}
=============

Here we gather together some definitions and results regarding row reduced bases, minimal approximant basis, and their shifted variants. For a matrix $A$ we denote by $A_{i,j}$ the entry in row $i$ and column $j$. For a matrix $A$ over ${{\mathsf{K}}}[x]$ we denote by ${{\textnormal{Row}}}(A)$ the ${{\mathsf{K}}}[x]$-linear row space of $A$.

Degrees and shifted degrees
---------------------------

The degree of a nonzero vector ${\bm{v}} \in {{\mathsf{K}}}[x]^{1 \times m}$ or matrix $A \in {{\mathsf{K}}}[x]^{n\times m}$ is denoted by $\deg {\bm{v}}$ or $\deg A$, and is the maximal degree of entries of ${\bm{v}}$ or $A$. If $A$ has no zero rows the [*row degrees*]{} of $A$, denoted by ${{\textnormal{rowdeg}}}\, A$, is the tuple $(d_1,\ldots,d_n)$ with $d_i = \deg
{{\textnormal{row}}}(A,i)$.

The (row-wise) [*leading matrix*]{} of $A$, denoted by ${\rm LM}(A) \in {{\mathsf{K}}}^{n \times m}$, has ${\rm LM}(A)_{i,j}$ equal to the coefficient of $x^{d_i}$ of $A_{i,j}$.

Next we recall [@barel_general_1992; @zhou_efficient_2012; @jeannerod_computation_2016] the shifted variants of the notion of degree, row degrees, and leading matrix. For a [*shift*]{} ${\bm{s}} =
(s_1,\ldots,s_n) \in {\mathbb Z\xspace}^n$, define the $n \times n$ diagonal matrix $x^{{\bm{s}}}$ by $$x^{{\bm{s}}} := \left [ \begin{array}{ccc} x^{s_1}
& & \\
 &  \ddots & \\ & & x^{s_n} \end{array} \right ].$$ Then the [*${{\bm{s}}}$-degree*]{} of $v$, the [*${{\bm{s}}}$-row degrees*]{} of $A$, and the [*${\bm{s}}$-leading matrix*]{} of $A$, are defined by $\deg_{{\bm{s}}} v := \deg v x^{{\bm{s}}}$, ${{\textnormal{rowdeg}}}_{{\bm{s}}} A := {{\textnormal{rowdeg}}}\, Ax^{{\bm{s}}}$, and ${\rm LM}_{{\bm{s}}}(A) := {\rm
LM}(Ax^{{\bm{s}}})$. Note that we pass over the ring of Laurent polynomials only for convenience; our algorithms will only compute with polynomials. As pointed out in [@jeannerod_computation_2016], up to negation the definition of ${{\bm{s}}}$-degree is equivalent to that used in [@BeckermannLabahnVillard06] and to the notion of [*defect*]{} in [@beckermann_uniform_1994].

For an instance $({\bm{S}}, {\bm{g}}, {\bm{N}})$ of \[prob:sim\_pade\], in the context of defining matrices, we will be using ${\bm{S}}$ and ${\bm{g}}$ as vectors, and by ${  \ifx\null\null  \Gamma_{{\bm{g}}}  \else  \Gamma_{{\bm{g}}_{\null}}  \fi}$ denote the diagonal matrix with the entries of ${\bm{g}}$ on its diagonal.

Row reduced
-----------

Although row reducedness can be defined for matrices of arbitrary shape and rank, it suffices here to consider the case of matrices of full row rank. A matrix $R \in {{\mathsf{K}}}[x]^{n \times m}$ is [*row reduced*]{} if ${\rm LM}(R)$ has full row rank, and [*${\bm{s}}$-row reduced*]{} if ${\rm LM}_{{\bm{s}}}(R)$ has full row rank. Every $A \in {{\mathsf{K}}}[x]^{n \times m}$ of full row rank is left equivalent to a matrix $R \in {{\mathsf{K}}}[x]^{n \times m}$ that is ${{\bm{s}}}$-row reduced. The rows of $R$ give a basis for ${{\textnormal{Row}}}(A)$ that is minimal in the following sense: the list of ${{\bm{s}}}$-degrees of the rows of $R$, when sorted in non-decreasing order, will be lexicographically minimal. An important feature of row reduced matrices is the so-called “predictable degree”-property [@kailath_linear_1980 Theorem 6.3-13]: for any ${\bm{v}} \in {{\mathsf{K}}}[x]^{1 \times n}$, we have $$\deg_{{\bm{s}}}({\bm{v}} R) = \max_{i=1,\ldots,n}( \deg_{{\bm{s}}} {\rm row}(R,i)
+ \deg v_i ) \ .$$

A canonical ${\bm{s}}$-reduced basis is provided by the ${{\bm{s}}}$-Popov form. Although an ${{\bm{s}}}$-Popov form can be defined for a matrix of arbitrary shape and rank, it suffices here to consider the case of a non-singular matrix. The following definition is equivalent to [@jeannerod_computation_2016 Definition 1.2].

\[def:popov\] A non-singular matrix $R \in {{\mathsf{K}}}[x]^{n\times n}$ is in ${{\bm{s}}}$-Popov form if ${\rm LM}_{{\bm{s}}}(R)$ is unit lower triangular and the degrees of off-diagonal entries of $R$ are strictly less than the degree of the diagonal entry in the same column.

Adjoints of row reduced matrices
--------------------------------

For a non-singular matrix $A$ recall that the adjoint of $A$, denoted by ${\rm adj}(A)$, is equal to $(\det A)A^{-1}$, and that entry ${{{\textnormal{adj}}}(A)}_{i,j}{^\top}$ is equal to $(-1)^{i+j}$ times the determinant of the $(n-1) \times (n-1)$ sub-matrix that is obtained from $A$ by deleting row $i$ and column $j$.

\[lem:adjointRowReduced\] Let $A \in {{\mathsf{K}}}[x]^{n \times n}$ be ${\bm{s}}$-row reduced. Then ${{\textnormal{adj}}}(A){^\top}$ is $(-{\bm{s}})$-row reduced with $${{\textnormal{rowdeg}}}_{(-{\bm{s}})} {{\textnormal{adj}}}(A){^\top}=(\eta - s - \eta_1,\ldots , \eta - s -\eta_n) \ ,$$ where ${\bm{\eta}} = {{\textnormal{rowdeg}}}_{{\bm{s}}} A$, $\eta = \sum_i \eta_i$ and $s = \sum_i s_i$.

Since $A$ is ${\bm{s}}$-row reduced then $A x^{{\bm{s}}}$ is row reduced. Note that ${{\textnormal{adj}}}(A x^{{\bm{s}}}){^\top}(A x^{{\bm{s}}}){^\top}= (\det A x^{{\bm{s}}}) I_m$ with $\deg
\det A x^{{\bm{s}}} = \eta$. It follows that row $i$ of ${{\textnormal{adj}}}(A x^{{\bm{s}}}){^\top}$ must have degree at least $\eta - \eta_i$ since $\eta_i$ is the degree of column $i$ of $(A x^{{\bm{s}}}){^\top}$. However, entries in row $i$ of ${{\textnormal{adj}}}(A x^{{\bm{s}}}){^\top}$ are minors of the matrix obtained from $A x^{{\bm{s}}}$ by removing row $i$, hence have degree at most $\eta
- \eta_i$. It follows that the (row-wise) leading coefficient matrix of ${{\textnormal{adj}}}(A x^{{\bm{s}}}){^\top}$ is non-singular, hence ${{\textnormal{adj}}}(A x^{{\bm{s}}}){^\top}$ is row reduced. Since ${{\textnormal{adj}}}(A x^{{\bm{s}}}){^\top}=  
(\det x^{{\bm{s}}}) {{\textnormal{adj}}}(A){^\top}x^{-{\bm{s}}}$ we conclude that ${{\textnormal{adj}}}(A){^\top}$ is $(-{\bm{s}})$-row reduced with ${{\textnormal{rowdeg}}}_{(-{\bm{s}})} {{\textnormal{adj}}}(A) =
(\eta - \eta_1 - s, \ldots, \eta - \eta_n - s)$.

Minimal approximant bases
-------------------------

We recall the standard notion of minimal approximant basis, sometimes known as order basis or $\sigma$-basis [@beckermann_uniform_1994]. For a matrix $A \in {{\mathsf{K}}}[x]^{n \times m}$ and order $d \in
{\mathbb Z\xspace}_{\geq 0}$, an *order $d$ approximant* is a vector ${\bm{p}} \in
{{\mathsf{K}}}[x]^{1 \times n}$ such that ${\bm{p}}A \equiv {\bm{0}} \mod x^d.$ An *approximant basis of order $d$* is then a matrix $F \in
{{\mathsf{K}}}[x]^{n \times n}$ which is a basis of all order $d$ approximants. Such a basis always exists and has full rank $n$. For a shift ${\bm{s}} \in {\mathbb Z\xspace}^n$, $F$ is then an *${\bm{s}}$-minimal approximant basis* if it is ${\bm{s}}$-row reduced.

Let ${{\ensuremath{\mathsf{MinBasis}}\xspace}}(d,A,{\bm{s}})$ be a function that returns $(F,{\bm{\delta}})$, where $F$ is an ${\bm{s}}$-minimal approximant basis of $A$ of order $d$, and ${\bm{\delta}}  = {{\textnormal{rowdeg}}}_{{\bm{s}}} F$. The next lemma recalls a well known method of constructing minimal approximant bases recursively. Although the output of ${{\ensuremath{\mathsf{MinBasis}}\xspace}}$ may not be unique, the lemma holds for *any* ${\bm{s}}$-minimal approximant basis that ${{\ensuremath{\mathsf{MinBasis}}\xspace}}$ might return.

\[lem:paderec\] Let $A = \left [ \begin{array}{c|c}
A_1 & A_2 \end{array} \right ]$ over ${{\mathsf{K}}}[x]$. If $(F_1, {\bm{\delta}}_1) = {{\ensuremath{\mathsf{MinBasis}}\xspace}}(d,A_1,{\bm{s}})$ and $(F_2,{\bm{\delta}}_2) =
{{\ensuremath{\mathsf{MinBasis}}\xspace}}(d,F_1A_2,{\bm{\delta}}_1)$, then $F_2F_1$ is an ${\bm{s}}$-minimal approximant basis of $A$ of order $d$ with ${\bm{\delta}}_2 = {{\textnormal{rowdeg}}}_{{\bm{s}}} F_2 F_1$.

Sometimes only the [*negative part*]{} of an ${\bm{s}}$-minimal approximant bases is required, the submatrix of the approximant bases consisting of rows with negative ${\bm{s}}$-degree. Let function ${{\ensuremath{\mathsf{NegMinBasis}}\xspace}}(d,A,{\bm{s}})$ have the same output as ${{\ensuremath{\mathsf{MinBasis}}\xspace}}$, but with $F$ restricted to the negative part.

\[lem:paderecprune\] \[lem:paderec\] still holds if ${{\ensuremath{\mathsf{MinBasis}}\xspace}}$ is replaced by ${{\ensuremath{\mathsf{NegMinBasis}}\xspace}}$, and “an ${\bm{s}}$-minimal” is replaced with “the negative part of an ${\bm{s}}$-minimal.”

Using for example the algorithm `M-Basis` of [@giorgi_complexity_2003], it is easy to show that any order $d$ approximant basis $G$ for an $A$ of column dimension $m$ has $\det G = x^D$ for some $D \in {\mathbb Z\xspace}_{\geq 0}$ with $D \leq md$.

Many problems of ${{\mathsf{K}}}[x]$ matrices or approximations reduce to the computation of (shifted) minimal approximant bases, see e.g. [@beckermann_uniform_1994; @giorgi_complexity_2003], often resulting in the best known asymptotic complexities for these problems.

Direct solving of Simultaneous [Padé]{}approximations {#sec:direct_solve}
-----------------------------------------------------

Let $({\bm{S}}, {\bm{g}}, {\bm{N}})$ be an instance of \[prob:sim\_pade\_basis\] of size $n$. We recall some known approaches for computing a solution specification using row reduction and minimal approximant basis computation.

### Via reduced basis {#sec:direct_reduced_basis}

Using the predictable degree property it is easy to show that if $R \in {{\mathsf{K}}}[x]^{(n+1)
\times (n+1)}$ is an $(-{\bm{N}})$-reduced basis of $$A = 
 \left [ \begin{array}{c|c} 
  1 & {\bm{S}} \\\hline
 & {  \ifx\null\null  \Gamma_{{\bm{g}}}  \else  \Gamma_{{\bm{g}}_{\null}}  \fi}\end{array} \right]
       \in {{\mathsf{K}}}[x]^{(n+1) \times (n+1)},$$ then the sub-matrix of $R$ comprised of the rows with negative $(-{\bm{N}})$-degree form a solution basis. A solution specification $({\bm{\lambda}}, {\bm{\delta}})$ is then a subvector ${\bm{\lambda}}$ of the first column of $R$, with ${\bm{\delta}}$ the corresponding subtuple ${\bm{\delta}}$ of ${{\textnormal{rowdeg}}}_{(- {\bm{N}})} R$.

Mulders and Storjohann [@mulders_lattice_2003] gave an iterative algorithm for performing row reduction by successive cancellation; it is similar to but faster than earlier algorithms [@kailath_linear_1980; @lenstra_factoring_1985]. Generically on input $F \in {{\mathsf{K}}}[x]^{m \times m}$ it has complexity $O(n^3 (\deg F)^2)$. Alekhnovich [@alekhnovich_linear_2005] gave what is essentially a Divide & Conquer variant of Mulders and Storjohann’s algorithm, with complexity ${O^{\scriptscriptstyle \sim}\!}(n^{\omega+1}\deg F)$. Nielsen remarked [@nielsen_generalised_2013] that these algorithms perform fewer iterations when applied to the matrix $A$ above, due to its low *orthogonality defect*: ${\rm OD}(F) = \sum{{\textnormal{rowdeg}}}F - \deg \det F$, resulting in $O(n^2(\deg A)^2)$ respectively ${O^{\scriptscriptstyle \sim}\!}(n^\omega \deg A)$. Nielsen also used the special shape of $A$ to give a variant of the Mulders–Storjohann algorithm that computes coefficients in the working matrix in a lazy manner with a resulting complexity $O(n \,\mathsf{P}(\deg A))$, where $\mathsf P(\deg A) = (\deg A)^2$ when the $g_i$ are all powers of $x$, and $\mathsf P(\deg A) = {{\mathsf{M}}}(\deg A)\deg A$ otherwise.

Giorgi, et al. [@giorgi_complexity_2003] gave a reduction for performing row reduction by computing a minimal approximant basis. For the special matrix $A$, this essentially boils down to the approach described in the following section.

When $n = 1$, the extended Euclidean algorithm on input $S_1$ and $g_1$ can solve the approximation problem by essentially computing the reduced basis of the $2 \times 2$ matrix $A$: each iteration corresponds to a reduced basis for a range of possible shifts [@sugiyama_further_1976; @justesen_complexity_1976; @gustavson_fast_1979]. The complexity of this is $O({{\mathsf{M}}}(\deg g_1) \log \deg g_1)$.

### Via minimal approximant basis

First consider the special case when all $g_i = x^d$ for the same $d$. An approximant ${\bm{v}} = (\lambda, \phi_1, \ldots, \phi_n)$ of order $d$ of $$\begin{aligned}
   A &=  \left [ \begin{array}{c}
 - {\bm{S}} \\
 I
            \end{array} \right ]
       \in {{\mathsf{K}}}[x]^{(n+1) \times n}\end{aligned}$$ clearly satisfies $\lambda S_i \equiv \phi_i \mod x^d$ for $i =
1,\ldots,n$; conversely, any such vector ${\bm{v}}$ satisfying these congruences must be an approximant of $A$ of order $d$. So the negative part of a $(-{\bm{N}})$-minimal approximant basis of $A$ of order $d$ is a solution basis.

In the general case we can reduce to a minimal approximant bases computation as shown by \[alg:simpadedirect\]. Correctness of the algorithm follows from the following result.

\[thm:simpadedirect\] Corresponding to an instance $({\bm{S}}, {\bm{g}}, {\bm{N}})$ of \[prob:sim\_pade\_basis\] of size $n$, define a shift ${\bm{h}}$ and order $d$:

-   ${\bm{h}} := -({\bm{N}} \mid N_0 -1, \ldots, N_0-1) \in {\mathbb Z\xspace}^{2n+1}$

-   $d := N_0 + \max_i \deg g_i -1$

If $G$ is the negative part of an ${\bm{h}}$-minimal approximant basis of $$H = \left [ \begin{array}{c} -{\bm{S}} \\
 I \\
 {  \ifx\null\null  \Gamma_{{\bm{g}}}  \else  \Gamma_{{\bm{g}}_{\null}}  \fi}\end{array} \right ] \in {{\mathsf{K}}}[x]^{(2n+1) \times n}$$ of order $d$, then the submatrix of $G$ comprised of the first $n+1$ columns is a solution basis to the problem instance.

An approximant ${\bm{v}} = (\lambda, \phi_1,\ldots,\phi_n,q_1,\ldots, q_n)$ of order $d$ of $H$ clearly satisfies $$\begin{aligned}
 \label{eq:lambda} 
\lambda S_i = \phi_i + q_ig_i \bmod
x^d\end{aligned}$$ for $i=1,\ldots,n$; conversely, any such vector ${\bm{v}}$ satisfying these congruences must be an approximant of $H$ of order $d$.

Now suppose ${\bm{v}}$ is an order $d$ approximant of $H$ with negative ${\bm{h}}$-degree, so $\deg \lambda \leq N_0-1$, $\deg \phi_i \leq
N_i-1$, and $\deg q_i \leq  N_0 - 2$. Since \[prob:sim\_pade\] specifies that $\deg S_i < \deg g_i$ and $N_i \leq \deg g_i$, both $\lambda S_i$ and $q_i g_i$ will have degree bounded by $N_0 + \deg g_i - 2$. Since \[prob:sim\_pade\] specifies that $N_0 \geq 1$, it follows that both the left and right hand sides of (\[eq:lambda\]) have degree bounded by $N_0+ \deg g_i -2$, which is strictly less than $d$. We conclude that $$\begin{aligned}
 \label{eq:lambda2}
\lambda S_i = \phi_i + q_i g_i\end{aligned}$$ for $i=1,\ldots,n$. It follows that ${\bm{v}} H = 0$ so ${\bm{v}}$ is in the left kernel of $H$. Moreover, restricting ${\bm{v}}$ to its first $n+1$ entries gives $\bar{{\bm{v}}} := (\lambda, \phi_1,\ldots,\phi_n)$, a solution to the simultaneous [Padé]{}problem with $\deg_{- {\bm{N}}} \bar{{\bm{v}}} = \deg_{{\bm{h}}} {\bm{v}}$. Conversely, if $\bar{{\bm{v}}} = (\lambda, \phi_1,\ldots,\phi_n)$ is a solution to the simultaneous [Padé]{}problem, then the extension ${\bm{v}} = (\lambda, \phi_1,\ldots,\phi_n,q_1,\ldots,q_n)$ with $q_i
= (\lambda S_i - \phi_i)/g_i \in {{\mathsf{K}}}[x]$ for $i=1,\ldots,n$ is an approximant of $H$ of order $d$ with $\deg_{{\bm{h}}} {\bm{v}} =
\deg_{-{\bm{N}}} \bar{{\bm{v}}}$.

Finally, consider that a left kernel basis for $H$ is given by $$K = \left[ \begin{array}{c|c} K_1 & K_2 \end{array} \right ] = 
\left [ \begin{array}{cc|c} 1 & {\bm{S}} & \\
 & {  \ifx\null\null  \Gamma_{{\bm{g}}}  \else  \Gamma_{{\bm{g}}_{\null}}  \fi}& -I \end{array} \right ].$$ We must have $G = M K$ for some polynomial matrix $M$ of full row rank. But then $M K_1$ also has full row rank with ${{\textnormal{rowdeg}}}_{-{\bm{N}}} MK_1
= {{\textnormal{rowdeg}}}_{{\bm{h}}} G$.

${\bm{h}} {\leftarrow}-( {\bm{N}} \mid  N_0-1,\ldots,N_0-1) \in {\mathbb Z\xspace}^{2n+1}$ $d {\leftarrow}N_0 + \max_i \deg g_i - 1$ $H = 
  \left[ \begin{array}{c}
       -{\bm{S}} \\
  I \\
  {  \ifx\null\null  \Gamma_{{\bm{g}}}  \else  \Gamma_{{\bm{g}}_{\null}}  \fi}\end{array} \right]$ $(\left [ \begin{array}{c|c} {\bm{\lambda}} & \ast \end{array} 
  \right ], {\bm{\delta}}) {\leftarrow}{{\ensuremath{\mathsf{NegMinBasis}}\xspace}}(d, H, {\bm{h}})$ ${\textbf{return }}({\bm{\lambda}}, {\bm{\delta}})$

[$\mathsf{DirectSimPade}$]{} can be performed in time ${O^{\scriptscriptstyle \sim}\!}(n^{\omega} \deg H) = {O^{\scriptscriptstyle \sim}\!}(n^\omega \max_i \deg g_i)$ using the minimal approximant basis algorithm by Jeannerod, et al. [@jeannerod_computation_2016], see \[sec:subroutines\].

A closely related alternative to [$\mathsf{DirectSimPade}$]{} is the recent algorithm by Neiger [@neiger_fast_2016] for computing solutions to modular equations with general moduli $g_i$. This would give the complexity ${O^{\scriptscriptstyle \sim}\!}(n^{\omega-1} \sum_i \deg g_i) \subset {O^{\scriptscriptstyle \sim}\!}(n^\omega \max_i \deg g_i)$.

All of the above solutions ignore the sparse, simple structure of the input matrices, which is why they do not obtain the improved complexity that we do here.

Computational tools {#sec:subroutines}
===================

The main computational tool we will use is the following very recent result from Jeannerod, Neiger, Schost and Villard [@jeannerod_computation_2016] on minimal approximant basis computation.

\[thm:orderbasis\] There exists an algorithm ${{\ensuremath{\mathsf{PopovBasis}}\xspace}}(d,A, {\bm{s}})$ where the input is an order $d \in {\mathbb Z\xspace}_+$, a polynomial matrix $A \in {{\mathsf{K}}}[x]^{n
  \times m}$ of degree at most $d$, and shift ${\bm{s}} \in {\mathbb Z\xspace}^n$, and which returns $(F, {\bm{\delta}})$, where $F$ is an ${\bm{s}}$-minimal approximant basis of $A$ of order $d$, $F$ is in ${\bm{s}}$-Popov form, and ${\bm{\delta}} = {{\textnormal{rowdeg}}}_{{\bm{s}}} F$. [[$\mathsf{PopovBasis}$]{}]{}has complexity $O(n^{\omega-1}\, {{\mathsf{M}}}(\sigma)\, (\log \sigma) \, (\log \sigma /n)^2)$ operations in ${{\mathsf{K}}}$, where $\sigma = md$.

Our next result says that we can quickly compute the first row of ${{\textnormal{adj}}}(F)$ if $F$ is a minimal approximant basis in Popov form. In particular, since $F$ is an approximant basis $\det F = x^D$ for some $D \leq \sigma$, where $\sigma = md$ from \[thm:orderbasis\].

\[thm:fastsolver\] Let $F \in {{\mathsf{K}}}[x]^{n \times n}$ be in Popov form and with $\det F = x^D$ for some $D \in {\mathbb Z\xspace}_{\geq 0}$. Then the first row of ${{\textnormal{adj}}}(F)$ can be computed in $O(n^{\omega-1}\, {{\mathsf{M}}}(D)\, (\log D) \, (\log D/n))$ operations in ${{\mathsf{K}}}$.

Because $F$ is in ${\bm{s}}$-Popov form, $D$ is the sum of the column degrees of $F$. We consider two cases: $D \geq n$ and $D < n$.

First suppose $D
\geq n$. Partial linearisation [@GuptaSarkarStorjohannValeriote11 Corollary 2] can produce from $F$, with no operations in ${{\mathsf{K}}}$, a new matrix $G \in {{\mathsf{K}}}[x]^{\bar n \times \bar n}$ with dimension $\bar{n} < 2n$, $\deg G \leq \lceil D/n\rceil$, $\det G = \det F$, and such that $F^{{-1}}$ is equal to the principal $n \times n$ sub-matrix of $G^{{-1}}$. Let ${\bm{v}} \in {{\mathsf{K}}}[x]^{1 \times \bar{n}}$ be the first row of $x^DI_{\bar{n}}$. Then the first row of ${{\textnormal{adj}}}(F)$ will be the first $n$ entries of the first row of ${\bm{v}}G^{-1}$. High-order $X$-adic lifting [@storjohann_high-order_2003 Algorithm 5] using the modulus $X=(x-1)^{\lceil
D/n \rceil}$ will compute ${\bm{v}}G^{-1}$ in $O\big(n^{\omega}\,
{{\mathsf{M}}}(\lceil D/n \rceil) \,(\log \lceil D/n \rceil)\big)$ operations in ${{\mathsf{K}}}$ [@storjohann_high-order_2003 Corollary 16]. Since $D \geq n$ this cost estimate remains valid if we replace $\lceil D/n \rceil$ with $D/n$. Finally, from the super-linearity assumption on ${{\mathsf{M}}}(\cdot)$ we have $M(D/n) \leq (1/n) {{\mathsf{M}}}(D)$, thus matching our target cost.

Now suppose $D < n$. In this case we can not directly appeal to the partial linearisation technique since the resulting $O(n^{\omega} \lceil D/n\rceil)$ may be asymptotically larger than our target cost. But $D < n$ means that $F$ has — possibly many — columns of degree 0; since $F$ is in Popov form, such columns have a 1 on the matrix’s diagonal and are 0 on the remaining entries. The following describes how to essentially ignore those columns. $D$ is then greater than or equal to the number of remaining columns, thus effectuating the gain from the partial linearisation.

If $n-k$ is the number of such columns in $F$ that means we can find a permutation matrix $P$ such that $$\hat{F} := PFP{^\top}= \left [ \begin{array}{c|c}
                             F_1 & \\\hline
                             F_2 & I_{n-k}
                           \end{array} \right ] \ ,$$ with each column of $F_1$ having degree strictly greater than zero. Let $i$ be the row index of the single 1 in the first column of $P{^\top}$. Since $F^{-1} = P{^\top}\hat{F}^{-1}P$, we have $$\label{first}
{\rm row}({{\textnormal{adj}}}(F),1)P^{-1} = x^D\, {\rm row}(\hat{F}^{-1},i).$$ Considering that $$\hat{F}^{-1} = \left [ \begin{array}{c|c} F_1^{-1} & \\\hline -F_2F_1^{-1} & I_{n-k} \end{array} \right ],$$ it will suffice to compute the first $k$ entries of the vector on the right hand side of (\[first\]). If $i \leq k$ then let ${\bm{v}} \in {{\mathsf{K}}}[x]^{1 \times k}$ be row $i$ of $x^{D}I_k$. Otherwise, if $i>k$ then let ${\bm{v}}$ be row $i-k$ of $-x^{D}F_2$. Then in both cases, ${\bm{v}}F_1^{-1}$ will be equal to the first $k$ entries of the vector on the right hand side of (\[first\]). Like before, high-order lifting combined with partial linearisation will compute this vector in $O\big(k^{\omega}\, {{\mathsf{M}}}(\lceil D/k \rceil)\,(\log \lceil D/k \rceil)
\big)$ operations in ${{\mathsf{K}}}$. Since $D\geq k$ the cost estimate remains valid if $\lceil D/k \rceil$ is replaced with $D/k$.

Reduction to Hermite [PADÉ]{} {#sec:dual}
=============================

In this section we present an algorithm for solving \[prob:sim\_pade\_basis\] when $g_1 = \ldots = g_n = x^d$ for some $d \in {\mathbb Z\xspace}_{\geq 0}$. The algorithm is based on the well-known duality between the Simultaneous [Padé]{}problem and the Hermite [Padé]{}problem, see for example [@beckermann_uniform_1992]. This duality, first observed in a special case [@Mahler68], and then later in the general case [@beckermann_recursiveness_1997], was exploited in [@beckermann_fraction-free_2009] to develop algorithms for the fraction free computation of Simultaneous [Padé]{}approximation. We begin with a technical lemma that is at the heart of this duality.

\[lem:duality\] Let $\hat A, \hat B \in {{\mathsf{K}}}[x]^{(n+1)\times(n+1)}$ be as follows. $$\begin{aligned}
    \hat A &= 
         \left [ \begin{array}{c|c}
                   x^d & -{\bm{S}} \\\hline
                       & I
                 \end{array} \right ]
    &&
    \hspace*{-1em}\hat B &= 
           \left [ \begin{array}{c|cccc}
                     1 & \\{  \hline  \\[\dimexpr-\normalbaselineskip+2pt]  }
                     {\bm{S}}{^\top}& x^d I
                   \end{array} \right ]
  \end{aligned}$$ Then $\hat B$ is the adjoint of $\hat A{^\top}$. Furthermore, $\hat A{^\top}$ is an approximant basis for $\hat
  B$ of order $d$, and $\hat B{^\top}$ is an approximant basis of $\hat A$ of order $d$.

Direct computation shows that $\hat A{^\top}\hat B = x^d I_m =
\det \hat A{^\top}I_m$, so $\hat B$ is the adjoint of $\hat
A{^\top}$.

Let now $G$ be an approximant basis of $\hat B$. By the above computation the row space of $\hat A{^\top}$ must be a subset of the row space of $G$. But since $G \hat B = (x^dI_m) R$ for some $R \in {{\mathsf{K}}}[x]^{(n+1)\times(n+1)}$, then $\det G = x^d \det R$. Thus $x^d \mid \det G$. But $\det \hat A{^\top}= x^d$, so the row space of $\hat A{^\top}$ can not be smaller than the row space of $G$. That is, $\hat A{^\top}$ is an approximant basis for $B$ of order $d$. Taking the transpose through the argument shows that $\hat
B{^\top}$ is an approximant basis of $\hat B$ of order $d$.

\[thm:dualityMinbasis\] Let $A$ and $B$ be as follows. $$\begin{aligned}
    A &= \left [ \begin{array}{cccc}
                   -{\bm{S}} \\{  \hline  \\[\dimexpr-\normalbaselineskip+1pt]  }
                   I
                 \end{array} \right ] \in {{\mathsf{K}}}[x]^{(n+1) \times (n+1)}
    &&&
    \hspace*{-1em}B  &= \left [ \begin{array}{c}
                    1 \\ {\bm{S}}
                  \end{array} \right] \in {{\mathsf{K}}}[x]^{(n+1) \times 1}\end{aligned}$$ If $G$ is an ${\bm{N}}$-minimal approximant basis of $B$ of order $d$ with shift ${\bm{N}} \in {\mathbb Z\xspace}_{\geq 0}^{n+1}$, then ${{\textnormal{adj}}}(G{^\top})$ is a $(-{\bm{N}})$-minimal approximant basis of $A$ of order $d$. Moreover, if ${\bm{\eta}} = {{\textnormal{rowdeg}}}_{{\bm{N}}} G$, then ${{\textnormal{rowdeg}}}_{-{\bm{N}}} {{\textnormal{adj}}}(G) =(\eta - N - \eta_1,\ldots , \eta - N
-\eta_{n+1})$, where $\eta = \sum_i \eta_i$ and $N = \sum_i N_i$.

Introduce $\hat A$ and $\hat B$ as in \[lem:duality\]. Clearly $G$ is also an ${\bm{N}}$-minimal approximant basis of $\hat B$ of order $d$. Likewise, $\hat A$ and $A$ have the same minimal approximant bases for given order and shift.

Assume, without loss of generality, that we have scaled $G$ such that $\det G$ is monic. Since $\hat A{^\top}$ is also an approximant basis for $\hat B$ of order $d$, then $\det G = \det \hat A{^\top}= x^d$. By definition $G\hat B = x^d R$ for some matrix $R \in
  {{\mathsf{K}}}[x]^{(n+1)\times(n+1)}$. That means $$\begin{aligned}
    x^{2d}((G\hat B){^\top}))^{{-1}}&= x^{2d}((x^d R){^\top})^{{-1}}\ , & \textrm{so} \\ 
    (x^d(G{^\top})^{{-1}})(x^d(\hat B{^\top})^{{-1}}) &= x^{d}(R{^\top})^{{-1}}\ , & \textrm{that is} \\ 
    {{\textnormal{adj}}}(G{^\top}) \hat A &= x^d (R{^\top})^{{-1}}\ .
  \end{aligned}$$ Now $\det R = 1$ since $(x^d)^{n+1} \det R = \det(G\hat B) =
x^{d+nd}$, so $(R{^\top})^{{-1}}= {{\textnormal{adj}}}(R{^\top}) \in {{\mathsf{K}}}[x]^{(n+1)
\times (n+1)}$. Therefore ${{\textnormal{adj}}}(G{^\top})$ is an approximant basis of $\hat A$ of order $d$. The theorem now follows from \[lem:adjointRowReduced\] by noting that $G$ is ${\bm{N}}$-row reduced.

We apply \[thm:dualityMinbasis\] to the problem of \[ex:simpade\] with shifts ${\bm{N}} = (5, 3, 4, 5)$. We have $$\begin{aligned}
    A &= 
      \left[\begin{array}{rrr}
      x^{4} + x^{2} + 1 & x^{4} + 1 & x^{4} + x^{3} + 1 \\
      1 &   &   \\
        & 1 &   \\
        &   & 1
      \end{array}\right]
    \\
    B &= \left[\begin{array}{r}
             1 \\ x^4 + x^2 + 1 \\ x^{4} + 1 \\ x^{4} + x^{3} + 1
         \end{array}\right]
  \end{aligned}$$ An ${\bm{N}}$-minimal approximant basis to order $d = 5$ of $B$ is $$\begin{aligned}
    G &= 
  \left[\begin{array}{rrrr}
  x & 0 & x & 0 \\
  1 & x^{2} + 1 & 0 & 0 \\
  0 & 1 & x^{2} + 1 & 0 \\
  0 & x & x + 1 & 1
  \end{array}\right] ,  \textrm{ and}
  \\
  {{\textnormal{adj}}}(G){^\top}&= 
    \left[\begin{array}{rrrr}
    x^{4} + 1 & x^{2} + 1 & 1 & x^{3} + 1 \\
    x & x^{3} + x & x & x^{4} + x \\
    x^{3} + x & x & x^{3} + x & x^{4} + x^{3} + x \\
    0 & 0 & 0 & x^{5}
    \end{array}\right]
    \ .
  \end{aligned}$$ ${{\textnormal{adj}}}(G){^\top}$ can be confirmed to be an $(-{\bm{N}})$-minimal approximant basis of $A$, since ${{\textnormal{adj}}}(G){^\top}A \equiv 0 \mod x^d$, and since the $(-{\bm{N}})$-leading coefficient matrix of ${{\textnormal{adj}}}(G){^\top}$ has full rank.

Algorithm \[alg:simpadedual\] uses \[thm:dualityMinbasis\] to solve a Simultaneous [Padé]{}approximation by computing a minimal approximant basis of $B$ in Popov form.

$B {\leftarrow}[ 1, S_1, \ldots, S_n ]^T \in {{\mathsf{K}}}[x]^{(n+1) \times 1}$ $G {\leftarrow}{{\ensuremath{\mathsf{PopovBasis}}\xspace}}(d, B, {\bm{N}})$ \[line:simpadedual:basis\] ${\bm{\eta}} {\leftarrow}{{\textnormal{rowdeg}}}_{{\bm{N}}} G$ $\hat{{\bm{\lambda}}} {\leftarrow}$ first column of ${{\textnormal{adj}}}(G{^\top})$ \[line:simpadedual:firstcol\] $\hat{ {\bm{\delta}}} {\leftarrow}(\eta - N -\eta_1, \ldots, \eta - N - \eta_{n+1})$, where $\eta = \sum_i \eta_i$ and $N = \sum_i N_i$ \[line:simpadedual:degrees\] $I {\leftarrow}\{ i \mid \hat{{\bm{\delta}}}_i < 0 \}$, and $k {\leftarrow}|I|$ $({\bm{\lambda}},\ {\bm{\delta}}) {\leftarrow}\big( \hat{{\bm{\lambda}}}_{i \in I},\ (\hat{{\bm{\delta}}}_i)_{i \in I} 
\big) \in {{\mathsf{K}}}[x]^{k \times 1} \times {\mathbb Z\xspace}^{k} $ ${\textbf{return }}({\bm{\lambda}}, {\bm{\delta}})$

\[alg:simpadedual\] is correct. The cost of the algorithm is $O(n^{\omega-1}\, {{\mathsf{M}}}(d) (\log d) (\log d/n)^2)$ operations in ${{\mathsf{K}}}$.

Correctness follows from \[thm:dualityMinbasis\]. The complexity estimate is achieved if the algorithms supporting \[thm:orderbasis\] and \[thm:fastsolver\] are used for the computation in lines 2 and 4, respectively.

A Divide & Conquer algorithm {#sec:intersect}
============================

Our second algorithm can handle the full generality of \[prob:sim\_pade\_basis\]. It works by first solving $n$ single [Padé]{}approximations, one for each of the $S_i$ individually, and then intersecting these solutions to form approximations of multiple $S_i$ simultaneously. The intersection is structured in a Divide & Conquer tree, and performed by computing minimal approximant bases. Let $({\bm{S}}, {\bm{g}}, {\bm{N}})$ be an instance of \[prob:sim\_pade\_basis\] of size $n$.

The idea of the intersection algorithm is the following: consider that we have solution specifications for two different Simultaneous [Padé]{}problems, $({\bm{\lambda}}_1, {\bm{\delta}}_1)$ and $({\bm{\lambda}}_2, {\bm{\delta}}_2)$. We then compute an approximant basis $G$ of the following matrix: $$\label{eqn:intersect_R}
    R =
    \left[\begin{array}{@{}c|c@{}}
            1 & 1 \\ \hline
            -{\bm{\lambda}}_1 &        \\\hline
                   & -{\bm{\lambda}}_2 \\
    \end{array}\right]$$ $G$ then encodes the *intersection* of the ${{\mathsf{K}}}[x]$-linear combinations of the ${\bm{\lambda}}_1$ with the ${{\mathsf{K}}}[x]$-linear combinations of the ${\bm{\lambda}}_2$: any $\lambda \in {{\mathsf{K}}}[x]$ residing in both sets of polynomials will appear as the first entry of a vector in the row space of $G$. We compute $G$ as an ${\bm{r}}$-minimal approximant basis to high enough order, where ${\bm{r}}$ is selected carefully such that the ${\bm{r}}$-degree of any $(\lambda \mid \ldots) \in {{\textnormal{Row}}}(G)$ will equal the $(-{\bm{N}})$-degree of the completion of $\lambda$ according to the combined Simultaneous [Padé]{}problem, whenever this degree is negative. From those rows of $G$ with negative ${\bm{r}}$-degree we then get a solution specification for the combined problem.

Consider again \[ex:simpade\]. We divide the problem into two sub-problems ${\bm{S}}_1 = (S_1, S_2)$, ${\bm{N}}_1 = (5, 3, 4)$, and ${\bm{S}}_2 = (S_3)$ and ${\bm{N}}_2 = (5, 5)$. Note that $N_{1,0} = N_{2,0} = 5$, since this is the degree bound on the sought $\lambda$ for the combined problem. The sub-problems have the following solution specifications and their completions: $$\begin{aligned}
    ({\bm{\lambda}}_1, {\bm{\delta}}_1) &= \big( [ x^{4} + 1,\ x^3 + x ]{^\top},\ ( -1, -1 ) \big)
    \\
      A_1 &=
        \left(\begin{array}{rrr}
        x^{4} + 1 & x^{2} + 1 & 1 \\
        x^{3} + x & x & x^{3} + x
        \end{array}\right)
    \\
    ({\bm{\lambda}}_2, {\bm{\delta}}_2) &= \big( [ x^2,\ x^3 + x + 1 ]{^\top},\ ( -3, -2 ) \big)
    \\
      A_2 &=
        \left(\begin{array}{rr}
        x^{2} & x^{2} \\
        x^{3} + x + 1 & x + 1
        \end{array}\right)
  \end{aligned}$$ We construct $R$ as in , and compute $G$, a minimal approximant basis of $R$ of order $7$ and with shifts ${\bm{r}}= (-5 \mid -1, -1 \mid -3, -2)$ (the $G$ below is actually in ${\bm{r}}$-Popov form): $$G =  \left(\begin{array}{rrrrr}
    x^{8} & 0 & 0 & 0 & 0 \\
    x^{3} + x + 1 & x^{4} + 1 & 1 & 0 & 1 \\
    x^{3} + x^{2} + x + 1 & 1 & x + 1 & 1 & 1 \\
    x^{4} + x^{3} + x + 1 & 1 & 1 & x^{2} & 1 \\
    x^{4} + 1 & 1 & 0 & x + 1 & x + 1
    \end{array}\right)$$ $G$ has ${\bm{r}}$-row degree $(3, 3, 0, -1, -1)$. Only rows 4 and 5 have negative ${\bm{r}}$-degree, and their first entries are the linearly independent solutions $x^4 + x^3 + x + 1$ and $x^4 + 1$. Both solutions complete into vectors with $(-{\bm{N}})$-degree -1.

To prove the correctness of the above intuition, we will use \[alg:simpadedirect\] ([$\mathsf{DirectSimPade}$]{}). The following lemma says that to solve two simultaneous [Padé]{}approximations, one can compute a minimal approximant basis of one big matrix $A$ constructed essentially from two of the matrices employed in [$\mathsf{DirectSimPade}$]{}. Afterwards, \[lem:recursive\_R\_solves\] uses this to show that a minimal approximant basis of $R$ in provides the crucial information in a minimal approximant basis of $A$.

${\textbf{return }}{\ensuremath{\mathsf{DirectSimPade}}\xspace}({\bm{S}}, {\bm{g}}, {\bm{N}})$ ${\bm{S}}_1, {\bm{g}}_1 {\leftarrow}$ the first ${\lceil n/2 \rceil}$ elements of ${\bm{S}}, {\bm{g}}$ ${\bm{S}}_2, {\bm{g}}_2 {\leftarrow}$ the last ${\lfloor n/2 \rfloor}$ elements of ${\bm{S}}, {\bm{g}}$ ${\bm{N}}_1
         {\leftarrow}(N_0,N_1,\ldots,N_{{\lceil n/2 \rceil}})$ ${\bm{N}}_2  
         {\leftarrow}(N_0,N_{{\lceil n/2 \rceil}+1},\ldots,N_n)$ $({\bm{\lambda}}_1, {\bm{\delta}}_1) {\leftarrow}{\ensuremath{\mathsf{RecursiveSimPade}}\xspace}\big({\bm{S}}_1, {\bm{g}}_1, {\bm{N}}_1)$ $({\bm{\lambda}}_2, {\bm{\delta}}_2) {\leftarrow}{\ensuremath{\mathsf{RecursiveSimPade}}\xspace}
   \big({\bm{S}}_2, {\bm{g}}_2, {\bm{N}}_2)$ ${\bm{r}} {\leftarrow}(-N_0 \mid {\bm{\delta}}_1 \mid {\bm{\delta}}_2)$ $d {\leftarrow}N_0 + \max_i \deg g_i - 1$ \[line:recursive:choosed\] \[lineH\] $R {\leftarrow}\left[\begin{array}{c|c}
            1 & 1\\ \hline
            -{\bm{\lambda}}_1 &  \\\hline
              & -{\bm{\lambda}}_2 
    \end{array}\right]$ $(\left [\begin{array}{c|c} {\bm{\lambda}} & \ast \end{array} \right ],
 {\bm{\delta}}) {\leftarrow}{{\ensuremath{\mathsf{NegMinBasis}}\xspace}}(d, R, {\bm{r}})$ \[calltoNegMin\] ${\textbf{return }}({\bm{\lambda}}, {\bm{\delta}})$

\[lem:recursive\_big\_matrix\] Let $({\bm{S}}_1, {\bm{g}}_1,{\bm{N}}_1)$ and $({\bm{S}}_2, {\bm{g}}_2,{\bm{N}}_2)$ be two instances of \[prob:sim\_pade\_basis\] of lengths $n_1, n_2$ respectively, and where ${\bm{N}}_1 = (N_0 \mid \grave{{\bm{N}}_1})$ and ${\bm{N}}_2 = (N_0 \mid \grave{{\bm{N}}_2})$. Let ${\bm{S}} = ({\bm{S}}_1 \mid {\bm{S}}_2)$, ${\bm{g}} = ({\bm{g}}_1 \mid {\bm{g}}_2)$ and ${\bm{N}} = (N_0 \mid \grave{{\bm{N}}_1} \mid \grave{{\bm{N}}_2})$ be the combined problem having length $n = n_1 + n_2$.

Let ${\bm{h}}_i = (-{\bm{N}}_i \mid N_0 -1 \ldots N_0 -1) \in {\mathbb Z\xspace}^{2n_i+1}$ for $i=1,2$. Let $(F, {\bm{\delta}}) = {{\ensuremath{\mathsf{NegMinBasis}}\xspace}}(d, A, {\bm{a}})$, where $A$ of dimension $(2n+3) \times (n+2)$ is given as: $$A =  \left [ \begin{array}{c|c}
                   A_1          & A_2
                 \end{array} \right ]
      = \left [ \begin{array}{cc|cc}
                                &                 & 1 & 1  \\
                -{\bm{S}}_1       &                 & -1     \\
                I               &                 &        \\
                {  \ifx\null1  \Gamma_{{\bm{g}}}  \else  \Gamma_{{\bm{g}}_{1}}  \fi} &                 &        \\
                                & -{\bm{S}}_2       &   & -1 \\
                                & I               &        \\
                                & {  \ifx\null2  \Gamma_{{\bm{g}}}  \else  \Gamma_{{\bm{g}}_{2}}  \fi} & 
                \end{array} \right ] ,$$ with ${\bm{a}} = (- N_0 \mid {\bm{h}}_1 \mid {\bm{h}}_2)$ and $d = N_0 + \max_i \deg g_i - 1$. Then $({\bm{\lambda}}, {\bm{\delta}})$ is a solution specification to $({\bm{S}}, {\bm{g}}, {\bm{N}})$, where ${\bm{\lambda}}$ is the first column of $F$.

Note that the matrix $A$ is right equivalent to the following matrix $B$: $$B := A 
      \left [ \begin{array}{cccc} 
           &    & I               &           \\
           &    &                 & I         \\
        1  &    & {\bm{S}}_1        &           \\
           & 1  &                 & {\bm{S}}_2
      \end{array} \right ] = 
    \left [ \begin{array}{cc|cc}
       1   & 1  & - {\bm{S}}_1      & -{\bm{S}}_2 \\
        -1 &    &                 &           \\
           &    & I               &           \\
           &    & {  \ifx\null1  \Gamma_{{\bm{g}}}  \else  \Gamma_{{\bm{g}}_{1}}  \fi} &           \\
           & -1 &                 &           \\
           &    &                 & I         \\
           &    &                 & {  \ifx\null2  \Gamma_{{\bm{g}}}  \else  \Gamma_{{\bm{g}}_{2}}  \fi}
    \end{array} \right ].$$ Since $F$ is an ${{\bm{a}}}$-minimal approximant of $A$ of order $d$, then it will also be one for $B$. Let $P$ be the permutation matrix that produces the following matrix $C := P B$: $$\setlength{\arraycolsep}{.8\arraycolsep}
  C = P B = 
    \left [ \begin{array}{cc|cc}
    1  & 1  & - {\bm{S}}_1      & -{\bm{S}}_2       \\
       &    & I               &                 \\ 
       &    &                 & I               \\
       &    & {  \ifx\null1  \Gamma_{{\bm{g}}}  \else  \Gamma_{{\bm{g}}_{1}}  \fi} &                 \\
       &    &                 & {  \ifx\null2  \Gamma_{{\bm{g}}}  \else  \Gamma_{{\bm{g}}_{2}}  \fi} \\\hline
    -1 &    &                 &                 \\
       & -1 &                 & 
    \end{array}  \right ]  =
    \left [ \begin{array}{cc|c}
    1  & 1  & - {\bm{S}}                          \\     &  & I \\ &  & {  \ifx\null\null  \Gamma_{{\bm{g}}}  \else  \Gamma_{{\bm{g}}_{\null}}  \fi}\\\hline  
    -1 &    &                                   \\
       & -1 & 
    \end{array}  \right ].$$ Define ${\bm{c}} := {\bm{a}} P^{-1}$, and note that ${\bm{c}} = ({\bm{h}} \mid -N_0, -N_0)$. Since $F = {{\ensuremath{\mathsf{NegMinBasis}}\xspace}}(d, A, {\bm{a}})$, then $(FP^{-1}, {\bm{\delta}})$ is a valid output of ${{\ensuremath{\mathsf{NegMinBasis}}\xspace}}(d, C, {\bm{c}})$. Furthermore, since the first column of $P$ is $(1, 0, \ldots, 0)$, the first column of $F$ will be equal to the first column of $FP^{{-1}}$.

We are therefore finished if we can show that if $(F', {\bm{\delta}}')$ is any valid output of ${{\ensuremath{\mathsf{NegMinBasis}}\xspace}}(d, C, {\bm{c}})$, then the first column of $F'$ together with ${\bm{\delta}}'$ form a solution specification to $({\bm{S}}, {\bm{g}}, {\bm{N}})$.

Consider therefore such an $(F', {\bm{\delta}}')$. By the first two columns of $C$, we must have $F'_{*,1} \equiv F'_{*,2n+2} \equiv F'_{*,2n+3} \mod x^d$, where $F'_{*,i}$ denotes the $i$’th column of $F'$. Since each row of $F'$ have negative ${\bm{c}}$-degree, and since $N_0 < d$, then the congruences must lift to equalities. We can therefore write $F = [ G \mid F'_{*,1} \mid F'_{*,1} ]$ for some $G \in {{\mathsf{K}}}[x]^{k \times (2n+1)}$ for some $k$, and we have ${{\textnormal{rowdeg}}}_{{\bm{h}}} G = {{\textnormal{rowdeg}}}_{{\bm{c}}} F' = {\bm{\delta}}'$.

By the last $n$ columns of $C$, we have $G H \equiv 0 \mod x^d$, where $$H =
    \left[\begin{array}{c}
            -{\bm{S}} \\
            I \\
            {  \ifx\null\null  \Gamma_{{\bm{g}}}  \else  \Gamma_{{\bm{g}}_{\null}}  \fi}\end{array}\right] \ .$$ In fact, $(G, {\bm{\delta}}')$ is a valid output for ${{\ensuremath{\mathsf{NegMinBasis}}\xspace}}(d, H, {\bm{h}})$: for $G$ has full row rank since $F'$ does; $G$ is ${\bm{h}}$-row reduced since $F'$ is ${\bm{c}}$-row reduced; and any negative ${\bm{h}}$-order $d$ approximant of $H$ must clearly be in the span of $G$ since $F'$ is a negative ${\bm{c}}$-minimal approximant basis of $C$.

By the choice of $d$, then \[thm:simpadedirect\] therefore implies that the first column of $G$ together with ${\bm{\delta}}'$ form a solution specification to the problem $({\bm{S}}, {\bm{g}}, {\bm{N}})$. Since the first column of $G$ is also the first column of $F'$, this finishes the proof.

\[lem:recursive\_R\_solves\] In the context of \[lem:recursive\_big\_matrix\], let $({\bm{\lambda}}_1, {\bm{\delta}}_1)$ and $({\bm{\lambda}}_2, {\bm{\delta}}_2)$ be solution specifications to the two sub-problems, and let ${\bm{r}} = (-N_0 \mid {\bm{\delta}}_1 \mid {\bm{\delta}}_2)$. If $([ {\bm{\lambda}} \mid * ], {\bm{\delta}}) = {{\ensuremath{\mathsf{NegMinBasis}}\xspace}}(d, R, {\bm{r}})$, where ${\bm{\lambda}}$ is a column vector and $$R = 
    \left[\begin{array}{c|c}
            1 & 1\\ \hline
            -{\bm{\lambda}}_1 &  \\\hline
              & -{\bm{\lambda}}_2 
          \end{array}\right] \ ,$$ then $({\bm{\lambda}}, {\bm{\delta}})$ is a solution specification for the combined problem.

We will prove the lemma by using \[lem:paderecprune\] to relate valid outputs of ${{\ensuremath{\mathsf{NegMinBasis}}\xspace}}(d, R, {\bm{r}})$ with valid outputs of ${{\ensuremath{\mathsf{NegMinBasis}}\xspace}}(d, A, {\bm{a}})$ from \[lem:recursive\_big\_matrix\].

For $i=1,2$, since $({\bm{\lambda}}_i, {\bm{\delta}}_i)$ is a solution specification to the $i$’th problem, then by \[thm:simpadedirect\] there is some $G_i \in {{\mathsf{K}}}[x]^{k_i \times 2n_i+1}$ whose first column is ${\bm{\lambda}}_i$ and such that $G_i$ is a valid output of ${{\ensuremath{\mathsf{NegMinBasis}}\xspace}}(d, H_i, {\bm{h}}_i)$, where $$H_i = \left [ \begin{array}{c}
                  -{\bm{S}}_i \\
                  I \\
                  {  \ifx\nulli  \Gamma_{{\bm{g}}}  \else  \Gamma_{{\bm{g}}_{i}}  \fi}
                \end{array} \right ] \in {{\mathsf{K}}}[x]^{(2n_i+1) \times n_i} ,$$ and ${\bm{h}}_i$ is as in \[lem:recursive\_big\_matrix\]. Note now that if $$F_1 := \left [\begin{array}{ccc}
                    1 & &   \\
                      & G_1 \\
                      & & G_2
                  \end{array} \right ] \in {{\mathsf{K}}}[x]^{(k_1+k_2+1) \times (2n_1 + 2n_2 + 3)} ,$$ then $(F_1, {\bm{r}})$ is a valid output of ${{\ensuremath{\mathsf{NegMinBasis}}\xspace}}(d, A_1, {\bm{a}})$: for ${{\textnormal{rowdeg}}}_{{\bm{a}}} F_1$ is clearly ${\bm{r}}$; $F_1$ has full row rank and is ${\bm{r}}$-row reduced; and the rows of $F_1$ must span all ${\bm{a}}$-order $d$ approximants of $A_1$, since the three column “parts” of $F_1$ correspond to the three row parts of $A_1$. .

Note now that $F_1 A_2 = R$. Thus by \[lem:paderecprune\], if $(F_2, {\bm{\delta}}) = {{\ensuremath{\mathsf{NegMinBasis}}\xspace}}(d, R, {\bm{r}})$, then $(F_2 F_1, {\bm{\delta}})$ is a valid output of ${{\ensuremath{\mathsf{NegMinBasis}}\xspace}}(d, A, {\bm{a}})$. Note that by the shape of $F_1$ then the first column ${\bm{\lambda}}$ of $F_2 F_1$ is the first column of $F_2$. Thus ${\bm{\lambda}}, {\bm{\delta}}$ are exactly as stated in the lemma, and by \[lem:recursive\_big\_matrix\] they must be a solution specification to the combined problem.

\[alg:recsimpade\] is correct. The cost of the algorithm is $O(n^{\omega-1}\, {{\mathsf{M}}}(d) (\log d) (\log d/n)^2)$, $d = \max_i \deg g_i$.

Correctness follows from \[lem:recursive\_R\_solves\]. For complexity, note that the choice of order in \[line:recursive:choosed\] is bounded by $2\max_i \deg g_i$, i.e. twice the value of $d$ of this theorem. So if $T(n)$ is the cost \[alg:recsimpade\] for given $n$ and where the order will be bounded by $O(d)$, then we have the following recursion: $$T(n) = \left \{\begin{array}{ll}
                     2T(n/2) + P(n) & \textrm{if } n > 1 \\
                     O({{\mathsf{M}}}(d)\log d) & \textrm{if } n = 1 \textrm{ (see {\Fref{sec:direct_reduced_basis}})}
                   \end{array}\right . \ ,$$ where $P(n)$ is the cost of line \[calltoNegMin\]. Using algorithm [$\mathsf{PopovBasis}$]{} for the computation of the negative part of the minimal approximant bases we can set $P(n)$ to the target cost. The recursion then implies $T(n) \in O(P(n))$.

[**Acknowledgements.**]{} The authors would like to thank George Labahn for valuable discussions, and for making us aware of the Hermite–Simultaneous [Padé]{}duality. We would also like to thank Vincent Neiger for making preprints of [@jeannerod_computation_2016] available to us. The first author would like to thank the Digiteo Foundation for funding the research visit at Waterloo, during which most of the ideas of this paper were developed.

[^1]: $\copyright$ Johan Rosenkilde, Arne Storjohann. This is the authors’ version of the work. It is posted here for your personal use. Not for redistribution. The definitive version was published in ISSAC ’16, http://dx.doi.org/10.1145/2930889.2930933.

[^2]: The notions $(-{\bm{N}})$-degree, $\deg_{-{\bm{N}}}$ and $(-{\bm{N}})$-row reduced are recalled in \[sec:preliminaries\].

[^3]: The restriction $N_i \leq \deg g_i$ in \[prob:sim\_pade\] ensures that for a given $\lambda$, the only possibilities for the $\phi_i$ in a solution are ${{\textnormal{rem}}}(\lambda S_i, \ g_i)$. In particular, if we allowed $N_i > \deg g_i$ then $(0,\ldots, 0, g_i, 0, \ldots,
      0)$ would be a solution which can not be directly reconstructed from its first element.
