Initial commit

Demo_test.m

+%% Info
+% Demo program to demonstrate the use of the a kernel model trained as in
+% Demo_training.m
+clear all, close all
+startup()
+
+
+%% Path of the model (created by Demo_training)
+modelFile = 'Models/ModelExample.mat'; 
+
+
+%% Load the kernel model
+modelPar = importdata(modelFile);
+
+
+%% Initialize the model
+d = 3; % input dimension
+X = rand(1, 3); % a random d-dimensional point 
+kernelModel(X, modelPar); % set modelPar
+
+
+%% Compute predictions
+% Now the model can be evaluated for a new vector of inputs from
+% everywhere (provided 'kernelMode' is in the path), without passing again
+% the parameters modelPar
+Xte = rand(1000, d); % example of test inputs
+ste = kernelModel(Xte); % kernel model evaluated on Xte

Demo_training.m

+%% Info
+% Demo program to demonstrate the use of the
+% Vectorial Kernel Orthogonal Greedy Algorithm (VKOGA)
+% to train a kernel model.
+% The code validates the RBF shape par. and the regularizatio par.
+clear all, 
+close all
+startup()
+
+
+%% Path of various files
+% Data file
+dataFile = 'Data/DataExample.mat';
+% Permutation file used to split the data into train/validation/test.
+% If it exists, this permutation is used. Otherwise a random permutation
+% is used and saved with this name
+permFile = 'Random/DataPermutationExample.mat';
+% File where to save the trained model. Leave empty or remove to avoid saving
+modelFile = 'Models/ModelExample.mat';
+% Print info during training
+verboseFlag = 1;
+% Plot results after training
+plotFlag = 1;
+
+
+%% Parameters for VKOGA training and validation
+% Data parameters
+dataPar.percTe = 0.1;            % fraction of data to use as test set
+dataPar.normalization = 'none';  % 'none' or 'uniform' (scale to [0, 1] ^ d)
+dataPar.k = 5;                   % k-fold cross validation
+
+
+% Training parameters
+trPar.rbf = getRbf('gauss');     % rbf kernel (see getRbf.m)
+trPar.tolF = 1e-10;              % tolerance on the residual
+trPar.tolP = 1e-20;              % tolerance on the power function
+trPar.maxExpSize = 100;          % max number of selected points
+trPar.gType = 'pgreedy';         % greedy selection rule (fgreedy, pgreedy, fpgreedy)
+
+
+% Validation parameters
+valPar.spaceType = 'log';        % 'log' or 'lin' spacing
+valPar.minEp = 1e-1;             % min value of shape parameter
+valPar.maxEp = 1e1;               % max value of shape parameter
+valPar.nEp = 5;                 % number of shape parameters
+valPar.minLambda = 1e-10;        % min value of regularization parameter
+valPar.maxLambda = 1e-2;         % max value of regularization parameter
+valPar.nLambda = 4;              % number of regularization parameters
+valPar.tolFEp = trPar.tolF;
+valPar.tolPEp = trPar.tolP;
+valPar.maxExpSizeEp = trPar.maxExpSize;
+
+
+%% Load data
+% It assumes dataFile contains 'input' and 'output', where input is
+% an Nxd matrix of N d-dimensional points and output is a Nxq matrix of
+% N q-dimensional points.
+% Any other way to produce X, f is fine, as long as they contain one point 
+% per row.
+data = load(dataFile);
+X = data.input;
+f = data.output;
+
+
+%% Create dataset
+% - normalize data
+% - permute them
+% - split into train/validation/test
+
+% Load or create a permutation
+if exist('permFile', 'var') && ~isempty(permFile) && exist(permFile, 'file')
+    dataPar.perm = importdata(permFile);
+else
+    warning('Using random permutation')
+    pp = randperm(size(X, 1));
+    dataPar.perm = pp;
+    save(permFile, 'pp')
+end
+
+% Create the actual dataset
+dataset = setData(X, f, dataPar);
+clear('X'), clear('f')
+
+
+%% Train and save the model
+fprintf('Training kernel model\n')
+[modelPar, trainResults, valResults, testResults] = trainModelReg(dataset, valPar, trPar, verboseFlag);
+
+
+if exist('modelFile', 'var') && ~isempty(modelFile)
+    save( modelFile, 'modelPar')
+end
+
+
+%% Plot training/validation/test errors
+if plotFlag
+    % size of model
+    n = size(modelPar.a, 1);
+    
+    % Training error
+    figure(1)
+    semilogy(1 : n, trainResults.fmax, '.-', 'linewidth', 1.5)
+    xlabel('Kernel expansion size')
+    ylabel('Training error')
+    set(gca, 'fontsize', 14)
+    grid on
+    
+    % Power function on training set
+    figure(2)
+    semilogy(1 : n, trainResults.pmax, '.-', 'linewidth', 1.5), hold on
+    xlabel('Kernel expansion size')
+    ylabel('Max of power function')
+    set(gca, 'fontsize', 14)
+    grid on
+    
+    % Validation error
+    if length(valResults.Lambda) > 1 && length(valResults.Eps) > 1
+    figure(3)
+    [xx, yy] = meshgrid(valResults.Lambda, valResults.Eps);
+    if strcmp(valPar.spaceType, 'log')
+        plot(log10(valResults.lambda), log10(valResults.ep), 'ro'), hold on
+        contour(log10(xx), log10(yy), reshape(log10(valResults.errVal), size(xx)), 20, 'linewidth', 1.5)
+        hold off
+    else
+        plot(valResults.lambda, valResults.ep, 'ro'), hold on
+        contour(xx, yy, valResults.errVal, 20, 'linewidth', 1.5)
+        hold off
+    end
+    legend([], 'Selected parameters')
+    xlabel('Regularization parameter')
+    ylabel('Shape parameter')
+    title('Validation error (log_{10})')
+    set(gca, 'fontsize', 14)
+    grid on
+    end
+    
+    if dataPar.percTe > 0
+        % Test error (max and mean)
+        figure(4)
+        semilogy(1 : n, testResults.maxErr, '.-', 1 : n, testResults.meanErr, '.-', 'linewidth', 1.5)
+        xlabel('Kernel expansion size')
+        ylabel('Test error')
+        set(gca, 'fontsize', 14)
+        legend('Max', 'Mean')
+        grid on
+    end
+end
+
+

