From: Bruce Johnson Subject: Re: rational approximation Date: 27 Oct 2000 13:18:48 GMT Newsgroups: sci.math.num-analysis Summary: [missing] Royer Fr?d?ric wrote: > Hi, > in the aim of modelizing intensity (Z) on a detector (X,Y), I would > like to approximate it by a 2D rational function: > Z = {\sum_i P(X,Y) \over \sum_j Q(X,Y)} > I know it can be found a rational Chebyshev approximation but I would > like not to interpolate in the plane (X,Y) and simply use the mesh at my > disposal which is uniformly distributed. > Does anyone know references of algorithm that could solve this problem > ? > -- > + __ . Frederic ROYER > * .'( `. Observatoire de Geneve > ------ * Chemin des Maillettes, 51 > * | | + CH-1290 SAUVERNY > _|__N_|_ Tel: +41 22 755 26 11 / Fax : +41 22 755 39 83 Well, this came up in the theoretical chemistry application of trying to find a continuous local representation for multivariable potential energy surfaces: R. E. Stevens, J. L. Kinsey, and B. R. Johnson, "Rational approximants as analytic polyatomic potential surfaces," J. Phys. Chem. 96, 7619 (1992). (I remember there also turned out to be a Computer Physics Communications article somewhere in the 91-94 neighborhood that pursued the same issues, I think in 1D, but don't have the reference. I'm sure the generality of the problem has made it appear in many other places, though.) In the above article can be found earlier references to work on algorithms for evaluating the coefficients for the "Canterbury approximants," ratios of two- or higher-variable polynomials obeying some restrictions. Those restrictions were loosened in our work, which gave some useful flexibility. The basic idea for solving all of these problems is to set up a system of equations in which the known samples of the function are set equal to multivariate approximant at those points. This will generally be an over- determined system if you have an excess of samples over coefficients, and least squares solution is natural. Unfortunately, a direct measure of error will lead to a nonlinear least squares problem. Instead, one can multiply each of the sample equations through by the denominator and try to solve the resulting linearized equations by least squares using some other error measure. This was turned into an iterative weighted least squares procedure which was reasonably automatable and would converge toward minimization of the original direct error measure. However, more or less the same problem occurs as does for 1D. For 1D functions which undergo a lot of variation in the region of interest, a zero of the denominator can appear between consecutive data points, or else complex zeros can appear near the real line (assuming the data to be real). These will be partially masked by nearly-coincident zeros in the numerator, so one can miss the areas of abnormal behavior unless sampling the calculated approximant sufficiently closely. That is, these rational functions can be very tricky to use as interpolants. If you do not want them for that purpose, then there may be no problem. We found that we could adjust weighting to eliminate specific nearly-coincident pole-zero combinations, but never developed a satisfactory scheme for that problem. Nor, to my knowledge, has anyone else. In 1D, I thought about the possibility of using restrictions to polynomials constrained to not have zeros in particular regions. This would diminish the approximation power of the approximants, but guarantee that no glitches happened between data points. However, one would like to accomplish this in a way that could be extended to the multivariate case, and that seems like a real challenge. Hope that is of some use (even if only to make you decide not to use the approximants!). Bruce ============================================================================== From: spellucci@mathematik.tu-darmstadt.de (Peter Spellucci) Subject: Re: rational approximation Date: 27 Oct 2000 15:39:41 GMT Newsgroups: sci.math.num-analysis In article <8tbvbo$cb$1@joe.rice.edu>, Bruce Johnson writes: |> Royer Fr?d?ric wrote: |> > Hi, |> |> > in the aim of modelizing intensity (Z) on a detector (X,Y), I would |> > like to approximate it by a 2D rational function: |> > Z = {\sum_i P(X,Y) \over \sum_j Q(X,Y)} snip. Requirement: use only given discrete data set (x_i,y_i,Z_i) for computing P and Q. hence a discrete approximation problem. solution: a) nonlinear least squares (proposed below by Bruce Johnson) b) mimimization of maximum error via nonlinear programming : minimize eps subject to the constraints -eps <= z_i -P(x_i,y_i)/Q(x_i,y_i) <= eps for all i unknowns of the problem: eps and the coeffs of P and Q. software: http://plato.la.asu.edu/guide.html |> The basic idea for solving all of these problems is to set up a system of |> equations in which the known samples of the function are set equal to |> multivariate approximant at those points. This will generally be an over- |> determined system if you have an excess of samples over coefficients, and |> least squares solution is natural. Unfortunately, a direct measure of |> error will lead to a nonlinear least squares problem. Why unfortunate? there are lots of good methods for doing that! |> Instead, one can |> multiply each of the sample equations through by the denominator and try |> to solve the resulting linearized equations by least squares using some |> other error measure. This was turned into an iterative weighted least |> squares procedure which was reasonably automatable and would converge |> toward minimization of the original direct error measure. this idea was published already in 1962 by Linde Wittmeyer in BIT. unfortunately, the method a)may not converge at all. this is the case if the deviations are large and b) if it converges it doesn't converge against the solution of the original error measure. (see two papers of mine from the early 70's) > |> However, more or less the same problem occurs as does for 1D. For |> 1D functions which undergo a lot of variation in the region of interest, |> a zero of the denominator can appear between consecutive data points, unfortunatly this is true (and already known to L. Wittmeyer) snip |> We |> found that we could adjust weighting to eliminate specific nearly-coincident |> pole-zero combinations, but never developed a satisfactory scheme for that |> problem. Nor, to my knowledge, has anyone else. false. T. Pomentale wrote a paper on it . the essential idea is to use some upper bound for the coefficients of the denominator first, obtaining upper bounds for the derivative (in your case gradient) of Q(x,y) on the region of interest. from this one you can compute a grid measure delta, such that the additional requirement Q(x_k,y_k) >=eps>0 for all points (x_k,y_k) in the grid (these are _not_ the measurement pints, these are (much more) points on a regular grid) guarantee that indeed Q(x,y)>=eps**2>0 (say) for all(x,y) in [x_anf,x_end]x[y_anf,y_end] Hence we get a nonlinear least squares problem with linear inequality constraints and there are good codes for that, see the guide. |> However, one would like to accomplish this in a way that |> could be extended to the multivariate case, and that seems like a real |> challenge. yes, this is a so called semiinfinite approximation problem. there are numerical schemes for dealing with that. what I proposed above is the so called discretization approach for the semiinfinite problem, but one can do even better and directly tackle it, since the feasible set of denominator coefficients is convex. see the relevant literature. hope this helps peter ============================================================================== From: Bruce Johnson Subject: Re: rational approximation Date: 27 Oct 2000 16:39:35 GMT Newsgroups: sci.math.num-analysis Peter Spellucci wrote: > Bruce Johnson writes: > snip. Requirement: use only given discrete data set (x_i,y_i,Z_i) for > computing P and Q. > hence a discrete approximation problem. > solution: > a) nonlinear least squares (proposed below by Bruce Johnson) > b) mimimization of maximum error via nonlinear programming : > minimize eps subject to the constraints > -eps <= z_i -P(x_i,y_i)/Q(x_i,y_i) <= eps for all i > unknowns of the problem: eps and the coeffs of P and Q. > software: > http://plato.la.asu.edu/guide.html > |> least squares solution is natural. Unfortunately, a direct measure of > |> error will lead to a nonlinear least squares problem. > Why unfortunate? there are lots of good methods for doing that! True enough. Times have changed. The first people trying to do such things for diatomic potential energy curves (1D) in fact did use nonlinear ls sq. Back in the 70's, these were very time-consuming processes, a fact now mitigated. If the approach which we reinvented had one merit, it was quite rapid. A problem still in existence, however, is the presence of local minima, particularly as the polynomial orders of numerator and denominator increase. This will affect both a) and b). I recall finding neighboring minima, in particular, a big pain even for global methods (though I don't know of anybody who has attempted those with rational approximants). > |> Instead, one can > |> multiply each of the sample equations through by the denominator and try > |> to solve the resulting linearized equations by least squares using some > |> other error measure. This was turned into an iterative weighted least > this idea was published already in 1962 by Linde Wittmeyer in BIT. > unfortunately, the method a)may not converge at all. this is the case if > the deviations are large and b) if it converges it doesn't > converge against the solution of the original error measure. > (see two papers of mine from the early 70's) Thanks! I always appreciate the earliest references. Can you provide your particular ones? Empirically, I found convergence with our iterative procedure in many cases, including to the original error measure, but not always and wanted to find out more about this. > false. T. Pomentale wrote a paper on it . > the essential idea is to use some upper bound for the coefficients of the > denominator first, obtaining upper bounds for the derivative (in your case > gradient) of Q(x,y) on the region of interest. from this one you > can compute a grid measure delta, such that the additional requirement > Q(x_k,y_k) >=eps>0 for all points (x_k,y_k) in the grid > (these are _not_ the measurement pints, these are (much more) points > on a regular grid) guarantee that indeed > Q(x,y)>=eps**2>0 (say) for all(x,y) in [x_anf,x_end]x[y_anf,y_end] I think I see. I will check this out sometime soon. I am guessing that what this means is that the original problem of a fixed detector grid with a particular number of X-Y data would then influence the values of eps that could be used? Bruce