4
$\begingroup$

The main question

I was trying to obtain the piecewise cubic-spline interpolating polynomials of a series of very long lists of (x, f(x)) pairs (O(100) of them with length ~350,000), i.e. to get InterpolatingPolynomial on tab[[1;;4]], tab[[2;;5]], ..., with tab[[i]]={xi,yi}. The xi for all those lists are the same, on the good side.

Given that I'm working on a cluster, I tried to use ParallelTable to accelerate the calculation, with the code

(* tabvr is a list of 100*O(350,000) of (x, f(x)) pairs. For toy data of tabvr, see below *)

totpoint=Length[tabvr[[1]]];
wp=100;

LaunchKernels[100];
SetSharedVariable[tabvr];
vrpp$poly=ParallelTable[SetPrecision[InterpolatingPolynomial[Take[tabvr[[kk]],{i-3,i}],x]//Expand,wp]
          ,{kk,1,100},{i,4,totpoint}];

I requested 120 cores (5 nodes of 24 cores) and 400GB of memory. However, it turned out to be not enough and my job was killed as it goes beyond the memory limit. Moreover, the ParallelTable didn't seem to run as fast as I expected, as it takes ~50s for a single-core machine to process a single list, while it took about 2 hours before my job was terminated.

So, did I do anything wrong in this parallel calculation, or is such situation something I should expect for a parallel calculation? If it's the first case, what should I do to make the code work?

Edit:

If anyone want to actually play with the data, you may find a slice of the data here (the vrpp86.dat therein). You may just copy it 100 times to make a "complete" tabvr.

A side question

I also realized that, though I requested 5*24 cores, when I check the processor number

Print[$ProcessorCount];

I got the answer as 24 instead of 120. Does this mean that MMA's parallelization cannot work across different nodes?


Update

As inspired by @Henrik Schumacher 's answer, I managed to shorten the running time of the code by a lot after using Compile

xgrid=tabvr[[1]];
ygrid=tabvr[[2]];
totpoint=Length[xgrid];
wp=100;

result0=SetPrecision[Table[CoefficientList[InterpolatingPolynomial[tabvr[[i-3;;i]],x],x],{i,4,totpoints}].{1,x,x^2,x^3}, wp]//AbsoluteTiming;
Print[result0[[1]]];

(*
to get the code working, mat in the following line should be replaced by the result of
D[CoefficientList[InterpolatingPolynomial[Array[{x[[#]], a[#]} &, 4], x], x], {Array[a, 4], 1}]//Simplify
*) 
calc$coeff2 = Compile[{{x, _Real, 1}, {y, _Real, 1}}, mat.y, 
CompilationTarget -> "C", RuntimeAttributes -> {Listable}, Parallelization -> True, RuntimeOptions -> "Speed"];

result2 = (Table[SetPrecision[calc$coeff2[tabxgrid[[i - 3 ;; i]], ygrid[[i - 3 ;; i]]], wp], {i, 4, totpoints}].{1, x, x^2, x^3}) // AbsoluteTiming;
Print[result2[[1]]];

And I get the result

(* result0 *)
19.95942
(* result1 *)
2.976688
$\endgroup$
3
  • $\begingroup$ Welcome to Mathematica.SE! I hope you will become a regular contributor. To get started, 1) take the introductory tour now, 2) when you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge, 3) remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign, and 4) give help too, by answering questions in your areas of expertise. $\endgroup$
    – bbgodfrey
    Aug 12, 2021 at 14:38
  • $\begingroup$ Without (toy) data and without complete definitions, this is impossible to answer. So I expect this thread being closed soon, unless you edit your post. $\endgroup$ Aug 12, 2021 at 15:01
  • $\begingroup$ @HenrikSchumacher Thanks for the reminder. I added in a link to the toy data and modified the code a little bit. Now I guess one could play with it. $\endgroup$
    – llxy7
    Aug 12, 2021 at 17:13

1 Answer 1

4
$\begingroup$

Just don't use a cluster. You can do the full job on a costumer computer within a couple of seconds, at least if 15 digits of precision are sufficient. (I really wonder why you ask for 100 digits of precision.)

This is what I did: I downloaded the file vrpp86.dat into my home directory:

tabvr = ToPack[Join[##, 2] & @@ ConstantArray[ToPack@Import["~/vrpp86.dat"], 50]];

n = 10000;
x =.;
result0 = Table[
     CoefficientList[
      InterpolatingPolynomial[tabvr[[kk, i - 3 ;; i]] , x], x],
     {kk, 1, n}, {i, 4, Dimensions[tabvr][[2]]}]; // AbsoluteTiming // First

38.6408

The combination of CoefficientList and InterpolatingPolynomial is actually a linear operator acting on 4-vectors. We can generate the according matrix as follows:

A = Developer`ToPackedArray@N@
    D[
     CoefficientList[InterpolatingPolynomial[Array[a, 4], x], x],
     {Array[a, 4], 1}
     ];

This allows us to avoid all symbolic calculations and leads to a first speed-up:

result1 = Table[
     A . tabvr[[kk, i - 3 ;; i]],
     {kk, 1, n}, {i, 4, Dimensions[tabvr][[2]]}]; // AbsoluteTiming // First

2.59673

Compiling this leads to a greater speed-up:

cf = Compile[{{tabvr, _Real, 2}, {A, _Real, 2}, {k, _Integer}},
   Block[{n, m},
    n = Dimensions[A][[2]] - 1;
    m = Dimensions[tabvr][[2]];
    Table[A . tabvr[[k, i - n ;; i]], {i, n + 1, m}]
    ],
   CompilationTarget -> "C",
   RuntimeAttributes -> {Listable},
   Parallelization -> True,
   RuntimeOptions -> "Speed"
   ];

result2 = cf[tabvr, A, Range[n]]; // AbsoluteTiming // First

0.153785

And with a bit more effort, we can squeeze even more optimizations out of the compiler:

cf2 = Block[{k, i, X, XX},
   XX = Table[ Compile`GetElement[X, k, l], {l, {i - 3, i - 2, i - 1, i}}];
   With[{code = A . XX},
    Compile[{{X, _Real, 2}, {k, _Integer}},
     Table[code, {i, 4, Dimensions[X][[2]]}],
     CompilationTarget -> "C",
     RuntimeAttributes -> {Listable},
     Parallelization -> True,
     RuntimeOptions -> "Speed"
     ]
    ]
   ];
result3 = cf2[tabvr, Range[n]]; // AbsoluteTiming // First

0.028772

All variants have 15 digits of precision:

Max[Abs[1 - result0/result1]]
Max[Abs[1 - result0/result2]]
Max[Abs[1 - result0/result3]]

4.44089*10^-16

4.44089*10^-16

4.44089*10^-16

The compiled code is automatically parallelized. (This won't work on distributed systems.)

Timings were taken on my Quad Core Haswell CPU.

The full job is completed in less than a second:

fullresult = cf2[tabvr, Range[Length[tabvr]]]; // AbsoluteTiming // First

0.861428

$\endgroup$
3
  • 1
    $\begingroup$ Thanks a lot for the answer. It's very helpful, but not yet exactly what I want, mainly due to that I didn't describe the data clear enough (question updated). The data I provided is not a one dimensional list but a list of (x, f(x)) pairs. In that case, is there a way to redefine the A matrix? I tried to promote it to a (N-3)*4*4 packed array with the x-grid, but then I got a warning saying there's not enough memory for the calculation. $\endgroup$
    – llxy7
    Aug 13, 2021 at 12:40
  • $\begingroup$ It's very likely that you can compile also that. (I do not understand what you mean by "promote". Do not make any superfluous copies is all I can say.) $\endgroup$ Aug 14, 2021 at 6:39
  • $\begingroup$ By saying "promote" I mean to increase the dimension of A. The current A matrix is 4-by-4 and generates the coefficient list of an interpolating polynomial on{{1, A[1]}, {2, A[2]}, {3, A[3]}, {4, A[4]}}. But to make it work on arbitrary {{x[i], y[i]}, {x[i+1], y[i+1]}, {x[i+2], y[i+2]}, {x[i+3], y[i+3]}} I think I have to work out the 4-by-4 A for any index i, which effective changes A to an (N-3)*4*4 array with N the length of my table. $\endgroup$
    – llxy7
    Aug 15, 2021 at 7:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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