function [sig,C,alpha,b] = ams(X,Y,method,sphere,opt)%   AMS  Automatic Model Selection for Kernel Methods with%            - RBF Gaussian kernel%            - L2 penalization of the errors%          For classification: Support Vector Machines%          For regression: Kernel Ridge Regression%%  Usage: [sig,C,alpha,b] = AMS(X,Y,method,[sphere],[opt])%%    X: is an n times d matrix containing n training points in d dimensions  %    Y: is a n dimensional vector containing either the class labels%       (-1 or 1) in classification or the targets values in regression.  %    method is either:%       'loo': minimizes the leave-one-out error. The LOO is an accurate%              estimate of the generalization error but might be difficult%              to minimize (not smooth). Also, do not use with if sphere = 0%              and large d (long time for gradient calculations).%       'rm':  radius / margin bound. Easier to minimize. That can be the%              first thing to try%       'val': minimizes a validation error. Good if there are a lot%              a training points. The validation points are the%              nv points of the training set [nv specified in the options]%       'ev':  evidence maximization.%    sphere: if 1, spherical RBF kernel = exp(-0.5*||x-y||^2 / sigma^2)%            if 0, one sigma per component:%               exp(-0.5*sum (x_i-y_i)^2/sigma_i^2)%            [default 1]%    opt: structure containing the options (in brackets default value)%       itermax: maximum number of conjugate gradient steps [100]    %       tolX: gradient norm stopping criterion [1e-5]%       hregul: regularization on the hyperparameters (to avoid extreme values) [0.01]%       nonorm: if 1, do no not normalize the inputs [0]%       verb: verbosity [1]%       eta: smoothing parameter used in classif by 'ev' and 'loo'%            larger is smoother [0.1]%       sigmoid: slope of the sigmoid used by 'val' and 'loo' in classif [5] %       nv: number of validation points for 'val' [n/2]%%%    sig: the sigma(s) found by the model selection. Scalar if %         sphere = 1, d dimensional vector otherwise%    C: the constant penalizing the training errors%    alpha,b: the parameters of the prediction function%      f(x) = sum alpha_i K(x,x_i) + b  %    Copyright Olivier Chapelle%              olivier.chapelle@tuebingen.mpg.de%              last modified 21.04.2005  % Initialization  [n,d] = size(X);  if (size(Y,1)~=n) | (size(Y,2)~=1)    error('Dimension error');  end;      if all(abs(Y)==1)    classif = 1;      % We are in classification  else    classif = 0;      % We are in regression  end;    if nargin < 5       % Assign the options to their default values    opt = [];  end;  if ~isfield(opt,'itermax'), opt.itermax = 100;                end;  if ~isfield(opt,'length'),  opt.length  = opt.itermax;        end;  if ~isfield(opt,'tolX'),    opt.tolX    = 1e-5;               end;  if ~isfield(opt,'hregul'),  opt.hregul  = 1e-2;               end;  if ~isfield(opt,'nonorm'),  opt.nonorm  = 0;                  end;  if ~isfield(opt,'verb'),    opt.verb    = 0;                  end;  if ~isfield(opt,'eta'),     opt.eta     = 0.1;                end;  if ~isfield(opt,'sigmoid'), opt.sigmoid = 5;                  end;  if ~isfield(opt,'nv'),      opt.nv      = round(n/2);         end;      if nargin<4    sphere = 1;  end;    % Normalization  if ~opt.nonorm & (norm(std(X)-1) > 1e-5*d)    X = X./repmat(std(X),size(X,1),1);    fprintf(['The data have been normalized (each component is divided by ' ...             'its standard deviation).\n Don''t forget to do the same for the ' ...             'test set.\n']);  end;    % Default values of the hyperparameters (in log scale)  C = 0;  sig = .3*var(X)';  if ~sphere,    sig = log(d*sig) / 2;  else    sig = log(sum(sig)) / 2;  end;    % Do the optimization on the hyperparameters  param = [sig; C];  default_param = param;%  check_all_grad(param,X,Y,classif,default_param,opt); return;%  plot_criterion(X,Y,param,[-3:0.1:3],classif,default_param,opt);%  return;  [param,fX] = minimize(param,@obj_fun,opt,X,Y,method,classif,default_param,opt);    % Prepare the outputs  [alpha, b] = learn(X,Y,param,classif);  C = exp(param(end));  sig = exp(param(1:end-1));function plot_criterion(X,Y,param,range,classif,opt)  for i=1:length(range)    fprintf('%f\r',range(i));    for j=1:2      param2 = param;      param2(end-j+1) = param2(end-j+1)+range(i);      obj(1,i,j) = obj_fun(param2,X,Y,'loo',classif,opt);      obj(2,i,j) = obj_fun(param2,X,Y,'rm', classif,opt);      obj(3,i,j) = obj_fun(param2,X,Y,'val',classif,opt);      obj(4,i,j) = exp(obj_fun(param2,X,Y,'ev', classif,opt));    end;  end;  figure; plot(obj(:,:,1)');  figure; plot(obj(:,:,2)');function  check_all_grad(param,X,Y,classif,opt)  param = randn(size(param));  checkgrad(@obj_fun,param,1e-6,X,Y,'loo',classif,opt)  checkgrad(@obj_fun,param,1e-6,X,Y,'rm', classif,opt)  checkgrad(@obj_fun,param,1e-6,X,Y,'val',classif,opt)  checkgrad(@obj_fun,param,1e-6,X,Y,'ev', classif,opt)  function [obj, grad] = obj_fun(param,X,Y,meth,classif,default_param,opt)% Compute the model selection criterion and its derivatives with% respect to param which is a vector containing in log scale sigma% and C  % Add a regularization on the hyperparameters to avoid that  % they take extreme values  obj0 = opt.hregul * mean((param-default_param).^2);  grad0 = 2*opt.hregul * (param-default_param) / length(param);  if obj0 > 1               % Extreme value of hyperparamaters. The rest    obj = obj0;             % might not even work because of numercial    grad = grad0;           % instabilities -> exit    return;  end      if strcmp(meth,'val')    % Remove the validation set for the learning    [alpha,b] = learn(X(1:end-opt.nv,:),Y(1:end-opt.nv),param,classif);  else    [alpha,b] = learn(X,Y,param,classif);  end;    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  % Compute the objective function %  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    sv = find(alpha~=0);            % Find the set of support vectors    switch lower(meth)    %    % Leave-one-out     %   case 'loo'    K = compute_kernel(X(sv,:),[],param);    K(end+1,:) = 1;                 % To take the bias into account    K(:,end+1) = 1;    K(end,end) = 0;    D = opt.eta./abs(alpha(sv));    invKD = inv(K+diag([D;0]));    invKD(:,end) = []; invKD(end,:) = [];    span = 1./diag(invKD) - D;      % "Span" of the support vectors    % Estimated change of the output of each point when it's    % removed from the training set.    change = alpha(sv).*span;    if classif      out = 1-Y(sv).*change;      obj = sum(tanh(-opt.sigmoid*out)+1)/2/length(Y);    else      obj = mean(change.^2);    end;    invK = inv(K); invK(:,end) = []; invK(end,:) = [];        %    % Radius margin    %   case 'rm'    K = compute_kernel(X,[],param);    w2 = sum(alpha.*Y);    if classif      R2 = mean(diag(K)) - mean(mean(K)); % Radius estimated by variance    else      K(end+1,:) = 1;                 % Radius = mean span      K(:,end+1) = 1;      K(end,end) = 0;      invK = inv(K);      span = 1./diag(invK);              R2 = mean(span(1:end-1));    end;    obj = R2*w2 / length(Y);        %    % Validation set    %   case 'val'    K = compute_kernel(X(1:end-opt.nv,:),X(end-opt.nv+1:end,:),param);    Yv = Y(end-opt.nv+1:end);    out = K(end-opt.nv+1:end,:)*alpha+b;   % Output on the validation points    if classif      % Goes through sigmoid because step function is not differentiable      obj = mean(tanh(-opt.sigmoid*out.*Yv)+1)/2;    else      obj = mean((out-Yv).^2);    end;    Kl = K(sv,sv);    Kl(end+1,:) = 1;                   Kl(:,end+1) = 1;    Kl(end,end) = 0;    invK = inv(Kl);         %    % Evidence maximization    %   case 'ev'    K = compute_kernel(X,[],param);    if classif      % Normally, K should be computed only on the support vectors. But      % then, the evidence wouldn't be continuous. We thus make a soft      % transition for small alphas.      weight = min(0, (alpha.*Y)/(mean(alpha.*Y))/opt.eta - 1);      weight2 = 1 - weight.^2;      K0 = K;      K = K.*(max(weight2*weight2',eye(length(weight2))));          Kl = K0(sv,sv);                          Kl(end+1,:) = 1;                     Kl(:,end+1) = 1;      Kl(end,end) = 0;      invKl = inv(Kl); invKl(:,end) = []; invKl(end,:) = [];    end;    [invK, logdetK] = inv_logdet_pd_(K);    obj = logdetK/length(K) + log(mean(alpha.*Y));           otherwise    error('Unknown model selection criterion');  end;    %%%%%%%%%%%%%%%%%%%%%%%%%  % Compute the gradients %  %%%%%%%%%%%%%%%%%%%%%%%%%    for i=1:length(param)        switch lower(meth)      %      % Leave-one-out       %     case 'loo'      dK = compute_kernel(X(sv,:),[],param,i);      dalpha = -invK * (dK*alpha(sv));      dD = - opt.eta * dalpha.*sign(alpha(sv))./alpha(sv).^2;      dspan = diag(invKD).^(-2).* ...              diag(invKD*(dK+diag(dD))*invKD) - dD;      dchange = dalpha.*span + alpha(sv).*dspan;      if classif        dout = -Y(sv).*dchange;        grad(i) = -0.5*opt.sigmoid * ...                  sum(cosh(-opt.sigmoid*out).^(-2).*dout)/length(Y);      else        grad(i) = 2*mean(change.*dchange);      end;          %      % Radius margin      %     case 'rm'      dK = compute_kernel(X,[],param,i);      dw2 = -alpha'*dK*alpha;      if classif        dR2 = mean(diag(dK)) - mean(mean(dK));      else        dspan = diag(invK).^(-2).*diag(invK(:,1:end-1)*dK*invK(1:end-1,:));        dR2 = mean(dspan(1:end-1));      end;      grad(i) = (dR2*w2 + R2*dw2) / length(Y);          %      % Validation set      %     case 'val'      dK = compute_kernel(X(1:end-opt.nv,:),X(end-opt.nv+1:end,:),param,i);      dalpha = -invK(:,1:end-1) * (dK(sv,:)*alpha);      db = dalpha(end);      dalpha = dalpha(1:end-1);      dout =  K(end-opt.nv+1:end,sv)*dalpha + db + ...             dK(end-opt.nv+1:end,:) * alpha;       if classif        grad(i) = -0.5*opt.sigmoid * ...                  mean(cosh(-opt.sigmoid*out.*Yv).^(-2).*dout.*Yv);      else        grad(i) = 2*mean((out-Yv).*dout);      end;          %      % Evidence maximization      %     case 'ev'      dK = compute_kernel(X,[],param,i);      dK0 = dK;      if classif         dalpha = zeros(size(alpha));        dalpha(sv) = -invKl * (dK(sv,:)*alpha);        dweight = -2/opt.eta * weight.* ...            (dalpha.*Y/mean(alpha.*Y) - (alpha.*Y) * mean(dalpha.*Y)/mean(alpha.*Y)^2);        dK = dK .* max(weight2*weight2',eye(length(weight))) + ...             (K0-diag(diag(K0))) .* (weight2*dweight' + dweight*weight2');        end;      grad(i) = invK(:)'*dK(:)/length(K) - ...                alpha'*dK0*alpha / sum(alpha.*Y);    end;  end;  obj = obj0 + obj;  grad = grad0 + grad';      function [alpha,b] = learn(X,Y,param,classif)% Do the actual learning (SVM or kernel ridge regression)  K = compute_kernel(X,[],param);  svset = 1:length(Y);  old_svset = [];  iter = 0;  itermax = 20;    % If the set of support vectors has changed, we need to  % reiterate. Note that for regression, all points are support  % vectors and there is only one iteration.  while ~isempty(setxor(svset,old_svset)) & (iter<itermax)    old_svset = svset;    H = K(svset,svset);    H(end+1,:) = 1;                 % To take the bias into account    H(:,end+1) = 1;    H(end,end) = 0;    % Find the parameters    par = H\[Y(svset);0];    alpha = zeros(length(Y),1);    alpha(svset) = par(1:end-1)';    b = par(end);    % Compute the new set of support vectors    if classif      % Compute the outputs by removing the ridge in K      out = K*alpha+b - alpha/exp(param(end));       svset = find(Y.*out < 1);    end;    iter = iter + 1;  end;  if iter == itermax    warning('Maximum number of Newton steps reached');  elseif classif & any(alpha.*Y<0)    warning('Bug in learn');  end;    function K = compute_kernel(X,Y,param,diff)% Compute the RBF kernel matrix between  [X;Y] and X.% if diff is given, compute the derivate of the kernel matrix%  with respect to the diff-th parameter of param.      global K0   % Keep the kernel matrix in memory; we'll need to              % compute the derivatives.    Z = [X; Y];  Y = X; X = Z; clear Z;  ridge = exp(-param(end));    if nargin<4    % Compute the kernel matrix        % Divide by the length scales    if length(param)==2      X = X * exp(-param(1));      Y = Y * exp(-param(1));    else      X = X .* repmat(exp(-param(1:end-1)'),size(X,1),1);      Y = Y .* repmat(exp(-param(1:end-1)'),size(Y,1),1);    end;        normX = sum(X.^2,2);    normY = sum(Y.^2,2);        K0 = exp(-0.5*(repmat(normX ,1,size(Y,1)) + ...                  repmat(normY',size(X,1),1) - 2*X*Y'));    K = K0 + eye(size(K0))*ridge;  else  % Compute the derivate    if diff == length(param)    % With respect to log(C)      K = - ridge*eye(size(K0));    elseif length(param) == 2   % With respect to spherical log(sig)      K = -2*K0.*log(K0+eps);    else                        % With respect to log(sig_i)      dist = repmat(X(:,diff),1,size(Y,1)) - ...             repmat(Y(:,diff)',size(X,1),1);      K = K0 .* (dist/exp(param(diff))).^2;    end;  end;function [invA, logdetA] = inv_logdet_pd_(A)  % For those who don't have the inv_logdet_pd written by Carl  if exist('inv_logdet_pd')    [invA, logdetA] = inv_logdet_pd(A);  else    invA = inv(A);    C = chol(A);    logdetA = 2*sum(log(diag(C)));  end;