Generally efficient numerical programming in Mathematica

Mathematica allows you to write simple programs which can work for a very general  class of inputs.  For numerical programs, a few simple principles allows you  to write programs which are efficient for both machine numbers and bignums.   In this talk, I will briefly describe the types of numbers Mathematica works with, then go on to show you two examples which demonstrate what I  consider to be the essential principles of efficient numerical programming in  Mathematica.

You should be able to find the notebook and packages for this talk at http://members.wri.com/rknapp/DeveloperConference99

Numbers in Mathematica

One of the important features of Mathematica is that it works with a variety of types of numbers, using a consistent  set of semantics to combine different types when necessary.

Exact Quantities

Pi  <  E

100!

NumericQ[Sin[E]]

Precision[Sin[E]]

Semantic note:  Combinations of exact expression with exact expressions 
generally comes out exact.  When combining with inexact expressions, the 
final precision is typically determined by the precisions of the inexact 
parts.

Machine real numbers

In Mathematica, when you use machine numbers, you typically get whatever your machine's  hardware and library functions give you.  For example:

Sin[2. ]

Computations with machine numbers can be done very quickly.  In Mathematica 4, PackedArrays give even better speed.  For example:

data  =  Table[Random[], {1000000}];
Timing[Sin[data];]

Semantic note:  Even though machine numbers technically have precision equal 
to $MachinePrecision, this precision is not guaranteed in any way. Thus 
combinations of machine numbers with other types of numbers attempts coercion 
to machine numbers.

Bignums

These are approximate real numbers too large to fit into machine hardware 
real numbers (on most machines this is IEEE double precision with a 53 bit 
(approx 16 digit mantissa) and an 8 bit exponent.)

This has too large an exponent, so it is given by a bignum:

E^1000.

The arithmetic for bignums is implemented in software.  Because of this, it 
is possible to keep track of the number of significant bits through 
significance arithmetic.  If you are interested in details on this, I suggest 
looking at some excellent notebooks written by Mark Sofroniou on the subject.

Here is an example where keeping track of precision makes a big difference:

arg  =  2^60  -  21

Even though the argument is large, the machine library sine funciton returns 
a value:

Sin[N[arg]]

However, this is identical to the value for sine of 2^60, which is clearly 
wrong.

%  ===  Sin[2. ^60]

If we use precision > $MachinePrecision, then the computation keeps track of 
precision:

Sin[N[arg, 17]]

The zero means we have no precision in the result!

Precision[%]

For precisions > $MachinePrecision, the command N is adaptive, so

N[Sin[arg], 17]

internally uses higher precision to get a result which is good to 17 
digits.

The disadvantage of using bignums is that software arithmetic is bound to be 
slower than hardware arithmetic.

data  =  
        SetPrecision[Take[data, 1000], 17];
Timing[Sin[data];]

In Mathematica 4, a lot of effort was put into making sure that with very high  precision, operations are as efficient as possible.

In my opinion, the improvements in Mathematica 4 are shown off well by the following question which I was asked by a  representative of the Japanese Mathematica users group in late 1995:

"Can Mathematica compute Pi to 5,000,000 digits in principle?  In practice?"

At the time the question was asked, the answer was basically: in prinicple, 
yes, but in practice, it could take days or weeks.  Though I would like to do 
it before your very eyes, it still takes a little long.  Nonetheless, the 
computation is practical to do in the time of a coffee break:

$MaxPrecision  =  Infinity;
Timing[pi50k  =  N[Pi, 5000000];]

Yes, part of the speed improvement has to do with the nice 333MhZ  processor I have on this machine, but that is actually a small part of the  story.  Mathematica 3.0 would take days to weeks on the same computation on this machine.

The reason I consider this a good example is that there are two distinct 
improvements which went into this functionality.

First, multiplication is done with asymptotically efficient algorithms.   This change effects any computation you do in Mathematica with more than a few hundred digits. You will see an example of how this  algorithm works in a few minutes.

Second, the algorithm for computing π was improved to take advantage  of the faster multiplication.  This has been done for many fundamental  functions in Mathematica, such as factorials, exponentials, trig functions, and the like.

The end result is that, using the same syntax and commands as in Mathematica 3.0, you may get results an order of magnitude faster with Mathematica 4.

Some programming principles

These are general principles.  They apply to all types of numbers, including 
those stored in packed arrays, where using these principles can make an even 
bigger difference.

