# Padé-Hermite Approximants

## Introduction

A given rational function $F(z)\in k(z)$ over a field $k$ can be written as a power series $f(z)$ in the completion $k((z))$ of $k(z)$ at the place $(z)$. It thus is a good approximation of $f(z)$ in the sense that $F(z)$ is equal to $f(z)$ up to infinite order.

Padé-Hermite approximants deal with the converse. Given a (formal) power series $f(z)\in k[[z]]$ and some non- negative integers $n_P$ and $n_Q$, find polynomials $P(z), Q(z)\in k[z]$ of degrees at most $n_P$ and $n_Q$, respectively, such that $P(z)*f(z)-Q(z) = O(z^{n_P+n_Q+1})$ holds. In other words, $P/Q(z)$ is an approximation for $f(z)$ for polynomials $P$ and $Q$ of limited degree in $z$. In this notation the pair $[P,Q]$ is known as the Padé-Hermite approximant for the power series tuple $(f,-1)$.

The idea of approximating one power series by two polynomials can be extended to approximating several power series at the same time as follows. Consider the vector $\underline{f}^T:= (f_1, f_2, \ldots, f_m)^T$ in the $m$-dimensional vector space $k((z))^m$ over a power series ring. Let $\underline{n}=(n_1, n_2,\ldots, n_m)$ be an m-tuple of non-negative integers. A *Padé-Hermite approximant* of $\underline{f}^T(z)$ of type $\underline{n}$ is a non-zero vector of polynomials $\underline{P}=(P_1,P_2,\ldots,P_m)$ in $k[z]^m$ such that $\underline{P}\cdot\underline{f}^T = O(z^{N})$ , $N=n_1+n_2+\cdots+n_m+m-1$, is satisfied. A non-trivial Padé-Hermite approximant always exists. The approximant is contained in the sub-space $V_{\underline{f},N}=\{ \underline{Q}\in k[z]^m : \underline{Q}\cdot\underline{f}^T = O(z^{N})\}$.

