1 files changed, 113 insertions, 0 deletions
diff --git a/R/l1tf.R b/R/l1tf.R
new file mode 100644
index 00000000..eab5208d
--- /dev/null
+++ b/R/l1tf.R
@@ -0,0 +1,113 @@
+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;
+  }
+