function [X,w] = update_X_adm_bd(n,w,Y,P,beta,tol)

% parameters
tol = min(1e-6,tol);
sigma = 1e-4;
k = 10;
mu = 1/beta;

% initialization
tmp = Y - eye(n) - P/beta;
tmp = (tmp + tmp')/2;
[V,D] = eig(tmp);
u = diag(D);
fval = 0.5*norm(max(u-w,0))^2;
grad = -max(u-w,0);
step_init = 1;

iter = 1;
f = -inf*ones(k,1);
f(1) = fval;
fmax = max(f);

while 1

  x = max(u - w,0);
  dval = mu*max(x) + 0.5*norm(x - u)^2 - 0.5*norm(u)^2;
  if abs(fval + dval)/max(1,abs(fval - 0.5*norm(u)^2)) < tol
    X = V*diag(max(u - w,0))/V + eye(n);
    break
  end
  
  d = simplex_proj(w - step_init*grad, mu) - w;
  Delta = grad'*d;
  
  stepsize = 1;
  wnew = w + stepsize*d;
  fval = 0.5*norm(max(u-wnew,0))^2;
  
  while fval > fmax + stepsize*sigma*Delta && stepsize > 1e-10
    
    stepsize = stepsize/2;
    wnew = w + stepsize*d;
    fval = 0.5*norm(max(u-wnew,0))^2;
    
  end
  
  if stepsize <= 1e-10
    X = V*diag(max(u - w,0))/V + eye(n);
    break
  end
  
  gradnew = -max(u-wnew,0);
  
  dw = wnew - w;
  dg = gradnew - grad;
  step_init = min(max(norm(dw)^2/(dw'*dg),1e-8),1);
  
  grad = gradnew;
  w = wnew;
  f(mod(iter,k) + 1) = fval;
  fmax = max(f);
%   
%   fprintf(' iter, fval = %g, %g\n', iter, fval)
%   pause
  
  iter = iter + 1;
  
end