2
$\begingroup$

I am observing something which I do not understand with the following function function:

kelem={{n+0.,-n},{-n,n}};
assembleMatrice[diagK_]:=Block[{matrice},
  matrice=Sum[ArrayPad[kelem,{(i-1),(n-i)}],{i,1,n}];
  matrice[[2;;-1,2;;-1]]]
function[n_]:=Block[{},
  klist = ConstantArray[n+0.,n];
  mlist = ConstantArray[1./n,n];
  k = assembleMatrice[klist];
  m = DiagonalMatrix[mlist];
  b = Inverse[MatrixPower[m,.5]].MatrixPower[k,.5];
  q = Transpose@Eigenvectors[b.Transpose@b];
  xp:=Sum[q[[n,i]],{i,1,n}];]

When I evaluate AbsoluteTiming[function[n]] for $n\in\{1,\dots,249\}$, the time returned evolves smoothly. But why does it takes 50 times longer when $n$ reaches $250$?

Illustration:

Table[Block[{}, function[n]; AbsoluteTiming[xp;]], {n, 248, 252}][[All, 1]]
(* output: {0.000205, 0.000229, 0.014314, 0.014379, 0.014414} *)

You see that the difference between the third and second term, corresponding to $n$ reaching $250$, is much higher than the other differences.

Note: This is a MWE. xp is actually much more complicated than this, so I am not asking for better ways to calculate this xp. I am intersted in optimised definitions for the other variables though, if any.


Edit I just stumbled upon an interesting (for me...) observation.

Consider the following (stupid) function:

x1 = Function[{n, T, t}, (Table[Cos[(Mod[t, T] - T/2)]/Sin[T/2.], {j, 1, n}])[[1]]];
times249 = Table[x1[249, 123, 3] // AbsoluteTiming, {i, 1, 100}][[All, 1]];
times250 = Table[x1[250, 123, 3] // AbsoluteTiming, {i, 1, 100}][[All, 1]];
ListPlot[{times249, times250}, PlotRange -> {{0, 100}, {0, 0.005}}]

Output:

enter image description here

It appears that this time, the function x1 takes way more time for values below 250. This can corrected using

SetSystemOptions["CompileOptions" -> {"TableCompileLength" -> 1}]

Note the the value has to be taken lower than 250, contrary to the previous example. The conclusion is, do not increase the values of CompileOptions without thinking, or at least trying.

$\endgroup$
  • 1
    $\begingroup$ (55242) $\endgroup$ – Sektor Sep 30 '15 at 20:43
  • $\begingroup$ Try and see if the solution actually works -- if it does it will just get flagged as a duplicate $\endgroup$ – Sektor Sep 30 '15 at 20:46
  • $\begingroup$ Both answers give good result, especially changing SumCompileLength. Using Total@Table instead of Sum results in a much faster computation than the OP but still shows a threshold at 250. $\endgroup$ – anderstood Sep 30 '15 at 20:55
  • $\begingroup$ Glad it works good :) $\endgroup$ – Sektor Sep 30 '15 at 20:56
  • 1
    $\begingroup$ @Sektor I asked the author of the accepted answer if there is any drawback of (reasonably) increasing SumCompileLength. I am mentionning it here to get a wider audience. $\endgroup$ – anderstood Sep 30 '15 at 20:58