login  home  contents  what's new  discussion  bug reports     help  links  subscribe  changes  refresh  edit
AxiomContributions
PolyMake ← ↑ → RealClosure

Rational Interpolation
last edited 1 year ago by page

Edit detail for Rational Interpolation revision 1 of 1
1
Editor: page
Time: 2007/11/12 23:59:22 GMT-8
Note: transferred from axiom-developer.org 
changed:
-
Introduction

  This file contains an implementation of rational interpolation, where the data
points are element of any integral domain.

Questions and Outlook

  - Maybe this file should be joined with pinterp.spad, where polynomial
    Lagrange interpolation is implemented. This version parallels the structure
    of pinterp.spad closely. This also answers comments and questions from
    wyscc. He remarked

    - Abbreviations for a constructor should be limited to 7 letters (not 8).
      The system occasionally adds the 8th character to a package for internal
      use.

    - Function names begin with a lower case, so RationalInterpolation should
      be rationalInterpolation, or better, rationalInterpolate.

  - Regarding the types I used for the values, wyscc remarked
   
    - If we are doing a rational interpolation, presumably the values are
      rational, so it does not make sense to require the $y$-coordinates of
      inputs be integral. On the other hand, as in the above example, if one
      uses 'FRAC INT', problems can arise when this package is combined with
      other packages that constructs the quotient field of the parameter domain
      'F' because Axiom does not like constructing 'FRAC FRAC INT'.
      
    Note however, that the package would rather construct the type 'FRAC SUP
    FRAC INT', so this problem should not occur. Moreover, there are situations
    - for example in the package [mantepse.spad2], where we want to interpolate values
    from an IntegralDomain. Of course we could first convert them to the
    quotient field, however, the current approach seems more natural to me.

  - Finally, wyscc asked:
    If <code>p(xx) = interpolate(lx, ly, m, k)</code>, what is the purpose of
    <code>elt(px, qx) = p(qx)</code>, the composition of <code>p(xx)</code> and
    <code>qx</code>, especially when <code>qx</code> is from <code>FRAC UP(xx,
    F)</code> instead of from just <code>F</code>? and why is this function
    (the composition) also called <code>interpolate</code>?

    I do not really know - apart from a very superficial level: Clearly, the
    second function was intended to let the user easily plug in values into the
    interpolated function. I don't find this sensible and I would be happy to
    change it. Indeed, this would also get rid of the first parameter to
    'RINTERP', which is quite a nuisance.

    I think we should agree on a general interface for interpolation
    algorithms, and mark 'PINTERP' as obsolete. By the way, it seems that
    'RINTERP' is faster, too.
    
  - There are probably better ways to implement rational interpolation. Maybe
    http://www.cs.ucsb.edu/~omer/personal/abstracts/rational.html
    contains something useful. In particular, in my package [mantepse.spad2], in 'guessRat'
    and 'guessExpRat' I generate interpolating polynomials for all possible degrees 
    of numerator and denominator. The above article contains an algorithm that does 
    this in time $O(n^2)$, which would be quite nice. Currently, I need $O(n^2)$ 
    operations for *each* degree!

  - For polynomial interpolation, there seems to be an algorithm that needs 
    only $O(n\log(n)^2\log\log(n))$ operations. It can be found in van zur Gathen's book
    "Modern computer algebra", chapter 10.

  - For those who speak german,
    http://www.num.math.uni-goettingen.de/schaback/teaching/numath.ps
    contains quite a bit of information.

  - This implementation of rational interpolation neither takes care of
    unattainable points, nor does it check whether the values of the
    $x$-coordinates are all distinct.

  - Comments welcome!


