In our imaginary image processing example, we would thus get a speedup, and instead of the poor user waiting over 3 minutes, they wait less than half a second.

If your function is sufficiently simple, you can also code it up using Cython.  When this works, it is blazingly fast, though you have the overhead of compiling the code in the first place.  Also, for optimized function evaluation, you have to use Cython like the C language, so to get good performance you must know the basics of C programming:

• C datatypes like float and int
• the C library, e.g., math.h

Damn, that was fast!!!!  This is over 1000 times faster than our original symbolic function f!

FTW: In our imaginary image processing example, we would thus get a speedup of a factor of 1000, and instead of the poor user waiting over 3 minutes, they wait 0.184 seconds.  

By the way, fast_float works with symbolic functions of any number of variables:

To figure out how to use "pow", I typed 

       man 3 pow

at the bash shell command line.

Problem:

1. Write symbolic code, pure Python code, fast_float, and Cython code, as above, to evaluate the expression $$f(x,y,z) = x\cos(y) - xyz + x^3 + y^{10}$$ efficiently.  
2. Compare timings of evaluation of this function using each method.

Root Finding: A basic numerical algorithm

Suppose $f(x)$ is a continuous real-valued function on an interval $[a,b]$ and that $f(a)f(b) < 0$, i.e., there is a sign change. Then by calculus (the Intermediate Value Theorem) $f$ has a zero on this interval.

The main goal is to understand various ways of numerically finding a point $c \in [a,b]$ such that $f(c)$ is very close to 0.

Bisection

Perhaps the simplest method is bisection. Given $\epsilon>0$, the following algorithm find a value $c\in [a,b]$ such that $|f(c)| < \epsilon$.

1. If $|f((a+b)/2)| < \epsilon$, output $c=(a+b)/2$.
2. Chop the interval $[a,b]$ in two parts $[a,c] \cup [c,b]$.
3. If $f$ has a sign change on $[a,c]$ replace $[a,b]$ by $[a,c]$. Otherwise, replace $[a,b]$ by $[c,b]$. Then go to step 2.

At every step in the algorithm there is a point $c$ in the interval so that $f(c) = 0$. Since the width of the interval halves each time and each interval contains a zero of $f$, the sequence of intervals will converge to some root of $f$. Since $f$ is continuous, after a finite number of steps we will find a $c$ such that $|f(c)| < \epsilon$.

Double Precision Root Finding Using Bisection

The bisect we implemented above is slower:

Scipy's bisect is a lot faster than ours, of course. It's optimized code written in C. Let's look.

1. Extract the scipy spkg:

   wstein> ls spkg/standard/scipy*
   spkg/standard/scipy-0.7.p4.spkg
   spkg/standard/scipy_sandbox-20071020.p4.spkg
   wstein> tar jxvf spkg/standard/scipy-0.7.p4.spkg  # may have to get from http://sagemath.org/packages/standard/

2. Track down the bisect code. This just takes guessing and knowing how to get around source code. I found the bisect code in

     scipy-0.7.p4/src/scipy/optimize/Zeros/bisect.c
   

   Here it is:

   /* Written by Charles Harris charles.harris@sdl.usu.edu */
   
   #include "zeros.h"
   
   double
   bisect(callback_type f, double xa, double xb, double xtol, double rtol, int iter, default_parameters *params)
   {
       int i;
       double dm,xm,fm,fa,fb,tol;
       
       tol = xtol + rtol*(fabs(xa) + fabs(xb));
       
       fa = (*f)(xa,params);
       fb = (*f)(xb,params);
       params->funcalls = 2;
       if (fa*fb > 0) {ERROR(params,SIGNERR,0.0);}
       if (fa == 0) return xa;
       if (fb == 0) return xb;
       dm = xb - xa;
       params->iterations = 0;
       for(i=0; iiterations++;
           dm *= .5;
           xm = xa + dm;
           fm = (*f)(xm,params);
           params->funcalls++;
           if (fm*fa >= 0) {
               xa = xm;
           }
           if (fm == 0 || fabs(dm) < tol)
               return xm;
       }
       ERROR(params,CONVERR,xa);
   }
   

We can implement this directly in Sage using Cython.

Problem: Click on each of the .html link above to see the autogenerated code (double click on a line of Cython to see the corresponding C code).

Note: Putting r after a number keeps the Sage preparser from turning it into a Sage arbitrary precision real number (instead of a standard Python float).

Newton's Method

