[Axiom-math] [ANN] guessing formulas for sequences

[Martin Rubey]

Stupid me, I forgot to attach the packages...

%-*-Latex-*-
\documentclass{article}
\usepackage{axiom,amsmath,url}
\begin{document}
\title{rinterp.spad}
\author{Martin Rubey}
\maketitle
\begin{abstract}
Rational Interpolation
\end{abstract}
\section{Introduction}
This file contains a crude na\"ive implementation of rational interpolation,
where the coefficients of the rational function are in any given field.

\section{Questions and Outlook}
\begin{itemize}
\item Maybe this file should be joined with pinterp.spad, where polynomial
  Lagrange interpolation is implemented. This version parallels the structure of
  pinterp.spad closely.
\item There are probably better ways to implement rational interpolation. Maybe
  \url{http://www.cs.ucsb.edu/~omer/personal/abstracts/rational.html} contains
  something useful, but I don't know.
\item For those who speak german,
  \url{http://www.num.math.uni-goettingen.de/schaback/teaching/numath.ps}
  contains quite a bit of information.
\item The message: \lq\lq Warning: unattainable points!\rq\rq\ is a bit
  misleading. There can be unattainable points even if the message does not
  appear! I did not yet have time to fix this.
\item Comments welcome!
\end{itemize}

<<*>>=
<<RINTERP Header>>
<<RINTERP Exports>>
<<RINTERP Implementation>>
@

\section{RationalInterpolation}

<<RINTERP Header>>=
)abbrev package RINTERPA RationalInterpolationAlgorithms
++ Description:
++ This package exports rational interpolation algorithms
RationalInterpolationAlgorithms(F, P): Cat == Body   where
    F: Field 
    P: UnivariatePolynomialCategory(F)
@

<<RINTERP Exports>>= 
    Cat == with
        RationalInterpolation: (List F, List F, NonNegativeInteger, 
                                NonNegativeInteger) -> Fraction P
@

The implementation sets up a system of linear equations and solves it. 
<<RINTERP Implementation>>=
    Body == add
        RationalInterpolation(xlist, ylist, m, k) ==
@

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.
<<RINTERP Implementation>>=
            #xlist ^= #ylist =>
                error "Different number of points and values."
            #xlist ^= m+k+1 =>
                error "wrong number of points"
@

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{align*}
 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{align*}

This can be written as
$$
\begin{pmatrix}
  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{pmatrix}
\begin{pmatrix}
  a_0\\a_1\\\vdots\\a_m\\b_0\\b_1\\\vdots\\b_k
\end{pmatrix}=\mathbf 0
$$
We generate this matrix columnwise:
<<RINTERP Implementation>>= 
            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)])
@
Now we can solve the system:
<<RINTERP Implementation>>=
            res: List Vector F := nullSpace((transpose matrix collist) _
                                            ::Matrix F)
@

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.

More importantly, note that there may be unattainable points even if the
message does not appear, for example
<<inputfile>>=
interpolate([q,q^2,q^3],[0,x^1,x^2],0,2)$RINTERP(qn,FRAC POLY INT)
@

silently gives zero\dots

<<RINTERP Implementation>>=
            if #res~=1 then output("Warning: unattainable points!" _
                                   ::OutputForm)$OutputPackage
@

In this situation, all the solutions will be equivalent, thus we can always
simply take the first one:

<<RINTERP Implementation>>=
            reslist: List List P := _
                      [[monomial((res.1).(i+1),i) for i in 0..m], _
                      [monomial((res.1).(i+m+2),i) for i in 0..k]] 
@
Finally, we generate the rational function:
<<RINTERP Implementation>>=
            reduce((_+),reslist.1)/reduce((_+),reslist.2)

<<*>>=
)abbrev package RINTERP RationalInterpolation
++ Description:
++ This package exports interpolation algorithms
RationalInterpolation(xx, F): Cat == Body   where
    xx: Symbol
    F:  Field
    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)
@
\eject
\begin{thebibliography}{99}
\bibitem{1} nothing
\end{thebibliography}
\end{document}