\begin{spad}
)abbrev package RINTERPA RationalInterpolationAlgorithms
++ Description:
++ This package exports rational interpolation algorithms
RationalInterpolationAlgorithms(F, P): Cat == Body   where
    F: IntegralDomain 
    P: UnivariatePolynomialCategory(F)
    Cat == with
        RationalInterpolation: (List F, List F, NonNegativeInteger,
                                NonNegativeInteger) 
                               -> Fraction P
        +++ We assume that the elements of the first list are all distinct.
        +++ If they are not, division by zero might occur.

    Body == add
        RationalInterpolation(xlist, ylist, m, k) ==
            #xlist ^= #ylist =>
                error "Different number of points and values."
            #xlist ^= m+k+1 =>
                error "wrong number of points"
            tempvec: List F := [1 for i in 1..(m+k+1)]

            collist: List List F := cons(tempvec, 
                                         [(tempvec := [tempvec.i * xlist.i _
                                                       for i in 1..(m+k+1)]) _
                                          for j in 1..max(m, k)])

            collist := append([collist.j for j in 1..(m+1)], _
                              [[- collist.j.i * ylist.i for i in 1..(m+k+1)] _
                               for j in 1..(k+1)])
            resspace: List Vector F := nullSpace((transpose matrix collist) _
                                                 ::Matrix F)
            reslist: List List P := _
                      [[monomial((resspace.1).(i+1), i) for i in 0..m], _
                      [monomial((resspace.1).(i+m+2), i) for i in 0..k]]

            reduce((_+), reslist.1)/reduce((_+), reslist.2)
\end{spad}

\begin{spad}
)abbrev package RINTERP RationalInterpolation
++ Description:
++ This package exports interpolation algorithms
RationalInterpolation(xx, F): Cat == Body   where
    xx: Symbol
    F:  IntegralDomain
    UP  ==> UnivariatePolynomial
    SUP ==> SparseUnivariatePolynomial
 
    Cat == with
        interpolate: (Fraction UP(xx, F), List F, List F, _
                      NonNegativeInteger, NonNegativeInteger) _
                      -> Fraction UP(xx, F)

        interpolate: (List F, List F, NonNegativeInteger, NonNegativeInteger) _
                      -> Fraction SUP F

    Body == add
        RIA ==> RationalInterpolationAlgorithms

        interpolate(qx, lx, ly, m, k) ==
            px := RationalInterpolation(lx, ly, m, k)$RIA(F, UP(xx, F))

            elt(px, qx)
 
        interpolate(lx, ly, m, k) ==
            RationalInterpolation(lx, ly, m, k)$RIA(F, SUP F)
\end{spad}


First we check whether we have the right number of points and values. Clearly
the number of points and the number of values must be identical. Note that we
want to determine the numerator and denominator polynomials only up to a
factor. Thus, we want to determine $m+k+1$ coefficients, where $m$ is the degree
of the polynomial in the numerator and $k$ is the degree of the polynomial in
the denominator.

In fact, we could also leave - for example - $k$ unspecified and determine it
as $k=\#xlist-m-1$: I don't know whether this would be better.

The next step is to set up the matrix. Suppose that our numerator polynomial is
$p(x)=a_0+a_1x+\dots+a_mx^m$ and that our denominator polynomial is
$q(x)=b_0+b_1x+\dots+b_mx^m$. Then we have the following equations, writing $n$
for $m+k+1$:

\begin{eqnarray*}
 p(x_1)-y_1q(x_1)&=a_0+a_1x_1+\dots +a_mx_1^m-y_1(b_0+b_1x_1+\dots +b_kx_1^k)=0\\
 p(x_2)-y_2q(x_2)&=a_0+a_1x_2+\dots +a_mx_2^m-y_2(b_0+b_1x_2+\dots +b_kx_2^k)=0\\
                 &\;\;\vdots\\                                                 
 p(x_n)-y_nq(x_n)&=a_0+a_1x_n+\dots +a_mx_n^m-y_n(b_0+b_1x_n+\dots +b_kx_n^k)=0
\end{eqnarray*}

This can be written as
\begin{equation*}
\begin{bmatrix}
  1&x_1&\dots&x_1^m&-y_1&-y_1x_1&\dots&-y_1x_1^k\\
  1&x_2&\dots&x_2^m&-y_2&-y_2x_2&\dots&-y_2x_2^k\\
