function [w,elg,plg,blg] = hecgd(w, D, C, be, W, opt)
% HECGD Hyper-Ellipsoidal Conjugate Gradient Descent
% 
%  Use this function to minimize a differentiable cost function subject
%  to one positive definite quadratic constraint and none or a number of
%  linear orthogonality constraints. The function is particularly suited
%  to perform sparse optimization for second order constrained
%  optimization problems
%
%  [w,elg,plg] = HECGD(w, C, D, be, W)
%
%  minimizes 
%       w'*D*w + be*sum(abs(w))^2/(w'*w)
%  s.t.
%       w'*C*w ==1
%       w'*C*W ==0
%
%  where 
%    w   is an Nx1 vector of initial weights
%    C   is a square, symmetric, and positive definite matrix, and 
%    D   is an arbitrary square matrix
%    W   is an NxM matrix with M<N  (otherwise there would be no
%        solution), or [] if there are no orthogonality constraints
%    be  >=0 weigths the strength of the regularization term 
%
%  use
%   [w,elg,plg] = HECGD(w, C, D, be, W, opt)
%
%   to pass a structure of user specified parameters. opt can have the
%   following fields:
%
%    field       admissible values   default  meaning
%    --------------------------------------------------------------------------
%    'Nit',        integer >0      100       maximum number of iterarions
%
%    'Nit_inner'   integer >0      500       maximum number of line search iterarions
%
%    'Nit_restart' integer >=0     0         conjugate gradient reset period (0=never)
%
%    'reduction'   boolean         false     clamp elements of w to zero that
%                                            come closer to zero than
%                                            tol(3). 
%    'tol'         1x3 
%                  double>=0       [1e-4 1e-4 1e-3] termination tolerance
%    values 
%                                            The optimization terminates when
%                                            the changes in w are below
%                                            tol(1) and the changes  in the
%                                            cost function are below
%                                            tol(2). The maximum number of
%                                            iterations is Nit.  If a weight
%                                            falls below tol(3) it is
%                                            discarded if reduction is activated
%
%
%    'verbose'     boolean         true      print status information
%
%    'symmetry'    boolean         true      specifies whether the error
%                                            function is symmetric on the
%                                            line search ellipse. This
%                                            affects the maximum search
%                                            interval, which is either 2*pi
%                                            for non-symmetric functions or
%                                            otherwise pi
%
%    'errfun'    function handle   @errsokm  user specified unconstrained
%                                            error/gradient function. See
%                                            errsokm for description of the interface
%
%
%  Example: kernel-PCA
%  -------------------
%  compute 10 sparse kernel principal components from an NxT data set X
% 
%  1. compute the TxT kernel matrix between all pairs of data points
%  2. W=[];
%     D = -K*K';
%     C = K;
%     for i=1:10
%       w = hecdg(randn(T,1), D, C, be, W);
%       W = [W w];
%     end
%
%
%   the algorithm is described in length in:
%
%   [1] "R. Vollgraf and K. Obermayer, Sparse Optimization for Second Order
%   Kernel Methods, in Proc. of 2006 IJCNN, Vancouver, 2006"
%
%
%   See also errsokm.

% Copyright (C) 2006  TU-Berlin, Roland Vollgraf
% 
% This program is free software; you can redistribute it and/or
% modify it under the terms of the GNU General Public License
% as published by the Free Software Foundation; either version 2
% of the License, or (at your option) any later version.
% 
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
% 
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
% 
% If you use this software, please cite 
% "R. Vollgraf and K. Obermayer, Sparse Optimization for Second Order Kernel
% Methods, in Proc. of 2006 IJCNN, Vancouver, 2006"
% in your publications




gold = 0.38197;



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% global initializations %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% default parameter values %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Nit         = 100;                       % maximum number of iterarions
Nit_inner   = 500;                       % maximum number of line search iterarions
Nit_restart = 0;                         % conjugate gradient reset
                                         % period (0=never)

reduction   = false;                     % clamp elements of w close to zero
tol         = [1e-4 1e-4 1e-3];          % termination tolerance values 
verbose     = true;                      % print status information

interval    = pi;                        % line search interval, pi or 2*pi

errfun = @(w,D,be) errsokm(w,D,be);      % unconstrained error/gradient function


% User specified parameter values %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if nargin >= 6
  for field = fieldnames(opt)'
    switch(char(field))
      case 'Nit'           
        Nit         = opt.Nit;
      case 'Nit_restart'   
        Nit_restart = opt.Nit_restart;
      case 'Nit_inner'     
        Nit_inner   = opt.Nit_inner;        
      case 'reduction'     
        reduction   = opt.reduction;
      case 'tol'           
        tol         = opt.tol;
      case 'verbose'       
        verbose     = opt.verbose;
      case 'symmetry'
        if opt.symmetry
          interval = pi;
        else
          interval = 2*pi;
        end
      case 'errfun'        
        errfun = opt.errfun;
      
      otherwise
        error(sprintf('Unknown option ''%s''',char(field)));
    end
  end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


N=length(w);


term_fw = 0; % termination criterion counter
term_w  = 0; % termination criterion counter


elg = zeros(Nit,1); % log of the error function values
plg = zeros(Nit,N); % log of the parameters
                    


