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
Do
. It'sAppendTo
, which copies its argument on every call. UseSow
/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. $\endgroup$ – Oleksandr R. Feb 26 '14 at 10:19b = Flatten[ Table[Table[(v2[[i]] - v2[[j]])/(v1[[i]] - v1[[j]]), {j, i - 1, 1, -1}], {i, n, 2, -1}]]; // Timing
. $\endgroup$ – b.gates.you.know.what Feb 26 '14 at 10:27