Functions/autoNdGrid.m

+function [Z, N, varargout] = autoNdGrid(N, d)
+% [Z, N] = AUTONDGRID(N, d) returns a uniform grid in [-1, 1] ^ d consisting 
+% of approximately N points. The actual number of points is returned.
+% If d == 2, the function optionally returns 
+%   varargout{1} = xx, 
+%   varargout{2} = yy,
+% with X = [xx(:), yy(:)] (useful for surf plotting).
+%
+% WARNING: the final N may be much larger than the input N if d is large.
+
+box = repmat([-1 1], d, 1);
+
+n = ceil(N ^ (1 / d));
+x = linspace(-1, 1, n)';
+
+
+ndArgs = cell(d, 1);
+for j = 1 : d
+    ndArgs{j} = x * (box(j, 2) - box(j, 1)) / 2 + (box(j, 2) + box(j, 1)) / 2;
+end
+
+[Z{1 : d}] = ndgrid(ndArgs{:});
+if nargout == 4 && d == 2
+    varargout{1} = Z{1}; % xx
+    varargout{2} = Z{2}; % yy
+end
+Z = cat(d + 1, Z{:});
+Z = reshape(Z, [], d);
+N = size(Z, 1);
+
+end
\ No newline at end of file

Functions/distanceMatrix.m