%-*-Latex-*-
\documentclass{article}
\usepackage{axiom,amsmath,amssymb,url}
\DeclareMathOperator{\op}{op}
\newcommand{\axiom}{{\bf axiom}}
\begin{document}
\title{guess.spad}
\author{Martin Rubey}
\maketitle
\begin{abstract}
  Guessing formulas for sequences
\end{abstract}
\section{Introduction}
This file contains a slightly improved implementation of Christian
Krattenthalers Mathematica program [[rate.m]]. Given the first few terms of a
sequence of numbers, the program guesses a closed form that might fit the
sequence. 

More precisely, the program guesses formulas of the form
$${}^{(0)}\!\!\op_{i_0=0}^{n}{}^{(1)}\!\!\op_{i_1=0}^{i_0}\dots
  {}^{(k)}\!\!\op_{i_k=0}^{i_{k-1}}r(i_k)$$

where ${}^{(i)}\!\!\op$ is either $\prod$ or $\sum$ and $r(i)$ is a rational
function in $i$.

The basic idea is the following: assuming our sequence of numbers stems from a
rational function, we can simply use rational interpolation to determine the
formula. Otherwise, we transform the original sequence 
$$
x_1,x_2,\dots,x_{n-1},x_n
$$
in two different ways:

\begin{itemize}
\item we generate the sequence of successive quotients:
  $$
  x_2/x_1,x_3/x_2,\dots,x_n/x_{n-1}.
  $$
  If we are lucky, the original sequence was a product over a rational
  function. In this case, the sequence of successive quotients is a rational
  function and can be guessed using rational interpolation.
\item we generate the sequence of successive differences:
  $$
  x_2-x_1,x_3-x_2,\dots,x_n-x_{n-1}.
  $$
  Similarly, if the original sequence was a sum over a rational function, the
  sequence of successive differences is a rational function and can be guessed
  using rational interpolation.
\end{itemize}

To determine whether a given sequence stems from a rational function, we can
only use a heuristic: we interpolate using all but the last data point. If the
interpolation function interpolates also this point, we accept it, otherwise we
discard it.

\section{The implementation}

<<*>>=
<<GUESS Header>>
<<GUESS Exports>>
<<GUESS Implementation>>

<<GUESSINT>>
@

<<GUESS Header>>=
)abbrev package GUESS Guess
++ Description:
++ This package implements guessing of sequences
Guess(xx, F, basis, R, EXPRR, retract, coerce): Exports == Implementation where
    xx: Symbol                                                        
    F: Field                             -- zB.: FRAC POLY PF 5
    R: Join(OrderedSet, IntegralDomain)  -- zB.: FRAC POLY INT
    EXPRR: Join(ExpressionSpace, IntegralDomain, 
                RetractableTo R, RetractableTo Symbol, 
                RetractableTo Integer, CombinatorialOpsCategory) with
              _/ : (%,%) -> %
                                         -- zB.: EXPR INT
    EXPRI ==> Expression Integer
    basis: EXPRI -> EXPRR   -- zB.: i+->q^i
    retract: R -> F                      -- zB.: i+->i
    coerce: F -> EXPRR                   -- zB.: i+->i

    SUPF ==> SparseUnivariatePolynomial F
    FSUPF ==> Fraction SUPF

@

Unfortunately, we have to work around some difficulties arising from the design
of \axiom. This is why the header looks a little complicated.

In principle, all we need for our guessing game to work is an infinite field
[[F]] and a mapping [[basis]] from the positive integers into the field, which
will be used to generate the independent values for the rational interpolation.
If the field is not infinite, there are not too many functions around, so there
is no guessing game.

However, \axiom\ does not accept expressions over arbitrary fields, it requires
that the set of coefficients is an [[OrderedSet]]. For example, [[Expression
PrimeField 5]] is not a valid type. So, we specify a second domain, [[R]],
which is an [[OrderedSet]] which extends [[F]] to an [[OrderedSet]], a mapping
[[retract]] from [[R]] to [[F]], which retracts an element from [[R]] to [[F]]
and a mapping [[coerce]] which coerces an element of [[F]] to an element of
[[EXPRR]]. Fortunately, \axiom 's interpreter is very clever, so we can specify
these to functions simply as [[i+->i]]\dots

