integrated_newton.m

一个matlab的工具包,里面包括一些分类器 例如 KNN KMEAN SVM NETLAB 等等有很多.

% This funciton provides an integrated Newton's Method for IIS,
% MIS, etc.
% function [root, fval, bound] = integrated_newton(method, mycoeff, mypow, a0, maxiter)
%
% method:	'IIS': 		performs the newton's method for IIS.
%		'MIS': 		performs the newton's method for MIS.
%		'IIS_ONCE':	performs the newton approximation for IIS
%		'IIS_ONCE_DEBUG':	 with debug info.
%		'MIS_ONCE':	performs the newton approximation for MIS
%		'MIS_ONCE_DEBUG':	 with debug info.
%
% NOTE: newton approximation means only apply one iteration of
% newton's method.
%
% mycoeff:	the coefficients of the polynomial function,
% 		excluding the const.
% mypow:	the powers of the polynomial function, excluding
% 		the "0" power.
% a0:		the const of the polynomial function.
% maxiter:	the maximum iteration of the newton's method, only
% 		applicable to 'IIS' and 'MIS'.

% root:		the root of the function.
% fval:		the function value at root.
% bound:	used for the computation of newton approximation


function [root, fval, bound] = integrated_newton(method, mycoeff, mypow, a0, maxiter)

% set the default value
bound = 0.0;

% method IIS_ONCE:
if (strcmp(method, 'IIS_ONCE') == 1)
  root = 1.0;
  
  % caculate the function value at root;	
  fval = - (sum(mycoeff) + a0);
  % caculate the derivative function value at root;	
  dfval = - sum(mycoeff .* mypow);
  root = root - fval / dfval;

  return;
end

