Something faster than Do loops with AppendTo

Question

I am an "old" programmer used to Fortran and Pascal. I can't get rid of For, Do and While loops, but I know Mathematica can do things much faster!

I am using the following code

SeedRandom[3]
n = 10;
v1 = Range[n];
v2 = RandomReal[250., n];

a = {};
Do[
   Do[
    AppendTo[a, (v2[[i]] - v2[[j]])/(v1[[i]] - v1[[j]])
     ],
    {j, i - 1, 1, -1}], {i, n, 2, -1}
   ]; // Timing

If n is small, it runs fast enough, but for bigger n it slows down. I usually deal with n > 600.

How can the code be made faster?

It isn't due to Do. It's AppendTo, which copies its argument on every call. Use Sow/Reap instead and you'll most likely find that your code runs fine as it is. Not to say that a functional alternative might not be faster and clearer, but if you find it easier to think procedurally, there are ways to improve performance without sacrificing that. — Oleksandr R., Feb 26 '14 at 10:19
Another suggestion is to just get the result as a jigged array and flatten it at the end : b = Flatten[ Table[Table[(v2[[i]] - v2[[j]])/(v1[[i]] - v1[[j]]), {j, i - 1, 1, -1}], {i, n, 2, -1}]]; // Timing. — b.gates.you.know.what, Feb 26 '14 at 10:27
You posted this just as I was writing mine. I didn't mean to rip you off! — WalkingRandomly, Feb 26 '14 at 10:30
Would anyone mind if we rename this Q&A: "The dangers of using Append(To) in a loop" and clean it up, to make it as recognisable and easy to find as possible? — Jacob Akkerboom, Feb 26 '14 at 12:24

Simon Woods · Accepted Answer · 2014-02-26 11:07:16Z

10

A more functional approach:

