PageRenderTime 63ms CodeModel.GetById 28ms app.highlight 17ms RepoModel.GetById 8ms app.codeStats 0ms

HTML | 535 lines | 535 code | 0 blank | 0 comment | 0 complexity | 5772150156029cc047e41ebd92723a87 MD5 | raw file
  3<meta http-equiv="Content-Type" content="text/html; charset=US-ASCII">
  4<title>Sample Article (The Remez Method)</title>
  5<link rel="stylesheet" href="../../../../doc/src/boostbook.css" type="text/css">
  6<meta name="generator" content="DocBook XSL Stylesheets V1.75.2">
  7<link rel="home" href="../index.html" title="Document To Test Formatting">
  8<link rel="up" href="../index.html" title="Document To Test Formatting">
  9<link rel="prev" href="test.html" title="test HTML4 symbols">
 10<link rel="next" href="array.html" title="Array Example Boostbook XML Documentation">
 12<body bgcolor="white" text="black" link="#0000FF" vlink="#840084" alink="#0000FF">
 13<table cellpadding="2" width="100%"><tr>
 14<td valign="top"><img alt="Boost C++ Libraries" width="277" height="86" src="../../../../boost.png"></td>
 15<td align="center"><a href="../../../../index.html">Home</a></td>
 16<td align="center"><a href="../../../../libs/libraries.htm">Libraries</a></td>
 17<td align="center"><a href="">People</a></td>
 18<td align="center"><a href="">FAQ</a></td>
 19<td align="center"><a href="../../../../more/index.htm">More</a></td>
 22<div class="spirit-nav">
 23<a accesskey="p" href="test.html"><img src="../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../index.html"><img src="../../../../doc/src/images/up.png" alt="Up"></a><a accesskey="h" href="../index.html"><img src="../../../../doc/src/images/home.png" alt="Home"></a><a accesskey="n" href="array.html"><img src="../../../../doc/src/images/next.png" alt="Next"></a>
 25<div class="section">
 26<div class="titlepage"><div><div><h2 class="title" style="clear: both">
 27<a name="document_to_test_formatting.remez"></a><a class="link" href="remez.html" title="Sample Article (The Remez Method)"> Sample Article (The
 28    Remez Method)</a>
 31      The <a href="" target="_top">Remez algorithm</a>
 32      is a methodology for locating the minimax rational approximation to a function.
 33      This short article gives a brief overview of the method, but it should not
 34      be regarded as a thorough theoretical treatment, for that you should consult
 35      your favorite textbook.
 36    </p>
 38      Imagine that you want to approximate some function f(x) by way of a rational
 39      function R(x), where R(x) may be either a polynomial P(x) or a ratio of two
 40      polynomials P(x)/Q(x) (a rational function). Initially we'll concentrate on
 41      the polynomial case, as it's by far the easier to deal with, later we'll extend
 42      to the full rational function case.
 43    </p>
 45      We want to find the "best" rational approximation, where "best"
 46      is defined to be the approximation that has the least deviation from f(x).
 47      We can measure the deviation by way of an error function:
 48    </p>
 50      E<sub>abs</sub>(x) = f(x) - R(x)
 51    </p>
 53      which is expressed in terms of absolute error, but we can equally use relative
 54      error:
 55    </p>
 57      E<sub>rel</sub>(x) = (f(x) - R(x)) / |f(x)|
 58    </p>
 60      And indeed in general we can scale the error function in any way we want, it
 61      makes no difference to the maths, although the two forms above cover almost
 62      every practical case that you're likely to encounter.
 63    </p>
 65      The minimax rational function R(x) is then defined to be the function that
 66      yields the smallest maximal value of the error function. Chebyshev showed that
 67      there is a unique minimax solution for R(x) that has the following properties:
 68    </p>
 69<div class="itemizedlist"><ul class="itemizedlist" type="disc">
 70<li class="listitem">
 71          If R(x) is a polynomial of degree N, then there are N+2 unknowns: the N+1
 72          coefficients of the polynomial, and maximal value of the error function.
 73        </li>
 74<li class="listitem">
 75          The error function has N+1 roots, and N+2 extrema (minima and maxima).
 76        </li>
 77<li class="listitem">
 78          The extrema alternate in sign, and all have the same magnitude.
 79        </li>
 82      That means that if we know the location of the extrema of the error function
 83      then we can write N+2 simultaneous equations:
 84    </p>
 86      R(x<sub>i</sub>) + (-1)<sup>i</sup>E = f(x<sub>i</sub>)
 87    </p>
 89      where E is the maximal error term, and x<sub>i</sub> are the abscissa values of the N+2
 90      extrema of the error function. It is then trivial to solve the simultaneous
 91      equations to obtain the polynomial coefficients and the error term.
 92    </p>
 94      <span class="emphasis"><em>Unfortunately we don't know where the extrema of the error function
 95      are located!</em></span>
 96    </p>
 97<a name="document_to_test_formatting.remez.the_remez_method"></a><h5>
 98<a name="id771060"></a>
 99      <a class="link" href="remez.html#document_to_test_formatting.remez.the_remez_method">The Remez
