MATLAB Examples

# 4D Examples

This example demonstrates the use of regularizeNd with a 4D input and 1D output dataset. Note, that selecting the smoothness isn't always obvious until you have gained experience using regularizeNd. Often the best choice is to change the smoothness by factors of 10 until you have reached the smooth hyper surface that you want or reduce the smoothness until you have reached the fidelity you want. Even if my first guess of the smoothness is adequate, I change the smoothness by factors of 10 until it isn't smooth to see the effect of the smoothness parameter.

## Contents

Note this example has 30,240 grid points. This means the linear system solved, x = A\y, has 250,505 rows and 30,240 columns with 2,539,252 nonzero elements in A. This solves on my laptop in 6 seconds.

```clc; clear; close all; ```

convert to double because I saved them as single to save space.

```load('dataset_4d_1.mat') % inputs. x1 = double(x1); y1 = double(y1); z1 = double(z1); t1 = double(t1); % ouput w1 = double(w1); ```

## Look at the Input and Output Vectors

The main point is there is a lot of points. The input space is well filled. I am not implying good fill is necessary. Sometimes it is necessary. Sometimes it isn't. The fill need in the input space depends on each case.

```figure; histogram(x1) xlabel('x'); figure; histogram(y1) xlabel('y'); figure; histogram(z1) xlabel('z'); figure; histogram(t1) xlabel('t'); figure histogram(w1) xlabel('w'); figure; scatter3(x1,y1,z1, [], t1, '.'); view(-87.2, 12.4); xlabel('x'); ylabel('y'); zlabel('z'); h = colorbar; ylabel(h, 't'); title('Input Space fill'); ```

## Creating the regularized hypersurface

Note with regularizeNd, the grid vectors must always span the input data. The reason for this is that calculating the weights for a given query point is well defined for interpolation and it not well defined for extrapolation. Understand though, that the fact that the grid vectors can far exceed the expanses of the fitted data and the fitted surface is guaranteed to change close to linearly far away from fitted data is of great advantage over other fitting technieques. It allows extrapolation away from the fitted data in a known way. In polynomial fitting, there is no guarantee how the polynomial will behave outside of the test data set.

```% setup the grid points xGrid = linspace(min(x1), max(x1), 20); yGrid = 0:10:110; zGrid = [-0.11 0:5:30 50:50:500]; tGrid = [50 65 75 80 90 100 115]; gridPoints = {xGrid, yGrid, zGrid, tGrid}; % Setup to loop smoothness = [1e-4, 1e-3, 1e-3, 1e-2]; % smoothness = [7e-3, 1e-2, 8e-3, 1e-1]; % better values whileFlag = true; tPlotLevels = [65 85]; yPlotLevels = [20 40 60 80]; gridForPlots = cell(4,1); [gridForPlots{:}] = ndgrid(xGrid, yPlotLevels, zGrid, tPlotLevels); surfaceHandles = nan(length(tPlotLevels), length(yPlotLevels)); while whileFlag display(smoothness); % Calculate regularized surface and time it tic; wGrid = regularizeNd([x1,y1,z1,t1], w1, gridPoints, smoothness, 'linear'); toc; % create the griddedInterpolant function wFunction = griddedInterpolant(gridPoints, wGrid, 'linear'); % calculate values for the plots wPlot = wFunction(gridForPlots{:}); % create the plot levels. Often I use a lot more than 8 levels. Often, % I will put them on the same figure/axes. However, for simplicity and % clarity, I left them as individual figures here. for iT = 1:length(tPlotLevels) for iY = 1:length(yPlotLevels) if ishghandle(surfaceHandles(iT,iY)) set(surfaceHandles(iT,iY), 'Zdata', squeeze(wPlot(:,iY,:,iT))'); else figure; surfaceHandles(iT,iY) = surf(xGrid, zGrid, squeeze(wPlot(:,iY,:,iT))'); xlabel('x');ylabel('z'); zlabel('w'); title(sprintf('Level y = %g, Level t = %g. Use to check smoothness.', yPlotLevels(iY), tPlotLevels(iT))); end end end % Ask for a new smoothness answer = input('Enter a new smoothness to recaculate surface.'); % parse the answer if isempty(answer) % Do nothing elseif any(answer == -1) whileFlag = false; elseif any(answer<0) warning('Negative values of smoothness not excepted. No changes made.'); else if any(length(answer) == [1 length(gridPoints)]) smoothness = answer; else warning('smoothness must be either scalar or have same number elements as the length of gridPoints. No change to smoothness was made.'); end end end ```
```smoothness = 0.0001 0.0010 0.0010 0.0100 Elapsed time is 6.184256 seconds. Enter a new smoothness to recaculate surface.[7e-3, 1e-2, 8e-3, 1e-1] smoothness = 0.0070 0.0100 0.0080 0.1000 Elapsed time is 6.050388 seconds. Enter a new smoothness to recaculate surface.-1 ```

## Look at the residuals

```% calculate the residuals residualsAbsolute = w1 - wFunction(x1,y1,z1,t1); figure; histogram(residualsAbsolute, 'FaceColor', 'auto', 'EdgeColor', 'auto') xlabel('Residual Size'); ```