Create a gridded lookup table of scattered data in n dimensions.
regularizeNd Fits a nD lookup table with smoothness to scattered data.
regularizeNd answers the question what is the best possible lookup table that the scattered data input x and output y in the least squares sense with smoothing? regularizeNd is meant to calculate a smooth lookup table given nD scattered data. regularizeNd supports extrapolation from a scattered data set.
The calculated lookup table, yGrid, is meant to be used with griddedInterpolant class with the conservative memory form. Call griddedInterpolant like
xGrid = cell array of grid vectors
smoothness = smoothness value or vector
yGrid = regularizeNd(xData, yData, xGrid, smoothness);
F = griddedInterpolant(xGrid, yGrid).
Desirable properties of regularizeNd:
Calculates a relationship between the input x and the output y without definition of the functional form of x to y.
Often the fit is superior to polynomial type fitting without the wiggles.
Extrapolation is possible from a scattered data set.
After creating the lookup table yGrid and using it with griddedInterpolant, as the query point moves away from the scattered data, the relationship between the input x and output y becomes more linear because of the smoothness equations and no nearby fidelity equations. The linear relationship is a good choice when the relationship between x and y is unknown in extrapolation.
regularizeNd can handle 1D, 2D, nD input data to 1D output data. RegularizeData3D and gridfit can only handle 2D input and 1D out (total 3D).
regularizeNd can handle setting the smoothness to 0 in any/some axis/dimension. This means no smoothing is applied in a particular axis/dimension and the data is just a least squares fit of a lookup table in that axis/dimension. Note this is not recommended and often can lead to an illconditioned fitting problem. However, I have found it useful so I left this as an option.
The source code is locate here:
https://github.com/jasonnicholson/regularizeNd
For an introduction on how regularization of a lookup table works, start here: https://mathformeremortals.wordpress.com/2013/01/29/introductiontoregularizingwith2ddatapart1of3/
Acknowledgement
Special thanks to Peter Goldstein, author of RegularizeData3D, for his coaching and help through writing regularizeNd.
1.8  Added iterative solvers, corresponding documentation, and an iterative example. 

1.7  Just some formatting changes to the description. 

1.7  Fixed a figure for example1D that didn't display well on the web. 

1.7  Added a 1D example. 

1.6  Some formatting was messed up on the gridfit example. I had to republish. 

1.6  For some reason the description was lost on the last update. I had to redo it. Ugh. 

1.6  Fixed a spelling mistake in the gridfit examples.


1.5  more updates to examples. 

1.5  Made some fixes to the formatted documentation. 

1.5  More updates to examples. 

1.5  More example updates. 

1.5  More updates to the examples. 

1.4  Added several examples.


1.3  Small update to the function description. 

1.3  Added support for cubic interpolation. This is a major step. All of the interpolation methods that I intended to use are implemented. 

1.2  Added a note about ndgrid format of yGrid. yGrid is in the ndgrid format at not the meshgrid format. There is no analogue of the meshgrid format for higher dimensions. In 2d, ndgrid is the transpose of the meshgrid format. 

1.2  Updating the description. 

1.2  Description formatting. 

1.2  Fixed two issues with the internal documentation and example.


1.1  Fixed a bug I introduced in the last week. 

1.0  Added a picture. 

1.0  Grammar fix to description. 

1.0  Minor update to the link in the description. 

1.0  Updated the description fixing some grammar errors. 

