Fitting a Globally Monotonically Increasing Polynomial to Data in Python

I have a dataset to which I want to fit a nth degree univariate polynomial in a least squares fashion. The resulting polynomial should be monotonically increasing for all real numbers, not just the range of the data. I am using Python.

I've tried using cvxpy to create a problem where I added constraints that the derivative should be non-negative at specific points. This way I can create a polynomial that is monotonically increasing in the range of data, but it won't guarantee that the polynomial is monotonically increasing globally.

I figured it out. First of all, the degree of the polynomial must be odd.

Here's how I solved a 5th degree polynomial case:

A polynomial will be monotonically increasing globally if its derivative is a perfect square. In this case the degree of the derivative is 4. A fourth degree perfect square has the form (ax**2 + bx + c)**2. We can get a definition for our 5th degree polynomial if we expand this expression and then calculate its integral.

We start with parameters a, b, c, and d. d is serves as a bias parameter. We can then determine the coefficients of the original polynomial like this (the polynomial will be of form A*x**5 + B*x**4 + C*x**3 + D*x**2 + E*x + F):

A = a**2 / 5
B = a * b / 2
C = (2 * a * c + b**2) / 3
D = b * c
E = c**2
F = d

With these coefficient definitions we can directly fit the polynomial without needing to use any constraints.

A similar question is addressed here, where a package in R to do exactly this is discussed (MonoPoly). I don't know of a python implementation, but the R implementation details are available in this paper.

Depending on your problem, you might also be able to brute force it and compute the derivative everywhere.