Try to avoid using individual parts of lists

To help see the differences, I will use an alternative way of computing  the time it takes an expression to evaluate

time[Table[Sin[data[[i]]], {i, Length[data]}];]

data  =  N[data];

time[Table[Sin[data[[i]]], {i, Length[data]}];]

time[Sin[data];]

f[x_]  :=  Sin[x];
time[Table[f[data[[i]]], {i, Length[data]}]]

time[f[data]]

Use functional operations which operate on large chunks of your data.

This is really the same as the one above, but it is so important, I give it 
twice!

Try to avoid copying

Just the other day I got a question through tech. support about how to 
replace submatrices efficiently.  They had been using ReplacePart.  Here is 
an example:

matrix  =  Table[Random[], {100}, {100}];

parts  =  
      Transpose[{Range[10, 90], 50  +  0 Range[10, 90]}]

time[matrix  =  ReplacePart[matrix, 0. ,  parts]]

Extract[matrix,  parts]

time[Extract[matrix,  parts]]

Here is a faster way:

time[matrix[[Range[10, 90],  50]]  =  0. ]

Part is typically faster than Extract for submatrices

time[matrix[[Range[10, 90], 50]]]

Don't use Append in a loop

..or Prepend or AppendTo or PrependTo or Join.

This is actually the previous principle repeated in a specific instance.

This example is a bit contrived, but it actually closely relates to solve 
differential equations with automatic time step control.  (Subtle 
advertisment for my next talk later today.)

randomlist[x_]  :=  
        Module[{a  =  {},  sum  =  0}, 
While[
            sum  <  x, 
r  =  
              Random[]; 
AppendTo[a,  
              r]; 
sum  +=  
              r]; 
a];

randomlist[5]

That seemed pretty fast, but when we compute over a wider range of 
values.  Here I use the precomputed data to save time

timetable  =  (Hold[
            Table[
SeedRandom[
                1234]; 
{n  =  2^k,  
                time[randomlist[n]]}, 
{k, 12}]]; {{2, 
              0.00044921875000000083`}, {4, 0.000722656249999995`}, {8, 
              0.0012500000000000011`}, {16, 0.0023437500000000056`}, {32, 
              0.004687500000000011`}, {64, 0.009531250000000102`}, {128, 
              0.020312499999999956`}, {256, 0.047812500000000036`}, {512, 
              0.13374999999999915`}, {1024, 0.427500000000002`}, {2048, 
              1.5900000000000034`}, {4096, 
              5.859999999999999`}});

... it eventually slows down a lot

ListPlot[timetable,  PlotJoined -> True,  
        PlotRange -> All];

Fit[timetable,  {1, x, x^2}, x]

Here is a better way to do it.

