Why Interpolation is slower than BSplineFunction?

Question

This my question is concerned with previous about Tube cross-section. Let's use the same pts as in that question.

pts = {{0.`, 0.`, 0.`}, {7.497665500157259`*^-14, 
   0.00036555534709634656`, 
   0.041887873232098444`}, {2.9988378071419946`*^-13, 
   0.0014621100333863195`, 
   0.08376298663421493`}, {6.746528662189823`*^-13, 
   0.0032893300277940163`, 
   0.12561258426325045`}, {1.1991697498398906`*^-12, 
   0.005846658724918654`, 
   0.16742391794868797`}, {1.8732746539195436`*^-12, 
   0.009133317114586991`, 
   0.20918425117592362`}, {2.696762233490476`*^-12, 
   0.013148304019154402`, 
   0.25088086296604745`}, {3.669381638827055`*^-12, 
   0.017890396398482316`, 
   0.2925010517508915`}, {4.790836591859004`*^-12, 
   0.02335814972249914`, 
   0.33403213924216535`}, {6.0607854764232425`*^-12, 
   0.029549898411231134`, 
   0.37546147429349974`}, {7.478841442326563`*^-12, 
   0.03646375634216925`, 
   0.41677643675422266`}, {9.044572523187444`*^-12, 
   0.04409761742481735`, 
   0.45796444131369324`}, {1.0757501768021084`*^-11, 
   0.052449156242246726`, 
   0.499012941335023`}, {1.26171073865276`*^-11, 0.06151582875946169`,
    0.5399094326770161`}, {1.4622822908039118`*^-11, 
   0.07129487309836018`, 
   0.5806414575031648`}, {1.6774037354077335`*^-11, 
   0.08178331037905351`, 0.6211966080765395`}, {1.9070095424469`*^-11,
    0.09297794562728895`, 
   0.6615625305394178`}, {2.1510297696962602`*^-11, 
   0.10487536874769866`, 
   0.7017269286765003`}, {2.4093900840285483`*^-11, 
   0.11747195556257856`, 
   0.7416775676605682`}, {2.6820117840576448`*^-11, 
   0.13076386891588063`, 
   0.7814022777794407`}, {2.9688118241124915`*^-11, 
   0.14474705984208247`, 0.8208889581430966`}}

Now I need to construct parametric form of the curve. One method is using Interpolation:

path[u_] = Interpolation[#, u] & /@ Transpose[pts]

The main advantage of Interpolation is that path is going through all the points of pts exactly. Then I use code from @SquareOne:

frenet[u_] = FrenetSerretSystem[path[u], u][[2]]; // AbsoluteTiming
transform[u_] := 
  Composition[TranslationTransform[pts[[u]]], 
   FindGeometricTransform[
     frenet[u], {{0, 0, 1}, {1, 0, 0}, {0, 1, 0}}][[2]]];

I included AbsoluteTiming here to measure time of operation. And I see that frenet[u_] takes significant time, about 0.27--0.3 sec on my machine.

Another approach is using BSplineFunction:

path1[u_] := BSplineFunction[#][u] & /@ (Transpose@pts);
frenet1[u_] = FrenetSerretSystem[path1[u], u][[2]]; // AbsoluteTiming

This is much faster, it takes about 0.09--0.1 sec on my machine, i.e. BSplineFunction is 2.5--3 times faster than Interpolation. Why is this difference in speed? Can I speed-up Interpolation if I want to use it because of (as I think, probably wrong) it gives more accurate approximation.

Thank you.

What happens if you use Method -> "Spline" in Interpolation[]? — J. M.'s missing motivation, Dec 19 '16 at 10:12
If I run path[u_] = Interpolation[#, u, Method -> "Spline"] & /@ Transpose[pts], then frenet[u_] takes about 0.32--0.35 sec, more or the same time or even longer. — Alx, Dec 19 '16 at 10:28
Just to be clear, you want to speed up the construction of the code generated by FrenetSerretSystem, which is done only once for each set of pts? Interpolation itself is fast, but it is FrenetSerretSystem that takes the extra time when analyzing (I surmise) InterpolatingFunction compared to analyzing BSplineFunction. — Michael E2, Dec 19 '16 at 18:33
@MichaelE2, right. I have many such sets of points, and I had a hope it is possibe to decrease time of these operations. — Alx, Dec 19 '16 at 23:21

Michael E2 · Accepted Answer · 2020-12-06T18:47:28.167

Perhaps this will help: FrenetSerretSystem seems to take a symbolic approach, which slows it down. A more straightforward numeric approach, which will be much faster, may be more satisfactory.

The following function works for all (twice differentiable) paths. The savings are two-fold. One doesn't have to compute a new Frenet frame function for each set of points; only an interpolating path needs to be computed. Second, as seen below, this function executes faster than the one created with FrenetSerretFrame[] in the OP.

ClearAll[frenet2];
fsfC = Compile[{{v, _Real, 1}, {a, _Real, 1}},
   Block[{nv, na},
    If[v == {0, 0, 0},                        (* singular point *)
     {{0, 0, 0}, {0, 0, 0}, {0, 0, 0}},
     nv = v * Internal`ReciprocalSqrt[v.v];   (* use v / Sqrt[v.v] pre V10 *)
     If[a == {0, 0, 0},                       (* no turning *)
      {nv, {0, 0, 0}, {0, 0, 0}},
      na = # * Internal`ReciprocalSqrt[#.#] &[a - (a.nv) nv];
      {nv, na, Cross[nv, na]}
      ]]
    ]];
