Series as a Computational Differentiation Tool

Series is generally the most efficient and accurate way to compute high-order derivative values.  We evaluate the use and misuse of Mathematica's Series, D, and Simplify functions for working with nth-order derivatives.

by Richard D. Neidinger  (8/10/99)

Which is faster, computing a 50th derivative or computing a 50th order Taylor series expansion of the same function?  The following comparison, timed on a Pentium III 500Mhz machine using Mathematica 3.0 under Windows NT, may be surprising.

Timing[ d[50] = D[Tan[Log[x]], {x, 50}];  ][[1
      ]]

Timing[ ser = Series[Tan[Log[x]], {x, 2.0, 50}];  ]
        [[1]]

Both computations can be used to find a value for the 50th derivative at 2.0, but the answers differ in the seventh digit.  Which is more precise?

d[50] /. x -> 2.0 // InputForm

50!*SeriesCoefficient[ser, 50]  // InputForm

Arbitrary precision computation shows that the value computed by Series is only off by 1 in the 15th digit, double the precision of evaluating the symbolic derivative.

The following two sections briefly describe how Series can generate coefficients efficiently without symbolic differentiation and how to make sure our calls to Series do, in fact, use efficient machine number calculations.  The remainder of the paper comparatively evaluates alternatives for finding high-order derivatives in Mathematica.  On first reading, it is possible to skip the rather technical details of the first sections to get to the punch lines found in the charts, beginning with section 3.  The focus is on the question:  Just when does Series have a significant time and/or precision advantage over D, and when does Simplify help?  Detailed answers depends on the specific function expression f(x), the order n of the derivative, and the number m of points a at which the derivative is to be evaluated.  The general answer is that Series is the most efficient and reliable method for finding high-order derivatives.  The advantage usually begins between derivative order 3 and 6 and rapidly grows with n.

Series avoids symbolic differentiation