% sphere the space that is spanned by the linear constraints
if numel(W)>0
  CW = C*W;
  CW = CW*(CW'*CW)^(-1/2);
end




% some initializations for the conjugate gradient
dw   = zeros(N,1); % initial value for the old search direction
dwdw = 1;          % initial value for its squared length



iw = 1:N;  % indices to the non-zero elements of w if reduction is used 


if verbose
    fprintf(['It      fval     fval(unc)  It(GSS)  delta_w,     delta_fw\n'...
             '----------------------------------------------------------\n'])
end


% the main loop %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for n=1:Nit

    % reduction %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % elements of w that are closer to zeros that tol(3) become clamped
    % to zero, the corresponding rows and columns of C and D are discarded
    if reduction
        i = find(abs(w/sqrt(w'*w))>tol(3));
  
        if length(i)<length(w) && length(i)>2

            iw = iw(i);
            w  = w(i);
            C  = C(i,i);
            D  = D(i,i);
    
            if numel(W)>0

                CW = CW(i,:);

                % we discard small components of CW to avoid numerical
                % instabilities
                [v,d] = eig(CW'*CW);
                j = find(diag(d)>1e-8);
                CW = CW*v(:,j)*d(j,j)^(-1/2);
            end

            
            % dw is the old search direction, 
            % reset the conjugate gradient procedure after reduction
            dw = 0*dw(i);
        end
    end

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  
  % make w orthogonal to W 
  if numel(W)>0
      w = w- CW*pinv(CW)*w;
  end
  wCw = w'*C*w;
  if wCw <=0
      %still a result of too small dw in the last iteration ?  
      %should not happen anymore because of (a)
      warning('\nBreak due to zero gradient (b)')
      wnew =w;
      break;  
  end
  
  % This should not happen anymore because of the sphering of CW
  if abs(wCw-1)>1e-3 && n>1
      warning('Numerical instabilities!')
  end
  w = w/sqrt(wCw);

  
  
  % restart the conjugate gradient all Nit_restore iterations 
  % or when the last gradient was too small
  if mod(n, Nit_restart)==0 || dwdw==0
      gamma = 0;
  else
      gamma = 1/dwdw*dw;
  end
  
  %dw = gra(w,D,be);
  %fw = err(w,D,be);
  
  [fw,dw] = errfun( w, D, be);
  
  elg(n)  = fw;
  plg(n,iw)= w(:)';
 
  
  
  % make the new gradient, dw, orthogonal to W and w
  if numel(W)>0
    CWw = [CW, C*w];
  else
    CWw = [C*w];
  end
  dw = dw - CWw*pinv(CWw) * dw;
  
  dwdw = dw'*C*dw;  % we need this value for the computation of gamma in
                    % the next iteration 
  
  if dw'*C*dw <1e-20
    warning('\nBreak due to to small gradient')
    wnew=w;
    break;
  end
  
  

  
  % dw becomes the new search direction
  dw = dw + dwdw*gamma;
  

  dwCdw = dw'*C*dw;     % we must remember the length of the new search
                        % direction for the transport into the new
                        % tangent space below
  if dwCdw <=0
    % obviously dw was too small
    fprintf('\nBreak due to zero gradient (a)')
    wnew =w;
    break;
  end
  dw = dw/sqrt(dwCdw);
  
  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Golden section search %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  % setup the initial interval %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  a  = 0;
  fa = errfun(w,D,be);

  b  = gold*pi;
  fb = errfun( cos(b)*w - sin(b)*dw, D,be);

  c  = pi;
  fc = errfun( cos(c)*w - sin(c)*dw, D,be);

  while fb>fa || fb>fc
    c  = b;
    fc = fb;
  
    b  = gold*(c-a);
    fb = errfun( cos(b)*w - sin(b)*dw, D,be);
  end
  
  
  % start iterations %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  i=1;
  while (c-a)>tol(1) || fa-fb>tol(2) || fc-fb>tol(2) 

    if (b-a)>(c-b) % split the first half of the interval
      d = a+gold*(c-a);
      fd = errfun( cos(d)*w-sin(d)*dw, D,be);
      if fd<fb
        c=b; fc=fd;
        b=d; fb=fd;
      else
        a=d; fa=fd;
      end
    
    else % split the second half
      d = c-gold*(c-a);
      fd = errfun( cos(d)*w-sin(d)*dw, D,be);
      if fd<fb
        a=b; fa=fd;
        b=d; fb=fd;
      else
        c=d; fc=fd;
      end
    end
    
    if fb>fa || fb>fc % this really should not happen
      warning('Bracket was lost!')
    end
    
    i=i+1;
    
    if i>Nit_inner
      warning('Maximum number of iterations exceeded!')
      break
    end
    
  end

  wnew  = cos(b)*w-sin(b)*dw;
  fwnew = fb;

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % the line search was finished %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  
  % stopping criterion for fw %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  delta_fw = abs(fwnew-fw);
  if delta_fw < tol(2)
    term_fw = term_fw+1;
  else
    term_fw = 0;
  end
  
  % stopping criterion for w %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  delta_w = sqrt((w-wnew)'*C*(w-wnew));  
  if delta_w < tol(1)
    term_w = term_w+1;
  else
    term_w = 0;
  end
  
  
  if  term_w>=5 && term_fw>=5
    w  = wnew;
    fw = fwnew;
    break;
  end
  
  
  % transport of the search direction into the new tangent space %%%%%%%%%%%%%%%%
  dw = (cos(b)*dw+sin(b)*w)*sqrt(dwCdw); % it is assigned the saved
                                         % length dwCdw

  w  = wnew;
  fw = fwnew;

  
  if verbose
      fprintf('It %3d: %f, %f, %4d,   %f, %f\n',...
              n, fw, errfun(w,D,0), i,delta_w, delta_fw);
  end
  
end
fprintf('\n');

w     = zeros(N,1);
w(iw) = wnew;


elg(n+1)  = fw;
plg(n+1,:)= w';

elg = elg(1:n+1);
plg = plg(1:n+1,:);


