1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
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;
}
|
