Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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?

share|improve this question
7  
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 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.gatessucks Feb 26 at 10:27
    
You posted this just as I was writing mine. I didn't mean to rip you off! –  WalkingRandomly Feb 26 at 10:30
    
@WalkingRandomly No problem of course - well done. –  b.gatessucks Feb 26 at 10:35
1  
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 at 12:24

5 Answers 5

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
share|improve this answer
4  
You can also replace v1[[i]]-v1[[j]] with i-j. There is no need for v1. –  Ymareth Feb 26 at 10:31
    
Good point. Corrected. –  WalkingRandomly Feb 26 at 16:41

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.

share|improve this answer

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]]]];
share|improve this answer
    
+1, nice work. See my comment re: later subtraction. –  rasher Feb 26 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 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 at 17:33

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

share|improve this answer

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

share|improve this answer
    
it seems the same to me. –  lalmei Feb 26 at 12:18
    
@lalmei see the edit :) –  Jacob Akkerboom Feb 26 at 12:20
    
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. –  rasher Feb 26 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 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 at 16:14

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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