Polynomial fitting where polynomial must be monotonically increasing

Given a set of monotonically increasing data points (in 2D), I want to fit a polynomial to the data which is monotonically increasing over the domain of the data. If the highest x value is 100, I don't care what the slope of the polynomial is at x=101.

I'd also like to keep the degree of the polynomial below 7. My data usually has about 20 (x,y) pairs. Least squares fitting is probably best, but I'm open to other techniques.

I was thinking that maybe I could create a b-spline, sample that b-spline (generating, say, another 100 data points), then perform a polynomial fit. Does that approach sound sound?

• In response to alex' questions, I have attempted to reduce the vagueness of my question by replacing the part about "minimize the degree of the polynomial, while fitting fairly tightly to the data" with "keep the degree of the polynomial below 7". – splicer Aug 30 '11 at 3:54
• I'm going to hold off on selecting my accepted answer until I fully understand all the answers and can judge them fairly. Thanks so much! – splicer Aug 30 '11 at 3:58

I think the problem is not precisely stated: its unclear what "minimize the degree of the polynomial, while fitting fairly tightly to the data" means (what is "fairly tightly?").

Anyway, for fixed degree, this may be formulated as a semidefinite program which can be solved approximately in polynomial time; there is a lot of software out there that can solve semidefinite programs very efficiently, e.g., sedumi.

Indeed, suppose your data points are $(x_1,y_1), (x_2,y_2), \ldots, (x_n,y_n)$. Let $y$ be the vector that stacks up the $y_i$ and let $V$ be the Vandermonde matrix $$V = \begin{pmatrix} 1 & x_1 & x_1^2 & \cdots& x_1^n \\ 1 & x_2 & x_2^2 & \cdots & x_2^n \\ \vdots & \vdots & \vdots & & \vdots \\ 1 & x_n & x_n^2 & \cdots & x_n^n \end{pmatrix}$$ Then, you'd like the minimize $||y-Vp||_2^2$ subject to the constraint that the entries of the vector $p$, which we will call $p_i$, correspond to a monotonic polynomial on your domain $[a,b]$, or, in other words, $p'(x) \geq 0$ on $[a,b]$: $$p_1 + 2p_1 x + 3 p_2^2 + \cdots n p_n x^{n-1} \geq 0 ~~~ \mbox{ for all } x \in [a,b].$$ This last constraint can be written as a semidefinite program,as outlined in these lecture notes. I will briefly outline the idea.

A univariate polynomial which is nonnegative has a representation as a sum of squares of polynomials. In particular, if $p'(x) \geq 0$ on $[a,b]$, and its degree $d$ is, say, even, then it can be written as $$p'(x) = s(x) + (x-a)(b-x) t(x),~~~~~~~~(1)$$ where $s(x), t(x)$ are sums of squares of polynomials (this is Theorem 6 in the above-linked lecture notes; a similar formula is available for odd degree). The condition that a polynomial $s(x)$ is a sum of squares is equivalent to saying there is a nonnegative definite matrix $Q$ such that $$s(x) = \begin{pmatrix} 1 & x & \cdots x^{d/2} \end{pmatrix} Q \begin{pmatrix} 1 \\ x \\ x^2 \\ \vdots \\ x^{d/2} \end{pmatrix}.$$ This is Lemma 3 in the same lecture notes.

Putting it all together, what we optimize are the entries of the matrices $Q_1,Q_2$ that make the polynomials $s(x)=x^T Q_1 x,t(x)=x^T Q_2 x$ sums of squares, which means imposing the constraints $Q_1 \geq 0, Q_2 \geq 0$. Then Eq. (1) is a linear constraint on the entries of the matrices $Q_1,Q_2$. This gives you you have a semidefinite program you can input to your SDP solver.

• Sorry, I'm an artist, not a mathematician. ;) By "fairly tightly", I meant "as close as possible to intersecting every point while maintaining the constraint that the polynomial must be monotonically increasing". – splicer Aug 30 '11 at 3:27
• @splicer - right, but the question is how do you trade that off with the degree. For example, suppose a degree $4$ polynomial is off by a max of $0.1$, but a degree $5$ polynomial is off by a max of $0.05$. Do you go for degree four or five? This kinds of information needs to be specified before a complete answer can be given. – morgan Aug 30 '11 at 3:30
• The reason for minimizing the degree is the reduce the time it would take a computer program to evaluate the polynomial for a given input. Any polynomial with a degree less than 7 should be fine. In your example, I would pick degree 5 since it's more accurate, yet less than degree 7. Sorry for being vague... – splicer Aug 30 '11 at 3:44
• I've implemented this method and found it less than practical. It turns out that the matrix Q is extremely ill-conditioned in many cases, thus preventing the solver from converging. – Victor May Mar 13 '12 at 17:08