The implementation of Padé-Hermite approximants is based on [[Derksen, 1994](../../references.md#cite-derksen94)] and [[Beckermann and Labahn, 1994](../../references.md#cite-bl94)]. They are implemented in Magma as sequences rather than vectors.

## Ordering of Sequences

A Padé-Hermite approximant to some sequence of power series does not have to be unique. They can be ordered according to their maximum degree and their type of sequence. A sequence $\underline{P}=[P_1, P_2,..., P_m]$ is of degree $i$ if the maximum of the weak degrees of $P_1, P_2, \ldots, P_m$ is $i$. An extension of the definition maximum degree is when a distortion of non-negative integers $\underline{d}=[d_1, d_2, \ldots, d_m]$ on the degrees is allowed. In this case the maximum degree is defined as the maximum of $\deg{P_i}-d_i$. The type of $\underline{P}$ identifies the last $P_i$ in $\underline{P}$ whose weak degree equals the maximum degree of $\underline{P}$. Again a distortion on the degrees is allowed.

### `MaximumDegree(f): SeqEnum -> RngIntElt`

```magma
Distortion: SeqEnum                    Default: []
```

MaximumDegree returns the degree of a sequence of polynomials or power series, defined as the maximum of the degrees of $f[i]-d[i]$, where $d$ is the distortion. The value -infinity is returned in the case that $f$ is weakly equal to the zero-sequence.

### `Example: Degree Of Sequence (ex-ea74e8)`

```magma
> S<u> := PowerSeriesRing(Rationals());
> f := [u+u^2, 2+u^2+u^3,0];
> MaximumDegree(f);
3
> MaximumDegree(f:Distortion:=[]);
3
> MaximumDegree(f:Distortion:=[0,2,1]);
2
> MaximumDegree([S|0,0]);
-Infinity
> MaximumDegree([O(u)]);
-Infinity

```

### `TypeOfSequence(f): SeqEnum -> RngIntElt, RngIntElt`

```magma
Distortion: SeqEnum                    Default: []
```

Returns the highest index $i$ for those $f[i]$ whose (distorted) degree is weakly equal to the maximum of the degrees of all entries. The second integer returned is the maximum degree of the sequence.

### `Example: Type Of Sequence (ex-8d75ca)`

```magma
> S<u> := PowerSeriesRing(Rationals());
> f := [u+u^2, 2+u^2+u^3,0];
> TypeOfSequence(f);
2 3
> TypeOfSequence(f:Distortion:=[]);
2 3
> TypeOfSequence(f:Distortion:=[0,2,1]);
1 2
> TypeOfSequence([S|0,0]);
2 -Infinity

```

The Padé-Hermite approximant of $\underline{f}$ can be seen as an element of the free module $k[z]^m$ of rank $m$, or as an element of the sub-space $V_{\underline{f},N}$.

A free sub-model $V\subset k[z]^m$ is generated by $m$ polynomial vectors $\underline{Q}_i$, $1=1,2,\ldots,m$, such that $\underline{Q}_1(z),  \underline{Q}_2(z), \ldots, \underline{Q}_m(z)$ form a minimal vector sequence of $V$. Such a sequence is defined as a sequence $S$ of $m$ vectors in $V$, such that $S[i]$ is a non-trivial polynomial vector in $V$ of minimal degree of type $i$, for $i=1,2,\ldots,m$. A minimal vector sequence is not unique.

Two variations on a minimal vector sequence are implemented. The first allows a distortion as an attribute. The sequence is the based on the distorted maximal degree. The second attribute sets the positive power $p$ of $z$ in each of the $\underline{Q}_i$ as follows. Instead of considering $\underline{Q}_i(z)$ for each $i$, one considers $\underline{Q}_i(z^p)$.

### `MinimalVectorSequence(f, n): SeqEnum, RngIntElt -> SeqEnum`

```magma
Distortion: SeqEnum                      Default: []
Power     : RngIntElt                    Default: 1
```

A minimal sequence of vectors $\underline{Q}_1, \underline{Q}_2, \ldots, \underline{Q}_m$ with respect to the sequence $f$ of length $m$ whose entries are polynomials or power series. The order of $\underline{Q}_i\cdot f$ is at least $m$.

### `Example: Minimal Vector Sequence (ex-67b928)`

```magma
> S<u> := PowerSeriesRing(Rationals());
> f := [u+u^2, 2+u^2+u^3];
> seq := MinimalVectorSequence(f, 2);
> seq;
[
    (u 0),
    (   -1 1/2*u)
]
> sums := [&+([Q[i]*f[i]: i in [1..#f]]) : Q in seq];
> sums;
[
    u^2 + u^3,
    -u^2 + 1/2*u^3 + 1/2*u^4
]
> seq := MinimalVectorSequence(f,3);
> #seq eq 2,  seq[1], seq[2];
true (u^2   0)
(-1 + u  1/2*u)

```

### `Example: The Next Example (ex-d1c2b6)`

```magma
> L<x> := PolynomialRing(Rationals());
> f := [1+x, 3-x^2, 5+x+x^3-x^5];
> seq := MinimalVectorSequence(f, 10);
> seq[1];
(-2*x^2 + 6   -2*x - 2          0)
> sums := [&+([Q[i]*f[i]: i in [1..#f]]) : Q in seq];
> sums;
[
    0,
    0,
    x^10
]

```

### `Example: Another Example (ex-328015)`

```magma
> S<u> := PowerSeriesRing(Rationals());
> f := [2*u^4,2+u^3+u^6];
> seq := MinimalVectorSequence(f, 10);
> seq;
[
    (1/2*u^6       0),
    (-1/2 - 1/4*u^3        1/2*u^4)
]
> sums := [&+([Q[i]*f[i]: i in [1..#f]]) : Q in seq];
> sums;
[
    u^10,
    1/2*u^10
]

```

### `Example: One More (ex-1d53c0)`

```magma
> S<u> := PowerSeriesRing(Rationals());
> f :=  [1+u-7*u^2, 6-3*u+1/2*u^2-u^3, 5-u+u^2];
> seq := MinimalVectorSequence(f, 5);
>  [&+([Q[i]*f[i]: i in [1..#f]]) : Q in seq];
[
    0,
    u^5,
    0
]
> seq := MinimalVectorSequence(f, 5:Distortion:=[2,0,1]);
> [&+([Q[i]*f[i]: i in [1..#f]]) : Q in seq];
[
    -7/2*u^5,
    u^5,
    0
]
> p:=2;
> seq := MinimalVectorSequence(f, 5:Distortion:=[2,0,1], Power:=p);
> sums := [&+([Q[i]*f[i]: i in [1..#f]]) : Q in seq];
> sums;
> mp:= map<S->S|  x :-> (IsWeaklyZero(x) select 0
>     else  &+([Coefficient(x,i)*(S.1)^(p*i) :  i in Exponents(x)]))
>       + (ISA(Type(v),RngIntElt) select O((S.1)^(p*v))
>   else S!0 where v := AbsolutePrecision(x))>;
> sums := [&+([mp(Q[i])*f[i]: i in [1..#f]]) : Q in seq];
> sums;
[
    u^6 + u^7 - 7*u^8,
    u^5 - 23/3*u^6,
    -5/3*u^5 + 35/3*u^6
]

```

## Approximants

The Padé-Hermite approximant of type $\underline{d}=[d_1,d_2,\ldots,d_m]$ with respect to the tuple $\underline{f}^T\in k((z))^m$ is an element of the space $V_{\underline{f},N}=\{ \underline{Q}\in k[z]^m : \underline{Q}\cdot\underline{f}^T = O(z^{N})\}$ for $N$ equal to $d_1+d_2+\cdots+d_m+m-1$. This space is generated by the vectors in the minimal vector sequence with respect to $\underline{f}$ with distortion $\underline{d}$. The routine [`PadeHermiteApproximant`](#function-padhermapp) returns one that is smallest with respect to the degree on sequences. The input sequence $\underline{f}$ must be a sequence of polynomial ring elements, or be a power series sequence. While the Padé Hermite approximants theoretically are polynomials, Magma returns them as elements of the same ring the entries of $\underline{f}$ are contained in.

### `PadeHermiteApproximant(f, d): SeqEnum, SeqEnum -> ModTupRngElt, SeqEnum, RngIntElt`

```magma
Power: RngIntElt                    Default: 1
```

Returns a Padé-Hermite form $\underline{P}$ of $f$ with distortion $d$, smallest with respect to the degree on sequences, and the corresponding minimal vector sequence. The third argument returned is the order in the order term of $\underline{P}\cdot f$.

### `Example: Pade Hermite Approximants (ex-05a6f0)`

This example can be found on page 813 in [[Beckermann and Labahn, 1994](../../references.md#cite-bl94)].

```magma
> S<u> := PowerSeriesRing(Rationals());
> f := [1,u,u/(1-u^4)+u^10+O(u^16),u/(1+u^4)+u^12+O(u^16)];
> pade, padebasis, ord := PadeHermiteApproximant(f,[2,2,2,2]);
> pade, ord;
( u -1  0  0)
11
> BaseRing(Parent(pade)) eq S;
true
> MinimalVectorSequence(f,10);
[
    ( u -1  0  0),
    (   0  u^4 -1/2  1/2),
    (         0    1 - u^4 -1/2 + u^4       -1/2),
    (    0    -u 1/2*u 1/2*u)
]
>
> p := 2;
> seq := MinimalVectorSequence(f,10: Distortion :=[2,2,2,2],Power := p);
> seq;
[
    (u^5   0   0   0),
    (   0  u^2 -1/2  1/2),
    (         0    1 - u^2 -1/2 + u^2       -1/2),
    (    0    -u 1/2*u 1/2*u)
]
> mp:= map<S->S|  x :-> (IsWeaklyZero(x) select 0
>     else  &+([Coefficient(x,i)*(S.1)^(p*i) :  i in Exponents(x)]))
>       + (ISA(Type(v),RngIntElt) select O((S.1)^(p*v))
>   else S!0 where v := AbsolutePrecision(x))>;
> sums := [&+([mp(Q[i])*f[i]: i in [1..#f]]) : Q in seq];
> [Valuation(v) : v in sums];
[ 10, 10, 10, 11 ]
> sums;
[
    u^10,
    -1/2*u^10 + 1/2*u^12 - u^13 + O(u^16),
    -1/2*u^10 - 1/2*u^12 + u^13 + u^14 + O(u^16),
    u^11 + 1/2*u^12 + 1/2*u^14 + O(u^18)
]

```

### `Example: Last Example (ex-e98d49)`

This example covers the example on page 815 in [[Beckermann and Labahn, 1994](../../references.md#cite-bl94)].

```magma
> S<u> := PowerSeriesRing(Rationals());
> f := [1,u, -1-u^4-2*u^8+u^10+u^11-u^12+O(u^16),-u-u^5-u^9-u^14-u^15+O(u^16)];
> dist:=[2,2,3,3];
> seq := MinimalVectorSequence(f,13:Distortion:=dist);
> pade, padebasis, ord :=  PadeHermiteApproximant(f,dist);
> pade, ord;
(-u  1  0  0)
13
> padebasis;
[
    (-u  1  0  0),
    (1/2*u -1/2*u^4 1/2*u + 1/2*u^3 + 1/2*u^4 -1/2*u^2 - 1/2*u^3 - u^4),
    (-1/2 1/2*u^3 -1/2 - 1/2*u^2 - 1/2*u^3 - u^4 1/2*u + 1/2*u^2 + 2*u^3),
    (      -u        0        0 -1 + u^4)
]
> padebasis eq seq;
true
> [[Valuation(w[i]): i in [1..Degree(w)]] : w in seq];
[
    [ 1, 0, Infinity, Infinity ],
    [ 1, 4, 1, 2 ],
    [ 0, 3, 0, 1 ],
    [ 1, Infinity, Infinity, 0 ]
]
> [[MaximumDegree([w[i]])-dist[i]: i in [1..Degree(w)] ] : w in seq];
[
    [ -1, -2, -Infinity, -Infinity ],
    [ -1, 2, 1, 1 ],
    [ -2, 1, 1, 0 ],
    [ -1, -Infinity, -Infinity, 1 ]
]
> p:=2;
> seq := MinimalVectorSequence(f,12:Distortion:=dist,Power:=p);
> seq;
[
    (  u - u^3         0 u - 2*u^3         0),
    (-1 - u + u^2 -u^3 -1 - u + 2*u^2 + u^3 -u^3),
    (u + u^2 - u^3 0 u + u^2 - 2*u^3 - u^4 0),
    (      0       1       0 1 - u^2)
]
%> [[MaximumDegree([w[i]])-dist[i]: i in [1..Degree(w)] ] : w in seq];
%[
%    [ 1, -Infinity, 0, -Infinity ],
%    [ 0, 1, 0, 0 ],
%    [ 1, -Infinity, 1, -Infinity ],
%    [ -Infinity, -2, -Infinity, -1 ]
%]
> seq[1]-seq[3];
(      -u^2          0 -u^2 + u^4          0)
> mp:= map<S->S|  x :-> (IsWeaklyZero(x) select 0
>     else  &+([Coefficient(x,i)*(S.1)^(p*i) :  i in Exponents(x)]))
>       + (ISA(Type(v),RngIntElt) select O((S.1)^(p*v))
>   else S!0 where v := AbsolutePrecision(x))>;
>  [Valuation(&+([mp(Q[i])*f[i]: i in [1..#f]])) : Q in seq];
[ 12, 12, 13, 13 ]

```

### `Example (ex-ea09f6)`

This example considers the example on page 816 in [[Beckermann and Labahn, 1994](../../references.md#cite-bl94)].

```magma
> S<u> := PowerSeriesRing(Rationals());
> f := [1,u,-1-u^4-2*u^8+u^10+O(u^12),-u-u^5-u^9+u^10+O(u^12)];
> dist := [2,2,3,3];
> seq := MinimalVectorSequence(f,12: Distortion := dist);
> [[MaximumDegree([w[i]])-dist[i]: i in [1..Degree(w)] ] : w in seq];
[
    [ -1, -2, -Infinity, -Infinity ],
    [ -2, 1, 0, 0 ],
    [ -Infinity, -Infinity, 1, 0 ],
    [ -1, -Infinity, 0, 1 ]
]
> [Valuation(&+([(Q[i])*f[i]: i in [1..#f]]) ) : Q in seq];
[ Infinity, 12, 12, 12 ]
> [MaximumDegree([ &+([(Q[i])*f[i]: i in [1..#f]]) ]) : Q in seq];
[ -Infinity, -Infinity, 14, -Infinity ]
> PadeHermiteApproximant(f,[2,2,3,3]);
% error
%
% redefine with more precision
%> g := [1,u,-1-u^4-2*u^8+u^10+O(u^13),-u-u^5-u^9+u^10+O(u^13)];
%> seq := MinimalVectorSequence(g,12:Distortion:=[2,2,3,3]);
%> PadeHermiteApproximant(g,[2,2,3,3]);
%(-u  1  0  0)
%[
%    (-u  1  0  0),
%    (1/4*u -1/4*u^4 1/4*u + 1/4*u^3 - 1/4*u^4 -1/4*u^2 + 1/4*u^3 - 1/2*u^4),
%    (-1/4 1/4*u^3 -1/4 - 1/4*u^2 + 1/4*u^3 - u^4 1/4*u - 1/4*u^2 + 3/2*u^3),
%    (1/2 - u -1/2*u^3 1/2 - 1/2*u^2 + 1/2*u^3 -1 + 1/2*u - 1/2*u^2 - u^3 + u^4)
%]
%13
> p := 2;
> seq := MinimalVectorSequence(f,12:Distortion:=[2,2,3,3],Power:=p);
> seq;
[
    (  1 - u^2        -1 1 - 2*u^2  -1 + u^2),
    (   0 -u^4    0 -u^4),
    (      -u       -1 -u + u^3 -1 + u^2),
    (      0       u       0 u - u^3)
]
> [[MaximumDegree([w[i]])-dist[i]: i in [1..Degree(w)] ] : w in seq];
[
    [ 0, -2, -1, -1 ],
    [ -Infinity, 2, -Infinity, 1 ],
    [ -1, -2, 0, -1 ],
    [ -Infinity, -1, -Infinity, 0 ]
]
> mp:= map<S->S|  x :-> (IsWeaklyZero(x) select 0
>     else  &+([Coefficient(x,i)*(S.1)^(p*i) :  i in Exponents(x)]))
>       + (ISA(Type(v),RngIntElt) select O((S.1)^(p*v))
>   else S!0 where v := AbsolutePrecision(x))>;
>
> [Valuation(&+([mp(Q[i])*f[i]: i in [1..#f]])) : Q in seq];
[ 12, 13, 12, 12 ]

```

A variant of the Padé-Hermite approximant is when the exponent in the order term is set rather that the type of the sequence. It is also possible to let $\underline{f}$ be a sequence such that its entries themselves are vectors of polynomials or power series.

### `PadeHermiteApproximant(f, m): SeqEnum, RngIntElt -> ModTupRngElt, SeqEnum`

```magma
Power: RngIntElt                    Default: 1
```

Returns a Padé-Hermite form of minimal degree in the corresponding minimal vector sequence, such that its inproduct with $f$ has order at least $m$. The second argument returned is the corresponding minimal vector sequence.

### `Example: Pade Hermite Approximants Vectors (ex-b5d830)`

This example can be found on page 813 in [[Beckermann and Labahn, 1994](../../references.md#cite-bl94)].

```magma
> S<u> := PowerSeriesRing(Rationals());
> f := [Vector([1]), Vector([u])];
> pade, seq := PadeHermiteApproximant(f,3);
Calculating the Pade'-Hermite approximant for the sequence [
    1,
    u
]
with order term 3 and power 1 .
> pade;
( u -1)
> seq;
[
    ( u -1),
    (  0 u^2)
]
> mat := Matrix([Eltseq(v): v in f]);
> pade*mat;
(0)
> PadeHermiteApproximant([1,u],5);
( u -1)
[
    ( u -1),
    (  0 u^4)
]
> PadeHermiteApproximant(f,3:Power:=2);
> g:= [Vector([1,0,0]), Vector([0,1,0]), Vector([1+u,2+u^2,u^3])];
> pade := PadeHermiteApproximant(g,5);
Calculating the Pade'-Hermite approximant for the sequence [
    1,
    u,
    1 + 2*u + u^3 + u^7 + u^11
]
with order term 15 and power 3 .
> pade;
(-u^3 - u^4     -2*u^3        u^3)
> pade*Matrix([Eltseq(v): v in g]);
(  0 u^5 u^6)

```

### `Example (ex-15b388)`

This example considers Padé-Hermite approximants for some series that have non-trivial power series expansions.

```magma
> S<u> := PowerSeriesRing(Rationals());
> f := [Sin(u), Cos(u), Exp(u)];
> [Valuation(f[i]) : i in [1..#f]], [Degree(f[i]) : i in [1..#f]];
[ 1, 0, 0 ]
[ 19, 20, 20 ]
> [AbsolutePrecision(f[i]) : i in [1..#f]];
[ 21, 22, 21 ]
> dist := [3,2,5];
> pade, seq, ord := PadeHermiteApproximant(f,dist);
> 1/420*pade;
(-1275 - 255*u + 45*u^2 + 5*u^3 120 + 495*u + 75*u^2 -120 + 900*u - 600*u^2 +
    160*u^3 - 20*u^4 + u^5)
> ord eq &+(dist)+#f-1, ord;
true 12
> [Degree(pade[i]) : i in [1..Degree(pade)]];
[ 3, 2, 5 ]
% The next one is in the Maple package nummapprox[hermite-pade].
> g:= [Cos(2*u)*(u+1)+3,Cos(u)^2+u*Cos(u)+1,Cos(2*u)+1,Cos(u)];
> pade, basis := PadeHermiteApproximant(g,20);
> 131/75880*pade;
(          2    -4 + 2*u        -3*u 4*u - 2*u^2)
% This is one from Maple's OrderBasis.
% It returns a different minimal vector sequence.
> h := [ 1+u^2-u^7+u^12, Sin(u), Exp(u) ];
> dist:=[3,1,2];
> seq := MinimalVectorSequence(h,8:Distortion := dist);
> sums := [&+([Q[i]*h[i]: i in [1..#f]]) : Q in seq];
> [Valuation(s) : s in sums];
[ 8, 8, 8 ]
> [[MaximumDegree([w[i]]): i in [1..Degree(w)] ] : w in seq];
[
    [ 4, 1, 2 ],
    [ 4, 2, 2 ],
    [ 3, 1, 2 ]
]
>  [[MaximumDegree([w[i]])-dist[i]: i in [1..Degree(w)] ] : w in seq];
[
    [ 1, 0, 0 ],
    [ 1, 1, 0 ],
    [ 0, 0, 0 ]
]

```