# Something faster than Do loops with AppendTo

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.gatessucks 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
@WalkingRandomly No problem of course - well done. – b.gatessucks Feb 26 '14 at 10:35
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

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 = InternalBag[];
n = Length@v;
Do[InternalStuffBag[result, (v[[i]] - v[[j]])/(i - j)], {i, n,
2, -1}, {j, i - 1, 1, -1}];
InternalBagPart[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

-
it seems the same to me. – lalmei Feb 26 '14 at 12:18
@lalmei see the edit :) – Jacob Akkerboom Feb 26 '14 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. – ciao Feb 26 '14 at 12:28
@JacobAkkerboom thanks for the inclusion :) though I wasn't going for speed. Nice answer btw! What is InternalStuffBag ? 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

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

-
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

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.

-

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

-

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]]]];

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