Furthermore, the argument domain of [[Expression]] must not contain
variables. For example, [[Expression Polynomial Integer]] is an invalid
type. Unfortunately, there is currently no way in \axiom\ to obtain the domain
of numbers of a given domain. Thus, [[Guess]] takes an argument [[EXPRR]], that
specifies the expression domain, over which we will have our result.

Finally, the modeline of summation and product are 
[[(%, SegmentBinding %%)->%]], i.e., the variable of summation and product
takes values in [[%]], rather than in [[Expression Integer]]. Thus we need
that [[EXPRR]] is [[RetractableTo Integer]].

<<GUESS Exports>>=
    Exports == with
      if not (F has Finite) then  
        guess: List F -> List EXPRR
@

The exported function takes a list of elements of [[F]] and returns a list of
expressions over [[R]] that seem to fit the data. For example,

<<inputfile>>=
guess([1,1+q,1+q+q^2])
$GUESS(n,FRAC POLY INT,i+->q^i,FRAC POLY INT, i+->i)
@

returns

[[      n 
       q  - 1
 (1)  [------]                                      
        q - 1                                       
             Type: List Expression Integer 
]]

Furthermore, if we type
<<inputfile>>=
guess([(i+q^i) for i in 0..3])
$GUESS(n,FRAC POLY INT,i+->q^i,FRAC POLY INT, EXPR INT, i+->i, i+->i)
@

we obtain

[[                       
                         s
          n - 1           4
           --+   (q - 1)q   + q
   (2)  [( >     --------------) + 1]
           --+          q
          s = 1
            4
                                      Type: List Expression Integer
]]

As a third example, we have
<<inputfile>>=
g:=guess([i+q^i for i in 1..4])
   $GUESS(m, FRAC POLY PF 5, i+->q^i, FRAC POLY INT, EXPR INT, i+->i, i+->i)
@

which gives

[[
           m - 1          s
            --+            4
   (3)  [( >     (q + 4)q   + 1) + q + 1]
            --+
           s = 1
            4
                                         Type: List Expression Integer
]]

Since \axiom\ knows no expressions over finite fields, it is up to the user to
interpret the result as such an expression.

\section{Helper functions}

<<GUESS Implementation>>=
    Implementation == add
      basis2: Integer -> F
      basis2 i == retract(retract(basis(i::EXPRI))@R)

@
This is used only in [[guessCons]] to generate the list of independent values
for the rational interpolation.

<<GUESS Implementation>>=
      guessCons: List F -> List FSUPF
      guessCons(list) ==
