# Creating a function incrementally

Let's say that I need a computer program to generate a function that fits a data set which has $(x,y)$ data points.

Let's say that I get the data set one point at a time, I don't know how many data points there are, and I have a finite amount of memory for storage as I build up this function.

Are there any good ways to fit a function to a data set as it comes in, with limited storage space? I think I'm looking for an "online" algorithm for fitting a function to a data set.

My best guess at an algorithm would be to decide on a fixed amount of storage to use, and have that define the degree of a polynomial. For instance, I might decide on a cubic curve and store the coefficients of the power series of $x$; $A,B,C,D$ in the equation $y = Ax^3+Bx^2+Cx+D$.

The next part of the problem would be to start the coefficients at some value, and then as data points come in, I would account for their influence in the equation some how.

Making polynomials from data points is fairly easy, using something like Lagrange Interpolation, but unfortunately I'm not sure how to turn that into an "online" algorithm where I don't know all of the data points in advance.

With fixed memory constraints, and a fixed degree polynomial (when i may get any number of points), I know that my resulting equation will only be an approximation, and won't actually necessarily be a true representation of my data set, but getting as close as reasonably possible is a desired result.

Does anyone know any techniques or methods for doing this?

• Assuming your data points are of the form $(x_i,y_i)$ do you know in advance the smallest and largest possible values of the $x_i$ so that you can restrict your function to some interval $[x_{min},x_{max}]$. If so, perhaps you could break that interval into a fixed number of $n$ sub-intervals and do something like a $2nd$ or third degree spline on each sub-interval. With each new data point, update the parameters of the $n$ different splines to adjust for the added data point. Dec 19, 2016 at 22:41
• Yeah, both x and y have a known min and max value so we could restrict both to $\in [0,1]$. Interesting idea. Dec 19, 2016 at 22:45
• You might actually want the number of sub intervals to change as the number of data points to increase also. For example with only one or two data points your function would be just a straight line so you could let $n=1$ until you got at least three data points. Assuming storage is the main problem, not processing power, the subintervals would not have to be of equal length but could also be adjusted to accomodate a roughly equal number of data points in each subinterval. Lots of possibilities. Dec 19, 2016 at 22:50
• A problem I see with those things are that I'd have to use storage space for knowing how many points I've seen, and also for the length of the intervals. Maybe something clever can be thought up, but this is a "per pixel" thing in a computer graphics situation, so literally if I could have as few as 4 values of storage, that'd be the most attractive to me! Dec 19, 2016 at 22:52
• So are you wanting something like a running average for each pixel? Dec 19, 2016 at 23:01

Suppose you fix a degree $d$ in advance for your polynomial, and ask for a linear least-squares fit to your polynomial. With data $(x_j, y_j)$, $j = 1\ldots n$, the coefficients $c_0, \ldots, c_d$ give the least-squares solution to the system $Y = A c$ where $Y = [y_1, \ldots, y_n]^T$, and $A_{jk} = x_j^k$, $j=1\ldots, n$, $k = 0\ldots d$. This least-squares solution is given by solving the normal equations $$A^T A c = A^T Y$$ assuming $A^T A$ has rank $d+1$, which it does as soon as there are at least $d+1$ distinct $x_j$, $c = (A^T A)^{-1} A^T Y$. Note that $(A^T A)_{ik} = \sum_j x_j^{i+k}$, so $A^T A$ is constant on antidiagonals and requires only $O(d)$ storage, and is easily updated. For each new data point $(x_j, y_j)$ you add $x_j^{i+k}$ to $(A^T A)_{ik}$, and $x_j^i y_j$ to $(A^T Y)_i$, $i=0 \ldots, d$, $k = 0 \ldots d$.

• Thank you for the answer. I'm parsing it, but am not sure what you mean by $A_{jk} = x_j^k$, can you explain? Sorry, it's likely that I'm just not familiar with that notation. I don't understand how $x$ can have both a subscript (array index) as well as a super script (what does this mean?). Thanks! Dec 19, 2016 at 23:23
• The superscript is just an exponent. Thus in the case $x_j = 2$, $A_{j0} = 1$, $A_{j1} = 2$, $A_{j2} = 4$, etc. Dec 19, 2016 at 23:32
• I've worked through this and it's very nice. I doubt you can get the same results using less storage. Thank you so much Robert! Dec 20, 2016 at 19:14
• For other folks reading this question and answer in the future, I've written up a blog post about this, including some simple working C++ code, an interactive HTML5 demo, and some details you'd likely find interesting if implementing this in software towards a specific goal: blog.demofox.org/2016/12/22/… Dec 23, 2016 at 16:30

Let's look at the case where the trial function is $y(x) = a_{0} + a_{1}x + a_{2}x^{2}$. Start with a sequence of measurements $\left\{ x_{k}, y_{k} \right\}_{k=1}^{m}$.

The column vector structure of the system matrix is $\mathbf{A}=\left[ \mathbf{1} \ x \ x^{2} \right]$. The normal equations are \begin{align} \mathbf{A}^{*} \mathbf{A} x &= \mathbf{A}^{*} y \\ % \left[ \begin{array}{rrr} \mathbf{1}\cdot\mathbf{1} & \mathbf{1}\cdot x & \mathbf{1}\cdot x^{2} \\ x\cdot\mathbf{1} & x\cdot x & x\cdot x^{2} \\ x^{2}\cdot\mathbf{1} & x^{2}\cdot x & x^{2}\cdot x^{2} \\ \end{array} \right] \left[ \begin{array}{c} a_{0} \\ a_{1} \\ a_{2} \end{array} \right] &= \left[ \begin{array}{r} \mathbf{1} \cdot y \\ x \cdot y \\ x^{2} \cdot y \end{array} \right] \end{align} Of course, the system has yet to be solved. While we can solve this quickly, the process bogs down rapidly as the fit order grows.

Proceeding, we add a new point $\left( x_{m+1}, y_{m+1}\right)$. The updates follow this pattern: \begin{align} \mathbf{1} \cdot \mathbf{1} &: m \to m+1 \\ \mathbf{1} \cdot x &: \sum_{k=1}^{m} x_{k} \to \sum_{k=1}^{m} x_{k}+x_{m+1} \\ \mathbf{1} \cdot x^{2} &: \sum_{k=1}^{m} x_{k}^{2} \to \sum_{k=1}^{m} x_{k}^{2}+x_{m+1}^{2} \\ \vdots \end{align}

Your question seems to be this: given a solution vector $a$, how do we perturbatively update this solution when the data is perturbatively updated? Find $\Delta a$ given a new data point. We can write updates for the lowest order fits, but not for higher orders.