betterrandomlist[x_]  :=  
        Module[{sum  =  0}, 
NestWhileList[Random[] &,  
            Random[],  ((sum  +=  #)  <  x) &]];

betterrandomlist[5]

Now check timings for this command (again precomputed so you don't ahve to wait)

timetable2  =  (Hold[
            Table[
SeedRandom[
                1234]; 
{n  =  2^k,  
                time[betterrandomlist[n]]}, 
{k, 
                12}]]; {{2, 0.0004736328125000029`}, {4, 
              0.0008007812500000003`}, {8, 0.0014453125000000039`}, {16, 
              0.0028124999999999956`}, {32, 0.005351562500000018`}, {64, 
              0.011250000000000093`}, {128, 0.021093750000000355`}, {256, 
              0.04187500000000011`}, {512, 0.08375000000000021`}, {1024, 
              0.16750000000000043`}, {2048, 0.332499999999996`}, {4096, 
              0.6749999999999972`}});

ListPlot[timetable2,  PlotJoined -> True,  
        PlotRange -> All];

Fit[timetable2,  {1, x, x^2}, x]

Compare the size of the quadratic coefficient.

Think carefully about how many things will have to be evaluated.

This is subtle, but important.  Here is an example which illustrates the idea

n  =  10000;

Timing[
      For[sum  =  0; i  =  1,  i  ≤  n,  i++,  
        sum  +=  i];  sum]

Timing[sum  =  0;  Do[sum  +=  i, {i, n}];  sum]

The example below is a nice trick for packed arrays...

time[Apply[Plus,  Range[n]]]

time[Tr[Range[n]]]

When you can't possibly avoid doing one of the first two, try to use the 
compiler.

Here is a cute example which shows how asymptotically faster multiplication 
works at the same time.

a  =  
        1839081381093812098309830918381238183908381038138;
 

b  =  
        8237489274892747624237592768267459874897489724987284972897489274728497  289478927489274982734987238497;

To simplify the discussion, I will use radix 10, though the code written 
below will work for any positive integer radix.

r  = 10;

A  =  Reverse[IntegerDigits[a, r]];
B  =  Reverse[IntegerDigits[b, r]];

len  =  2*Max[Length[A],  Length[B]]

Here is a way to multiply the two numbers by using the discrete convolution 
theorem on the list of digits.

AB  =  
      Round[Take[
          InverseFourier[
            Fourier[PadRight[A, len]] Fourier[PadRight[B, len]],  
            FourierParameters -> {-1, 1}], 
          Length[A]  +  Length[B]]]

...but there is a problem in that you need to deal with the carries.

PropagateCarries[x_,   r_]  :=  
      Module[{car  =  0,  res  =  x,  
          t}, 
Do[
t  =  
            res[[i]]  +  car; 
car  =  
            Quotient[t,  r]; 
res[[i]]   =  
            t  -  r*car, 
{i, 
            Length[x]}]; 
If[Last[res]  ==  0,  
          Drop[res,  -1], res]]

PropagateCarries[AB,  r]

Now check the result...

FromDigits[Reverse[%], 10]  -  a*b

Mathematica 4 provides a more convenient (and generally faster way) to do this via  the discrete convolution command ListConvolve

ListConvolve[A,  B,  {1, Length[A] + 1},  0]

using this, make a function to do the job

ListMultiply[r_][x_,  y_]  :=  
 
      PropagateCarries[ListConvolve[x,  y,  {1, Length[x]  +  1},  0],  
        r]

ListMultiply[r][A, B]

or with another radix:

ListMultiply[27][Reverse[IntegerDigits[a, 27]], 
      Reverse[IntegerDigits[b, 27]]]

BaseForm[a*b, 27]

a*b  -  FromDigits[Reverse[%%],  27]

Now lets try this with a bigger number

a  =  
        Random[Integer,  {10^99999,  
            10^100000}];
b  =   Random[
          Integer,  {10^99999,  10^100000}];
Timing[A  =  Reverse[IntegerDigits[a, r]]; 
      B  =  Reverse[IntegerDigits[b, r]];]

Timing[ListMultiply[r][A, B];]

Timing[AB  =  
          ListConvolve[A, B, {1, Length[A]  +  1}, 0];]

Here is an alternate way to do the carries:

ListPropagateCarries[x_,  
        r_]  :=  
Module[{car}, 
FixedPoint[
          Function[
            car  =  Quotient[#, r];  
#  -  
              r*car  +  RotateRight[car]],  x,  10]]

Timing[ListPropagateCarries[AB,  r];]

If we can bear to restrict ourselves to reaonably small radix, then we can  use the Mathematica  compiler.

CompiledPropagateCarries  =  
        Compile[{{x,  _Integer,  1},  {r,  _Integer}}, 
          Module[{car  =  0,  res  =  x,  
              t  =  0}, 
Do[
t  =  
                res[[i]]  +  car; 
car  =  
                Quotient[t,  r]; 
res[[i]]   =  
                t  -  r*car, 
{i, 
                Length[x]}]; 
If[
              Last[res]  ==  0,  Drop[res,  -1], 
              res]]];

Timing[CompiledPropagateCarries[AB,  r];]

Clear[PropagateCarries];
PropagateCarries  =  CompiledPropagateCarries;

Even with this work, it is still slower than the internal Mathematica Times:

Timing[ListMultiply[r][A, B];]

Timing[a*b;]

which is using basically the same algorithm.  If we turn it off, you get  the Mathematica 3.0 speed.

Developer`SetSystemOptions["TimesCrossoverDigits" -> 
          Infinity];

Timing[a*b;]

More involved Examples

There are many examples I could show you.  I have picked two particular  ones.

The first example came from a user I helped to speed up his program  to compute solutions of the Schrodinger equation using the Crank-Nicholson method.  Analysis of the program led to a huge speed increase.

The second example is from the demonstrative package  LinearAlgebra`GaussianElimination.  With some significant modifications, the core routine of this package can be made nearly as fast as the built  in code for bignums.