Compact Rational Function Fitting

Question

I have some data from a study that I am trying to fit with a rational function of arbitrary form:

$y(x) = \frac{\sum_{i=n}^m A_ix^{i/q}}{\sum_{j=o}^p B_jx^{j/r}}$

where $i$, $j$, $k$, $l$, $m$, $n$, $o$, $p$, and $q$ are all integers. The goal is to produce a function of this form with as few terms as possible that still provides a reasonable match to the data. For instance, I have some data that looks like this:

where the black line is a fit to the data, and the fit is a function of the following form:

$y(x) = \frac{-0.31+1.2 \sqrt{x}-1.2 x }{1.0+1.3 \sqrt{x}-4.2 x}$

My current method of generating these functions is to produce a rational function of a high order, and then to successively eliminate terms if their coefficients are small compared to other terms in the expression. My algorithm is pretty ad-hoc, and the cutoffs, maximum order to use, and fractional powers to include are all parameters that one can pass as inputs to the function. I was curious if anyone had a better way to produce a compact rational fit to an arbitrary set of data that perhaps didn't rely on such an ad-hoc approach.

Here is the code I currently use to generate these fits:

compactRationalFit[data_, x_, coeffthresh_: 0.1, minorder_: 0, 
   maxorder_: 4, powerfrac_: 1] := 
  Module[{lcoeffthresh, bestfit, cld, cln, csc, i, dd, nn, norderlist,
     dorderlist, ndelpos, ddelpos, result, totdel, norigdelems, 
    norignelems},
   dorderlist = norderlist = Range[minorder, maxorder, 1/powerfrac];
   norignelems = Length[dorderlist];
   norigdelems = Length[dorderlist];
   totdel = 1;
   cld = ConstantArray[1, Length[dorderlist]];
   cln = ConstantArray[1, Length[norderlist]];
   While[totdel != 
      0 && ! (Length[dorderlist] == 1 && Length[norderlist] == 1),
lminorder = Min[norderlist~Join~dorderlist];
norderlist -= lminorder;
dorderlist -= lminorder;
lncoeffthresh = coeffthresh*Length[norderlist]/norignelems;
ldcoeffthresh = coeffthresh*Length[dorderlist]/norigdelems;
lncoeffthresh = coeffthresh;
ldcoeffthresh = coeffthresh;
cld[[Flatten[Position[cld, x_ /; x < 0]]]] = 1;

bestfit = 
 Normal[NonlinearModelFit[SetPrecision[data, 20],(*{*)
   Total[Table[
       ToExpression[
        "nn" <> 
         ToString[PaddedForm[i, 2, NumberPadding -> "0"]]], {i, 
        Length[norderlist]}] Table[x^i, {i, norderlist}]]/
    Total[Table[
       ToExpression[
        "dd" <> 
         ToString[PaddedForm[i, 2, NumberPadding -> "0"]]], {i, 
        Length[dorderlist]}] Table[x^i, {i, dorderlist}]](*,
   Apply[And,Table[ToExpression["dd"<>ToString[PaddedForm[i,2,
   NumberPadding->"0"]]]>0,{i,Length[dorderlist]}]]}*), 
   Table[{ToExpression[
       "dd" <> ToString[PaddedForm[i, 2, NumberPadding -> "0"]]], 
      cld[[i]]}, {i, Length[dorderlist]}]~Join~
    Table[{ToExpression[
       "nn" <> ToString[PaddedForm[i, 2, NumberPadding -> "0"]]], 
      cln[[i]]}, {i, Length[norderlist]}], x, 
   AccuracyGoal -> 3, PrecisionGoal -> 3, MaxIterations -> 1000, 
   ConfidenceLevel -> .95, 
   VarianceEstimatorFunction -> (Mean[#^2] &)]];
    cld = 
     Table[Coefficient[Expand[Denominator[bestfit]], x, i], {i, 
       dorderlist}]; 
    cln = Table[
      Coefficient[Expand[Numerator[bestfit]], x, i], {i, 
       norderlist}];
    cln = Delete[cln, Position[cln, x_ /; x == 0]]; 
    cld = Delete[cld, Position[cld, x_ /; x == 0]];
    ndelpos = 
     If[Length[cln] == 1, {}, 
      Position[cln, x_ /; Abs[x] < lncoeffthresh Max[Abs[cln]]]];
    ddelpos = 
     If[Length[cld] == 1, {}, 
      Position[cld, x_ /; Abs[x] < ldcoeffthresh Max[Abs[cld]]]];
    totdel = Length[Flatten[ndelpos]] + Length[Flatten[ddelpos]];
    norderlist = Delete[norderlist, ndelpos]; 
    dorderlist = Delete[dorderlist, ddelpos];
    cln = Delete[cln, ndelpos]; cld = Delete[cld, ddelpos];
    csc = First[cld];
    result = 
     SetPrecision[Total[cln/csc*Table[x^i, {i, norderlist}]]/
      Total[cld/csc*Table[x^i, {i, dorderlist}]], 2];
    ];
   Return[result];
   ];

Here, I have made the simplifying assumption that the numerator and the denominator have the same limits/fractional order, i.e. $n=o$, $p=m$, and $q=r$. A call to this function (using the default parameters) looks like this:

compactRationalFit[myData, x]

with the restriction that the number of data points + 1 should not exceed the total number of terms in the numerator and denominator combined.

Some of the other ideas I have been thinking of involve trying all possible combinations simultaneously, and selecting the one with the least amount of terms that satisfies an error threshold, but this involves finding fits to potentially thousands of possible combinations of power series. My current method is relatively fast, as it only requires at most $q^2(m-n)(m-n-1)/2$ iterations, but often times the fits it finds are not perfect, as rational functions are prone to zeroes in the denominator that can produce sharp discontinuities in-between data points:

If anyone could come up with a more robust solution that results in better, more compact fitting formulae, I would be grateful.

Answer 1

Try using a complex independent variable and fitting your function to the magnitude of the rational polynomial. This avoids the sharp discontinuities between data points due to poles in the denominator. With integer exponents the coefficients will still be real; I don't know what happens with fractional exponents.

myData = Function[x, {x, (-0.31 + 1.2*Sqrt[x] - 1.2*x)/(1. + 1.3*Sqrt[x] - 
    4.2*x)}] /@ Flatten[{Range[0.5, 1.5, 0.03], 2., 2.5}]; 
ratpoly[n_, m_] := Sum[a[k]*$s^(n - k), {k, 0, n}]/Sum[b[l]*$s^(m - l), 
  {l, 0, m}] /. b[m] -> 1; 
fitfnc = ratpoly[5, 5] /. $s -> I*w; 
vars = Cases[fitfnc, (a | b)[_], Infinity]; 
fit = fitfnc /. FindFit[myData, fitfnc, vars, w, NormFunction -> 
  (Abs[Norm[#1]] & ), Method -> "NMinimize"]; 
Show[{ListPlot[myData, PlotRange -> All, PlotStyle -> 
 Directive[PointSize[0.015], Opacity[0.5], ColorData[1, 3]]], 
Plot[Abs[fit], {w, 0.5, 2.5}, PlotRange -> All]}]