l1tf <- function(y,lambda){

# PARAMETERS
ALPHA     = 0.01;   # backtracking linesearch parameter (0,0.5]
BETA      = 0.5;    # backtracking linesearch parameter (0,1)
MU        = 2;      # IPM parameter: t update
MAXITER   = 40;     # IPM parameter: max iteration of IPM
MAXLSITER = 20;     # IPM parameter: max iteration of line search
TOL       = 1e-4;   # IPM parameter: tolerance

# DIMENSIONS
n   = length(y);    # length of signal x
m   = n-2;          # length of Dx


# OPERATOR MATRICES
I2  = diag(n-2);
O2  = matrix(0,n-2,1);
D   = cbind(I2,O2,O2) + cbind(O2-2*I2,O2) + cbind(O2,O2,I2);

DDT = D%*%t(D);
Dy  = D%*%y;

# VARIABLES
z   = rep(1,m);   % dual variable
mu1 = rep(1,m);    % dual of dual variable
mu2 = rep(1,m);    % dual of dual variable

t    = 1e-10; 
pobj =  Inf;
dobj =  0;
step =  Inf;
f1   =  z-lambda;
f2   = -z-lambda;

for(iters in 0:MAXITER){

    DTz  = t(t(z)%*%D);
    DDTz = D%*%DTz;
    w    = Dy-(mu1-mu2);
    
    pobj1 = 0.5*t(w)%*%solve(DDT,w)+lambda*sum(mu1+mu2);
    pobj2 = 0.5*t(DTz)%*%DTz+lambda*sum(abs(Dy-DDTz));
    pobj = min(pobj1,pobj2);
    dobj = -0.5*t(DTz)%*%DTz+t(Dy)%*%z;

    gap  =  pobj - dobj;

    # STOPPING CRITERION
    if (gap <= TOL){
        x = y-t(D)%*%z;
        return;
      }

    if (step >= 0.2){
        t =max(2*m*MU/gap, 1.2*t);
      }

    # CALCULATE NEWTON STEP
    
    rz      =  DDTz - w;
    S       =  DDT-sparse(1:m,1:m,mu1/f1+mu2/f2);
    r       = -DDTz + Dy + (1/t)/f1 - (1/t)/f2;
    dz      =  solve(S,r);
    dmu1    = -(mu1+((1/t)+dz*mu1)/f1);
    dmu2    = -(mu2+((1/t)-dz*mu2)/f2);

    resDual = rz;
    resCent = rbind(-mu1*f1-1/t, -mu2*f2-1/t);
    residual= rbind(resDual,resCent);

    # BACKTRACKING LINESEARCH
    negIdx1 = (dmu1 < 0); 
    negIdx2 = (dmu2 < 0);
    step = 1;
    if any(negIdx1){
        step = min( step, 0.99*min(-mu1[negIdx1)/dmu1[negIdx1]) );
      }
    if (any(negIdx2){
        step = min( step, 0.99*min(-mu2[negIdx2]/dmu2[negIdx2]) );
      }

    for( liter in 1:MAXLSITER){
        newz    =  z  + step*dz;
        newmu1  =  mu1 + step*dmu1;
        newmu2  =  mu2 + step*dmu2;
        newf1   =  newz - lambda;
        newf2   = -newz - lambda;

        # UPDATE RESIDUAL
        newResDual  = DDT*newz - Dy + newmu1 - newmu2;
        newResCent  = rbind(-newmu1.*newf1-1/t, -newmu2*newf2-1/t);
        newResidual = rbind(newResDual, newResCent);
        
        if ( max(max(newf1),max(newf2)) < 0 && ...
            norm(newResidual) <= (1-ALPHA*step)*norm(residual) )
            break;
      }
        step = BETA*step;
      }
    # UPDATE PRIMAL AND DUAL VARIABLES
    z  = newz; mu1 = newmu1; mu2 = newmu2; f1 = newf1; f2 = newf2;
  }

# The solution may be close at this point, but does not meet the stopping
# criterion (in terms of duality gap).
x = y-t(D)%*%z;
if (iters >= MAXITER){
    status = "maxiter exceeded";
    cat(status,"\n");
    return;
  }