100      Method</a>
101    </h5>
103      The Remez method is an iterative technique which, given a broad range of assumptions,
104      will converge on the extrema of the error function, and therefore the minimax
105      solution.
106    </p>
108      In the following discussion we'll use a concrete example to illustrate the
109      Remez method: an approximation to the function e<sup>x</sup> over the range [-1, 1].
110    </p>
112      Before we can begin the Remez method, we must obtain an initial value for the
113      location of the extrema of the error function. We could "guess" these,
114      but a much closer first approximation can be obtained by first constructing
115      an interpolated polynomial approximation to f(x).
116    </p>
118      In order to obtain the N+1 coefficients of the interpolated polynomial we need
119      N+1 points (x<sub>0</sub>...x<sub>N</sub>): with our interpolated form passing through each of those
120      points that yields N+1 simultaneous equations:
121    </p>
123      f(x<sub>i</sub>) = P(x<sub>i</sub>) = c<sub>0</sub> + c<sub>1</sub>x<sub>i</sub> ... + c<sub>N</sub>x<sub>i</sub><sup>N</sup>
124    </p>
126      Which can be solved for the coefficients c<sub>0</sub>...c<sub>N</sub> in P(x).
127    </p>
129      Obviously this is not a minimax solution, indeed our only guarantee is that
130      f(x) and P(x) touch at N+1 locations, away from those points the error may
131      be arbitrarily large. However, we would clearly like this initial approximation
132      to be as close to f(x) as possible, and it turns out that using the zeros of
133      an orthogonal polynomial as the initial interpolation points is a good choice.
134      In our example we'll use the zeros of a Chebyshev polynomial as these are particularly
135      easy to calculate, interpolating for a polynomial of degree 4, and measuring
136      <span class="emphasis"><em>relative error</em></span> we get the following error function:
137    </p>
139      <span class="inlinemediaobject"><img src="../images/remez-2.png" alt="remez-2"></span>
140    </p>
142      Which has a peak relative error of 1.2x10<sup>-3</sup>.
143    </p>
145      While this is a pretty good approximation already, judging by the shape of
146      the error function we can clearly do better. Before starting on the Remez method
147      propper, we have one more step to perform: locate all the extrema of the error
148      function, and store these locations as our initial <span class="emphasis"><em>Chebyshev control
149      points</em></span>.
150    </p>
151<div class="note"><table border="0" summary="Note">
153<td rowspan="2" align="center" valign="top" width="25"><img alt="[Note]" src="../../../../doc/src/images/note.png"></td>
154<th align="left">Note</th>
156<tr><td align="left" valign="top">
158        In the simple case of a polynomial approximation, by interpolating through
159        the roots of a Chebyshev polynomial we have in fact created a <span class="emphasis"><em>Chebyshev
160        approximation</em></span> to the function: in terms of <span class="emphasis"><em>absolute
161        error</em></span> this is the best a priori choice for the interpolated form
162        we can achieve, and typically is very close to the minimax solution.
163      </p>
165        However, if we want to optimise for <span class="emphasis"><em>relative error</em></span>,
166        or if the approximation is a rational function, then the initial Chebyshev
167        solution can be quite far from the ideal minimax solution.
168      </p>
170        A more technical discussion of the theory involved can be found in this
171        <a href="" target="_top">online
172        course</a>.
173      </p>
176<a name="document_to_test_formatting.remez.remez_step_1"></a><h5>
177<a name="id771248"></a>
178      <a class="link" href="remez.html#document_to_test_formatting.remez.remez_step_1">Remez Step 1</a>
179    </h5>
181      The first step in the Remez method, given our current set of N+2 Chebyshev
182      control points x<sub>i</sub>, is to solve the N+2 simultaneous equations:
183    </p>
185      P(x<sub>i</sub>) + (-1)<sup>i</sup>E = f(x<sub>i</sub>)
186    </p>
188      To obtain the error term E, and the coefficients of the polynomial P(x).
189    </p>
191      This gives us a new approximation to f(x) that has the same error <span class="emphasis"><em>E</em></span>
192      at each of the control points, and whose error function <span class="emphasis"><em>alternates
193      in sign</em></span> at the control points. This is still not necessarily the
194      minimax solution though: since the control points may not be at the extrema
195      of the error function. After this first step here's what our approximation's
196      error function looks like:
197    </p>
199      <span class="inlinemediaobject"><img src="../images/remez-3.png" alt="remez-3"></span>
200    </p>
202      Clearly this is still not the minimax solution since the control points are
203      not located at the extrema, but the maximum relative error has now dropped
204      to 5.6x10<sup>-4</sup>.
205    </p>
206<a name="document_to_test_formatting.remez.remez_step_2"></a><h5>
207<a name="id771342"></a>
208      <a class="link" href="remez.html#document_to_test_formatting.remez.remez_step_2">Remez Step 2</a>
209    </h5>
211      The second step is to locate the extrema of the new approximation, which we
212      do in two stages: first, since the error function changes sign at each control
213      point, we must have N+1 roots of the error function located between each pair
214      of N+2 control points. Once these roots are found by standard root finding
215      techniques, we know that N extrema are bracketed between each pair of roots,
216      plus two more between the endpoints of the range and the first and last roots.
217      The N+2 extrema can then be found using standard function minimisation techniques.
218    </p>
220      We now have a choice: multi-point exchange, or single point exchange.
221    </p>
223      In single point exchange, we move the control point nearest to the largest
224      extrema to the absissa value of the extrema.
225    </p>
227      In multi-point exchange we swap all the current control points, for the locations
228      of the extrema.
229    </p>
231      In our example we perform multi-point exchange.
232    </p>
233<a name="document_to_test_formatting.remez.iteration"></a><h5>
234<a name="id771387"></a>
235      <a class="link" href="remez.html#document_to_test_formatting.remez.iteration">Iteration</a>
236    </h5>
238      The Remez method then performs steps 1 and 2 above iteratively until the control
239      points are located at the extrema of the error function: this is then the minimax
240      solution.
241    </p>
243      For our current example, two more iterations converges on a minimax solution
244      with a peak relative error of 5x10<sup>-4</sup> and an error function that looks like:
245    </p>
247      <span class="inlinemediaobject"><img src="../images/remez-4.png" alt="remez-4"></span>
248    </p>
249<a name="document_to_test_formatting.remez.rational_approximations"></a><h5>
250<a name="id771441"></a>
251      <a class="link" href="remez.html#document_to_test_formatting.remez.rational_approximations">Rational
252      Approximations</a>
253    </h5>
255      If we wish to extend the Remez method to a rational approximation of the form
256    </p>
258      f(x) = R(x) = P(x) / Q(x)
259    </p>
261      where P(x) and Q(x) are polynomials, then we proceed as before, except that
262      now we have N+M+2 unknowns if P(x) is of order N and Q(x) is of order M. This
263      assumes that Q(x) is normalised so that it's leading coefficient is 1, giving
264      N+M+1 polynomial coefficients in total, plus the error term E.
265    </p>
267      The simultaneous equations to be solved are now:
268    </p>
270      P(x<sub>i</sub>) / Q(x<sub>i</sub>) + (-1)<sup>i</sup>E = f(x<sub>i</sub>)
271    </p>
273      Evaluated at the N+M+2 control points x<sub>i</sub>.
274    </p>
276      Unfortunately these equations are non-linear in the error term E: we can only
277      solve them if we know E, and yet E is one of the unknowns!
278    </p>
280      The method usually adopted to solve these equations is an iterative one: we
281      guess the value of E, solve the equations to obtain a new value for E (as well
282      as the polynomial coefficients), then use the new value of E as the next guess.
283      The method is repeated until E converges on a stable value.
284    </p>
286      These complications extend the running time required for the development of
287      rational approximations quite considerably. It is often desirable to obtain
288      a rational rather than polynomial approximation none the less: rational approximations
289      will often match more difficult to approximate functions, to greater accuracy,
290      and with greater efficiency, than their polynomial alternatives. For example,
291      if we takes our previous example of an approximation to e<sup>x</sup>, we obtained 5x10<sup>-4</sup> accuracy
292      with an order 4 polynomial. If we move two of the unknowns into the denominator
293      to give a pair of order 2 polynomials, and re-minimise, then the peak relative
294      error drops to 8.7x10<sup>-5</sup>. That's a 5 fold increase in accuracy, for the same
295      number of terms overall.
296    </p>
297<a name="document_to_test_formatting.remez.practical_considerations"></a><h5>
298<a name="id771550"></a>
299      <a class="link" href="remez.html#document_to_test_formatting.remez.practical_considerations">Practical
300      Considerations</a>
301    </h5>
303      Most treatises on approximation theory stop at this point. However, from a
304      practical point of view, most of the work involves finding the right approximating
305      form, and then persuading the Remez method to converge on a solution.
306    </p>
308      So far we have used a direct approximation:
309    </p>
311      f(x) = R(x)
312    </p>
314      But this will converge to a useful approximation only if f(x) is smooth. In
315      addition round-off errors when evaluating the rational form mean that this
316      will never get closer than within a few epsilon of machine precision. Therefore
317      this form of direct approximation is often reserved for situations where we
318      want efficiency, rather than accuracy.
319    </p>
321      The first step in improving the situation is generally to split f(x) into a
322      dominant part that we can compute accurately by another method, and a slowly
323      changing remainder which can be approximated by a rational approximation. We
324      might be tempted to write:
325    </p>
327      f(x) = g(x) + R(x)
328    </p>
330      where g(x) is the dominant part of f(x), but if f(x)/g(x) is approximately
331      constant over the interval of interest then:
332    </p>
334      f(x) = g(x)(c + R(x))
335    </p>
337      Will yield a much better solution: here <span class="emphasis"><em>c</em></span> is a constant
338      that is the approximate value of f(x)/g(x) and R(x) is typically tiny compared
339      to <span class="emphasis"><em>c</em></span>. In this situation if R(x) is optimised for absolute
340      error, then as long as its error is small compared to the constant <span class="emphasis"><em>c</em></span>,
341      that error will effectively get wiped out when R(x) is added to <span class="emphasis"><em>c</em></span>.
342    </p>
344      The difficult part is obviously finding the right g(x) to extract from your
345      function: often the asymptotic behaviour of the function will give a clue,
346      so for example the function __erfc becomes proportional to e<sup>-x<sup>2</sup></sup>/x as x becomes
347      large. Therefore using:
348    </p>
350      erfc(z) = (C + R(x)) e<sup>-x<sup>2</sup></sup>/x
351    </p>
353      as the approximating form seems like an obvious thing to try, and does indeed
354      yield a useful approximation.
355    </p>
357      However, the difficulty then becomes one of converging the minimax solution.
358      Unfortunately, it is known that for some functions the Remez method can lead
359      to divergent behaviour, even when the initial starting approximation is quite
360      good. Furthermore, it is not uncommon for the solution obtained in the first
361      Remez step above to be a bad one: the equations to be solved are generally
362      "stiff", often very close to being singular, and assuming a solution
363      is found at all, round-off errors and a rapidly changing error function, can
364      lead to a situation where the error function does not in fact change sign at
365      each control point as required. If this occurs, it is fatal to the Remez method.
366      It is also possible to obtain solutions that are perfectly valid mathematically,
367      but which are quite useless computationally: either because there is an unavoidable
368      amount of roundoff error in the computation of the rational function, or because
369      the denominator has one or more roots over the interval of the approximation.
370      In the latter case while the approximation may have the correct limiting value
371      at the roots, the approximation is nonetheless useless.
372    </p>
374      Assuming that the approximation does not have any fatal errors, and that the
375      only issue is converging adequately on the minimax solution, the aim is to
376      get as close as possible to the minimax solution before beginning the Remez
377      method. Using the zeros of a Chebyshev polynomial for the initial interpolation
378      is a good start, but may not be ideal when dealing with relative errors and/or
379      rational (rather than polynomial) approximations. One approach is to skew the
380      initial interpolation points to one end: for example if we raise the roots
381      of the Chebyshev polynomial to a positive power greater than 1 then the roots
382      will be skewed towards the middle of the [-1,1] interval, while a positive
383      power less than one will skew them towards either end. More usefully, if we
384      initially rescale the points over [0,1] and then raise to a positive power,
385      we can skew them to the left or right. Returning to our example of e<sup>x</sup> over [-1,1],
386      the initial interpolated form was some way from the minimax solution:
387    </p>
389      <span class="inlinemediaobject"><img src="../images/remez-2.png" alt="remez-2"></span>
390    </p>
392      However, if we first skew the interpolation points to the left (rescale them
393      to [0, 1], raise to the power 1.3, and then rescale back to [-1,1]) we reduce
394      the error from 1.3x10<sup>-3</sup>to 6x10<sup>-4</sup>:
395    </p>
397      <span class="inlinemediaobject"><img src="../images/remez-5.png" alt="remez-5"></span>
398    </p>
400      It's clearly still not ideal, but it is only a few percent away from our desired
401      minimax solution (5x10<sup>-4</sup>).
402    </p>
403<a name="document_to_test_formatting.remez.remez_method_checklist"></a><h5>
404<a name="id771737"></a>
405      <a class="link" href="remez.html#document_to_test_formatting.remez.remez_method_checklist">Remez
406      Method Checklist</a>
407    </h5>
409      The following lists some of the things to check if the Remez method goes wrong,
410      it is by no means an exhaustive list, but is provided in the hopes that it
411      will prove useful.
412    </p>
413<div class="itemizedlist"><ul class="itemizedlist" type="disc">
414<li class="listitem">
415          Is the function smooth enough? Can it be better separated into a rapidly
416          changing part, and an asymptotic part?
417        </li>
418<li class="listitem">
419          Does the function being approximated have any "blips" in it?
420          Check for problems as the function changes computation method, or if a
421          root, or an infinity has been divided out. The telltale sign is if there
422          is a narrow region where the Remez method will not converge.
423        </li>
424<li class="listitem">
425          Check you have enough accuracy in your calculations: remember that the
426          Remez method works on the difference between the approximation and the
427          function being approximated: so you must have more digits of precision
428          available than the precision of the approximation being constructed. So
429          for example at double precision, you shouldn't expect to be able to get
430          better than a float precision approximation.
431        </li>
432<li class="listitem">
433          Try skewing the initial interpolated approximation to minimise the error
434          before you begin the Remez steps.
435        </li>
436<li class="listitem">
437          If the approximation won't converge or is ill-conditioned from one starting
438          location, try starting from a different location.
439        </li>
440<li class="listitem">
441          If a rational function won't converge, one can minimise a polynomial (which
442          presents no problems), then rotate one term from the numerator to the denominator
443          and minimise again. In theory one can continue moving terms one at a time
444          from numerator to denominator, and then re-minimising, retaining the last
445          set of control points at each stage.
446        </li>
447<li class="listitem">
448          Try using a smaller interval. It may also be possible to optimise over
449          one (small) interval, rescale the control points over a larger interval,
450          and then re-minimise.
451        </li>
452<li class="listitem">
453          Keep absissa values small: use a change of variable to keep the abscissa
454          over, say [0, b], for some smallish value <span class="emphasis"><em>b</em></span>.
455        </li>
457<a name="document_to_test_formatting.remez.references"></a><h5>
458<a name="id771857"></a>
459      <a class="link" href="remez.html#document_to_test_formatting.remez.references">References</a>
460    </h5>
462      The original references for the Remez Method and it's extension to rational
463      functions are unfortunately in Russian:
464    </p>
466      Remez, E.Ya., <span class="emphasis"><em>Fundamentals of numerical methods for Chebyshev approximations</em></span>,
467      "Naukova Dumka", Kiev, 1969.
468    </p>
470      Remez, E.Ya., Gavrilyuk, V.T., <span class="emphasis"><em>Computer development of certain approaches
471      to the approximate construction of solutions of Chebyshev problems nonlinearly
472      depending on parameters</em></span>, Ukr. Mat. Zh. 12 (1960), 324-338.
473    </p>
475      Gavrilyuk, V.T., <span class="emphasis"><em>Generalization of the first polynomial algorithm
476      of E.Ya.Remez for the problem of constructing rational-fractional Chebyshev
477      approximations</em></span>, Ukr. Mat. Zh. 16 (1961), 575-585.
478    </p>
480      Some English language sources include:
481    </p>
483      Fraser, W., Hart, J.F., <span class="emphasis"><em>On the computation of rational approximations
484      to continuous functions</em></span>, Comm. of the ACM 5 (1962), 401-403, 414.
485    </p>
487      Ralston, A., <span class="emphasis"><em>Rational Chebyshev approximation by Remes' algorithms</em></span>,
488      Numer.Math. 7 (1965), no. 4, 322-330.
489    </p>
491      A. Ralston, <span class="emphasis"><em>Rational Chebyshev approximation, Mathematical Methods
492      for Digital Computers v. 2</em></span> (Ralston A., Wilf H., eds.), Wiley, New
493      York, 1967, pp. 264-284.
494    </p>
496      Hart, J.F. e.a., <span class="emphasis"><em>Computer approximations</em></span>, Wiley, New York
497      a.o., 1968.
498    </p>
500      Cody, W.J., Fraser, W., Hart, J.F., <span class="emphasis"><em>Rational Chebyshev approximation
501      using linear equations</em></span>, Numer.Math. 12 (1968), 242-251.
502    </p>
504      Cody, W.J., <span class="emphasis"><em>A survey of practical rational and polynomial approximation
505      of functions</em></span>, SIAM Review 12 (1970), no. 3, 400-423.
506    </p>
508      Barrar, R.B., Loeb, H.J., <span class="emphasis"><em>On the Remez algorithm for non-linear families</em></span>,
509      Numer.Math. 15 (1970), 382-391.
510    </p>
512      Dunham, Ch.B., <span class="emphasis"><em>Convergence of the Fraser-Hart algorithm for rational
513      Chebyshev approximation</em></span>, Math. Comp. 29 (1975), no. 132, 1078-1082.
514    </p>
516      G. L. Litvinov, <span class="emphasis"><em>Approximate construction of rational approximations
517      and the effect of error autocorrection</em></span>, Russian Journal of Mathematical
518      Physics, vol.1, No. 3, 1994.
519    </p>
521<table xmlns:rev="" width="100%"><tr>
522<td align="left"></td>
523<td align="right"><div class="copyright-footer">Copyright &#169; 2007 John Maddock, Joel de Guzman, Eric Niebler and Matias
524      Capeletto<p>
525        Distributed under the Boost Software License, Version 1.0. (See accompanying
526        file LICENSE_1_0.txt or copy at <a href="" target="_top"></a>)
527      </p>
531<div class="spirit-nav">
532<a accesskey="p" href="test.html"><img src="../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../index.html"><img src="../../../../doc/src/images/up.png" alt="Up"></a><a accesskey="h" href="../index.html"><img src="../../../../doc/src/images/home.png" alt="Home"></a><a accesskey="n" href="array.html"><img src="../../../../doc/src/images/next.png" alt="Next"></a>