diff options
Diffstat (limited to 'l1tf.R')
| -rw-r--r-- | l1tf.R | 113 |
1 files changed, 0 insertions, 113 deletions
diff --git a/l1tf.R b/l1tf.R deleted file mode 100644 index eab5208d..00000000 --- a/l1tf.R +++ /dev/null @@ -1,113 +0,0 @@ -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;
- }
-
|