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
Closely related: professionalize/optimize my code that calls a function — Jacob Akkerboom, Feb 26 '14 at 12:13
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
@Ymareth correctly points out that v1 is not useful, as v1[[i]] == i for all i. I would suggest simply correcting this flaw in the question to make it as clean as possible. — Jacob Akkerboom, Feb 26 '14 at 12:55
I concur with @JacobAkkerboom; having a canonical pitfalls of AppendTo and some suggested alternatives would be helpful. The Sow/Reap suggestion is not intuitive (IMHO). — bobthechemist, Feb 26 '14 at 13:48
I think I will make a new Q&A, so that we can make that one canonical. This will remove distractions like the existence of v1 and comments. I hope we will close this question as a duplicate of the new one even if this sounds strange. — Jacob Akkerboom, Feb 26 '14 at 18:33

Simon Woods · Answer 1 · 2014-02-26T11:07:16.227

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

84,945
8
175
324

+1, nice work. See my comment re: later subtraction. – ciao Feb 26 '14 at 12:29
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 :(. – Jacob Akkerboom Feb 26 '14 at 15:28
@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. – Simon Woods Feb 26 '14 at 17:33

Jacob Akkerboom · Answer 2 · 2014-02-26T18:27:09.120

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 · Answer 3 · 2014-02-26T12:26:54.150

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 · Answer 4 · 2014-02-26T16:41:24.373

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

4,107
1
20
36

4

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

score 3 · Answer 5 · answered Feb 26 '14 at 11:48

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}}}

Something faster than Do loops with AppendTo

5 Answers5

Linked