4
$\begingroup$

I have a pretty large dataset ($x,y$) consisting of a few million points. There is a lot of noise in the data. I want to find a smooth but simple approximation/representation for this dataset, so that for a given value of $x$ I have a decent estimate of $y$. I expect the data to have some kind of a cubic/non-linear relationship (hence simple linear regression won't work and since I want a simple representation, kernels are probably ruled out). I want to try doing this using cubic splines with a few knots (say $5$).

In other words, I want to find the cubic spline approximation with a specified number of knots that minimizes the l1 or l2 norm with respect to the input data points. Do you have any ideas about how to do this? Any pointers will be pretty useful. I found a couple of pointers online for approximating using a cubic spline, but none of them seemed scalable for the amount of data I have.

Thanks a lot! Gaurav

PS: Both $x$ and $y$ are scalar i.e. I have a simple 2D dataset, and am just trying to learn a simple representation for the non linear relationship between $x$ and $y$.

$\endgroup$
  • $\begingroup$ ...Interpolating a huge amount of data points into $C^1$ polynomial space efficiently is already very challenging...and now you are asking minimizing $L^1$- and $L^2$-norm simultaneously?? Oh my... $\endgroup$ – Shuhao Cao Jun 18 '13 at 4:55
  • $\begingroup$ Wait do you want to minimize L1 and L2 simultaneously or do you want only (any) one of them minimized? Also a plot of the data (or at least some of the data showing most of its important features for example) would help. In addition, how many points would you be interpolating at...like needing only one point is very different than interpolating 100,000 points? Are the points to be interpolating all bunched up together or are they spread over the entire domain of the data? $\endgroup$ – Fixed Point Jun 18 '13 at 8:01
  • $\begingroup$ Sorry in case I wasn't clear. I meant I want to minimize either the L1 or the L2 norm, not both of them simultaneously. I would be happy with either, whichever is easier. My guess is that minimizing L2 should be easier? The points are all bunched up together. So even though I have a large number of points, they are not too far from each other in the domain (in fact, since I have noise, I might end up having a bunch of points with the same x value and same or different y values). I was hoping to be able to interpolate very few points, not 100000. For example, just have 5/6 knots. $\endgroup$ – user10 Jun 18 '13 at 12:27
  • $\begingroup$ Maybe what I should do is to first have a grid (say distance between the points is uniform) over the x values, and compute the value of y at these grid points by minimizing L2 norm. I could do this to get 200-500 points. Then once I have only 200-500 points, I can just do a spline interpolation on these. Does this sound reasonable? $\endgroup$ – user10 Jun 18 '13 at 12:44
  • $\begingroup$ Splines or polynomials? $\endgroup$ – AnilB Jun 19 '13 at 20:00
1
$\begingroup$

If you're using Python, scipy.interpolate.LSQUnivariateSpline can easily handle 1M points:

""" LSQUnivariateSpline( x, y, user knots ) """
# http://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.LSQUnivariateSpline.html
# cf https://github.com/scipy/scipy/issues/2611  many knots -> all NaN

from __future__ import division
import sys
from time import time
import numpy as np
import pylab as pl
from scipy.interpolate import LSQUnivariateSpline

__date__ = "2013-08-02 aug denis"

    # LSQSpline nx points, interpolate at nsample --
nx = 1e6
nsample = 1000
ncycle = 1
noise = .5
nknot = 5
plot = 0
seed = 0

exec( "\n".join( sys.argv[1:] ))  # run this.py nx= ...  from sh or ipython
np.set_printoptions( 1, threshold=100, edgeitems=10, suppress=True )
np.random.seed(seed)
nx = int(nx)

title = "scipy.interpolate.LSQUnivariateSpline  nx %d  noise %.2g  %d knots" % (
    nx, noise, nknot )
print 80 * "-"
print title

def func( x ):
    return np.sin( 2*np.pi * ncycle * x )

x = np.sort( np.random.uniform( size=nx ))
xs = np.linspace( 0, 1, nsample )
y = func(x) + noise * np.random.randn(nx)

#...............................................................................
knots = np.linspace( xs[1], xs[-2], nknot )
    # knots = np.linspace( xs[0], xs[-1], nknot )
    # ValueError: Interior knots t must satisfy Schoenberg-Whitney conditions
t0 = time()

spline = LSQUnivariateSpline( x, y, knots, k=3 )
ys = spline(xs)  # interpolate at xs
print "%.2g sec  setup + interpolate %d points" % (time() - t0, nsample)
print "y interpolated:", ys
yexact = func(xs)
averr = np.fabs( ys - yexact ).mean()
print "av |yspline - yexact| %.2g" % averr

if plot:
    pl.title( title )
    pl.plot( ys, label="spline" )
    pl.plot( yexact, label="exact" )
    # pl.plot( ys - func(xs), label="diff" )
    pl.legend()
    pl.show()

Without Python (what language are you using ?), start with knots [ 0 .2 .4 .6 .8 1] and fit 5 separate cubics in each interval -- least-squares fit [1 x x^2 x^3].
If the ends are close enough at .2 .4 .6 .8, you're done.
If not, set up a global least-squares with 5 * 4 columns in all, plus dummy points at and near the knots with weights 1000000 to force continuity and smoothness there.

If the points x are very non-uniformly spaced, use non-uniform intervals, e.g.
np.percentile( x, [0,20,40,60,80,100] ) .

If you really want L1 approximation, I believe a cheap near-L1 can be done by iterative L2 with weighting; ask a new question on that.

Hope this helps. You might ask on stackoverflow, see
https://stackoverflow.com/search?q=[spline]+knots ,
and of course google "spline (intro|tutorial) user-knots" .

$\endgroup$
1
$\begingroup$

You can think of spline interpolation as linear regression with a special feature map.

Consider linear regression with coefficients $\beta_m$ and feature transformation $h_m(X)$: $$f(X) = \sum_{m=1}^{M} \beta_m h_m(X) = \beta^\top h(X)$$

The feature map $h_m(X)$ for a cubic spline with knots at $\xi_1$ and $\xi_2$ can be achieved by:

$$h_1(X) = 1, \quad h_2(X) = X \\ h_3(X) = X^2 \quad h_4(X) = X^3 \\ h_5(X) = (X - \xi_1)_+^3 \quad h_6(X)_ = (X- \xi_2)_+^3$$

with $t_+$ denoting the positive part. By this you would map your data from a vector $\mathbb{R}^{N\times1}$ to a matrix of $\mathbb{R}^{N\times6}$.

This is a purely cubic spline approach, but often times one would prefer to have linear segments at the borders. This is commonly called natural cubic splines. See [1] for further details. Often knots are chosen to be equally spaced, or placed at the quantiles of the data.

With that feature map at hand you can perform linear regression as usual. That is, whenever you feel your data set is to large for the analytic solution, switch to the stochastic (mini batch) gradient descent.

That is the error of $$E(\beta) = \frac{1}{2} \sum_{n=1}^{N} \left( t_n - \beta^\top h(x_n) \right)^2$$ Results in a gradient $$\nabla E(\beta) = \sum_{n=1}^{N} \left( t_n - \beta^\top h(x_n) \right) h(x_n)^\top $$ And we descent on individual points (or subsets) $$\beta^{(t+1)} = \beta^{(t)} + \eta \nabla E_n $$ $$\beta^{(t+1)} = \beta^{(t)} + \eta (t_n - \beta^{(t)^\top} h(x_n) ) h(x_n) $$

You can use the sklearn implementation: https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.SGDRegressor.html It also has an option for l2 and l1 regularization, or even both (elasticnet).

[1] Elements of Statistical Learning, Hastie, Tibshirani, Friedman - Chapter 5

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.