This is exactly what I tried to say in my first post. Using few samples per point (thus large error), but more points (which can capture the gradients more easily) with spatial filtering should be better than exact interpolation. For example like this (blue line is the exact function, dots are sampled estimates and the dotted line is the smoothed curve across the points)

MWE
using DIVAnd
using PyPlot
realf(x)=x^2
NSamples=200
noiselevel=1.0
cgrid=collect(-1.0:0.01:1.0)
csample=-1.0 .+ 2.0.*rand(Float64,NSamples)
fsample=realf.(csample) .+ noiselevel*randn(Float64,NSamples)
figure()
plot(cgrid,realf.(cgrid),"-",csample,fsample,".")
?DIVAndfun
apprfun=DIVAndfun((csample,),fsample;xi=(cgrid,),len=1.0,epsilon2=2.0)
figure()
plot(cgrid,realf.(cgrid),"-",csample,fsample,".",cgrid,apprfun.(cgrid),":")