performance issue with evaluating after interpolation
oguzumut-cnrs opened this issue · 7 comments
Hello,
I have a very complicated/long 12-th order (analytical) polynomial with 3 variables f(x,y,z) that I use in an iterative process. I create a discrete data using it and I interpolate successfully with SPLINTER. However, the computational cost of evaluating bsplines.eval(x) is much higher than calculating f(x,y,z) itself. I wonder whether it is normal? I was expecting bsplines.eval(x) performs faster because it has a much simpler cubic form.
Thanks
Hello,
I would need more information to be able to predict the computational cost of the two functions; what do you mean by much higher computational cost? It would help if posted an example showing the difference in evaluation time.
A tricubic B-spline, in k=3 variables, all degrees equal to p=3, and with all knot vectors having size m, should require approximately O(k(p^2 + m) + 2*(p+1)^k) computations (take this bound with a grain of salt).
Remember that a great advantage of the B-spline, over the high-order polynomial function, is that it is computationally stable.
Hello,
OK sorry, I wasn't precise enough.
The analytical function is very long and can be found here: https://link.springer.com/article/10.1007/s00205-004-0311-z
see below the piece of code that I am using to calculate the function.
For one call, I obtain the following result:
SPLINES :
Wall Time = 6.1988831e-05
CPU Time = 3.8000000e-05
Cubic B-spline at x: -6.5891931e-02
ANALYTICS :
Wall Time = 9.5367432e-07
CPU Time = 0.0000000e+00
Analytical function at x: -6.5891931e-02
c++ code for the analytical function (c11;c22;c12) = (x,y,z) :
energy = energy = beta_Var*(pow(pow(c11-c22,2.0)*(c11-c12*4.0+c22)-pow(c11-c12*4.0+c22,3.0)*(1.0/9.0),2.0)*9.46969696969697E-4-pow(pow(c11-c12*4.0+c22,2.0)*(1.0/1.2E1)+pow(c11-c22,2.0)*(1.0/4.0),3.0)*(4.1E1/9.9E1)+(pow(c11-c12*4.0+c22,2.0)*(1.0/1.2E1)+pow(c11-c22,2.0)*(1.0/4.0))*pow(c11*(1.0/3.0)-c12*(1.0/3.0)+c22*(1.0/3.0),4.0)+(pow(c11-c22,2.0)*(c11-c12*4.0+c22)-pow(c11-c12*4.0+c22,3.0)*(1.0/9.0))*(pow(c11-c12*4.0+c22,2.0)*(1.0/1.2E1)+pow(c11-c22,2.0)*(1.0/4.0))*(c11*(7.0/3.0)-c12*(7.0/3.0)+c22*(7.0/3.0))*(1.0/6.6E1))+pow(pow(c11-c22,2.0)*(c11-c12*4.0+c22)-pow(c11-c12*4.0+c22,3.0)*(1.0/9.0),2.0)*(1.7E1/5.28E2)+(pow(c11-c22,2.0)*(c11-c12*4.0+c22)-pow(c11-c12*4.0+c22,3.0)*(1.0/9.0))*pow(c11*(1.0/3.0)-c12*(1.0/3.0)+c22*(1.0/3.0),3.0)+pow(pow(c11-c12*4.0+c22,2.0)*(1.0/1.2E1)+pow(c11-c22,2.0)*(1.0/4.0),3.0)*(4.0/1.1E1)-(pow(c11-c22,2.0)*(c11-c12*4.0+c22)-pow(c11-c12*4.0+c22,3.0)*(1.0/9.0))*(pow(c11-c12*4.0+c22,2.0)*(1.0/1.2E1)+pow(c11-c22,2.0)*(1.0/4.0))*(c11*(8.0/3.0)-c12*(8.0/3.0)+c22*(8.0/3.0))*(1.0/1.1E1) +gamma_Var*pow(burgers-detc,4.0)+theta_Var*pow(burgers-detc,2.0);
That is one impressive polynomial 😮
Still, I am not surprised by the run times. Even if it is slower, the run time for evaluating the B-spline is still very low. I would not be too worried, unless you plan on evaluating it millions of times, in which case you would have to wait several seconds.
Let us know how you proceed, and if the run time in fact is a problem for your application!
Edit: I will keep these run times in mind the next time I work on the B-spline evaluation code!
Hello,
Actually, it is a compiler related issue. I compared Intel icpc and c++ (clang Apple LLVM version). Intel optimizes very efficiently the analytical function that makes the calculation faster than SPLINTER. With clang++, SPLINTER performs better.
In my applications, I work with a numerical grid of, for instance, 1000X1000. In each grid node, I evaluate the analytical function and its derivative and I monitor the time evolution of the system. So the run time is crucial. But you are right, SPLINTER is fast enough to be used in the calculations.
I attached an image of a stress-field related to dislocations calculated with SPLINES. I think it is a world premier:-)
Thanks
ANALYTICS with intel:
Function at x: -0.0013437
CPU Time = 3.0000000e-06
SPLINES with intel:
Cubic B-spline at x: -1.3437012e-03
CPU Time = 8.0000000e-06
ANALYTICS with clang++:
Function at x: -0.0013437
CPU Time = 1.0000000e-05
SPLINES with clang++:
Cubic B-spline at x: -1.3437012e-03
CPU Time = 8.0000000e-06
Cool stuff. Happy you found a use for SPLINTER.
PS: In the current version, the derivative of the B-spline is computed point-wise. In a future release, we plan to also provide the derivative of a B-spline as a B-spline of one lesser degree. This would speed up the computations for use cases in which the derivative must be evaluated many times.
Dear Grimstad,
Sorry to repeat the same issue but the spline.eval is extremely slow when compared to the evaluation time of analytical functions. Now, we have finalized some serious testing. Do you know if we can obtain the parameters of the splines in the intervals on which they are defined? So that we can save them inside arrays that we can access whenever needed. This is what we do when we interpolate in Matlab and import the data into our c++ code; which is extremely fast. Our functions are 1D. So we basically need to have the s_i such that and x1_i
F = s0 +s1*(x-x1) + s2*(x-x1)^2 + s3*(x-x1)^3.
Thanks for your help
Hi again,
Yes, it is possible to to extract the piecewise polynomials of the B-spline. The procedure would be inline with the following for the univariate case. Assume that you have a cubic spline defined by the knot vector t = [0, 0, 0, 0, 1, 2, 3, 3, 3, 3].
- Decompose the B-spline by calling BSpline::decomposeToBezierForm(). This will increase the knot multiplicity to degree + 1 for all knots. The knot vector will become [0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3].
- Now, in each of the knot spans ([0, 1), [1, 2), [2, 3]), the B-spline will be defined by 4 basis functions. So, the first four control points belong to [0, 1], the next four to [1, 2], etc. You may store these control points along with the knot spans in a look-up table (or other data structure).
- To evaluate at x the B-spline function using analytical polynomials you must first find the corresponding knot span. Lets say that x is in [a, b), and that the corresponding control points are [c_1, c_2, c_3, c_4].
- Then, compute u = (x - a) / (b - a). u will be in [0, 1).
- Evaluate the four Bernstein basis functions of degree 3 at u. Lets us denote them [b_1, b_2, b_3, b_4]
- The function value of the B-spline at x will then be equal to y = c_1 * b_1 + ... + c_4 * b_4.
Hope this is useful.
PS: If you want to be extra clever, and speed up your computations even more, you may want to consider using the power basis instead of the Bernstein basis. This will require that you perform a change of basis by applying a linear map to the coefficients (unfortunately, I do not have time to explain the details now). You can do this in advance when you create the look-up table. This will allow you to evaluate the function as y = d_1 + d_2 * x + d_3 * x^2 + d_4 * x^3, where the ds are the new coefficients. And, you may even compute the powers recursively to save some additional computations. Good luck :)