1. Choose a starting value $c$.
2. Let $c = c - f(c)/f'(c)$.
3. If $|f(c)| < \epsilon$ output $c$ and terminate. Otherwise, go to 2.
   /// Often f is differentiable and its values can help us give an even more efficient root finding algorithm called Newton's Method.
   
   The equation of the tangent line to the graph of f(x) at the point (c,f(c)) is
   y = f(c) + f'(c)(x-c).
   This is because the above line goes through the point (c,f(c)) and is a line of slow f'(c).
   
   Newton's method: Approximate f(x) by the line above, find a root of that line, then replace c by that root. Iterate. In particular, by easy algebra, the root of the linear function is
   c - \frac{f(c)}{f'(c)}.
   
   
   As an algorithm:
   1. Choose a starting value c.
   2. Let c = c - f(c)/f'(c).
   3. If |f(c)| < \epsilon output c and terminate. Otherwise, go to 2.
      }}} {{{id=73| x = var('x') @interact def _(f = cos(x) - x, c = float(1/2), a=float(0), b=float(1)): show(plot(f,a,b) + point((c,f(c)),rgbcolor='black',pointsize=20) + \ plot(f(c) + f.derivative()(c)*(x-c), a,b,color='red'), a,b) /// }}} {{{id=26| def newton_method(f, c, eps, maxiter=100): x = f.variables()[0] fprime = f.derivative(x) try: g = f._fast_float_(x) gprime = fprime._fast_float_(x) except AttributeError: g = f; gprime = fprime iterates = [c] for i in xrange(maxiter): fc = g(c) if abs(fc) < eps: return c, iterates c = c - fc/gprime(c) iterates.append(c) return c, iterates html("

      Double Precision Root Finding Using Newton's Method

      ") var('x') @interact def _(f = x^2 - 2, c = float(0.5), eps=(-3,(-16..-1)), interval=float(0.5)): eps = 10^(eps) print "eps = %s"%float(eps) time z, iterates = newton_method(f, c, eps) print "root =", z print "f(c) = %r"%f(z) n = len(iterates) print "iterations =", n html(iterates) P = plot(f, z-interval, z+interval, rgbcolor='blue') h = P.ymax(); j = P.ymin() L = sum(point((w,(n-1-float(i))/n*h), rgbcolor=(float(i)/n,0.2,0.3), pointsize=10) + \ line([(w,h),(w,j)],rgbcolor='black',thickness=0.2) for i,w in enumerate(iterates)) show(P + L, xmin=z-interval, xmax=z+interval) /// }}} {{{id=25| float(sqrt(2)) /// 1.4142135623730951 }}} {{{id=93| /// }}} {{{id=117| scipy.optimize.newton? /// }}} {{{id=78| x = var('x') f = x^2 - 2 g = f._fast_float_(x) gprime = f.derivative()._fast_float_(x) scipy.optimize.newton(g, 1, gprime) /// 1.4142135623730951 }}}

      Interestingly, newton in scipy is written in PURE PYTHON, so might (!?) be easy to beat using Cython.

      {{{id=79| scipy.optimize.newton?? /// }}}

      Write our own Newton method in Cython:

      {{{id=23| %cython def newton_cython(f, double c, fprime, double eps, int maxiter=100): cdef double fc cdef int i for i from 0 <= i < maxiter: fc = f(c) absfc = -fc if fc < 0 else fc if absfc < eps: return c c = c - fc/fprime(c) return c /// }}}

      Test them all

      {{{id=80| x = var('x') f = (x^2 - cos(x) + sin(x^2))^3 g = f._fast_float_(x) gprime = f.derivative()._fast_float_(x) timeit('scipy.optimize.newton(g, 0.5r, gprime)') timeit('newton_cython(g, 0.5r, gprime, 1e-13)') /// 625 loops, best of 3: 127 µs per loop 625 loops, best of 3: 37.6 µs per loop }}} {{{id=82| newton_cython(g, 0.5r, gprime, 1e-13) /// 0.63810381041535857 }}} {{{id=6| plot(g, 0,1) /// }}} {{{id=81| /// }}}

      Easy wrappers

      Sage has code that uses Scipy, fast_float, etc., behind the scenes.  With your help it could have a lot more!

      {{{id=86| var('x') f = cos(x)-x /// }}} {{{id=85| f.find_root(0,1) # fast /// 0.7390851332151559 }}} {{{id=109| /// }}}

      Problem: Find a root of something else on another interval.

      {{{id=83| /// }}}

      Problem: Type f.find_root??  to read the source code.  The function just calls into another module.  Type sage.numerical.optimize.find_root?? to read the source code of that, which does involve calling scipy.

      {{{id=114| /// }}} {{{id=16| /// }}} {{{id=115| /// }}} {{{id=75| /// }}} {{{id=125| /// }}} {{{id=124| /// }}} {{{id=123| /// }}} {{{id=122| /// }}} {{{id=77| /// }}}