`PiecewisePolynomial` vectorisation and speed issues

[duanyutong]

Before you post:

• [x] Have you reviewed the documentation on Getting Help?
• [x] Have you reviewed existing GitHub issues and StackOverflow posts?
• [x] Have you rerun install_prereqs.sh on your current commit of drake?

When using PiecewisePolynomial in pydrake, the evaluation of the function via .value(x) does not take vector input. This is likely the reason why drake is roughly a factor of 200 slower than the same functionality in scipy. It would be useful to support numpy ndarray as x.

All piecewise cubic splines have degree 3 and order 4 by definition, and automatically guarantees continuity in first and second derivatives by the choice of basis. The naming of the function CubicWithContinuousSecondDerivatives is very misleading, as it implies that there is a version of cubic spline in mathematics or an alternative algorithmic implementation that does not offer continuous second derivatives.

[SeanCurtis-TRI]

re: continuity

Generally, piecwise polynomials don't make any guarantees about continuity between the pieces. Thus, while a general cubic piecwise polynomial has C2 continuity inside each interval, it won't necessarily be continuous at the transition. CubicWithContinuousSecondDerivatives makes the strong promise that it will be continuous over the breaks (by making sure the first and second derivatives are shared as the break is approached from either side.) So, the name, while a bit unwieldy, is appropriate and does distinguish it from cases which lack that promise of continuity.

re: vectorised input. Yep, vectorized evaluation (such that the iteration happens in C++ instead of python) would definitely provide a pydrake speed boost.

[duanyutong]

A spline of order M on interval [a, b] by definition has (M-2) continuous derivatives on the entire interval [a, b], including at the knots. There are various texts/lecture notes that state this. Maybe in practice there are some bad spline solvers that aren't implemented properly and lack the continuity constraints? I understand that you'd like to emphasise the continuity which is fine.

To solve for cubic splines, even with the 2nd-order derivative continuity constraint there are only (4M-2) constraints, and one has to seek two extra constraints from some boundary condition (e.g. forcing naturalness at endpoints a and b). I don't know if it's possible for any other cubic spline solves to not guarantee 2nd-order derivative continuity.

[SeanCurtis-TRI]

True. But PiecewisePolynomail ≠ spline.

[duanyutong]

Understood. It leaves the average user wondering if CubicWithContinuousSecondDerivatives == CubicSpline or not. If so, maybe PiecewisePolynomial.CubicSpline would be a better name for the function to be consistent with the conventional name in literature, and the docstring can emphasise the continuity constraints.

[SeanCurtis-TRI]

I can certainly see that. Understanding the implications of the one function does depend on understanding the whole of PiecewisePolynomial.

However, it's not clear to me what the best way forward would be. Just calling it CubicSpline is ambiguous. Is it a Bezier spline? A bspline? Something else? Those two examples have significantly different semantics (and the latter can be given knots with high multiplicity that break continuity over the "breaks"). So, it may well be that any name would confuse someone, depending on their background.

The first general solution is to make sure that the documentation is clear. If the answer cannot be found in the documentation, that is a clear defect. Even we can find the perfect name, so much better. But even the perfect name doesn't excuse poor documentation.

[sherm1]

Assigning this to you for now @ToffeeAlbina-TRI for disposition.

[RussTedrake]

Just saw this.

We already have the vectorized method: https://drake.mit.edu/doxygen_cxx/classdrake_1_1trajectories_1_1_trajectory.html#a7607cc5c56514557ed9441479ae70e71

PiecewisePolynomial is absolutely a more general class than trajectories with continuous derivatives. The names we have chosen, though verbose, were chosen precisely to prevent the confusion you expressed. Indeed, there are plenty of ways to construct cubic piecewise polynomials that do not guarantee second-order continuity -- many of them are implemented. E.g.: https://drake.mit.edu/doxygen_cxx/classdrake_1_1trajectories_1_1_piecewise_polynomial.html#a5412c17187dde2f8984dcfb3eb67ca1d https://drake.mit.edu/doxygen_cxx/classdrake_1_1trajectories_1_1_piecewise_polynomial.html#aa360cbd2285db12eb014dc6e34868814

Barring a specific suggestion, I do not believe there is any action needed here.

[duanyutong]

Thanks for pointing to the vector_values function.

It is still over 10 times slower than scipy. Also note how my samples below are of shape m x n but I had to take its transpose. This is in contradiction to the doc which states that breaks are m x 1 and samples are m x n. Here in the n=1 case samples have to be 1 x m.