% method IIS_ONCE_DEBUG:
if (strcmp(method, 'IIS_ONCE_DEBUG') == 1)
  root = 1.0;
  
  % caculate the function value at root;	
  fval = - (sum(mycoeff) + a0);
  % caculate the derivative function value at root;	
  dfval = - sum(mycoeff .* mypow);
  root = root - fval / dfval;
  
  % compute the bound
  B = - sum((mycoeff .* mypow) * (power(root, (mypow-1))'));
  bound = B / dfval;
%  fprintf('root = %f, bound = %f (B=%f, dfval=%f)\n', root, bound, B, dfval);
  
  bound = (0.5 * bound - 1.0) * fval * fval / dfval;

%  if (bound < 0.0)
%    fprintf('fval = %f, dfval = %f\n', fval, dfval);
%    fprintf('The bound should greater than 0: %f\n', bound);
%  end
  return;
end


% method IIS:
if (strcmp(method, 'IIS') == 1)
  % initial function and derivative function's value
  fval = 100.0;
  dfval = 0.0;
  dxold = 1000000;
  dx = 1000000;

  % high & low bound of the root 
  highbound = 1000000;
  lowbound = 0;

  epsilon = 1.0e-8; % stop criteria
  root = 1.0001; % initial guess of the root

  
  % main loop of newton's method 
  for newton_iter = 1 : maxiter
    % caculate the function value at root;	
    fval = sum(mycoeff * (power(root, mypow)')) + a0;
    % caculate the derivative function value at root;	
    dfval = sum((mycoeff .* mypow) * (power(root, (mypow-1))'));
    
    if (fval < 0)
      lowbound = root;
    else
      highbound = root;
    end
    dxold = dx;

    % if the fval is infinity, bisect the region, otherwise use newton still.
    if (fval >= realmax) | (2 * abs(fval) > abs(dxold * dfval)) | (fval ...
						  ~= fval) | (dfval ...
						  ~= dfval)
      dx = (highbound - lowbound) / 2;
      root = lowbound + dx;
    else
      if (abs(fval) <= eps) & (abs(dfval) <= eps)
	dx = 0.0;	
      elseif (abs(dfval) <= eps)
	dx = 0.0;
	fprintf('Divided by zero!\n');
      else
	dx = fval / dfval;
      end
      root = root - dx;

      % if newton thinks to go outside the region, bisect the region
      if (root < lowbound) | (root > highbound)
	dx = (highbound - lowbound) / 2.0;
	root = lowbound + dx;
	%      fprintf('low = %f, high = %f, come here, dx = %f\n', lowbound, highbound, dx);
      end
    end
    
    % check converge
    if abs(fval) < epsilon
      break;
    end
  end
  if abs(fval) > epsilon
    fprintf('not converged, fval = %f, root = %f\n', fval, root);
  end
  
  return;
end






% method MIS_ONCE:
if (strcmp(method, 'MIS_ONCE') == 1)

  % rescale the pow
  mycoeff = [mycoeff, a0];
  mypow = [mypow, 0];

  maxpow = max(mypow);
  mypow = mypow / (maxpow/1000.0);

  % some parameters
  % old values
  dxold = 1000000;
  dx = 1000000;

  % initial function and derivative function's value
  fval = 100.0;
  dfval = 0.0;

  % high & low bound of the root 
  highbound = 1000000;
  lowbound = 0;

  epsilon = 1.0e-8; % stop criteria
  root = 1.0000; % initial guess of the root, in fact, for IIS, the

  % First loop of newton's method 
    % caculate the function value at root;	
    fval = sum(mycoeff);
    % caculate the derivative function value at root;	
    dfval = sum(mycoeff .* mypow);
    
    if (fval < 0)
      lowbound = root;
    else
      highbound = root;
    end
    dxold = dx;  

    % if the fval is infinity, bisect the region, otherwise use newton still.
    if (fval >= realmax) | (2 * abs(fval) > abs(dxold * dfval)) | (fval ...
						  ~= fval) | (dfval ...
						  ~= dfval)
      dx = (highbound - lowbound) / 2;
      root = lowbound + dx;
    else
      if (abs(fval) <= eps) & (abs(dfval) <= eps)
			dx = 0.0;	
      elseif (abs(dfval) <= eps)
			dx = 0.0;	
			fprintf('Divided by zero!\n');
      else
			dx = fval / dfval;
      end
      root = root - dx;

      % if newton thinks to go outside the region, bisect the region
      if (root < lowbound) | (root > highbound)
	dx = (highbound - lowbound) / 2.0;
	root = lowbound + dx;
      end
    end
    
    % root = log(root) / maxpow * 1000.0;  
    root = - dx /maxpow * 1000.0;
    return;
end






% method MIS_ONCE_DEBUG:
if (strcmp(method, 'MIS_ONCE_DEBUG') == 1)

  % rescale the pow
  mycoeff = [mycoeff, a0];
  mypow = [mypow, 0];

  maxpow = max(mypow);
  mypow = mypow / (maxpow/1000.0);

  % some parameters
  % old values
  dxold = 1000000;
  dx = 1000000;

  % initial function and derivative function's value
  fval = 100.0;
  dfval = 0.0;

  % high & low bound of the root 
  highbound = 1000000;
  lowbound = 0;

  epsilon = 1.0e-8; % stop criteria
  root = 1.0000; % initial guess of the root, in fact, for IIS, the

  % First loop of newton's method 
    % caculate the function value at root;	
    fval = sum(mycoeff);
    % caculate the derivative function value at root;	
    dfval = sum(mycoeff .* mypow);
    
    if (fval < 0)
      lowbound = root;
    else
      highbound = root;
    end
    dxold = dx;  

    % if the fval is infinity, bisect the region, otherwise use newton still.
    if (fval >= realmax) | (2 * abs(fval) > abs(dxold * dfval)) | (fval ...
						  ~= fval) | (dfval ...
						  ~= dfval)
      dx = (highbound - lowbound) / 2;
      root = lowbound + dx;
    else
      if (abs(fval) <= eps) & (abs(dfval) <= eps)
	dx = 0.0;	
      elseif (abs(dfval) <= eps)
	dx = 0.0;	
	fprintf('Divided by zero!\n');
      else
	dx = fval / dfval;
      end
      root = root - dx;

      % if newton thinks to go outside the region, bisect the region
      if (root < lowbound) | (root > highbound)
	dx = (highbound - lowbound) / 2.0;
	root = lowbound + dx;
      end
    end
    
  % compute the bound
  B = - sum((mycoeff .* mypow) * (power(root, (mypow-1))'));
  
  if (dfval ~=0.0)
    bound = (0.5 * B/dfval - 1.0) * fval * fval / dfval;
    if (bound < 0.0)
%      fprintf('fval = %f, dfval = %f\n', fval, dfval);
%      fprintf('The bound should greater than 0: %f\n', bound);
    else
      root = log(root) / maxpow * 1000.0;  
%      disp('return from once newton approximation\n');
      return;
    end
  end
end


% method MIS:
if (strcmp(method, 'MIS') == 1)

  % rescale the pow
  mycoeff = [mycoeff, a0];
  mypow = [mypow, 0];

  maxpow = max(mypow);
  mypow = mypow / (maxpow/1000.0);

  % some parameters
  % old values
  dxold = 1000000;
  dx = 1000000;

  % initial function and derivative function's value
  fval = 100.0;
  dfval = 0.0;

  % high & low bound of the root 
  highbound = 1000000;
  lowbound = 0;

  epsilon = 1.0e-6; % stop criteria
  root = 1.0001; % initial guess of the root, in fact, for IIS, the

  % main loop of newton's method 
  for newton_iter = 1 : maxiter
    % caculate the function value at root;	
    fval = sum(mycoeff * (power(root, mypow)'));
    % caculate the derivative function value at root;	
    dfval = sum((mycoeff .* mypow) * (power(root, (mypow-1))'));
    
    if (fval < 0)
      lowbound = root;
    else
      highbound = root;
    end
    dxold = dx;  

    % if the fval is infinity, bisect the region, otherwise use newton still.
    if (fval >= realmax) | (2 * abs(fval) > abs(dxold * dfval)) | (fval ...
						  ~= fval) | (dfval ...
						  ~= dfval)
      dx = (highbound - lowbound) / 2;
      root = lowbound + dx;
    else
      if (abs(fval) <= eps) & (abs(dfval) <= eps)
	dx = 0.0;	
      elseif (abs(dfval) <= eps)
	dx = 0.0;	
	fprintf('Divided by zero!\n');
      else
	dx = fval / dfval;
      end
      root = root - dx;

      % if newton thinks to go outside the region, bisect the region
      if (root < lowbound) | (root > highbound)
	dx = (highbound - lowbound) / 2.0;
	root = lowbound + dx;
      end
    end
    
    % check converge
    if abs(fval) < epsilon
      break;
    end
  end
  if abs(fval) > epsilon
    fprintf('not converged, fval = %f\n', fval);
  end
%  root = power(root, 1.0/(maxpow/1000.0));
  root = log(root) / maxpow * 1000.0;    
  return;
end