a = With[{f = Subtract @@@ Subsets[Reverse@#, {2}] &}, f[v2]/f[v1]]

For a bit more speed you could do this:

ii = Join @@ Table[ConstantArray[i, i - 1], {i, n, 2, -1}];
jj = Join @@ Table[Range[j, 1, -1], {j, n - 1, 1, -1}];
a = Divide[Subtract[v2[[ii]], v2[[jj]]], Subtract[v1[[ii]], v1[[jj]]]];

edited Feb 26 '14 at 11:07

answered Feb 26 '14 at 10:40

Simon Woods

79.5k8 gold badges159 silver badges299 bronze badges

$\begingroup$ +1, nice work. See my comment re: later subtraction. $\endgroup$ – ciao Feb 26 '14 at 12:29
$\begingroup$ At first I only looked at your second method, as you said this was faster. But now that I looked at the first method as well, it turns out that one is very slow :(. $\endgroup$ – Jacob Akkerboom Feb 26 '14 at 15:28
$\begingroup$ @rasher, I assumed that the v1 given in the question was an example and it would not generally be Range[n], but it's a good point worth mentioning. $\endgroup$ – Simon Woods Feb 26 '14 at 17:33

add a comment |

Jacob Akkerboom · Accepted Answer · 2014-02-26 18:27:09Z

As mentioned, the main issue is using AppendTo in loops like this. In this answer, I want to show that using Compile can make procedural code very fast. Below is a comparison of timings of all the answers, as well as the OPs code.

Here is a slight modification of the code by the OP. I have modified it because I wanted to focus on AppendTo, so I have cleaned it up a bit.

questionCode :=
  (
   qCRes = {};
   Do[AppendTo[qCRes, (v[[i]] - v[[j]])/(i - j)],
    {i, n, 2, -1}, {j, i - 1, 1, -1}]
   );

Alternatives

My code uses a compiled function, where code that is similar to that of the OP is compiled to C. If you do not have a C compiler, simply remove CompilationTarget->"C". My code also uses undocumented functions, notably Internal'Bag. This is basically an implementation of a linked list structure, which is especially useful inside Compile.

jacobCfu =
  Compile[
   {{v, _Real, 1}},

   Block[
    {result, n},
    result = Internal`Bag[];
    n = Length@v;
    Do[Internal`StuffBag[result, (v[[i]] - v[[j]])/(i - j)], {i, n, 
      2, -1}, {j, i - 1, 1, -1}];
    Internal`BagPart[result, All]], CompilationTarget -> "C"
   ];

jacobCode :=
  (
   jacobRes = jacobCfu[v]
   );

Other answerers' code. This includes a slightly modified version of Simon Woods (SW) code I made (sWCodeE)

rasherCode1 :=
  (
   rasherRes1 = 
    Table[(v2[[i]] - v2[[j]])/(v1[[i]] - v1[[j]]), {i, n, 2, -1}, {j, 
       i - 1, 1, -1}] // Flatten
   );

rasherCode2 :=
  (
   rasherRes2 =
    Block[
     {s, i1, i2},
     s = Subsets[Range[n, 1, -1], {2}];
     {i1, i2} = {s[[All, 1]], s[[All, 2]]}; 
     Divide[Subtract[v2[[i1]], v2[[i2]]], i1 - i2]
     ]
   );

wRCode :=
  (
   wRRes = 
    Flatten[Table[(v2[[i]] - v2[[j]])/(v1[[i]] - v1[[j]]), {i, n, 
       2, -1}, {j, i - 1, 1, -1}]]
   );

sWCode1 :=
  (sWRes1 = 
    With[{f = Subtract @@@ Subsets[Reverse@#, {2}] &}, f[v2]/f[v1]]);

sWCode2 :=
  (
   sWRes2 =
    Block[
     {ii, jj},

     ii = Join @@ Table[ConstantArray[i, i - 1], {i, n, 2, -1}];
     jj = Join @@ Table[Range[j, 1, -1], {j, n - 1, 1, -1}];
     Divide[Subtract[v2[[ii]], v2[[jj]]], Subtract[v1[[ii]], v1[[jj]]]]
     ]
   );

sWCodeE :=
  (
   sWResE =
    Block[
     {ii, jj},
     ii = Join @@ Table[ConstantArray[i, i - 1], {i, n, 2, -1}];
     jj = Join @@ Table[Range[j, 1, -1], {j, n - 1, 1, -1}];
     Divide[Subtract[v2[[ii]], v2[[jj]]], Subtract[ii, jj]]
     ]
   );

lalmeiCode :=
 (
  (lalmeiRes = 
    First@Last@
      Reap[Scan[
        Function[{x}, 
         Sow[(v2[[x[[1]]]] - v2[[x[[2]]]])/(v1[[x[[1]]]] - 
              v1[[x[[2]]]])];], 
        Table[{i, j}, {i, n, 2, -1}, {j, i - 1, 1, -1}], {2}]])
  )

Timing comparison functions

timing = Function[Null, First@Timing@#, HoldAll];
timingAndName = 
  Function[Null, {timing@#, ToString@Unevaluated@#}, HoldAll];
timingsAndNamesTable = 
  Function[Null, TableForm[timingAndName /@ Unevaluated[{##}]], 
   HoldAll];

formattedTimingsAndComparison =
  Function[Null,
   Block[{timingTable},
    timingTable = timingsAndNamesTable@##;

    Column[
     {
      StringForm["Comparison for n = ``", n]
      ,
      timingTable
      ,
      If[
       SameQ @@ resultNames
       ,
       "results are equal"
       ,
       Row[{"results are ", Style["not ", Bold], "equal"}]
       ]
      }
     ,
     Spacings -> 2 ]
    ]
   ,
   HoldAll
   ];

Initialisation

initialize[nn_] :=
 (
  SeedRandom[3];
  n = nn;
  v = RandomReal[250., n];
  v1 = Range[n];
  v2 = v;
  )

Timing comparison

initialize[200];
resultNames = 
  Hold[qCRes, jacobRes, sWResE, sWRes2, wRRes, rasherRes1, 
   rasherRes1, rasherRes2, lalmeiRes];

formattedTimingsAndComparison[
 questionCode,
 jacobCode,
 sWCodeE,
 sWCode2,
 rasherCode2,
 wRCode,
 rasherCode1,
 lalmeiCode
 ]

Gives

Comparison for n = 200

 1.168769   questionCode
 0.000674   jacobCode
 0.001093   sWCodeE
 0.001200   sWCode2
 0.004740   rasherCode2
 0.070625   wRCode
 0.069372   rasherCode1
 0.176819   lalmeiCode

results are equal

Let's look at a larger value of n as well

initialize[1000]
resultNames = Hold[jacobRes, sWResE, sWRes2, rasherRes2];

formattedTimingsAndComparison[
 jacobCode,
 sWCodeE,
 sWCode2,
 rasherCode2
 ]

Gives

Comparison for n = 1000
0.020356    jacobCode
0.052107    sWCodeE
0.110584    sWCode2
0.526762    rasherCode2
results are equal

C beats mathematica? Who would have thought it... Simon's fine idea can be further improved by removing the last subtraction, probably unbeatable using pure MMA and non-compiled, documented functions. Well done, SW. — ciao, Feb 26 '14 at 12:28
@JacobAkkerboom thanks for the inclusion :) though I wasn't going for speed. Nice answer btw! What is Internal`StuffBag ? Is it an undocumented function ? — lalmei, Feb 26 '14 at 15:58
@lalmei no problemo! Yeah, it is an undocumented function, as are Internal'Bag and Internal'BagPart, or any function in the Context Internal for that matter. It is basically an implementation of a linkedlist data structure, which can (even) be used in Compile. I quite like it :). — Jacob Akkerboom, Feb 26 '14 at 16:14

ciao · Accepted Answer · 2014-02-26 12:26:54Z

Just changing it to something like:

Table[(v2[[i]] - v2[[j]])/(v1[[i]] - v1[[j]]), {i, n, 2, -1}, {j, i - 1, 1, -1}] // Flatten

Should net you a nice boost. Edit - oops, ninja'd

About twice as fast as any so far on large N:

s = Subsets[Range[n, 1, -1], {2}];
{i1, i2} = {s[[All, 1]], s[[All, 2]]};
result = Divide[Subtract[v2[[i1]], v2[[i2]]], i1 - i2];

Edit - Simon beat me to the faster way to create indexes. His can be improved by removing the whole latter subtraction, as in mine, netting his a 25% boost in my tests.

WalkingRandomly · Accepted Answer · 2014-02-26 16:41:24Z

4

Your version:

SeedRandom[3]
n = 250;
v1 = Range[n];
v2 = RandomReal[250., n];

a = {};
Do[
   Do[
    AppendTo[a,
     (v2[[i]] - v2[[j]])/(v1[[i]] - v1[[j]])],
    {j, i - 1, 1, -1}], {i, n, 2, -1}]; // AbsoluteTiming

on my machine this takes 1.85 seconds

This version

AbsoluteTiming[A = Flatten[
  Table[(v2[[i]] - v2[[j]])/(i-j), {i, n, 2, -1}, {j, 
    i - 1, 1, -1}]] 
    ]

Takes 0.112 seconds

They give the same result.

In[75]:= a == A

Out[75]= True

edited Feb 26 '14 at 16:41

answered Feb 26 '14 at 10:29

WalkingRandomly

3,88715 silver badges36 bronze badges

4

$\begingroup$ You can also replace v1[[i]]-v1[[j]] with i-j. There is no need for v1. $\endgroup$ – Ymareth Feb 26 '14 at 10:31
$\begingroup$ Good point. Corrected. $\endgroup$ – WalkingRandomly Feb 26 '14 at 16:41

add a comment |

lalmei · Accepted Answer · 2014-02-26 11:48:58Z

Here is an example using Reap and Sow.

It is a bit overkill in this case, but if you need to append only some of the values it should go much faster with Reap and Sow, or if you need to perform other functions while appending.

Last@Reap[
Scan[Function[{x}, 
     Sow[(v2[[x[[1]]]] - v2[[x[[2]]]])/(v1[[x[[1]]]] -v1[[x[[2]]]])]; 
     ],
     Table[{i, j}, {i, n, 2, -1}, {j, i - 1, 1, -1}], 
     {2}]]//Timing

{0.00062, {{124.857, 17.9221, 27.095, 23.4391, 36.1562, 31.8519, 19.882, 27.9697, 11.8103, -89.0129, -21.7861, -10.3669, 13.9809, 13.2509, 2.38619, 14.1286, -2.3205, 45.4407, 28.956, 48.3122, 38.8168, 20.666, 31.3189, 10.0641, 12.4714, 49.748, 36.6089, 14.4723, 28.4945, 4.16804, 87.0246, 48.6777, 15.1393, 32.5003, 2.50737, 10.3307, -20.8033, 14.3255, -18.6219, -51.9374, 16.3229, -28.2728, 84.5831, -16.4406, -117.464}}}

Stack Exchange Network

current community

your communities

more stack exchange communities

Something faster than Do loops with AppendTo

5 Answers 5

Your Answer

Not the answer you're looking for? Browse other questions tagged performance-tuning or ask your own question.

Linked

Hot Network Questions

Something faster than Do loops with AppendTo

5 Answers 5

Your Answer

Sign up or log in

Post as a guest

Not the answer you're looking for? Browse other questions tagged performance-tuning or ask your own question.

Linked

Related

Hot Network Questions