Wolfram Language: Tips for Writing Fast Code

Wolfram Research
training@wolfram.com

Use Floating-Point Numbers (and Use Them Early)

N[Det[Table[1/(1+Abs[i-j]),{i,1,150},{j,1,150}]]]//AbsoluteTiming

N[Det[Table[1/(1.+Abs[i-j]),{i,1.,150.},{j,1.,150.}]]]//AbsoluteTiming

Nest[(Sqrt[1+#])&,2,1000]//N//AbsoluteTiming

Nest[(Sqrt[1+#])&,2.,1000]//N//AbsoluteTiming

Use the Special Data Types

PackedArray

data=RandomReal[1,{1000000}];
ByteCount[data]

ByteCount[Append[data,1]]

Developer`PackedArrayQ[data]

Developer`PackedArrayQ[Append[data,1]]

data2=Developer`ToPackedArray[Append[data,1],Real];

Developer`PackedArrayQ[data2]

ByteCount[data2]

SparseArray

data=ConstantArray[0,{1000,1000}];
Do[data[[RandomInteger[{1,1000}],RandomInteger[{1,1000}]]]=RandomReal[100],{400}];

sparsedata=SparseArray[{},{1000,1000}];
Do[sparsedata[[RandomInteger[{1,1000}],RandomInteger[{1,1000}]]]=RandomReal[100],{400}];

ByteCount[data]

ByteCount[sparsedata]

AbsoluteTiming[data.data;]

AbsoluteTiming[sparsedata.sparsedata;]

AbsoluteTiming[sparsedata.data;]

SparseArray[data]

StructuredArray

Compile Low-Level Code

data=Range[ -50.,50, 0.25];

fn=Function[{x},
Block[{sum=1.0,inc=1.0},Do[inc=inc*x/i;sum=sum+inc,{i,10000}];sum]];
Map[fn,data];//AbsoluteTiming

cfn=Compile[{x},
Block[{sum=1.0,inc=1.0},Do[inc=inc*x/i;sum=sum+inc,{i,10000}];sum]];
Map[cfn,data];//AbsoluteTiming

cfn2=Compile[{x},
Block[{sum=1.0,inc=1.0},Do[inc=inc*x/i;sum=sum+inc,{i,10000}];sum],
RuntimeAttributes->{Listable},Parallelization->True];
cfn2[data];//AbsoluteTiming

cfn2C=Compile[{x},
Block[{sum=1.0,inc=1.0},Do[inc=inc*x/i;sum=sum+inc,{i,10000}];sum],
RuntimeAttributes->{Listable},Parallelization->True,CompilationTarget->"C"];
cfn2C[data];//AbsoluteTiming

Use Built-in Functions

data=RandomReal[1,{1000000,2,2}];

Map[Flatten,data];//AbsoluteTiming

Flatten[data,{{1},{2,3}}];//AbsoluteTiming

Block[{temp=0},
Do[temp=temp+1/i^2,{i,1.,100000}];
temp];//AbsoluteTiming

NSum[1/i^2,{i,1.,100000}];//AbsoluteTiming

Understand Your Code with the Profiler and Debugger

fib[n_]:=fib[n-1]+fib[n-2];
fib[1]=1;
fib[2]=1;
fib[35];//AbsoluteTiming

Fibonacci[35];//AbsoluteTiming

Do More Than One Thing at a Time

Automatic for many numerical operations

data=RandomReal[1,{1000,1000}];

While[True,Det[data]]

But mostly not automatic for everything else

While[True,Expand[(1+x)^1000]]

Parallelize

Parallelize your own code

Table[PrimeQ[x],{x,10^1000,10^1000+5000}];//AbsoluteTiming

LaunchKernels[]

ParallelTable[PrimeQ[x],{x,10^1000,10^1000+5000}];//AbsoluteTiming

...but it is not always a good thing

data=RandomReal[1,{500,100,100}];

Map[Transpose,data];//AbsoluteTiming

ParallelMap[Transpose,data];//AbsoluteTiming

Guidelines for Parallelizing and Using CUDA

Parallelize when:

The question is small

The answer is small

Computation is expensive

The questions are independent

Or use the GPU when:

The questions are completely independent

The question is simple

Question or answer data can be copied over just once

Try to Work Things Out Only Once

Remember values for reuse

f[x_]:=f[x]=(*What the function does*)

fib[n_]:=fib[n]=fib[n-1]+fib[n-2];

fib[35]//AbsoluteTiming

Try to Work Things Out Only Once

Watch out for functions with nonstandard evaluation

Attributes[Plot]

Plot[Integrate[1/(1+x)^17,{x,0,n}],{n,1,30}]//AbsoluteTiming

Plot[Evaluate[Integrate[1/(1+x)^17,{x,0,n}]],{n,1,30}]//AbsoluteTiming

With[{arg=Integrate[1/(1+x)^17,{x,0,n}]},
Plot[arg,{n,1,30}]]//AbsoluteTiming

data=Table[x^2+5000RandomReal[],{x,1,100}];
fit=LinearModelFit[data,{x,x^2},x]

Show[ListPlot[data],
Plot[fit["MeanPredictionBands",ConfidenceLevel->0.9999],{x,0,100}]]//AbsoluteTiming

Show[ListPlot[data],
Plot[Evaluate[fit["MeanPredictionBands",ConfidenceLevel->0.9999]],{x,0,100}]]//AbsoluteTiming

f[x_]:=Nest[Sin[#+Sin[2#]]&,x,20]

FindRoot[f[x],{x,0}]//Timing

g[x_?NumericQ]:=Nest[Sin[#+Sin[2#]]&,x,20]

FindRoot[g[x],{x,0}]//Timing

Try to Work Things Out Only Once

Use Dynamic when refreshing small parts of graphics

ptdata=RandomReal[1,{10000,2}];

Manipulate[Graphics[
{Point/@ptdata,PointSize[0.05],Red,Point[{0.5,t}]}],
{t,0.05,0.95},SaveDefinitions->True]

Manipulate[Graphics[
{Point/@ptdata,PointSize[0.05],Red,Dynamic[Point[{0.5,t}]]}],
{t,0.05,0.95},SaveDefinitions->True]

Use Sow and Reap to Accumulate Data (Not AppendTo)

data={};Do[AppendTo[data,RandomReal[x]],{x,0,40000}];//AbsoluteTiming

data=Reap[Do[Sow[RandomReal[x]],{x,0,40000}]][[2]];//AbsoluteTiming

Go Easy on Pattern Matching

data=RandomReal[1,{200}];
data//.{a___,b_,c_,d___}/;b>c->{a,c,b,d};//AbsoluteTiming

(flag=True;
While[TrueQ[flag],
flag=False;
Do[If[
data[[i]]>data[[i+1]],
temp=data[[i]];
data[[i]]=data[[i+1]];
data[[i+1]]=temp;
flag=True],{i,1,Length[data]-1}]];
data);//AbsoluteTiming

Sort[RandomReal[1,{200}]];//AbsoluteTiming

Use Block or With rather than Module

Do[Module[{x=2.},1/x],{1000000}];//AbsoluteTiming

Do[Block[{x=2.},1/x],{1000000}];//AbsoluteTiming

Do[With[{x=2.},1/x],{1000000}];//AbsoluteTiming

Try Doing Things Differently!

temp=1;Do[temp=temp i,{i,2^16}];//AbsoluteTiming

Apply[Times,Range[2^16]];//AbsoluteTiming

AbsoluteTiming[65536!;]

Partition[Flatten[data],4];//AbsoluteTiming

Summary

Do not be more accurate than you need to be (use reals)

Use (the fewest) built-in commands

Do things once and in parallel

Compile

Experiment