+function D = distanceMatrix(X, Y)
+% D = DISTANCEMATRIX(X, Y) computes the distance matrix D_ij = ||X_i - Y_j||
+
+% Try to use built-in distance matrix
+if exist('pdist2', 'file')
+    D = pdist2(X, Y);
+else
+    D = sqrt(bsxfun(@plus, sum(X .^ 2, 2), sum(Y .^ 2, 2)') - 2 * X * Y');
+end
+
+end
\ No newline at end of file

Functions/getRbf.m

+function rbf = getRbf(type, varargin)
+% GETRBF Returns a strictly positive definite RBF kernel.
+%   rbf  =  GETRBF(type) returns a function handle representing the radial
+%   kernel of type 'type'.
+%   The handled function is in the form
+%       rbf(ep, r)
+%   where ep is the shape parameter and r is the radial variable.
+%   Additional parameters (d, k) are required for the Wendland kernels. 
+%
+%   See also getWen
+
+
+switch type
+    case 'gauss'
+        rbf = @(ep, r) exp(-(ep*r).^2);
+    case 'imq'
+        rbf = @(ep, r) 1./sqrt(1+(ep*r).^2);
+    case 'wen'
+        if nargin ~= 3
+            error('Wendland kernels: getRbf("wen", d, k)')
+        else
+            rbf = getWen(varargin{1}, varargin{2});
+        end
+    otherwise
+        error(['RBF "' type '" not available'])
+end
\ No newline at end of file

Functions/getWen.m

+function rbf = getWen(d, k)
+% GETWEN returns a Wendland RBF.
+%   rbf = GETWEN(d, k) returns a function handle representing the Wendland 
+%   kernel of smoothness k in dimension d, where k is an integer in [0, 4].
+%   The handled function is in the form
+%       rbf(ep, r)
+%   where ep > 0 is the shape parameter and r > 0 is the radial variable.
+%
+%   See also getRbf
+
+
+switch k
+    case 0
+        p = '1';
+    case 1
+        p = '(l + 1) * r + 1';
+    case 2
+        p = '(l + 3) * (l + 1) * r .^ 2 + 3 * (l + 2) * r + 3';
+    case 3
+        p = '(l + 5) * (l + 3) * (l + 1) * r .^ 3 + (45 + 6 * l * (l + 6)) * r .^ 2 + (15 * (l + 3)) * r + 15';
+    case 4
+        p = '(l + 7) * (l + 5) * (l + 3) * (l + 1) * r .^ 4 + (5 * (l + 4) * (21 + 2 * l * (8 + l))) * r .^ 3 + (45 * (14 + l * (l + 8))) * r .^ 2 + (105 * (l + 4)) * r + 105';
+    otherwise
+        error('Order between 0 and 4!')
+end
+
+l = floor(d / 2) + k + 1;
+c = factorial(l + 2 * k) / factorial(l);
+e = l + k;
+
+p = ['max(1 - r, 0) .^ e .* (', p, ') / c']; 
+
+p = strrep(p, 'l', num2str(l));
+p = strrep(p, 'c', num2str(c));
+p = strrep(p, 'e', num2str(e));
+p = strrep(p, 'r', '(ep * r)');
+
+p = ['@(ep, r) ', p]; 
+
+rbf = str2func(p);
+
+end
+
+
+
+
+
+
+
+

Functions/kernelModel.m

+function f = kernelModel(x, varargin)
+% f = KERNELMODEL(x, varargin) evaluates a kernel model.
+%   f = KERNELMODEL(x, par) evaluates the kernel model defined by the
+%   structure par on the Nxd dimensional point x and returns the Nxq
+%   dimensional result f (one point per row).
+%   The par structure needs to be passed as input only at the first call of
+%   the function.
+
+persistent par
+if nargin > 1
+    par = varargin{1};
+    par.rbf = str2func(par.rbf);
+end
+
+x = bsxfun(@rdivide, bsxfun(@minus, x, par.translX), par.scaleX);
+f = par.rbf(par.ep, distanceMatrix(x, par.X)) * par.a;
+f = bsxfun(@plus, bsxfun(@times, f, par.scaleF), par.translF);
+
+end
\ No newline at end of file

Functions/newtonYYDataMF.m

+function [a, indI, n, c, Cut, fmax, pmax, Vx] = newtonYYDataMF(Afun, f, tolF, tolP, maxIter, gType, diagonal)
+% NEWTONYYDATAMF runs a matrix-free VKOGA.
+%   [a, indI, n] = NEWTONYYDATAMF(Afun, f, tolF, tolP, maxIter, gType)
+%   returns the indexes indI of the selected points, the coefficient
+%   matrix a of the corresponding interpolant, and the number n of selected
+%   points.
+%   The points are selected from a training set Xtr, which is implicitly
+%   defined via the kernel matrix function Afun, which accepts two indexes
+%   vectors as input, e.g.,
+%   Afun((1 : 10)', (2 : 20)')
+%   and returns the corresponding submatrix of the kernel matrix of Xtr.
+%   The interpolant can then be evaluated in the test points Xte as
+%   ste = ker(ep, Xte, Xtr(indI, :)) * a;
+%   where ker = ker(ep, x, y) is a SPD kernel.
+%   The matrix f contains the evaluations of the unknown function on the
+%   training points, one point per row.
+%   The algorithm is terminated by tolerances tolF on the square of the
+%   residual, tolP on the square of the power function, or when maxIter
+%   points are selected.
+%   The greedy selection is made by the criterion gType, which can be one
+%   among fgreedy, pgreedy, fpgreedy.
+%   The implementation is based on the Matrix-Free implementation of the
+%   Newton basis by M. Pazouki and R. Schaback. This means that only the
+%   rows/columns of the kernel matrix corresponding to the selected points
+%   are computed, while the full kernel matrix is never stored nor
+%   computed.
+
+
+[N, q] = size(f);               % number of points and dimension of output
+maxIter = min(maxIter, N);      % reduces the max number of points if > N
+indI = zeros(maxIter, 1);       % selected indexes
+notIndI = (1 : N)';             % not yet selected indexes
+Vx = zeros(N, maxIter);         % values of the Newton basis on the points
+c = zeros(maxIter, q);          % coefficients of the interpolant/expansion
+fmax = zeros(maxIter, 1);       % maximal residuals on training set
+pmax = zeros(maxIter, 1);       % maximal power function on training set
+if ~exist('diagonal', 'var') || isempty(diagonal)
+    p = Afun(1, 1) * ones(N, 1);    % power function squared (assuming transl. inv. ker.)
+else
+    p = diagonal;
+end
+Cut = zeros(maxIter);           % matrix of change of basis
+
+% Iterative selection of new points
+for n = 1 : maxIter
+    % select the current index
+    if strcmp(gType, 'fgreedy')
+        [fmax(n), i] = max(sum(f(notIndI, :) .^ 2, 2));
+        pmax(n) = max(p(notIndI));
+    elseif strcmp(gType, 'fpgreedy')
+        [~, i] = max(sum(f(notIndI, :) .^ 2, 2) ./ p(notIndI));
+        fmax(n) = max(sum(f(notIndI, :) .^ 2, 2));
+        pmax(n) = max(p(notIndI));
+    elseif strcmp(gType, 'pgreedy')
+        [pmax(n), i] = max(p(notIndI));
+        fmax(n) = max(sum(f(notIndI, :) .^ 2, 2));
+    else
+        error('Not a valid greedy method')
+    end
+    % add the current index
+    indI(n) = notIndI(i);
+    
+    if fmax(n) >= tolF && pmax(n) >= tolP
+        % compute the nth basis
+        Vx(notIndI, n) = Afun(notIndI, indI(n)) - Vx(notIndI, 1 : n - 1) * Vx(indI(n), 1 : n - 1)';
+        % normalize the nth basis
+        Vx(notIndI, n) = Vx(notIndI, n) / sqrt(p(indI(n)));
+        % update the change of basis
+        Cut(n, 1 : n) = [-Vx(indI(n), 1 : n - 1) * Cut(1 : n - 1, 1 : n - 1),  1] / Vx(indI(n), n);
+        % compute the nth coefficient
+        c(n, :) = f(indI(n), :) ./ sqrt(p(indI(n)));
+        % update the power function
+        p(notIndI) = p(notIndI) - Vx(notIndI, n) .^ 2;
+        % update the residual
+        f(notIndI, :) = f(notIndI, :) - bsxfun(@times, Vx(notIndI, n), c(n, :));
+        % remove the nth index from the dictionary
+        notIndI = notIndI([1 : i - 1, i + 1 : end]);
+    else
+        n = n - 1;
+        break
+    end
+end
+
+% Final resize to match the actual number of selected points
+c = c(1 : n, :);
+Cut = Cut(1 : n, 1 : n);
+a = Cut' * c;
+indI = indI(1 : n);
+fmax = sqrt(fmax(1 : n));
+pmax = sqrt(pmax(1 : n));
+
+end
\ No newline at end of file

Functions/plotPoints.m

+function plotPoints(X, color, titleString)
+% PLOTPOINTS Plots points in arbitrary dimension.
+%   PLOTPOINTS(X) plots the N d-dimensional points contained in the Nxd
+%   matrix X. If d>3, a PCA is applied to the matrix and only the first
+%   three dimensions are plotted.
+%   Color specifications and a title can be passed as inputs.
+
+if ~exist('color', 'var')
+    color = '.';
+end
+if ~exist('titleString', 'var')
+    titleString = [];
+end
+
+d = size(X, 2);
+
+if d == 1
+    plot(X, 0 * X, color, 'markersize', 12)
+elseif d == 2
+    plot(X(:, 1), X(:, 2), color, 'markersize', 12)
+elseif d == 3
+    plot3(X(:, 1), X(:, 2), X(:, 3), color, 'markersize', 12)
+else
+    [U, ~] = svd(X' * X);
+    X = X * U(:, 1 : 3);
+    plot3(X(:, 1), X(:, 2), X(:, 3), color, 'markersize', 12)
+end
+title(titleString)
+axis equal
+
+end
+

Functions/setData.m

+function dataset = setData(X, f, par)
+% dataset = SETDATA(X, f, par) creates a train/val/test dataset.
+
+
+% Decide sizes
+M = size(X, 1);
+Mte = floor(M * par.percTe);
+Mtr = M - Mte;
+
+% Permute data
+X = X(par.perm, :);
+f = f(par.perm, :);
+
+% Training set (including validation)
+Xtr = X(1 : Mtr, :);
+ftr = f(1 : Mtr, :);
+
+% Indexes for k-fold cross validation
+par.k = min(par.k, Mtr); % move to LOOCV is dataset is small
+MSubsample = floor(Mtr / par.k);
+indVal = zeros(par.k, 2);
+indVal(:, 1) = MSubsample * (0 : par.k - 1) + 1;
+indVal(:, 2) = MSubsample * ( 1 : par.k);
+indVal(end, 2) = Mtr;
+
+
+% Compute normalization variables
+if strcmp(par.normalization, 'uniform')
+    translXtr = min(Xtr);
+    scaleXtr =  max(Xtr) - min(Xtr);
+    translFtr = min(ftr);
+    scaleFtr = max(ftr) - min(ftr);
+elseif strcmp(par.normalization, 'none')
+    translXtr = 0;
+    scaleXtr = 1;
+    translFtr = 0;
+    scaleFtr = 1;
+else
+    error('Specify a valid normalization type')
+end
+
+% Avoid dividing by zero
+scaleXtr(scaleXtr == 0) = 1;
+scaleFtr(scaleFtr == 0) = 1;
+
+% Normalize
+Xtr = bsxfun(@rdivide, bsxfun(@minus, Xtr, translXtr), scaleXtr);
+ftr = bsxfun(@rdivide, bsxfun(@minus, ftr, translFtr), scaleFtr);
+
+
+% Set up output
+dataset.Xtr = Xtr;
+dataset.ftr = ftr;
+dataset.translXtr = translXtr;
+dataset.scaleXtr = scaleXtr;
+dataset.translFtr = translFtr;
+dataset.scaleFtr = scaleFtr;
+dataset.indVal = indVal;
+dataset.Xte = X(end - Mte + 1: end, :);
+dataset.fte = f(end - Mte + 1: end, :);
+
+end
\ No newline at end of file

Functions/trainModelReg.m

+function [modelPar, trainResults, valResults, testResults] = trainModelReg(dataset, valPar, trPar, verboseFlag)
+% TRAINMODELREG trains a kernel model using VKOGA, with optimization of the
+% kernel shape parameter and the regularizatio parameter via k-fold cross
+% validation.
+
+
+% Kernel
+ker = @(ep, x, y) trPar.rbf(ep, distanceMatrix(x, y));
+
+% Parameter optimization
+errFun = @(s, f) sqrt(sum((f - s) .^ 2, 2)); 
+if ~isfield(valPar, 'spaceType') || strcmp(valPar.spaceType, 'lin')
+    Eps = linspace(valPar.minEp, valPar.maxEp, valPar.nEp);
+    Lambda = linspace(valPar.minLambda, valPar.maxLambda, valPar.nLambda);
+elseif strcmp(valPar.spaceType, 'log')
+    Eps = logspace(log10(valPar.minEp), log10(valPar.maxEp), valPar.nEp);
+    Lambda = logspace(log10(valPar.minLambda), log10(valPar.maxLambda), valPar.nLambda);
+end
+
+[ee, ll] = ndgrid(Eps, Lambda);
+parameters = [ee(:) ll(:)];
+nPar = size(parameters, 1);
+
+errVal = zeros(nPar, 1);
+fprintf('Parameter validation:\n')
+
+parfor k = 1 : nPar
+        tmp_err = zeros(size(dataset.indVal, 1), 1);
+        for l = 1 : size(dataset.indVal, 1)
+            Xtr = dataset.Xtr(setdiff(1 : size(dataset.Xtr, 1), dataset.indVal(l, 1) : dataset.indVal(l, 2)), :);
+            ftr = dataset.ftr(setdiff(1 : size(dataset.Xtr, 1), dataset.indVal(l, 1) : dataset.indVal(l, 2)), :);
+            Xval = dataset.Xtr(dataset.indVal(l, 1) : dataset.indVal(l, 2), :);
+            fval = dataset.ftr(dataset.indVal(l, 1) : dataset.indVal(l, 2), :);
+            Afun = @(i, j) ker(parameters(k, 1), Xtr(i, :), Xtr(j, :)) + parameters(k, 2) * double(bsxfun(@eq, i(:), j(:)'));
+            [a, indX, n] = newtonYYDataMF(Afun, ftr, valPar.tolFEp, valPar.tolPEp, valPar.maxExpSizeEp, trPar.gType);
+            s = ker(parameters(k, 1), Xval, Xtr(indX, :)) * a;
+            tmp_err(l) = max(errFun(s, fval));
+        end
+        errVal(k) = mean(tmp_err);
+        if verboseFlag
+            fprintf('ep = %2.2e, lambda = %2.2e, error = %2.2e, n = %2d\n', parameters(k, 1), parameters(k, 2), errVal(k), n)
+        end
+end
+
+[~, minErrParameter] = min(errVal);
+epOpt = parameters(minErrParameter, 1);
+lambdaOpt = parameters(minErrParameter, 2);
+fprintf('Selected paramers: ep = %2.2e, lambda = %2.2e\n', epOpt, lambdaOpt)
+
+% Model training
+Afun = @(i, j) ker(epOpt, dataset.Xtr(i, :), dataset.Xtr(j, :)) + lambdaOpt * double(bsxfun(@eq, i(:), j(:)'));
+[a, indX, n, c, Cut, fmax, pmax] = newtonYYDataMF(Afun, dataset.ftr, trPar.tolF, trPar.tolP, trPar.maxExpSize, trPar.gType);
+
+% Define model parameters
+modelPar.ep = epOpt;
+modelPar.rbf = func2str(trPar.rbf);
+modelPar.translX = dataset.translXtr;
+modelPar.scaleX = dataset.scaleXtr;
+modelPar.translF = dataset.translFtr;
+modelPar.scaleF = dataset.scaleFtr;
+
+% Test error computation
+fprintf('Test error computation:\n')
+if size(dataset.Xte, 1)
+    errFun = @(s) sqrt(sum((dataset.fte - s) .^ 2, 2)); % ./ sqrt(sum((dataset.fte) .^ 2, 2)); 
+    
+    meanErr = zeros(n, 1);
+    maxErr = zeros(n, 1);
+    for k = 1 : n
+        ak = Cut(1 : k, 1 : k)' * c(1 : k, :);
+        modelPar.X = dataset.Xtr(indX(1 : k), :);
+        modelPar.a = ak;
+        s = kernelModel(dataset.Xte, modelPar);
+        pointwiseErr =  errFun(s);
+        meanErr(k) = mean(pointwiseErr);
+        maxErr(k) = max(pointwiseErr);
+    end
+else
+    meanErr = [];
+    maxErr = [];
+end
+
+% Update model parameters
+modelPar.X = dataset.Xtr(indX, :);
+modelPar.a = a;
+
+% Set up outputs
+trainResults.fmax = fmax;
+trainResults.pmax = pmax;
+
+valResults.ep = modelPar.ep;
+valResults.lambda = lambdaOpt;
+valResults.Eps = Eps;
+valResults.Lambda = Lambda;
+valResults.errVal = errVal;
+
+testResults.meanErr = meanErr;
+testResults.maxErr = maxErr;
+
+end
+
+

startup.m

+%% Set up before running VKOGA
+
+% Add relevant paths
+addpath('Functions/')
+
+% Start parpool
+if isempty(gcp('nocreate'))
+    gcp();
+end
\ No newline at end of file