1.0  I just uploaded the single m file instead of attaching the repository like when it is connected to github. 
Inspired by: Surface Fitting using gridfit, RegularizeData3D
Michal Kvasnicka (view profile)
very good function ... thanks
Jason Nicholson (view profile)
To answer Blago's question about 'pchip', 'cubic', and 'spline' interpolation methods (found in griddedInterpolant class) as they apply to regularizeNd:
Given a set of x values, 'pchip', and 'spline' produce a nonlinear set of equations when solving for the lookup values, y. In interpolation, then nonlinearity is inconsequential because the y values are already known. Therefore, I would have to write a iterative solution to find the lookup table that uses the 'pchip' and 'spline' type interpolation methods. For now, I don't really want to deal with that issue.
Bhartendu (view profile)
Abhijit Das (view profile)
Very helpful submission. Thanks for the contribution, Jason.
Jason Nicholson (view profile)
Thanks for the feedback Blago. You have a good point. If I knew or learn how to implement pchip or spline, I would/will added it for sure.
Sorry that I just saw your comments. :)
Blagovest Mihaylov (view profile)
To clarify, the first statement in my previous comment depends on your input data and query points. In general, it is a regression rather than simple interpolation.
Also one could still use griddedInterpolant to increase the point density for most types of data. I'm just being pedantic. I still think the function is great and it saves me a lot of time. So thank you again.
Blagovest Mihaylov (view profile)
Hi Jason,
If the smoothing parameters is set to zero it becomes an interpolant itself. (Well you would only use the fidelity part for computational reasons)
My suggestion is that the usefulness of the function would improve if it matches griddedInterpolant, because you would be able to use regularizeNd to get a smooth look up table at some relativelyquicktocompute grid resolution and then use the interpolant to possibly increase the grid density of this smooth hypersurface. This is currently the case for linear, but this is more important for cubic anyway.
Unfortunately this is beyond me as well, but I do find it interesting.I'll let you know if I happen to implement cubic convolution or spline.
Jason Nicholson (view profile)
Blago,
Thanks for the feedback. :)
I am aware that regularizeNd is not the same as 'cubic' in griddedInterpolant. The intent of regularizeNd is a fitting algorithm of a smooth lookup table. It is not intended as an interpolant itself. It is nice that the linear case that you get back the same values as you put into regularizeNd but this is not a requirement nor was I trying to intend that this would happen. I intended the regularizeNd cubic option to do 3rd order Lagrange polynomial interpolation. The cubic option in regularizeNd is just intended to allow better curvature and shape preservation with fewer grid points. The regularizer helps keep the Lagrange polynomials well behaved on very unequally spaced grids. The regularization is why Lagrange polynomial interpolation is sufficient.
Implementing 'pchip', 'spline', and griddedInterpolant's 'cubic' is beyond the scope of what I know how to do. Maybe someday I can implement these interpolation types but I would have quite a bit of learning to do to get these other options working. If you know how, please by all means do it. If you want to talk more, we can figure out how to get in contact with each other. I am hesitant to post my email here.
Blagovest Mihaylov (view profile)
Dear Jason,
I find this function very useful and the linear interpolation is great. Thank you for you contribution.
I however think there is a small problem with the cubic interpolation, possibly being slightly different algorithm than the one used in Matlab in particular in griddedinterpolant. Matlab I think uses Cubic Convolution Interpolation and this code is based on Lagrange Cubic Interpolation. Frankly, I haven't looked into the details of these algorithms. I think calling it cubic implies it's consistent and that could be a problem. It was confusing in my case.
Say we use
[yGrid residuals]=regularizeNd(x,y,xGrid,'cubic');
Where the residuals are calculated within the function (residuals=A*yGrid y) and we only take the fidelity part, i.e. the first length(y) elements.
We then compare to the interpolated points at the original locations (x) via:
F=griddedinterpolant(xGrid,yGrid,'cubic');
y2=F(x);
The difference y2y is not the same as the residuals.This can be a small or a large problem depending on the use.
There is a however a perfect match when linear interpolation is used.
Finally,I would also suggest to implement spline or pchip interpolation if possible.
Once again thank you for your contribution.
Blago
BRYAN LIAO (view profile)
Thank you, Jason.
I was not aware of this awkward behavior meshgrid, ndgrid, as well as the associated awkward behavior of surf, mesh, etc.
After clearing that up, the results from regularizeNd match fairly well to other fitting functions griddata and gridfit.
Looks good so far.
Jason Nicholson (view profile)
Bryan,
you must use the ndgrid format. ndgrid format is the transpose of the meshgrid format in 2d. There is no analogue in higher dimensions to meshgrid and therefore you must use ndgrid.
x=100:2:100;
y=x;
[X,Y]=ndgrid(x,y);
Z=regularizeNd([data(:,2),data(:,3)],data(:,4),{x,y},.0001);
Z2=griddata(data(:,2),data(:,3),data(:,4),X,Y);
Z3=gridfit(data(:,2),data(:,3),data(:,4),x,y);
BRYAN LIAO (view profile)
Hello,
I have tried to use this function for 2D input data with 1D output. IE, position as X,Y and a function of data as Z. In this case, I could compare the data to that of gridfit and griddata.
Perhaps I am using the syntax incorrectly, but somehow the data for this function is flipped in X/Y compared to that of gridfit and griddata (as well as the raw data).
I am curious to use this function for higher dimensions, but until this issue is cleared up, I am hesitant to do so. Any ideas?
example code (for some data with xposition in col 2, yposition in col 3, and zdata in col 4)
x=100:2:100;
y=x;
[X,Y]=meshgrid(x,y);
Z=regularizeNd([data(:,2),data(:,3)],data(:,4),{x,y},.0001);
Z2=griddata(data(:,2),data(:,3),data(:,4),X,Y);
Z3=gridfit(data(:,2),data(:,3),data(:,4),x,y);