function [w,b, invH]= logist_X(X,Y,alpha, inind,w0)
  % X is of size n x k
  % Y is in 1..C
  % hence alpha is of length k x C-1
  % inind is of length k (not k x C)
  % w0 is length k x C-1
  % output:
  % w is length k x C-1
  % b is length C-1
  
  % check inputs and fill in defaults
  [n,k]=size(X);
  Y=Y';
  C=max(Y);
  assert(length(Y)==n);

  if (~exist('alpha','var'))
    alpha=zeros(k*(C-1),1);
  end
  if (exist('alpha','var') & length(alpha)==1 & (k*(C-1)>1))
    alpha=alpha*ones(k*(C-1),1);
  end
  if (~exist('inind','var'))
    inind=true(k,1);
  end
  if (~exist('w0','var')),
    w0=zeros(size(alpha));
  end
  inind=inind(:);
  assert(islogical(inind));

  % augment constant to input and prune input according to inind
  X=[X(:,inind) ones(n,1)];
  k_in=sum(inind);
  b_ind=(k_in+1)*[1:C-1]';
  w_ind=setdiff([1:(k_in+1)*(C-1)], b_ind);
  
  % set optimization constants
  opts=optimset('MaxIter',100);
  opts=optimset(opts,'Display','iter');
  alpha_const=1E-5;  % the alpha associated with the constant feature
  
  % initialize w and alpha 
  wb0=zeros((k_in+1)*(C-1),1);
  in_kc = repmat(inind,C-1,1);
  wb0(w_ind)=w0(in_kc);
  alpha_wb=zeros(size(wb0));
  alpha_wb(w_ind)=alpha(in_kc);
  alpha_wb(b_ind)=alpha_const;

%  [wb,H]=newton(@(w) logist_X_fgH(w,X,Y,alpha_wb), wb0, opts);
  [wb,H]=newton(@logist_X_fgH, wb0, opts, X,Y,alpha_wb);
  
  % return values
  invH=inv(H(w_ind,w_ind));
  w=w0;
  w(in_kc)=wb(w_ind);
  b=wb(b_ind);