The algorithms that perform D and Series are as different as nature of the objects that are produced.  The introductory application of D produced a 50th derivative (f^((50)))(x) which is a huge symbolic expression.  The application of Series produced a 50th order series about 2.0 which is, essentially, just a list of 50 numerical coefficients.  The last series coefficient (Subscript[a, 50]) could be found by an evaluation of the symbolic derivative, since \(TraditionalForm
      \`a\_50\  == \ \(\(f\^\((50)\)\)(2.0)\)\/\(50!\)\).  However, the above timing shows that this is definitely not how it's done.  In fact, it is faster and more precise to find (f^((50)))(2.0) by computing (50!) (Subscript[a, 50]).  The Mathematica Book [Wolfram 1996, A.9.5] states that "Series works by recursively composing series expansions of functions with series expansions of their arguments."  For example, consider the series for sin((e^x)) about 2.0.  A calculus student could develop algorithms for the series coefficients of (e^x) about 2.0 and for sin(x) about the constant (e^2.0).  Now, the problem is reduced to composing the series, a task in numerical list manipulation.  While D performs symbolic manipulations on a growing (as order n increases) expression, Series does floating point computations on fixed-length (n) coefficient lists.  For tan(log(x)), the series task is nearly the same as for sin((e^x)) but the symbolic derivatives grow much more rapidly.  Consequently, the series computation is faster, uses less memory, and results in better accuracy.

The series composition idea belongs to a field known as automatic differentiation or computational differentiation.  Within this field, Mathematica's approach is a forward method of arbitrary-order Taylor series propagation using operator overloading (see the introductory article by Rall and Corliss in [Berz 1996], p. 1-18).  Operator overloading refers to fact that Mathematica uses a special data type (or object) for series so that arithmetic operators or scientific functions perform special operations on this data type.  The following section discusses this SeriesData type and shows the specific example Tan of Log of a SeriesData object.  (Readers interested in Mathematica programming techniques will find overloading of operators to be a great example of upvalues.)  At the core, the method depends on a specific algorithm for each operation that will algebraically combine series and an algorithm for each function that will be composed with a series.  Although the specific algorithms used by Mathematica are proprietary information, there are a number of published possibilities.  The previous paragraph described the use of a generic series composition algorithm as in [Nijenhuis & Wilf 1978].  There are more efficient algorithms for specific compositions such as sin(u(x)) where the series for u(x) is already known (computed).  Such algorithms have been around for many years, see [Moore 1979]  p. 24-29 and historical notes on p. 31.  A slightly more complete listing of such formulas is found in [Waldvogel 1984].  An efficient algorithm for tan(u(x)) is given in [Neidinger 1995], where the algorithms are structured to handle multivariable series.  Mathematica's Series function can handle multivariable series but we only study single variable series here.  The simple introductory article [Neidinger 1989] describes similar, though less efficient, algorithms to directly produce lists of derivative values through order n, as opposed to series coefficients.

SeriesData and insuring floating point computation

The output of Series is essentially a list of coefficients.  The default display of the output of Series is a bit misleading, as if it were similar to a symbolic polynomial.  Let's look at a third order series.

serF  =  Series[Tan[Log[x]], {x, 2.0, 3}]

It is really a special data structure called SeriesData where the list of coefficients is the third part.

FullForm[serF]

?SeriesData

The function SeriesCoefficient is a reliable way to access items in the list in SeriesData.

SeriesCoefficient[serF, 3]

Given the structure of SeriesData, one could directly extract the list of coefficients by serF[[3]] and this particular coefficient by serF[[3,4]].  However, this is risky since leading or trailing zeros (including 0.0) are not actually on the list but are accounted for in the parameters nmin and nmax that follow the list in SeriesData.  As an aside, the Normal function makes a more fundamental change than just removing the last + O((((x - 2. ))^4)).  Normal takes a SeriesData object and coverts it to an actual polynomial in x; a purpose for the only non-numeric argument x of the above SeriesData.  On the flip side, input +O[x]^4 does more than just add a term to a polynomial, it converts it to a SeriesData object about zero.

The coefficient lists are certainly not restricted to machine numbers, although machine numbers account for most of the efficiency.  Indeed, a symbolic point of evaluation will cause series composition algorithms to generate a derivative expression; a reminder that we are really creating derivative values, not Taylor polynomial approximations.

Clear[a]3!*SeriesCoefficient[Series[Tan[Log[x]], {x, a, 3}], 3]  // 
      Simplify

D[Tan[Log[x]], {x, 3}]  // Simplify

The difference in these expressions before simplification indicates that the algorithms were very different.  For such a symbolic result, the D is more efficient.  The coefficients in the SeriesData are unsimplified symbolic expressions that grow exponentially as a function of the list index.  Exact numeric expressions can be just as cumbersome; look at the 2nd-order coefficient if we evaluate about integer 2, rather than machine number 2.0.

SeriesCoefficient[Series[Tan[Log[x]], {x, 2, 3}], 2]

Clearly, this destroys the advantage of Series over D.

We need a reliable way to force machine numbers and floating point computations in all of the series manipulations.  It turns out that evaluating about a machine number doesn't always work.  Series about 0.0 usually returns exact expressions, as does Series[Sin[x-2.],{x,2.,5}], as explained in the next paragraph.  Our fix is based on the principle that "Series works by recursively composing series expansions of functions with series expansions of their arguments."  If we take this to the innermost level, any computation of Series begins with a series expansion of x.  Indeed, f[Series[x,{x,a,n}]] produces the same result as Series[f[x],{x,a,n}], in nearly identical time, for any explicit function f.  We will force machine numbers at this innermost level by computing f[N[Series[x,{x,a,n}]]].  As a fringe benefit, this method emphasizes the use of composition algorithms.  Let's follow this schema for Tan[Log[x]].

(serNx  =  N[Series[x, {x, 2, 3}]])  // FullForm

The N makes sure the entire list of coefficients are machine numbers.

Log[serNx]  // FullForm

Tan[Log[serNx]]  // FullForm

These compositions show the nature of the task that is accomplished by series composition algorithms.

The above method avoids the problem that can occur when N only applies to the point a.  Specifically, Series[x,{x,N[2],3}]] creates a SeriesData object with coefficient list {2.0, 1}, where the second coefficient (the first-order term) is the exact integer one.  For the expression x - 2, the constant term 0.0 is dropped leaving a list containing only the exact integer one.

Series[x - 2, {x, N[2], 3}]  // FullForm

This is the result whether subtracting 2 or 2.0, either inside or outside the Series function.  Thus, a naive application of Series to any f(x - 2) about 2.0 will result in exact numeric expressions for the coefficients.  This problem will occur for any Series about 0.0.  Our solution only has one flaw:  in any method that uses SeriesData, it is possible for an exact numeric expression to sneak into the constant term (after any step that has a zero constant term).  This is annoying but insignificant in terms of efficiency.  The above method reliably uses floating point computations for all other coefficients.

To summarize, we define functions that will always use floating point computations and Mathematica's series composition algorithms.  The first computes series and the second computes an nth-order derivative value.

numSeries[expr_, {x_, a_, n_}]  :=  
      expr  /.  x -> N[Series[x, {x, a, n}]]

\(seriesCoefD[expr_, {x_, a_, n_}]\  := 
	
      \(n!\)*SeriesCoefficient[numSeries[expr, {x, a, n}], n]\)

The function numSeries should not perform differently from the built-in Series, except in those cases where the built-in function results in exact numeric expressions for the coefficients.  On the other hand, direct derivative functions use symbolic manipulation that is very different from seriesCoefD.

Alternatives for high-order derivatives

We evaluate three basic alternatives for finding nth-order derivative values.   The alternative that we'll call "evaluating D" is given by:

(d[n] = D[f[x], {x, n}]) /. x -> N[a]

The following alternative is referred to as "evaluating after Simplify":

(sd[n] = Simplify[d[n]]) /. x -> N[a]

Finally, the function name "seriesCoefD" will encapsulate the process of finding the nth-order series, using numerical coefficients, about a and multiplying the nth coefficient by n!.  A call to this function, defined at the end of the preceding section, has the form:

seriesCoefD[f[x], {x, a, n}]

Usually, the intrinsic expression n!*SeriesCoefficient[Series[f[x],{x,N[a],n}],n]
will perform the same as a call to seriesCoefD.  In a few examples (notably a= 0 or 0.0), our modification of this syntax is required to trigger the machine number computations.

While these alternatives give the basic strategies, we will break them down and consider modification as appropriate.  Section 6 considers evaluation at multiple points, as in plotting, where the derivative is thought of more as a function than as a value.  There, the one-time symbolic D and Simplify steps must be separated from the repeated point evaluations (substitutions).  Section 7 considers the idea of iterating D together with Simplify to produce the sequence of derivative expressions.

All of the comparisons involve timing and precision measurements that are very much dependent on hardware, as well as the specific version of Mathematica.  To be consistent, the measurements must be performed on one platform and the discussion based on one set of results.  Thus, we will simply report or chart the results for one machine.  The actual code to create the charts is found in the accompanying notebook SeriesFigures.nb.  To encourage testing on other platforms and other examples, the accompanying notebook SeriesTest.nb will take one function and perform all of the comparisons discussed herein, controlled by time and space parameters to do what is possible within a reasonable amount of time (several minutes) and RAM.  The absolute measurements and specific details will change between machines and as hardware improves.  However, we expect similar relative and qualitative phenomena, such as the shape of chart curves and order-of-magnitude comparisons.  New versions of Mathematica, with improved implementation of Series, D or Simplify, could change some conclusions, especially for specific  f(x) cases.  Still, symbolic manipulation is inherently different from series composition algorithms on lists of machine numbers.  Comparisons are intended to highlight this inherent difference and how it can be exploited, rather than an evaluation of hardware or software.

Two technical notes on timings:  The very first evaluation of Series during a session will take a little longer while Series routines are loaded into the session Kernel.  Second, we use an average timing function, discussed in Section 5, to measure time intervals smaller than the minimum time interval recorded on the computer system.

Precision of derivative values

Precision should be studied first since the computations are pointless if the result is not reliable.  Often, all three ways of computing nth-order derivative values result in 13 to 16 digits of precision, i.e., the number of significant decimal digits.  Such results are as good as can be expected, given the massive number of operations in high-order derivatives and assuming $MachinePrecision is 16.  Differences within this range might be curious but are not important.  However, for some functions and points of evaluation, the resulting precision is much worse.  In all tested examples where precision of at least one method falls below 13, seriesCoefD is the best, evaluating D is the worst, and evaluating after Simplify can range between the two.

In the following figure, all data points represent the precision in the computed value of the nth-derivative of ArcTan evaluated at 2.3.  Each curve shows the precision of one method as a function of n.  To compute the digits of precision we first find the relative error of each value, relErr= Abs[(value-trueValue)/trueValue], which might be around the order of (10^(-14)).  The digits of precision are then given by -Log[10,relErr].  To find trueValue, we use arbitrary precision (starting with 32 or 64 digits) so that the computed precision of the result is more than 16 digits.  The full code to reproduce all plots is available with the electronic version in SeriesFigures.nb.

The most dramatic observation is that the result of evaluating D has lost all precision by order n=45; it is totally useless.  To gain some insight, consider what is happening for the specific function ArcTan.  This example is a, relatively well-behaved, rational function in disguise since all derivatives are a linear combinations of powers of x over powers of ((x^2) + 1).  The coefficients of these terms grow to huge integers that alternate in sign, a likely source of loss of significance.  Simplify gets a common denominator and factors the numerator.  The huge integers are avoided by seriesCoefD; it works with machine number series coefficients which are much smaller (decreasing steadily to around (10^(-20))).  Finally, seriesCoefD multiplies by n! to find the derivative value.  The true value of the 50th derivative at 2.3 is 1.89578538494699202`*^45.  It is very common for high-order derivatives values of any function to be huge.  In particular, any expression involving division is converted into Power[#,-1] in it's internal FullForm and consecutive derivatives of (x^(-1)) produce a factorial.

Simplification does not always improve the precision of the result as seen in the nth-derivative of Exp[-x^4] evaluated at 2.3.  Derivatives of this are polynomials times (e^(-(x^4))).  Actually, the result of D leaves (e^(-(x^4))) distributed into every term and Simplify factors it out, along with a constant and/or x.  This does reduce the evaluation time (since (e^(-(x^4))) is only evaluated once) but doesn't help with the imprecision.  Again, the polynomial coefficients grow to huge integers that alternate in sign.  The alternation is critical since derivatives of Exp[x^4] are identical except for the alternating sign and result in constant precision of 14 for all three methods.

The digits of precision do depend on the point of evaluation.  For a<1, derivative values of both ArcTan and Exp[-x^4] avoid the loss of significance as the powers of a counteract the huge integer coefficients.  Both functions also show a curious, though unimportant, difference between a=2.00 and a=2.01 in the curve for evaluating after Simplify.  These curves are very similar after they drop below precision 13, but the curve for 2.01 gradually descends to precision 13 while the curve for 2.00 is nearly perfect precision until it abruptly drops.  Yet, the simplified expressions didn't change in any way, just the point of evaluation!  At a=1, the true values of many derivatives of ArcTan are zero, making it impossible to measure digits of precision.

The function (Sin[x](1-Cos[x]))/(Exp[x^3]-1), found in a calculus text exercise, is a rare example where even seriesCoefD struggles with digits of precision.  Still, it results in about 4 meaningful digits in the 50th derivative when just evaluating D loses all precision at the 14th derivative.  Simplification is not a practical alternative since Simplify[d[3]] took over 30 seconds on my computer, with the same precision in the resulting value, and Simplify[d[4]] took over 30 minutes after crashing the kernel on the previous trial.

One should not conclude that this is typical of complicated functions.  For the monstrosity (x*Csc[x])/Exp[ArcTan[Log[x]]] at 2.0, both evaluating D and seriesCoefD result in over 13 digits of precision in the 50th derivative.  However, seriesCoefD uses much less time and space to arrive at the result, as we see in the next section.  The cost of Simplify is prohibitive for high-order derivatives of most complicated functions.

The bottom line is that seriesCoefD results in much more reliable precision than the more straight-forward alternatives.  In the vast majority of tests, seriesCoefD correctly computes the 50th derivative with over 13 of the 16 possible digits.  In a worst case scenario, such as the preceeding figure, seriesCoefD is the only practical alternative.

Timings for one high-order derivative value

Setting the issue of precision aside, how does the time-cost of seriesCoefD compare to the alternatives for finding the value of a high-order derivative at a point?  For most functions, seriesCoefD is faster for order n>5 with a dramatically growing gap.  For computing one value, the large time-cost of Simplify can never be recouped by the savings in the smaller evaluation time.  Thus, we compare only the time-cost of seriesCoefD[f[x],{x,a,n}] with (d[n]=D[f[x],{x,n}])/.x->a.  Timings do not depend on the specific point a and the dominance of seriesCoefD is seen for all except simple (non-growing derivative) expressions, unlike precision where there is a delicate dependence on both a and f(x).

For simple function expressions, such that consecutive derivative expressions do not grow in size, symbolic manipulation is very effective for all orders of derivatives.  This includes functions with simple derivative patterns such as polynomials, Sin[x], Exp[-3x], and Log[x].  Even some non-trivial expressions such as Exp[-2x]*Sin[3x] and (Sin[x])^2 can be seen to have simple one or two term expressions for any order derivative.  It is possible, though uncommon, for seriesCoefD to have a little advantage in some of these cases.  Still, evaluating D is not a problem, since growing symbolic derivative expressions are the primary source of the time, space, and precision problems in high-order derivatives.  One interesting observation is that seriesCoefD appears much less efficient for Sin[x] than for Exp[x].

A typical function, with growing consecutive derivative expressions, is Sin[x^2].  The following graph comparing timings of the alternatives shows an unmistakable difference in the growth rate.  A zoom into the lower-order derivatives, shows that evaluating D starts out a little faster than seriesCoefD and the curves cross around order n=5.  This is a common but conservative result.  The function Exp[-x^2] has an even more dramatic gap with seriesCoefD timings nearly flat at 0.001 sec while evaluating D grows to 0.38 sec.

Concern about the reliability of the timings motivates a few comments about how our graphs were generated.  The Timing function cannot be directly used for small time intervals that fall under or even near $TimeUnit, the minimum time interval in seconds recorded on your computer system.  Hence, we use Timing on a Do loop of repetitions, over a time interval that must exceed 50*$TimeUnit seconds.  Finally, our average timing function divides by the number of repetitions (code found in SeriesFigures.nb).  Repeated trials of the average timing function are consistent enough to preserve the curves, though individual data points do move around a bit.  The timings are not affected by numeric caching (see [Sofroniou 1996]) which saves values of the N function; specifically, N[NumericQ,Clear] was inserted between every evaluation within the loop with no effect.

For complicated functions, the timing difference becomes a serious issue in practice.

Using D to compute the 50th derivative of (x*Csc[x])/Exp[ArcTan[Log[x]]] took about 2 hours of CPU time and used 250 MB of RAM.  (With 384 MB of RAM, I don't think virtual memory slowed the calculation.)  In comparison, seriesCoefD computed the 50th derivative value in 0.06 seconds CPU time and the numerical SeriesData object has a ByteCount of about 1 KB.  At least, both methods for finding this 50th derivative value, and that of Sin[x^2] from the previous chart, result in over 13 digits of precision.

Evaluation at multiple points

Instead of just seeking a value at one point, we now consider the goal to be a derivative function, for evaluation at multiple points.  The three alternatives methods result in different definitions of nth derivative function for a specifically defined function f[x]; here n=50 for clarity.

\(df50[x_] = D[f[x], {x, 50}]; 
df50Sim[x_] = Simplify[df50[x]]; 

    df50Ser[a_] := seriesCoefD[f[x], {x, a, 50}]; \)

The use of Set (=), in the first two definitions, performs the symbolic computation once and for all.  Subsequent calls, such as df50[2.3], only involve the cost of substituting for x and evaluating the expression.  In contrast, every call to df50Ser generates the entire series "from scratch," so that df50Ser[2.3] and df50Ser[2.4] do not use any common result.  Nevertheless, for machine number input, these functions all seek the same numerical result.  The moral of the Section 4  is that df50Ser is the most numerically reliable of these alternatives.  In terms of time-efficiency, the superior method is sensitive to the specific function expression f(x).

We consider three rough categories for f(x):  (1) Complicated functions where series is the only practical alternative; (2) borderline functions where any of the alternatives can be the most efficient, depending on the order n of the derivative and the number m of point evaluations; and (3) primarily trigonometric or trivial functions where evaluating D is more efficient for multiple point evaluations.

Complicated functions where series is the only practical alternative.

For complicated function expressions, use seriesCoefD to handle high-order derivative functions.  Substitution into complicated symbolic derivatives (just df50[2.3]) can actually take longer than the series computation (df50Ser[2.3]), making the symbolic operation D pointless from a time-saving perspective.  Simplify quickly becomes impossible as the order rises above five.  For example, using seriesCoefD is the only practical way to produce the following interesting plot of a 20th derivative on my computer.  (The 50th derivative on interval [1,2] is just as successful and interesting, though it takes five times as long.)

\(f[x_] = \((Sin[x] \((1 - Cos[x])\))\)/\((Exp[x^3] - 1)\); 

    df20Ser[a_] := seriesCoefD[f[x], {x, a, 20}]; 

    Timing[Plot[df20Ser[x], {x,  .5, 2.5}]; ]\)

Even if you have the time to wait for computing and evaluating D, the precision problem (shown in the last chart of Section 4) would make the plot meaningless.  It is interesting to see that the plots agree for n=10, although using D takes longer than seriesCoefD.  Then, for n=15, the plots start out the same on the left but the plot using D degrades in precision as x increases to the right.

Borderline functions where any of the alternatives can be the most efficient, depending on order n and number of point evaluations m

For some simpler expressions f(x), either evaluating D, evaluating after Simplify, or seriesCoefD can be the most efficient, depending on the order n of the derivative and the number m of point evaluations.  We focus on the example f(x) = (e^(-(x^2))), where the internal implementations of D, Simplify, and Series are all especially efficient.  For this example, the symbolic operations have similar time costs with the rare occurrence that the cost of Simplify actually falls below the cost of D.

These symbolic time-costs are at a different order of magnitude from seriesCoefD which is very comparable to the evaluation (substitution) times.

The timing for substitution into d[n] nears the time-cost of the series computation around n=50, making the evaluating D option ineffective for any number of point evaluations if n>50.  For lower values of n, there is a more delicate balance of up-front symbolic cost and subsequent time savings in evaluations (substitutions).

Using these timings, we can create a partition of the (n,m)-plane into regions where each of the alternatives are superior.

To explain the curves in the this chart, we name the five different time-costs.

({{TIME IN SECONDS, FOR COMPUTATION OF}, {dTime[[n]], d[n] = D[f[x], {x, n}]}, {simpTime[[n]], sd[n] = Simplify[d[n]]}, {dvalTime[[n]], d[n] /. x -> a}, {sdvalTime[[n]], sd[n] /. x -> a}, {serTime[[n]], seriesCoefD[f[x], {x, a, n}]}})

For seriesCoefD to be the most time efficient, requires that the cost  m*serTime < dTime + m*dvalTime  and  m*serTime < dTime + simpTime + m*sdvalTime.  Solving for m, this means

yielding the break-even curve for seriesCoefD.  Similarly, evaluating after Simplify is the most efficient when  \(TraditionalForm
      \`dTime + simpTime + m*sdvalTime < dTime + m*dvalTime\)  and  dTime + simpTime + m*sdvalTime < m*serTime.  Solving for m, yields the break-even curve for Simplify,

These two curves divide the (n,m)-plane into the three regions.

Other borderline functions, that behave similarly to E^(-x^2), include simple rational functions, inverse trigonometric functions, and other simple functions that are not primarily trigonometric.  The functions ArcTan and ArcSin have smaller (n,m)-plane regions where seriesCoefD and Simplify are superior.  For both, the break-even curve for seriesCoefD only rises to about 125 at n=50, creating a horizontal strip around m=200 where evaluating D is the most efficient for all orders n<50.

Primarily trigonometric or trivial functions where evaluating D is more efficient for multiple point evaluations.

Evaluating D for multiple point evaluation can be relatively efficient, trivially when derivative expressions do not grow in size and, more interestingly, because Series and Simplify are less efficient for trigonometric functions.  Consider primarily trigonometric functions such that consecutive derivative expressions do grow but not so much that they become impractical.  For the example Sin[x^2], both Series and Simplify are more than an order of magnitude slower than for Exp[-x^2], while D is a little faster.  The charts for E^(-x^2), just above, show that the series computation timing is on a comparable scale with the evaluation timings, way below the time-cost of symbolic operations.  In contrast, the chart for Sin[x^2] in Section 5 shows a much larger series computation time that is on the same scale as the symbolic differentiation time.  For the 50th derivative of Sin[x^2], the series computation time is about one-sixth of the of the symbolic differentiation time.  Since the evaluation time is relatively insignificant, the cost of the symbolic 50th derivative is justified after just six point evaluations; so evaluate D for m>6.  Simplification is another matter, it would take thousands of point evaluations to justify the high time-cost of Simplify.

Often, the time-cost of simplification can be greatly reduced by tailoring it to the problem.  The high cost of simplification for the 50th derivative of Sin[x^2] is cut by a factor of 10 (from around 4.0 sec to 0.4 sec) by the Simplify option Trig->False.  Surprisingly, the resulting simplified expressions are identical, primarily factoring out Sin[x^2] and Cos[x^2]; the extra time was wasted in looking for applicable trigonometric identities that were never found.  Still, it would take about 1000 point evaluations to justify the lower simplification time of 0.4, since each simplified evaluation saves only about 0.0004 sec.  In a different example, Simplify does successfully use trigonometric identities on derivatives of Tan[x] to lower the evaluation time for order n<30.  However, around the 30th derivative of Tan[x], the simplification time makes a dramatic jump up while at the same time losing all of the advantage in evaluation time.  This problem is "fixed" by setting the Simplify option ComplexityFunction->LeafCount as discussed in The Mathematica Book [Wolfram 1996, end of 3.3.9].  Still, it takes over a thousand point evaluations to make it worthwhile.  There is a similar break-down in Simplify at the 18th derivative of the rational function ((x^3) + x + 1)/((x^4) + (x^2) + 2), even with ComplexityFunction->LeafCount.  The solution for the rational function is to use Together instead of Simplify.  The relatively low cost of Together does make it more efficient for only 20 point evaluations at n=6, growing to around 1000 point evaluations at n=30.  We conclude that customizing Simplify (or alternatives) may be helpful but it takes a lot of foresight (or, more likely, hindsight) and may still be difficult to justify in term of time savings in evaluation.

Derivative sequences, iterative D and iterative Simplify

Our previous analysis of finding the 50th derivative is essentially equivalent to finding the entire sequence of derivatives, from 1 to 50, a goal that is perhaps more desirable.  However, the above notation should not directly be used to produce this sequence of results.  A naive repetition of above notation would compute the sequence of 50 symbolic derivatives by Do[ d[n]=D[f[x],{x,n}], {n,50}] and the list of 50 derivative values by Table[ seriesCoefD[f[x],{x,a,n}], {n,50}].  More appropriately, each task can be accomplished at the time-cost of only the very last computation in each of these repetitions.  In contrast, the great expense of Simplify cannot usually improve upon Do[ sd[n]=Simplify[d[n]], {n,50}].  While iterative use of D is helpful, mixing simplification into the iteration does not produce savings and can backfire.

The bulk of every call to seriesCoefD is a call numSeries to produce the entire list of series coefficients.  Thus, the list of derivative values should be produced from one call to numSeries (recall that this is equivalent to a call to Series that insists on using machine numbers).  If ser=numSeries[f[x],{x,a,50}], then the list of derivative values are given by Map[(#!*SeriesCoefficient[ser,#])&, Range[50]].  The time-cost of this latter step is fixed (around 0.0017 seconds on my computer), independent of f[x], and is only relevant for the very fastest series computations.  There are ways to speed this latter step but, in the overall picture, it is not very significant.  The cost of seriesCoefD for order n is roughly equivalent to the cost of computing the list of derivative values.

The sequence of symbolic derivatives should be computed with an iterative application of a first derivative, using either NestList or a Do loop as seen in this example.

\(d[0] = Exp[\(-x^2\)]; 
Timing[Do[d[n] = D[d[n - 1], x], {n, 50}]]\)

The time-cost of this iteration (with either Do or NestList) is essentially the same as just computing the 50th derivative directly, found in the Section 6 graph of D timings for E^(-x^2).  The direct and iterative 50th derivative expressions are identical.  Thus, our above charts showing the cost of the nth derivative can be interpreted as the cost of derivatives up through n.

Although iterative application of D works well, mixing Simplify into each step does not.  It may seem that the idea would help, since many problems of inefficiency and imprecision do stem from the huge size of symbolic derivatives.  Rather than waiting to Simplify at the end, after iterative D generates huge expressions, it is natural to try using Simplify to reduce at each step.  This process, shown below, will be referred to as "iterative Simplify."

\(sd[0]\  = \ Exp[\(-x^2\)]; 

    Timing[Do[sd[n] = Simplify[D[sd[n - 1], x]], {n, 50}]; ]\)

This timing is an order of magnitude slower than Simplify[D[Exp[-x^2],{x,50}] and is about the same cost as just using Simplify on each symbolic derivative d[n].

Timing[Do[sd[n] = Simplify[d[n]], {n, 50}]; ]

For most functions, iterative Simplify has about the same cost as iterative D plus using Simplify on each result; a cost that is higher than anything else proposed in this paper.  The only examples where iterative Simplify produces real savings were when the original expression had cancellations that could (not just rewrite but) eliminate much of the expression.  For many complicated expressions, the idea actually backfires and produces larger expressions!

To see how Simplify can actually complicate subsequent derivatives, consider the expressions created in one step of iterative process.  The size of the expressions are measured by LeafCount.  Let g denote the symbolic fifth derivative of Exp[-x^2], called d[5] above.

g = D[Exp[-x^2],{x,5}]

Let h denote this expression after simplification, called sd[5] above.

h  =  Simplify[g]

Which is better to differentiate?  Let's look at the sequence of sizes given by the original expression, it's derivative, and the simplification of the derivative.

Map[LeafCount, {g, D[g, x], Simplify[D[g, x]]}]

Map[LeafCount, {h, D[h, x], Simplify[D[h, x]]}]

Although h starts out smaller, it's derivative is actually larger than the derivative of g.  In the end, the Simplify function returns to the same simplified sixth derivative expression by either method.  Using h can make the D operation more efficient but forces the Simplify function to work on a larger expression.  Thus, it is not helpful to mix Simplify with the iteration of D.  Exp[-x^2] is not a special function; many functions (without cancellation) appear to create a larger derivative after Simplify.

A more serious problem can occur for larger expressions where the above phenomena is magnified and Simplify does not succeed in returning to the same simplified derivative.  In this case, iterative Simplify can be disastrous, actually producing larger simplified derivatives and making everything less efficient.  This usually appears when you start pushing the boundary of what Simplify can handle, where the options become very important and my kernel tends to crash.  Ironically, the problem can be exasperated by the option ComplexityFunction->LeafCount.  (Try g = D[(x*Csc[x])/Log[ArcTan[Exp[x]]],{x,2}] if you want to experiment.)

For simpler functions, Simplify usually succeeds in producing the same expressions whether iterated or not, with about the same time-cost whether iterated or not.  It is rare for iterative Simplify to have any advantage.

Conclusion

In all of the widely varying tests for producing high-order derivatives, only one method always performs reasonably in terms of precision, time and space -- that method uses Series with lists of coefficients that are machine numbers.  Usually,  Series[f[x],{x,N[a],n}] results in machine number coefficients, though some cases require f[N[Series[x,{x,a,n}]]].  To produce a derivative value, multiply n! times the nth coefficient from the result (as in seriesCoefD).  Series uses composition algorithms based on lists rather than symbolic derivatives to produce the results.  By using numerical lists, computations for Series stay reasonable for any function and order, whereas symbolic computations using D or Simplify vary widely in complexity.  For complicated functions, Simplify can become impractical for derivatives as low as order 3.  Just evaluating the symbolic results of D (or Simplify) at a machine number can produce precision problems, especially over order 10.  For a function that is not obviously trivial, we suggest using Series to find derivative values for any order over 10.  Tests were performed up to order 50 but Series makes derivatives of several hundred

Of course, all of the methods do have a purpose and the choice depends on the goal.  Certainly, one should not use symbolic derivatives to find Taylor series expansions of functions.  This is an important conclusion for anyone implementing a series operation in any computing environment.  At the other extreme, the goal may be to study patterns or operations in lists of symbolic derivatives.  This goal clearly calls for symbolic differentiation and, probably, simplification.  Iterative use of Simplify does not help produce a simple list of derivatives.  Just simplify the expressions from applying D iteratively, if that is the goal, and beware that simplification is a costly operation.  If  the visible expression is not really required, Series can be used to efficiently and accurately form high-order derivative functions that return machine number values.  A single high-order derivative value is computed by a Series about point a.  This "function of a" (as in Section 6) defines the true derivative function, not a series approximation, that can be plotted.  Series also naturally produces (as in Section 7) the sequence of derivative values of order 0 through n, if that is the goal.

The specific relative efficiency of using Series, D or Simplify depends greatly on the function, order of derivative, and number of points of evaluation.  The conclusions roughly correspond to four categories of functions.  For functions where the consecutive derivative expressions do not grow in size, just evaluating D works fine; there is little point in using Series or Simplify.  At the other extreme are very complicated functions, where Series is the only practical method.  The time-cost of Simplify can be prohibitive and generating the Series about point a can actually beat the time-cost and precision of just substituting a into such complicated symbolic derivatives.  Thus, for complicated functions, Series is more efficient for any number of point evaluations, usually starting at a derivative of order between 3 and 6.  This leaves a middle group of functions that are fairly simple but have consecutive derivatives that do grow; for such functions, Series is the more efficient method to accomplish the goal of a high-order derivative value at a single point.  However, the category splits when considering the most efficient method for multiple point evaluations.  Series is less effective for primarily trigonometric functions in this middle group, so that computing symbolic D is time-efficient for more than just a handful of evaluations.  Coincidentally, the same functions usually take thousands of point evaluations to justify the high cost of Simplify.  The remaining category, of middle group functions that are not primarily trigonometric, has a delicate balance of efficiencies for each method of multiple point evaluations.  Series can be effective for dozens or even hundreds of point evaluations as the order n increases.  When several hundred point evaluations are desired, Simplify can be time-efficient for such functions.

In later versions of Mathematica, and with later hardware, all of the methods should become more efficient.  There are indications that the series manipulation algorithms are not the most efficient.  For example, algorithms as in [Waldvogel 1984] show that sin of a series should take roughly twice as long as exp of a series, yet it takes well over 20 times as long in Mathematica 3.0.  Thus, improvements may make Series even more effective for multiple point evaluations; eliminating the "primarily trigonometric" category above.  Still, the broad conclusions should remain the same because of the inherent difference between series methods that manipulate lists and symbolic methods that manipulate function expressions.

References

Wolfram, Stephen.  1996.  The Mathematica Book. Third edition.  Wolfram Media–Cambridge University Press, NewYork.

Sofroniou, Mark.  1996.  Numerics in Mathematica 3.0.  The Mathematica Journal, 6(4): 64-73.

Moore, Ramon E.  1979.  Methods and Applications of Interval Analysis.  SIAM Studies in Applied Mathematics 2.  SIAM, Philadelphia, PA.  p. 24-31.

Waldvogel, J.  1984.  Der Tayloralgorithmus.  Journal of Applied Mathathematics and Physics (ZAMP), 35(November): 780-189.

Neidinger, Richard D.  1989.  Automatic Differentiation and APL.  College Mathematics Journal, 20(May): 238-251.

Neidinger, Richard D.  1995.  Computing Multivariable Taylor Series to Arbitrary Order.  ACM APL95 Proceedings, APL Quote Quad, 25(4): 134-144.

Berz, Martin, et. al.  1996.  Computational Differentiation: Techniques, Applications, and Tools.  SIAM Proceedings Series.  SIAM, Philadelphia, PA.

Nijenhuis, A. and Wilf, H.S.  1978.  Composition of Power Series (POWSER) in Combinatorial Algorithms.  Academic Press, New York.

About the author

Electronic subscriptions