frenet2[path_, {u_, u0_?NumericQ}] :=
  With[{v = D[path, u]}, With[{a = D[v, u]}, Block[{u = u0},
     fsfC[v, a]
     ]]];

Using a single InterpolatingFunction instead of three seems to speed things up a little, too:

path2[u_] = Interpolation[MapIndexed[{#2, #1} &, pts]][u]

Timings. Note that frenet2 works on the OP's paths. The frame is slightly different for path1 because the B-spline path is slightly different from the interpolated path.

uu0 = 12.5;
frenet2[path[u], {u, uu0}] // RepeatedTiming
frenet2[path1[u], {u, Rescale[uu0, {1, 21}]}] // RepeatedTiming
frenet2[path2[u], {u, uu0}] // RepeatedTiming
(*
  {0.00013,
    {{4.08916*10^-11, 0.199371, 0.979924},
     {2.00986*10^-10, 0.979924, -0.199371},
     {-1., 2.05103*10^-10, 1.9387*10^-26}}}
{0.00016,
    {{4.0365310^-11, 0.196805, 0.980443},
     {2.0109210^-10, 0.980443, -0.196805},
     {-1., 2.0510310^-10, 9.6935210^-26}}}
{0.000077,
    {{4.0891610^-11, 0.199371, 0.979924},
     {2.0098610^-10, 0.979924, -0.199371},
     {-1., 2.0510310^-10, 1.938710^-26}}}
*)

The first and third frenet2 frames agree with the frame computed by FrenetSerretFrame used in the OP's frenet; they're also much faster than frenet:

frenet[uu0] // RepeatedTiming
(*
  {0.0011,
    {{4.08916*10^-11, 0.199371, 0.979924},
     {2.00986*10^-10, 0.979924, -0.199371},
     {-1., 2.05103*10^-10, 1.28793*10^-26}}
*)

sorry for late question, had no time. As I undestand, in your fsfC parameters nv and na are normalized. So we can write nv = Normalize[v], and na = Normalize[a - (a.nv) nv]. This works. Next step is replace (a.nv) nv with Projection[a, nv], and this gives me errors that parameters have different rank and uncompiled function will be used. What is wrong? I read many times that one shoud preferably use built-in functions, they are extremely optimized. — Alx, Dec 22 '16 at 09:17
@Alx: Yes, built-in are often (should be always but in practice not) to be preferred. If we're using Compile, then we should use functions built into the compiler. Execute the following: Needs["CompiledFunctionTools`"]; then CompilePrint@Compile[{{v, _Real, 1}, {w, _Real, 1}}, Projection[v, w]] and CompilePrint@Compile[{{v, _Real, 1}}, Normalize[v]]. They show MainEvaluate[..] called on Projection and Normalize, which means the functions are not "compilable"..... — Michael E2, Dec 22 '16 at 14:28
@Alx ...Further Compile has no (internal) declared return type for Normalize[], so it assumes it returns a real number. You can declare the return type by added a declaration at the end: fsfC = Compile[{{v, _Real, 1}, {a, _Real, 1}}, Block[{nv, na},,,], {{Normalize[_], _Real, 1}}]. See "Compiled Function Operation" in the compiler manual. Just to be clear, the use of Projection and Normalize inside Compile should be slower. — Michael E2, Dec 22 '16 at 14:37
thank you for clarification. I've read topic on List of compilable functions and the main conclusion (as I understood) is that (in most cases) simple ("low-level") functions such as Sin or Dot can be used in Compile, but not "high-level" already optimized ones as e.g. NIntegrate or Projection. — Alx, Dec 22 '16 at 15:07

Why Interpolation is slower than BSplineFunction?

1 Answers1