Ignoring the question of minimizing degree, it is certainly possible.

One can find a positive function $f$ for which, if $p'$ is close to $f$ on the whole interval then $p = \int^x f(t) dt$ (with integration constant chosen so that $p(x_1)$ is close to $x_1$), then $p(x)$ will approximate the data points to given accuracy. Choose some positivity-preserving approximation scheme such as Bernstein polynomials to get a positive polynomial $p'$ within any desired distance of $f'$ on an interval containing all $x_i$. Integrate to produce $p$.

• Having written $p$ as a positive linear combination of, say, antiderivatives of Bernstein polynomials, can't one directly minimize $\sum \lvert p(x_i) - y_i \rvert^2$, subject to positivity of the basis coefficients? I don't see the need to introduce an intermediate function $f$. – Rahul Aug 30 '11 at 4:54
• $f$ is used to show that (Bernstein or other) polynomials of degree N suffice for some finite N that can be extracted from the proof. Given a vector space known to contain solutions, such as degree N polynomials, one can search for an optimal one by least-squares or linear programming. But how do you bound N without first going to a larger space of functions where there is more room to move, and then showing that polynomials can approximate that (positively)? – zyx Aug 30 '11 at 7:11

Well, it is definitely not the best answer from the practical point of view, however, it does two things. Firstly, it proves an existence of such polynomial. Secondly, by working out some constants, one can try to get an upper bound on the degree.

Suppose, that we are given $x_0<x_2<...x_n$ and $y_0<y_2<...y_n$ and we are looking for an increasing polynomial such that $P(x_i)=y_i,$ $0\le i\le n.$ It easy to construct a $C^1[x_0,x_n]$ function $f,$ such that $f(x_i)=y_i$ and $f'(x)\ge\delta>0$ for some real $\delta.$ For each $\varepsilon > 0$ there exists a polynomial $P_{\varepsilon}$ such that $$\|f'-P_{\varepsilon}\|\le \varepsilon.$$ Consider $$Q_{\varepsilon}(x)=y_0+\int_{x_0}^{x}P_{\varepsilon}(t)dt.$$ Then, $\|f-Q_{\varepsilon}\|\le \varepsilon (x_n-x_0).$ Let $R_{\varepsilon}$ be the Lagrange polynomial, which interpolates $R_{\varepsilon}(x_i)=y_i-Q_{\varepsilon}(x_i),$ $0\le i\le n.$ Consider a liner operator $A: \mathbb{R}^{n+1}\to \mathbb{R}^m,$ such that $A(z_0,z_1,...z_n)=R'(z_0,...z_n),$ where $R'(z_0,...z_n)$ is a derivative of Lagrange polynomial $R,$ such that $R(x_i)=z_i.$ Since this operator is a linear operator between two finite dimensional spaces, there is a constant $C,$ such that $\|R'\|\le C\max |z_k|.$ Choose $\varepsilon > 0$ such that $C(x_n-x_0)\varepsilon+\varepsilon \le \delta.$ For $R_{\varepsilon}=R_{\varepsilon}(y_0-Q_{\varepsilon}(x_0)...y_n-Q_{\varepsilon}(x_n)),$ we have $\|R'\|\le \delta-\varepsilon$ and it is easy to see that $R_{\varepsilon}+Q_{\varepsilon}$ satisfies all the requirements.

• I'm going to have to polish-up on my notation to be able to fully understand your answer; it's been a few years ;) – splicer Aug 30 '11 at 3:06

I'm presuming that since you said "fit", the polynomial doesn't have to pass through your given points.

I would say that LS indeed is good enough for your purposes. You can do things stepwise (if you fit using orthogonal polynomials, this is quite easily done) and either monitor your polynomials graphically (see plots of your polynomial and the original points) or monitor quantities like the (adjusted) $R^2$ as you increase the degree of your fitting polynomial.

• Correct, the polynomial does not have to pass through the data points. – splicer Aug 30 '11 at 2:33
• Does "LS" stand for "linear spline"? If so, how should I go about fitting the polynomial? Do I simply sample an arbitrary number of points from the linear spline, then perform the polynomial regression, or is there another way? – splicer Aug 30 '11 at 2:40
• I meant "least squares"; somehow, I don't think you'd want to do "connect-the-dots" with straight lines on your data. :) Maybe start with a quadratic fit and see if the fits get better as the degree goes higher. – J. M. is a poor mathematician Aug 30 '11 at 2:45
• But a straightforward least-squares fitting of polynomials won't ensure that the polynomial is monotonically increasing over the domain of the data, will it? – splicer Aug 30 '11 at 2:59
• @splicer: indeed, hence I recommended monitoring your fits graphically. – J. M. is a poor mathematician Aug 30 '11 at 3:56