remez.html 30 KB

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