File Exchange

jasonnicholson/regu​larizeNd

version 1.8 (5.61 MB) by

Create a gridded lookup table of scattered data in n dimensions.

Updated

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 n-D 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 ill-conditioned 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/introduction-to-regularizing-with-2d-data-part-1-of-3/

Acknowledgement
Special thanks to Peter Goldstein, author of RegularizeData3D, for his coaching and help through writing regularizeNd.

Michal Kvasnicka

Michal Kvasnicka (view profile)

very good function ... thanks

Jason Nicholson

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 non-linearity 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

Abhijit Das

Abhijit Das (view profile)

Very helpful submission. Thanks for the contribution, Jason.

Jason Nicholson

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.

Blagovest Mihaylov

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

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 relatively-quick-to-compute 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

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

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 y2-y 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

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

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

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);

MATLAB Release
MATLAB 9.1 (R2016b)
Acknowledgements

Inspired by: Surface Fitting using gridfit, RegularizeData3D

Play today