\vdots\\
  1&x_n&\dots&x_n^m&-y_n&-y_nx_n&\dots&-y_nx_2^k
\end{bmatrix}
\begin{bmatrix}
  a_0\\a_1\\\vdots\\a_m\\b_0\\b_1\\\vdots\\b_k
\end{bmatrix}=\mathbf 0
\end{equation*}

We generate this matrix columnwise, then we can solve the system using 'nullSpace'.

Note that it may happen that the system has several solutions. In this case,
some of the data points may not be interpolated correctly. However, the
solution is often still useful, thus we do not signal an error.

Since all the solutions of 'nullSpace' will be equivalent, we can always
simply take the first one. Finally, we return the rational function.

Examples

  To conclude we present some examples. To begin with, the following interpolation 
illustrates the concept of unattainable points:

\begin{axiom}
interpolate([q,q^2,q^3],[0,x^1,x^2],0,2)$RINTERP(qn, FRAC POLY INT)
\end{axiom}

\begin{axiom}
f(x) == (x^3+5*x-3)/(x^2-3)
xlist := [1/2, 4, 1/6, 8, 1/10, 12]
ylist := [f(x) for x in xlist]

interpolate(xlist, ylist, 3, 2)$RINTERP('x, FRAC INT)
interpolate(1/6::FRAC UP(x,FRAC INT), xlist, ylist, 3, 2)$RINTERP('x,FRAC INT)
\end{axiom}

A harder example:

\begin{axiom}
dom := DMP([z],INT);
g: FRAC dom -> FRAC dom;
g(x) == (x^3*z+5*z^2*x -3*z^3)/(z*x^2 - 3)
xxlist: List FRAC dom := [1/(2*z), 4*z, 1/(6*z), 8*z, 1/(10*z), 12*z]
yylist := [g(x) for x in xxlist]
interpolate(xxlist, yylist, 3, 2)$RINTERP('x, FRAC dom)
interpolate(4*z::FRAC UP(x,dom), xxlist, yylist, 3, 2)$RINTERP('x, FRAC dom)
\end{axiom}

Introduction

This file contains an implementation of rational interpolation, where the data points are element of any integral domain.

Questions and Outlook

  • Maybe this file should be joined with pinterp.spad, where polynomial Lagrange interpolation is implemented. This version parallels the structure of pinterp.spad closely. This also answers comments and questions from wyscc. He remarked
    • Abbreviations for a constructor should be limited to 7 letters (not 8). The system occasionally adds the 8th character to a package for internal use.
    • Function names begin with a lower case, so RationalInterpolation should be rationalInterpolation, or better, rationalInterpolate.
  • Regarding the types I used for the values, wyscc remarked
    • If we are doing a rational interpolation, presumably the values are rational, so it does not make sense to require the LatexWiki Image-coordinates of inputs be integral. On the other hand, as in the above example, if one uses FRAC INT, problems can arise when this package is combined with other packages that constructs the quotient field of the parameter domain F because Axiom does not like constructing FRAC FRAC INT.

    Note however, that the package would rather construct the type FRAC SUP FRAC INT, so this problem should not occur. Moreover, there are situations - for example in the package mantepse.spad2, where we want to interpolate values from an IntegralDomain?. Of course we could first convert them to the quotient field, however, the current approach seems more natural to me.

  • Finally, wyscc asked: If p(xx) = interpolate(lx, ly, m, k), what is the purpose of elt(px, qx) = p(qx), the composition of p(xx) and qx, especially when qx is from FRAC UP(xx, F) instead of from just F? and why is this function (the composition) also called interpolate?

    I do not really know - apart from a very superficial level: Clearly, the second function was intended to let the user easily plug in values into the interpolated function. I don't find this sensible and I would be happy to change it. Indeed, this would also get rid of the first parameter to RINTERP, which is quite a nuisance.

    I think we should agree on a general interface for interpolation algorithms, and mark PINTERP as obsolete. By the way, it seems that RINTERP is faster, too.

  • There are probably better ways to implement rational interpolation. Maybe http://www.cs.ucsb.edu/~omer/personal/abstracts/rational.html contains something useful. In particular, in my package mantepse.spad2, in guessRat and guessExpRat I generate interpolating polynomials for all possible degrees of numerator and denominator. The above article contains an algorithm that does this in time LatexWiki Image, which would be quite nice. Currently, I need LatexWiki Image operations for each degree!
  • For polynomial interpolation, there seems to be an algorithm that needs only LatexWiki Image operations. It can be found in van zur Gathen's book "Modern computer algebra", chapter 10.
  • For those who speak german, http://www.num.math.uni-goettingen.de/schaback/teaching/numath.ps contains quite a bit of information.
  • This implementation of rational interpolation neither takes care of unattainable points, nor does it check whether the values of the LatexWiki Image-coordinates are all distinct.
  • Comments welcome!

spad
)abbrev package RINTERPA RationalInterpolationAlgorithms
++ Description:
++ This package exports rational interpolation algorithms
RationalInterpolationAlgorithms(F, P): Cat == Body   where
    F: IntegralDomain 
    P: UnivariatePolynomialCategory(F)
    Cat == with
        RationalInterpolation: (List F, List F, NonNegativeInteger,
                                NonNegativeInteger) 
                               -> Fraction P
        +++ We assume that the elements of the first list are all distinct.
        +++ If they are not, division by zero might occur.

    Body == add
        RationalInterpolation(xlist, ylist, m, k) ==
            #xlist ^= #ylist =>
                error "Different number of points and values."
            #xlist ^= m+k+1 =>
                error "wrong number of points"
            tempvec: List F := [1 for i in 1..(m+k+1)]

            collist: List List F := cons(tempvec, 
                                         [(tempvec := [tempvec.i * xlist.i _
                                                       for i in 1..(m+k+1)]) _
                                          for j in 1..max(m, k)])

            collist := append([collist.j for j in 1..(m+1)], _
                              [[- collist.j.i * ylist.i for i in 1..(m+k+1)] _
                               for j in 1..(k+1)])
            resspace: List Vector F := nullSpace((transpose matrix collist) _
                                                 ::Matrix F)
            reslist: List List P := _
                      [[monomial((resspace.1).(i+1), i) for i in 0..m], _
                      [monomial((resspace.1).(i+m+2), i) for i in 0..k]]

            reduce((_+), reslist.1)/reduce((_+), reslist.2)
spad
   Compiling FriCAS source code from file 
      /var/zope2/var/LatexWiki/959208282017174297-25px001.spad using 
      old system compiler.
   RINTERPA abbreviates package RationalInterpolationAlgorithms 
------------------------------------------------------------------------
   initializing NRLIB RINTERPA for RationalInterpolationAlgorithms 
   compiling into NRLIB RINTERPA 
   compiling exported RationalInterpolation : (List F,List F,NonNegativeInteger,NonNegativeInteger) -> Fraction P
Time: 0.83 SEC.

(time taken in buildFunctor:  0)

;;;     ***       |RationalInterpolationAlgorithms| REDEFINED

;;;     ***       |RationalInterpolationAlgorithms| REDEFINED
Time: 0 SEC.

   Cumulative Statistics for Constructor RationalInterpolationAlgorithms
      Time: 0.83 seconds

   finalizing NRLIB RINTERPA 
   Processing RationalInterpolationAlgorithms for Browser database:
--------(RationalInterpolation ((Fraction P) (List F) (List F) (NonNegativeInteger) (NonNegativeInteger)))---------
--->-->RationalInterpolationAlgorithms((RationalInterpolation ((Fraction P) (List F) (List F) (NonNegativeInteger) (NonNegativeInteger)))): Improper first word in comments: +
"+ We assume that the elements of the first list are all distinct. + If they are not,{} division by zero might occur."
--------constructor---------
------------------------------------------------------------------------
   RationalInterpolationAlgorithms is now explicitly exposed in frame 
      initial 
   RationalInterpolationAlgorithms will be automatically loaded when 
      needed from /var/zope2/var/LatexWiki/RINTERPA.NRLIB/code

spad
)abbrev package RINTERP RationalInterpolation
++ Description:
++ This package exports interpolation algorithms
RationalInterpolation(xx, F): Cat == Body   where
    xx: Symbol
    F:  IntegralDomain
    UP  ==> UnivariatePolynomial
    SUP ==> SparseUnivariatePolynomial

    Cat == with
        interpolate: (Fraction UP(xx, F), List F, List F, _
                      NonNegativeInteger, NonNegativeInteger) _
                      -> Fraction UP(xx, F)

        interpolate: (List F, List F, NonNegativeInteger, NonNegativeInteger) _
                      -> Fraction SUP F

    Body == add
        RIA ==> RationalInterpolationAlgorithms

        interpolate(qx, lx, ly, m, k) ==
            px := RationalInterpolation(lx, ly, m, k)$RIA(F, UP(xx, F))

            elt(px, qx)

        interpolate(lx, ly, m, k) ==
            RationalInterpolation(lx, ly, m, k)$RIA(F, SUP F)
spad
   Compiling FriCAS source code from file 
      /var/zope2/var/LatexWiki/6898989268314666018-25px002.spad using 
      old system compiler.
   RINTERP abbreviates package RationalInterpolation 
   processing macro definition UP ==> UnivariatePolynomial 
   processing macro definition SUP ==> SparseUnivariatePolynomial 
------------------------------------------------------------------------
   initializing NRLIB RINTERP for RationalInterpolation 
   compiling into NRLIB RINTERP 
   processing macro definition RIA ==> RationalInterpolationAlgorithms 
   compiling exported interpolate : (Fraction UnivariatePolynomial(xx,F),List F,List F,NonNegativeInteger,NonNegativeInteger) -> Fraction UnivariatePolynomial(xx,F)
Time: 0.04 SEC.

   compiling exported interpolate : (List F,List F,NonNegativeInteger,NonNegativeInteger) -> Fraction SparseUnivariatePolynomial F
Time: 0.01 SEC.

(time taken in buildFunctor:  0)

;;;     ***       |RationalInterpolation| REDEFINED

;;;     ***       |RationalInterpolation| REDEFINED
Time: 0 SEC.

   Warnings: 
      [1] interpolate: not known that (IntegralDomain) is of mode (CATEGORY package (SIGNATURE interpolate ((Fraction (UnivariatePolynomial xx F)) (Fraction (UnivariatePolynomial xx F)) (List F) (List F) (NonNegativeInteger) (NonNegativeInteger))) (SIGNATURE interpolate ((Fraction (SparseUnivariatePolynomial F)) (List F) (List F) (NonNegativeInteger) (NonNegativeInteger))))

   Cumulative Statistics for Constructor RationalInterpolation
      Time: 0.05 seconds

   finalizing NRLIB RINTERP 
   Processing RationalInterpolation for Browser database:
--->/usr/local/lib/axiom/target/x86_64-unknown-linux/../../src/algebra/RINTERP.spad-->RationalInterpolation((interpolate ((Fraction (UP xx F)) (Fraction (UP xx F)) (List F) (List F) (NonNegativeInteger) (NonNegativeInteger)))): Not documented!!!!
--->/usr/local/lib/axiom/target/x86_64-unknown-linux/../../src/algebra/RINTERP.spad-->RationalInterpolation((interpolate ((Fraction (SUP F)) (List F) (List F) (NonNegativeInteger) (NonNegativeInteger)))): Not documented!!!!
--------constructor---------
------------------------------------------------------------------------
   RationalInterpolation is now explicitly exposed in frame initial 
   RationalInterpolation will be automatically loaded when needed from 
      /var/zope2/var/LatexWiki/RINTERP.NRLIB/code

First we check whether we have the right number of points and values. Clearly the number of points and the number of values must be identical. Note that we want to determine the numerator and denominator polynomials only up to a factor. Thus, we want to determine LatexWiki Image coefficients, where LatexWiki Image is the degree of the polynomial in the numerator and LatexWiki Image is the degree of the polynomial in the denominator.

In fact, we could also leave - for example - LatexWiki Image unspecified and determine it as LatexWiki Image: I don't know whether this would be better.

The next step is to set up the matrix. Suppose that our numerator polynomial is LatexWiki Image and that our denominator polynomial is LatexWiki Image. Then we have the following equations, writing LatexWiki Image for LatexWiki Image:

LatexWiki Image 

This can be written as

LatexWiki Image 

We generate this matrix columnwise, then we can solve the system using nullSpace.

Note that it may happen that the system has several solutions. In this case, some of the data points may not be interpolated correctly. However, the solution is often still useful, thus we do not signal an error.

Since all the solutions of nullSpace will be equivalent, we can always simply take the first one. Finally, we return the rational function.

Examples

To conclude we present some examples. To begin with, the following interpolation illustrates the concept of unattainable points:

axiom
interpolate([q,q^2,q^3],[0,x^1,x^2],0,2)$RINTERP(qn, FRAC POLY INT)
LatexWiki Image(1)
Type: Fraction SparseUnivariatePolynomial? Fraction Polynomial Integer

axiom
f(x) == (x^3+5*x-3)/(x^2-3)
Type: Void
axiom
xlist := [1/2, 4, 1/6, 8, 1/10, 12]
LatexWiki Image(2)
Type: List Fraction Integer
axiom
ylist := [f(x) for x in xlist]
axiom
Compiling function f with type Fraction Integer -> Fraction Integer
LatexWiki Image(3)
Type: List Fraction Integer
axiom
interpolate(xlist, ylist, 3, 2)$RINTERP('x, FRAC INT)
LatexWiki Image(4)
Type: Fraction SparseUnivariatePolynomial? Fraction Integer
axiom
interpolate(1/6::FRAC UP(x,FRAC INT), xlist, ylist, 3, 2)$RINTERP('x,FRAC INT)
LatexWiki Image(5)
Type: Fraction UnivariatePolynomial(x,Fraction Integer)

A harder example:

axiom
dom := DMP([z],INT);
Type: Domain
axiom
g: FRAC dom -> FRAC dom;
Type: Void
axiom
g(x) == (x^3*z+5*z^2*x -3*z^3)/(z*x^2 - 3)
Type: Void
axiom
xxlist: List FRAC dom := [1/(2*z), 4*z, 1/(6*z), 8*z, 1/(10*z),
12*z]
LatexWiki Image(6)
Type: List Fraction DistributedMultivariatePolynomial?([z]?,Integer)
axiom
yylist := [g(x) for x in xxlist]
axiom
Compiling function g with type Fraction 
      DistributedMultivariatePolynomial([z],Integer) -> Fraction 
      DistributedMultivariatePolynomial([z],Integer)
LatexWiki Image(7)
Type: List Fraction DistributedMultivariatePolynomial?([z]?,Integer)
axiom
interpolate(xxlist, yylist, 3, 2)$RINTERP('x, FRAC dom)
LatexWiki Image(8)
Type: Fraction SparseUnivariatePolynomial? Fraction DistributedMultivariatePolynomial?([z]?,Integer)
axiom
interpolate(4*z::FRAC UP(x,dom), xxlist, yylist, 3, 2)$RINTERP('x, FRAC
dom)
LatexWiki Image(9)
Type: Fraction UnivariatePolynomial(x,Fraction DistributedMultivariatePolynomial?([z]?,Integer))