-- this function applies rationalInterpolation to the list of data points
-- [(list.i,basis.i) for i in 1..#list-1]. The last data point is used to test
-- for plausibility.

        len := (#list-1)::PositiveInteger 
        xlist := [basis2(i) for i in 1..len]
        x:F := basis2(len+1)
        ylist: List F := first(list,len)
        y:F := last(list)
        res: List FSUPF := []

        for i in 0..(len-1) repeat
          ri: FSUPF := interpolate(xlist, ylist, _
               (len-1-i)::NonNegativeInteger, i)$RationalInterpolation(xx, F)

@
We try to interpolate with the numerator and denominator polynomials having all
the possible degrees.

<<GUESS Implementation>>=
          de: F := elt(denom(ri), x)
          if (de ~= 0) and (elt(numer(ri), x) = y*de) _
             and not member?(ri, res) then
            res := cons(ri,res)  
@

If the denominator vanishes for the last data point [[x]], this point is
unattainable in the interpolation problem and we discard the result. Likewise,
we discard the result if the last point [[x]] does not have the right value
[[y]]. Finally we check whether the interpolating function [[ri]] is already in
the result list [[res]], if not, we add it.

<<GUESS Implementation>>=
        res 

@
We are done and return the list of interpolating functions.

<<GUESS Implementation>>=
      polyexpr: SUPF -> EXPRR
      polyexpr p == 
        zero? p => 0
        (coerce(leadingCoefficient p))::EXPRR * (xx::EXPRR)**degree p
           + polyexpr reductum p

@

This is just a tiny helper function that coerces a 
[[SparseUnivariatePolynomial F]] into an [[EXPRR]].

<<GUESS Implementation>>=
      expr: FSUPF -> EXPRR
      expr p == 
        (polyexpr numer p) / (polyexpr denom p)

@
This turns a [[Fraction SparseUnivariatePolynomial F]] into an [[EXPRR]].


\section{The main procedure}
<<GUESS Implementation>>=
      if not (F has Finite) then  
        guess(list) ==
@
Here comes the main procedure. We work recursively.

<<GUESS Implementation>>=
          res: List EXPRR := []
          len:= #list :: PositiveInteger
          if len > 1 then
@
If there is only one data point, we cannot do much, since we need at least a
second point to test our hypothesism, so we return.

<<GUESS Implementation>>=
            tempres: List FSUPF := guessCons(list)

            if not null(tempres) then 
              res := [eval(expr(tempres.i), kernel(xx), basis(xx::EXPRI)) _
                      for i in 1..#tempres]
@
If [[guessCons]] returns a result, we can add it to our list of results
immediately.

<<GUESS Implementation>>=
            if not member?(0$F, list) then
              prodList: List F := [(list.(i+1))/(list.i) for i in 1..(len-1)]
              if not every?(one?, prodList) then
                var: Symbol := subscript('p, [len::OutputForm])
                prodGuess := 
                    map(coerce(last(list) _
                               / retract(retract(eval(#1, kernel(xx), _
                                                      len::EXPRR))@R)) _
                        * #1, _
                        map(product(#1, equation(var, 1..xx::EXPRR-1)), _
                            guess(prodList)$Guess(var, F, basis, _
                                                  R, EXPRR, retract, coerce)))

                res := append(prodGuess, res)

@
If the sequence does not contain any zeros, we can generate the list of
successive quotients and recurse. The \lq\lq p\rq\rq\ in the subscript is for
\lq\lq product\rq\rq.

If [[guess]] returns something for the transformed sequence, note that this is
correct only up to a factor, since it fits the sequence of successive
quotients. In the very first version, we simply multiplied the result with the
first value of the sequence. However, this can lead to very wrong results in
special circumstances, that is, when the first point is an untainable point of
the interpolating problem. Therefore we use the last data point. Since this is
the point we tested against, it is certain to be correct. Note that even if
there are unattainable points, the result of [[guess]] can be useful.
Sometimes, for example, it is not clear how to define the first value of a
sequence.

<<GUESS Implementation>>=
            sumList: List F := [(list.(i+1))-(list.i) for i in 1..(len-1)]

            if not every?(#1=sumList.1, sumList) then
              var: Symbol := subscript('s, [len::OutputForm])
              sumGuess := 
                  map(coerce(last(list) _
                             - retract(retract(eval(#1, kernel(xx), _
                                                    len::EXPRR))@R)) _
                      + #1, _
                      map(summation(#1, equation(var, 1..xx::EXPRR-1)), _
                          guess(sumList)$Guess(var, F, basis, _
                                               R, EXPRR, retract, coerce)))

              res := append(sumGuess, res)

@

If the sequence of successive differences is constant, the original sequence is
already a rational function. Thus we procede only otherwise.

Again, if [[guess]] returns something for the transformed sequence, this is
correct with respect to the original sequence only up to a additive term. Thus
we have to do the same thing as in the former case.

<<GUESS Implementation>>=
          res
@
Done.

To ease the use of the package, I provide some packages that cover the most
common cases. Unfortunately, they do not seem to work right away. Lets see\dots

<<GUESSINT>>=
)abbrev package GUESSINT GuessInteger
++ Description:
++ This package exports guessing of integer sequences
GuessInteger(xx, basis): Exports == Implementation where
    xx: Symbol                                                        
    F ==> Fraction Integer
    R ==> Fraction Integer
    EXPRR ==> Expression Integer
    EXPRI ==> Expression Integer
    basis: EXPRI -> EXPRR

    Exports == with
      guess: List F -> List EXPRR

    Implementation == add
      retract: R -> F
      retract i == i::F

      coerce: F -> EXPRR
      coerce i == i::EXPRR

      guess(list) == 
        guess(list)$Guess(xx, F, basis, R, EXPRR, retract, coerce)
@
\end{document}

Excuse me,

Martin