Initial import from https://statweb.stanford.edu/~candes/software/l1magic/v1.11

Additional Clean up of Mac dirs and tex generated files

12 files changed, 1860 insertions, 0 deletions

diff --git a/Optimization/cgsolve.m b/Optimization/cgsolve.m
new file mode 100644
index 0000000..dbbc4d8
--- /dev/null
+++ b/Optimization/cgsolve.m
@@ -0,0 +1,76 @@
+% cgsolve.m
+%
+% Solve a symmetric positive definite system Ax = b via conjugate gradients.
+%
+% Usage: [x, res, iter] = cgsolve(A, b, tol, maxiter, verbose)
+%
+% A - Either an NxN matrix, or a function handle.
+%
+% b - N vector
+%
+% tol - Desired precision.  Algorithm terminates when 
+%    norm(Ax-b)/norm(b) < tol .
+%
+% maxiter - Maximum number of iterations.
+%
+% verbose - If 0, do not print out progress messages.
+%    If and integer greater than 0, print out progress every 'verbose' iters.
+%
+% Written by: Justin Romberg, Caltech
+% Email: jrom@acm.caltech.edu
+% Created: October 2005
+%
+
+function [x, res, iter] = cgsolve(A, b, tol, maxiter, verbose)
+
+if (nargin < 5), verbose = 1; end
+
+implicit = isa(A,'function_handle');
+
+x = zeros(length(b),1);
+r = b;
+d = r;
+delta = r'*r;
+delta0 = b'*b;
+numiter = 0;
+bestx = x;
+bestres = sqrt(delta/delta0); 
+while ((numiter < maxiter) && (delta > tol^2*delta0))
+
+  % q = A*d
+  if (implicit), q = A(d);  else  q = A*d;  end
+ 
+  alpha = delta/(d'*q);
+  x = x + alpha*d;
+  
+  if (mod(numiter+1,50) == 0)
+    % r = b - Aux*x
+    if (implicit), r = b - A(x);  else  r = b - A*x;  end
+  else
+    r = r - alpha*q;
+  end
+  
+  deltaold = delta;
+  delta = r'*r;
+  beta = delta/deltaold;
+  d = r + beta*d;
+  numiter = numiter + 1;
+  if (sqrt(delta/delta0) < bestres)
+    bestx = x;
+    bestres = sqrt(delta/delta0);
+  end    
+  
+  if ((verbose) && (mod(numiter,verbose)==0))
+    disp(sprintf('cg: Iter = %d, Best residual = %8.3e, Current residual = %8.3e', ...
+      numiter, bestres, sqrt(delta/delta0)));
+  end
+  
+end
+
+if (verbose)
+  disp(sprintf('cg: Iterations = %d, best residual = %14.8e', numiter, bestres));
+end
+x = bestx;
+res = bestres;
+iter = numiter;
+
diff --git a/Optimization/l1dantzig_pd.m b/Optimization/l1dantzig_pd.m
new file mode 100644
index 0000000..6a57dea
--- /dev/null
+++ b/Optimization/l1dantzig_pd.m
@@ -0,0 +1,234 @@
+% l1dantzig_pd.m
+%
+% Solves
+% min_x  ||x||_1  subject to  ||A'(Ax-b)||_\infty <= epsilon
+%
+% Recast as linear program
+% min_{x,u}  sum(u)  s.t.  x - u <= 0
+%                         -x - u <= 0
+%             A'(Ax-b) - epsilon <= 0
+%            -A'(Ax-b) - epsilon <= 0
+% and use primal-dual interior point method.
+%
+% Usage: xp = l1dantzig_pd(x0, A, At, b, epsilon, pdtol, pdmaxiter, cgtol, cgmaxiter)
+%
+% x0 - Nx1 vector, initial point.
+%
+% A - Either a handle to a function that takes a N vector and returns a K 
+%     vector , or a KxN matrix.  If A is a function handle, the algorithm
+%     operates in "largescale" mode, solving the Newton systems via the
+%     Conjugate Gradients algorithm.
+%
+% At - Handle to a function that takes a K vector and returns an N vector.
+%      If A is a KxN matrix, At is ignored.
+%
+% b - Kx1 vector of observations.
+%
+% epsilon - scalar or Nx1 vector of correlation constraints
+%
+% pdtol - Tolerance for primal-dual algorithm (algorithm terminates if
+%     the duality gap is less than pdtol).  
+%     Default = 1e-3.
+%
+% pdmaxiter - Maximum number of primal-dual iterations.  
+%     Default = 50.
+%
+% cgtol - Tolerance for Conjugate Gradients; ignored if A is a matrix.
+%     Default = 1e-8.
+%
+% cgmaxiter - Maximum number of iterations for Conjugate Gradients; ignored
+%     if A is a matrix.
+%     Default = 200.
+%
+% Written by: Justin Romberg, Caltech
+% Email: jrom@acm.caltech.edu
+% Created: October 2005
+%
+
+function xp = l1dantzig_pd(x0, A, At, b, epsilon, pdtol, pdmaxiter, cgtol, cgmaxiter)
+
+largescale = isa(A,'function_handle');
+
+if (nargin < 6), pdtol = 1e-3;  end
+if (nargin < 7), pdmaxiter = 50;  end
+if (nargin < 8), cgtol = 1e-8;  end
+if (nargin < 9), cgmaxiter = 200;  end
+
+N = length(x0);
+
+alpha = 0.01;
+beta = 0.5;
+mu = 10;
+
+gradf0 = [zeros(N,1); ones(N,1)];
+
+
+% starting point --- make sure that it is feasible
+if (largescale)
+  if (max( abs(At(A(x0) - b)) - epsilon ) > 0)
+    disp('Starting point infeasible; using x0 = At*inv(AAt)*y.');
+    AAt = @(z) A(At(z));
+    [w, cgres] = cgsolve(AAt, b, cgtol, cgmaxiter, 0);
+    if (cgres > 1/2)
+      disp('A*At is ill-conditioned: cannot find starting point');
+      xp = x0;
+      return;
+    end
+    x0 = At(w);
+  end
+else
+  if (max(abs(A'*(A*x0 - b)) - epsilon ) > 0)
+    disp('Starting point infeasible; using x0 = At*inv(AAt)*y.');
+    opts.POSDEF = true; opts.SYM = true;
+    [w, hcond] = linsolve(A*A', b, opts);
+    if (hcond < 1e-14)
+      disp('A*At is ill-conditioned: cannot find starting point');
+      xp = x0;
+      return;
+    end
+    x0 = A'*w;
+  end  
+end
+x = x0;
+u = (0.95)*abs(x0) + (0.10)*max(abs(x0));
+
+% set up for the first iteration
+if (largescale)
+  Atr = At(A(x) - b);
+else
+  Atr = A'*(A*x - b);
+end
+fu1 = x - u;
+fu2 = -x - u;
+fe1 = Atr - epsilon;
+fe2 = -Atr - epsilon;
+lamu1 = -(1./fu1);
+lamu2 = -(1./fu2);
+lame1 = -(1./fe1);
+lame2 = -(1./fe2);
+if (largescale)
+  AtAv = At(A(lame1-lame2));
+else
+  AtAv = A'*(A*(lame1-lame2));
+end
+
+% sdg = surrogate duality gap
+sdg = -[fu1; fu2; fe1; fe2]'*[lamu1; lamu2; lame1; lame2];
+tau = mu*(4*N)/sdg;
+
+% residuals
+rdual = gradf0 + [lamu1-lamu2 + AtAv; -lamu1-lamu2];
+rcent = -[lamu1.*fu1; lamu2.*fu2; lame1.*fe1; lame2.*fe2] - (1/tau);
+resnorm = norm([rdual; rcent]);
+
+% iterations
+pditer = 0;
+done = (sdg < pdtol) | (pditer >= pdmaxiter);
+while (~done)
+
+  % solve for step direction
+  w2 = - 1 - (1/tau)*(1./fu1 + 1./fu2);
+  
+  sig11 = -lamu1./fu1 - lamu2./fu2;
+  sig12 = lamu1./fu1 - lamu2./fu2;
+  siga = -(lame1./fe1 + lame2./fe2);
+  sigx = sig11 - sig12.^2./sig11;
+  
+  if (largescale)
+    w1 = -(1/tau)*( At(A(1./fe2-1./fe1)) + 1./fu2 - 1./fu1 );
+    w1p = w1 - (sig12./sig11).*w2;
+    hpfun = @(z) At(A(siga.*At(A(z)))) + sigx.*z;
+    [dx, cgres, cgiter] = cgsolve(hpfun, w1p, cgtol, cgmaxiter, 0);
+    if (cgres > 1/2)
+      disp('Cannot solve system.  Returning previous iterate.  (See Section 4 of notes for more information.)');
+      xp = x;
+      return
+    end
+    AtAdx = At(A(dx));
+  else
+    w1 = -(1/tau)*( A'*(A*(1./fe2-1./fe1)) + 1./fu2 - 1./fu1 );
+    w1p = w1 - (sig12./sig11).*w2;
+    Hp = A'*(A*sparse(diag(siga))*A')*A + diag(sigx);
+    opts.POSDEF = true; opts.SYM = true;
+    [dx, hcond] = linsolve(Hp, w1p,opts);
+    if (hcond < 1e-14)
+      disp('Matrix ill-conditioned.  Returning previous iterate.  (See Section 4 of notes for more information.)');
+      xp = x;
+      return
+    end
+    AtAdx = A'*(A*dx);
+  end
+  du = w2./sig11 - (sig12./sig11).*dx;
+  
+  dlamu1 = -(lamu1./fu1).*(dx-du) - lamu1 - (1/tau)*1./fu1;
+  dlamu2 = -(lamu2./fu2).*(-dx-du) - lamu2 - (1/tau)*1./fu2;
+  dlame1 = -(lame1./fe1).*(AtAdx) - lame1 - (1/tau)*1./fe1;
+  dlame2 = -(lame2./fe2).*(-AtAdx) - lame2 - (1/tau)*1./fe2;
+  if (largescale)  
+    AtAdv = At(A(dlame1-dlame2));  
+  else
+    AtAdv = A'*(A*(dlame1-dlame2));  
+  end
+	
+  
+  % find minimal step size that keeps ineq functions < 0, dual vars > 0
+  iu1 = find(dlamu1 < 0); iu2 = find(dlamu2 < 0); 
+  ie1 = find(dlame1 < 0); ie2 = find(dlame2 < 0);
+  ifu1 = find((dx-du) > 0); ifu2 = find((-dx-du) > 0); 
+  ife1 = find(AtAdx > 0); ife2 = find(-AtAdx > 0); 
+  smax = min(1,min([...
+    -lamu1(iu1)./dlamu1(iu1); -lamu2(iu2)./dlamu2(iu2); ...
+    -lame1(ie1)./dlame1(ie1); -lame2(ie2)./dlame2(ie2); ...
+    -fu1(ifu1)./(dx(ifu1)-du(ifu1)); -fu2(ifu2)./(-dx(ifu2)-du(ifu2)); ...
+    -fe1(ife1)./AtAdx(ife1); -fe2(ife2)./(-AtAdx(ife2)) ]));
+  s = 0.99*smax;
+  
+  % backtracking line search
+  suffdec = 0;
+  backiter = 0;
+  while (~suffdec)
+    xp = x + s*dx;  up = u + s*du;
+    Atrp = Atr + s*AtAdx;  AtAvp = AtAv + s*AtAdv;
+    fu1p = fu1 + s*(dx-du);  fu2p = fu2 + s*(-dx-du);
+    fe1p = fe1 + s*AtAdx;  fe2p = fe2 + s*(-AtAdx);
+    lamu1p = lamu1 + s*dlamu1;  lamu2p = lamu2 + s*dlamu2;
+    lame1p = lame1 + s*dlame1; lame2p = lame2 + s*dlame2;
+    rdp = gradf0 + [lamu1p-lamu2p + AtAvp; -lamu1p-lamu2p];
+    rcp = -[lamu1p.*fu1p; lamu2p.*fu2p; lame1p.*fe1p; lame2p.*fe2p] - (1/tau);
+    suffdec = (norm([rdp; rcp]) <= (1-alpha*s)*resnorm);
+    s = beta*s;
+    backiter = backiter+1;
+    if (backiter > 32)
+      disp('Stuck backtracking, returning last iterate.  (See Section 4 of notes for more information.)')
+      xp = x;
+      return
+    end
+  end
+    
+  % setup for next iteration
+  x = xp;  u = up;
+  Atr = Atrp;  AtAv = AtAvp;
+  fu1 = fu1p; fu2 = fu2p; 
+  fe1 = fe1p; fe2 = fe2p;
+  lamu1 = lamu1p; lamu2 = lamu2p; 
+  lame1 = lame1p; lame2 = lame2p;
+  
+  sdg = -[fu1; fu2; fe1; fe2]'*[lamu1; lamu2; lame1; lame2];
+  tau = mu*(4*N)/sdg;
+
+  rdual = rdp;
+  rcent = -[lamu1.*fu1; lamu2.*fu2; lame1.*fe1; lame2.*fe2] - (1/tau);
+  resnorm = norm([rdual; rcent]);
+  
+  pditer = pditer+1;
+  done = (sdg < pdtol) | (pditer >= pdmaxiter);
+  
+  disp(sprintf('Iteration = %d, tau = %8.3e, Primal = %8.3e, PDGap = %8.3e, Dual res = %8.3e',...
+    pditer, tau, sum(u), sdg, norm(rdual)));
+  if (largescale)
+    disp(sprintf('                CG Res = %8.3e, CG Iter = %d', cgres, cgiter));
+  else
+    disp(sprintf('                  H11p condition number = %8.3e', hcond));
+  end
+  
+end
diff --git a/Optimization/l1decode_pd.m b/Optimization/l1decode_pd.m
new file mode 100644
index 0000000..2757404
--- /dev/null
+++ b/Optimization/l1decode_pd.m
@@ -0,0 +1,175 @@
+% l1decode_pd.m
+%
+% Decoding via linear programming.
+% Solve
+% min_x  ||b-Ax||_1 .
+%
+% Recast as the linear program
+% min_{x,u} sum(u)  s.t.  -Ax - u + y <= 0
+%                          Ax - u - y <= 0
+% and solve using primal-dual interior point method.
+%
+% Usage: xp = l1decode_pd(x0, A, At, y, pdtol, pdmaxiter, cgtol, cgmaxiter)
+%
+% x0 - Nx1 vector, initial point.
+%
+% A - Either a handle to a function that takes a N vector and returns a M 
+%     vector, or a MxN matrix.  If A is a function handle, the algorithm
+%     operates in "largescale" mode, solving the Newton systems via the
+%     Conjugate Gradients algorithm.
+%
+% At - Handle to a function that takes an M vector and returns an N vector.
+%      If A is a matrix, At is ignored.
+%
+% y - Mx1 observed code (M > N).
+%
+% pdtol - Tolerance for primal-dual algorithm (algorithm terminates if
+%     the duality gap is less than pdtol).  
+%     Default = 1e-3.
+%
+% pdmaxiter - Maximum number of primal-dual iterations.  
+%     Default = 50.
+%
+% cgtol - Tolerance for Conjugate Gradients; ignored if A is a matrix.
+%     Default = 1e-8.
+%
+% cgmaxiter - Maximum number of iterations for Conjugate Gradients; ignored
+%     if A is a matrix.
+%     Default = 200.
+%
+% Written by: Justin Romberg, Caltech
+% Email: jrom@acm.caltech.edu
+% Created: October 2005
+%
+
+function xp = l1decode_pd(x0, A, At, y, pdtol, pdmaxiter, cgtol, cgmaxiter)
+
+largescale = isa(A,'function_handle');   
+
+if (nargin < 5), pdtol = 1e-3; end
+if (nargin < 6), pdmaxiter = 50; end
+if (nargin < 7), cgtol = 1e-8; end
+if (nargin < 8), cgmaxiter = 200; end
+
+N = length(x0);
+M = length(y);
+
+alpha = 0.01;
+beta = 0.5;
+mu = 10;
+
+gradf0 = [zeros(N,1); ones(M,1)];
+
+x = x0;
+if (largescale),  Ax = A(x);  else  Ax = A*x;  end
+u = (0.95)*abs(y-Ax) + (0.10)*max(abs(y-Ax));
+
+fu1 = Ax - y - u;
+fu2 = -Ax + y - u;
+
+lamu1 = -1./fu1;
+lamu2 = -1./fu2;
+
+if (largescale), Atv = At(lamu1-lamu2);  else Atv = A'*(lamu1-lamu2);  end
+
+sdg = -(fu1'*lamu1 + fu2'*lamu2);
+tau = mu*2*M/sdg;
+
+rcent = [-lamu1.*fu1; -lamu2.*fu2] - (1/tau);
+rdual = gradf0 + [Atv; -lamu1-lamu2];
+resnorm = norm([rdual; rcent]);
+
+pditer = 0;
+done = (sdg < pdtol)| (pditer >= pdmaxiter);
+while (~done)
+  
+  pditer = pditer + 1;
+  
+  w2 = -1 - 1/tau*(1./fu1 + 1./fu2);
+  
+  sig1 = -lamu1./fu1 - lamu2./fu2;
+  sig2 = lamu1./fu1 - lamu2./fu2;
+  sigx = sig1 - sig2.^2./sig1;
+  
+  if (largescale)
+    w1 = -1/tau*(At(-1./fu1 + 1./fu2));
+    w1p = w1 - At((sig2./sig1).*w2);
+    h11pfun = @(z) At(sigx.*A(z));
+    [dx, cgres, cgiter] = cgsolve(h11pfun, w1p, cgtol, cgmaxiter, 0); 
+    if (cgres > 1/2)
+      disp('Cannot solve system.  Returning previous iterate.  (See Section 4 of notes for more information.)');
+      xp = x;
+      return
+    end
+    Adx = A(dx);
+  else
+    w1 = -1/tau*(A'*(-1./fu1 + 1./fu2));
+    w1p = w1 - A'*((sig2./sig1).*w2);
+    H11p = A'*(sparse(diag(sigx))*A);
+    opts.POSDEF = true; opts.SYM = true;
+    [dx, hcond] = linsolve(H11p, w1p,opts);
+    if (hcond < 1e-14)
+      disp('Matrix ill-conditioned.  Returning previous iterate.  (See Section 4 of notes for more information.)');
+      xp = x;
+      return
+    end
+    Adx = A*dx;
+  end
+  
+  du = (w2 - sig2.*Adx)./sig1;
+  
+  dlamu1 = -(lamu1./fu1).*(Adx-du) - lamu1 - (1/tau)*1./fu1;
+  dlamu2 = (lamu2./fu2).*(Adx + du) -lamu2 - (1/tau)*1./fu2;
+  if (largescale), Atdv = At(dlamu1-dlamu2);  else  Atdv = A'*(dlamu1-dlamu2); end
+  
+  % make sure that the step is feasible: keeps lamu1,lamu2 > 0, fu1,fu2 < 0
+  indl = find(dlamu1 < 0);  indu = find(dlamu2 < 0);
+  s = min([1; -lamu1(indl)./dlamu1(indl); -lamu2(indu)./dlamu2(indu)]);
+  indl = find((Adx-du) > 0);  indu = find((-Adx-du) > 0);
+  s = (0.99)*min([s; -fu1(indl)./(Adx(indl)-du(indl)); -fu2(indu)./(-Adx(indu)-du(indu))]);
+  
+  % backtrack 
+  suffdec = 0;
+  backiter = 0;
+  while(~suffdec)
+    xp = x + s*dx;  up = u + s*du;
+    Axp = Ax + s*Adx;  Atvp = Atv + s*Atdv;
+    lamu1p = lamu1 + s*dlamu1;  lamu2p = lamu2 + s*dlamu2;
+    fu1p = Axp - y - up;  fu2p = -Axp + y - up;
+    rdp = gradf0 + [Atvp; -lamu1p-lamu2p];
+    rcp = [-lamu1p.*fu1p; -lamu2p.*fu2p] - (1/tau);
+    suffdec = (norm([rdp; rcp]) <= (1-alpha*s)*resnorm);
+    s = beta*s;
+    backiter = backiter + 1;
+    if (backiter > 32)
+      disp('Stuck backtracking, returning last iterate.  (See Section 4 of notes for more information.)')
+      xp = x;
+      return
+    end
+  end
+  
+  % next iteration
+  x = xp;  u = up;
+  Ax = Axp;  Atv = Atvp;
+  lamu1 = lamu1p;  lamu2 = lamu2p;
+  fu1 = fu1p;  fu2 = fu2p;
+  
+  % surrogate duality gap
+  sdg = -(fu1'*lamu1 + fu2'*lamu2);
+  tau = mu*2*M/sdg;
+  rcent = [-lamu1.*fu1; -lamu2.*fu2] - (1/tau);
+  rdual = rdp;
+  resnorm = norm([rdual; rcent]);
+  
+  done = (sdg < pdtol) | (pditer >= pdmaxiter);
+  
+  disp(sprintf('Iteration = %d, tau = %8.3e, Primal = %8.3e, PDGap = %8.3e, Dual res = %8.3e',...
+    pditer, tau, sum(u), sdg, norm(rdual)));
+  if (largescale)
+    disp(sprintf('                CG Res = %8.3e, CG Iter = %d', cgres, cgiter));
+  else
+    disp(sprintf('                  H11p condition number = %8.3e', hcond));
+  end  
+  
+end
+
diff --git a/Optimization/l1eq_pd.m b/Optimization/l1eq_pd.m
new file mode 100644
index 0000000..80c058e
--- /dev/null
+++ b/Optimization/l1eq_pd.m
@@ -0,0 +1,209 @@
+% l1eq_pd.m
+%
+% Solve
+% min_x ||x||_1  s.t.  Ax = b
+%
+% Recast as linear program
+% min_{x,u} sum(u)  s.t.  -u <= x <= u,  Ax=b
+% and use primal-dual interior point method
+%
+% Usage: xp = l1eq_pd(x0, A, At, b, pdtol, pdmaxiter, cgtol, cgmaxiter)
+%
+% x0 - Nx1 vector, initial point.
+%
+% A - Either a handle to a function that takes a N vector and returns a K 
+%     vector , or a KxN matrix.  If A is a function handle, the algorithm
+%     operates in "largescale" mode, solving the Newton systems via the
+%     Conjugate Gradients algorithm.
+%
+% At - Handle to a function that takes a K vector and returns an N vector.
+%      If A is a KxN matrix, At is ignored.
+%
+% b - Kx1 vector of observations.
+%
+% pdtol - Tolerance for primal-dual algorithm (algorithm terminates if
+%     the duality gap is less than pdtol).  
+%     Default = 1e-3.
+%
+% pdmaxiter - Maximum number of primal-dual iterations.  
+%     Default = 50.
+%
+% cgtol - Tolerance for Conjugate Gradients; ignored if A is a matrix.
+%     Default = 1e-8.
+%
+% cgmaxiter - Maximum number of iterations for Conjugate Gradients; ignored
+%     if A is a matrix.
+%     Default = 200.
+%
+% Written by: Justin Romberg, Caltech
+% Email: jrom@acm.caltech.edu
+% Created: October 2005
+%
+
+function xp = l1eq_pd(x0, A, At, b, pdtol, pdmaxiter, cgtol, cgmaxiter)
+
+largescale = isa(A,'function_handle');
+
+if (nargin < 5), pdtol = 1e-3;  end
+if (nargin < 6), pdmaxiter = 50;  end
+if (nargin < 7), cgtol = 1e-8;  end
+if (nargin < 8), cgmaxiter = 200;  end
+
+N = length(x0);
+
+alpha = 0.01;
+beta = 0.5;
+mu = 10;
+
+gradf0 = [zeros(N,1); ones(N,1)];
+
+% starting point --- make sure that it is feasible
+if (largescale)
+  if (norm(A(x0)-b)/norm(b) > cgtol)
+    disp('Starting point infeasible; using x0 = At*inv(AAt)*y.');
+    AAt = @(z) A(At(z));
+    [w, cgres, cgiter] = cgsolve(AAt, b, cgtol, cgmaxiter, 0);
+    if (cgres > 1/2)
+      disp('A*At is ill-conditioned: cannot find starting point');
+      xp = x0;
+      return;
+    end
+    x0 = At(w);
+  end
+else
+  if (norm(A*x0-b)/norm(b) > cgtol)
+    disp('Starting point infeasible; using x0 = At*inv(AAt)*y.');
+    opts.POSDEF = true; opts.SYM = true;
+    [w, hcond] = linsolve(A*A', b, opts);
+    if (hcond < 1e-14)
+      disp('A*At is ill-conditioned: cannot find starting point');
+      xp = x0;
+      return;
+    end
+    x0 = A'*w;
+  end  
+end
+x = x0;
+u = (0.95)*abs(x0) + (0.10)*max(abs(x0));
+
+% set up for the first iteration
+fu1 = x - u;
+fu2 = -x - u;
+lamu1 = -1./fu1;
+lamu2 = -1./fu2;
+if (largescale)
+  v = -A(lamu1-lamu2);
+  Atv = At(v);
+  rpri = A(x) - b;
+else
+  v = -A*(lamu1-lamu2);
+  Atv = A'*v;
+  rpri = A*x - b;
+end
+
+sdg = -(fu1'*lamu1 + fu2'*lamu2);
+tau = mu*2*N/sdg;
+
+rcent = [-lamu1.*fu1; -lamu2.*fu2] - (1/tau);
+rdual = gradf0 + [lamu1-lamu2; -lamu1-lamu2] + [Atv; zeros(N,1)];
+resnorm = norm([rdual; rcent; rpri]);
+
+pditer = 0;
+done = (sdg < pdtol) | (pditer >= pdmaxiter);
+while (~done)
+  
+  pditer = pditer + 1;
+  
+  w1 = -1/tau*(-1./fu1 + 1./fu2) - Atv;
+  w2 = -1 - 1/tau*(1./fu1 + 1./fu2);
+  w3 = -rpri;
+  
+  sig1 = -lamu1./fu1 - lamu2./fu2;
+  sig2 = lamu1./fu1 - lamu2./fu2;
+  sigx = sig1 - sig2.^2./sig1;
+  
+  if (largescale)
+    w1p = w3 - A(w1./sigx - w2.*sig2./(sigx.*sig1));
+    h11pfun = @(z) -A(1./sigx.*At(z));
+    [dv, cgres, cgiter] = cgsolve(h11pfun, w1p, cgtol, cgmaxiter, 0);
+    if (cgres > 1/2)
+      disp('Cannot solve system.  Returning previous iterate.  (See Section 4 of notes for more information.)');
+      xp = x;
+      return
+    end
+    dx = (w1 - w2.*sig2./sig1 - At(dv))./sigx;
+    Adx = A(dx);
+    Atdv = At(dv);
+  else
+    w1p = -(w3 - A*(w1./sigx - w2.*sig2./(sigx.*sig1)));
+    H11p = A*(sparse(diag(1./sigx))*A');
+    opts.POSDEF = true; opts.SYM = true;
+    [dv,hcond] = linsolve(H11p, w1p, opts);
+    if (hcond < 1e-14)
+      disp('Matrix ill-conditioned.  Returning previous iterate.  (See Section 4 of notes for more information.)');
+      xp = x;
+      return
+    end
+    dx = (w1 - w2.*sig2./sig1 - A'*dv)./sigx;
+    Adx = A*dx;
+    Atdv = A'*dv;
+  end
+  
+  du = (w2 - sig2.*dx)./sig1;
+  
+  dlamu1 = (lamu1./fu1).*(-dx+du) - lamu1 - (1/tau)*1./fu1;
+  dlamu2 = (lamu2./fu2).*(dx+du) - lamu2 - 1/tau*1./fu2;
+  
+  % make sure that the step is feasible: keeps lamu1,lamu2 > 0, fu1,fu2 < 0
+  indp = find(dlamu1 < 0);  indn = find(dlamu2 < 0);
+  s = min([1; -lamu1(indp)./dlamu1(indp); -lamu2(indn)./dlamu2(indn)]);
+  indp = find((dx-du) > 0);  indn = find((-dx-du) > 0);
+  s = (0.99)*min([s; -fu1(indp)./(dx(indp)-du(indp)); -fu2(indn)./(-dx(indn)-du(indn))]);
+  
+  % backtracking line search
+  suffdec = 0;
+  backiter = 0;
+  while (~suffdec)
+    xp = x + s*dx;  up = u + s*du; 
+    vp = v + s*dv;  Atvp = Atv + s*Atdv; 
+    lamu1p = lamu1 + s*dlamu1;  lamu2p = lamu2 + s*dlamu2;
+    fu1p = xp - up;  fu2p = -xp - up;  
+    rdp = gradf0 + [lamu1p-lamu2p; -lamu1p-lamu2p] + [Atvp; zeros(N,1)];
+    rcp = [-lamu1p.*fu1p; -lamu2p.*fu2p] - (1/tau);
+    rpp = rpri + s*Adx;
+    suffdec = (norm([rdp; rcp; rpp]) <= (1-alpha*s)*resnorm);
+    s = beta*s;
+    backiter = backiter + 1;
+    if (backiter > 32)
+      disp('Stuck backtracking, returning last iterate.  (See Section 4 of notes for more information.)')
+      xp = x;
+      return
+    end
+  end
+  
+  
+  % next iteration
+  x = xp;  u = up;
+  v = vp;  Atv = Atvp; 
+  lamu1 = lamu1p;  lamu2 = lamu2p;
+  fu1 = fu1p;  fu2 = fu2p;
+  
+  % surrogate duality gap
+  sdg = -(fu1'*lamu1 + fu2'*lamu2);
+  tau = mu*2*N/sdg;
+  rpri = rpp;
+  rcent = [-lamu1.*fu1; -lamu2.*fu2] - (1/tau);
+  rdual = gradf0 + [lamu1-lamu2; -lamu1-lamu2] + [Atv; zeros(N,1)];
+  resnorm = norm([rdual; rcent; rpri]);
+  
+  done = (sdg < pdtol) | (pditer >= pdmaxiter);
+  
+  disp(sprintf('Iteration = %d, tau = %8.3e, Primal = %8.3e, PDGap = %8.3e, Dual res = %8.3e, Primal res = %8.3e',...
+    pditer, tau, sum(u), sdg, norm(rdual), norm(rpri)));
+  if (largescale)
+    disp(sprintf('                  CG Res = %8.3e, CG Iter = %d', cgres, cgiter));
+  else
+    disp(sprintf('                  H11p condition number = %8.3e', hcond));
+  end
+  
+end
diff --git a/Optimization/l1qc_logbarrier.m b/Optimization/l1qc_logbarrier.m
new file mode 100644
index 0000000..388529e
--- /dev/null
+++ b/Optimization/l1qc_logbarrier.m
@@ -0,0 +1,116 @@
+% l1qc_logbarrier.m
+%
+% Solve quadratically constrained l1 minimization:
+% min ||x||_1   s.t.  ||Ax - b||_2 <= \epsilon
+%
+% Reformulate as the second-order cone program
+% min_{x,u}  sum(u)   s.t.    x - u <= 0,
+%                            -x - u <= 0,
+%      1/2(||Ax-b||^2 - \epsilon^2) <= 0
+% and use a log barrier algorithm.
+%
+% Usage:  xp = l1qc_logbarrier(x0, A, At, b, epsilon, lbtol, mu, cgtol, cgmaxiter)
+%
+% x0 - Nx1 vector, initial point.
+%
+% A - Either a handle to a function that takes a N vector and returns a K 
+%     vector , or a KxN matrix.  If A is a function handle, the algorithm
+%     operates in "largescale" mode, solving the Newton systems via the
+%     Conjugate Gradients algorithm.
+%
+% At - Handle to a function that takes a K vector and returns an N vector.
+%      If A is a KxN matrix, At is ignored.
+%
+% b - Kx1 vector of observations.
+%
+% epsilon - scalar, constraint relaxation parameter
+%
+% lbtol - The log barrier algorithm terminates when the duality gap <= lbtol.
+%         Also, the number of log barrier iterations is completely
+%         determined by lbtol.
+%         Default = 1e-3.
+%
+% mu - Factor by which to increase the barrier constant at each iteration.
+%      Default = 10.
+%
+% cgtol - Tolerance for Conjugate Gradients; ignored if A is a matrix.
+%     Default = 1e-8.
+%
+% cgmaxiter - Maximum number of iterations for Conjugate Gradients; ignored
+%     if A is a matrix.
+%     Default = 200.
+%
+% Written by: Justin Romberg, Caltech
+% Email: jrom@acm.caltech.edu
+% Created: October 2005
+%
+
+function xp = l1qc_logbarrier(x0, A, At, b, epsilon, lbtol, mu, cgtol, cgmaxiter)  
+
+largescale = isa(A,'function_handle');
+
+if (nargin < 6), lbtol = 1e-3; end
+if (nargin < 7), mu = 10; end
+if (nargin < 8), cgtol = 1e-8; end
+if (nargin < 9), cgmaxiter = 200; end
+
+newtontol = lbtol;
+newtonmaxiter = 50;
+
+N = length(x0);
+
+% starting point --- make sure that it is feasible
+if (largescale)
+  if (norm(A(x0)-b) > epsilon)
+    disp('Starting point infeasible; using x0 = At*inv(AAt)*y.');
+    AAt = @(z) A(At(z));
+    [w, cgres] = cgsolve(AAt, b, cgtol, cgmaxiter, 0);
+    if (cgres > 1/2)
+      disp('A*At is ill-conditioned: cannot find starting point');
+      xp = x0;
+      return;
+    end
+    x0 = At(w);
+  end
+else
+  if (norm(A*x0-b) > epsilon)
+    disp('Starting point infeasible; using x0 = At*inv(AAt)*y.');
+    opts.POSDEF = true; opts.SYM = true;
+    [w, hcond] = linsolve(A*A', b, opts);
+    if (hcond < 1e-14)
+      disp('A*At is ill-conditioned: cannot find starting point');
+      xp = x0;
+      return;
+    end
+    x0 = A'*w;
+  end  
+end
+x = x0;
+u = (0.95)*abs(x0) + (0.10)*max(abs(x0));
+
+disp(sprintf('Original l1 norm = %.3f, original functional = %.3f', sum(abs(x0)), sum(u)));
+
+% choose initial value of tau so that the duality gap after the first
+% step will be about the origial norm
+tau = max((2*N+1)/sum(abs(x0)), 1);
+                                                                                                                          
+lbiter = ceil((log(2*N+1)-log(lbtol)-log(tau))/log(mu));
+disp(sprintf('Number of log barrier iterations = %d\n', lbiter));
+
+totaliter = 0;
+
+for ii = 1:lbiter
+
+  [xp, up, ntiter] = l1qc_newton(x, u, A, At, b, epsilon, tau, newtontol, newtonmaxiter, cgtol, cgmaxiter);
+  totaliter = totaliter + ntiter;
+  
+  disp(sprintf('\nLog barrier iter = %d, l1 = %.3f, functional = %8.3f, tau = %8.3e, total newton iter = %d\n', ...
+    ii, sum(abs(xp)), sum(up), tau, totaliter));
+  
+  x = xp;
+  u = up;
+ 
+  tau = mu*tau;
+  
+end
+                   
diff --git a/Optimization/l1qc_newton.m b/Optimization/l1qc_newton.m
new file mode 100644
index 0000000..8a25cd2
--- /dev/null
+++ b/Optimization/l1qc_newton.m
@@ -0,0 +1,150 @@
+% l1qc_newton.m
+%
+% Newton algorithm for log-barrier subproblems for l1 minimization
+% with quadratic constraints.
+%
+% Usage: 
+% [xp,up,niter] = l1qc_newton(x0, u0, A, At, b, epsilon, tau, 
+%                             newtontol, newtonmaxiter, cgtol, cgmaxiter)
+%
+% x0,u0 - starting points
+%
+% A - Either a handle to a function that takes a N vector and returns a K 
+%     vector , or a KxN matrix.  If A is a function handle, the algorithm
+%     operates in "largescale" mode, solving the Newton systems via the
+%     Conjugate Gradients algorithm.
+%
+% At - Handle to a function that takes a K vector and returns an N vector.
+%      If A is a KxN matrix, At is ignored.
+%
+% b - Kx1 vector of observations.
+%
+% epsilon - scalar, constraint relaxation parameter
+%
+% tau - Log barrier parameter.
+%
+% newtontol - Terminate when the Newton decrement is <= newtontol.
+%         Default = 1e-3.
+%
+% newtonmaxiter - Maximum number of iterations.
+%         Default = 50.
+%
+% cgtol - Tolerance for Conjugate Gradients; ignored if A is a matrix.
+%     Default = 1e-8.
+%
+% cgmaxiter - Maximum number of iterations for Conjugate Gradients; ignored
+%     if A is a matrix.
+%     Default = 200.
+%
+% Written by: Justin Romberg, Caltech
+% Email: jrom@acm.caltech.edu
+% Created: October 2005
+%
+
+
+function [xp, up, niter] = l1qc_newton(x0, u0, A, At, b, epsilon, tau, newtontol, newtonmaxiter, cgtol, cgmaxiter) 
+
+% check if the matrix A is implicit or explicit
+largescale = isa(A,'function_handle');
+
+% line search parameters
+alpha = 0.01;
+beta = 0.5;  
+
+if (~largescale), AtA = A'*A; end
+
+% initial point
+x = x0;
+u = u0;
+if (largescale), r = A(x) - b; else  r = A*x - b; end
+fu1 = x - u;
+fu2 = -x - u;
+fe = 1/2*(r'*r - epsilon^2);
+f = sum(u) - (1/tau)*(sum(log(-fu1)) + sum(log(-fu2)) + log(-fe));
+
+niter = 0;
+done = 0;
+while (~done)
+  
+  if (largescale), atr = At(r); else  atr = A'*r; end
+  
+  ntgz = 1./fu1 - 1./fu2 + 1/fe*atr;
+  ntgu = -tau - 1./fu1 - 1./fu2;
+  gradf = -(1/tau)*[ntgz; ntgu];
+  
+  sig11 = 1./fu1.^2 + 1./fu2.^2;
+  sig12 = -1./fu1.^2 + 1./fu2.^2;
+  sigx = sig11 - sig12.^2./sig11;
+    
+  w1p = ntgz - sig12./sig11.*ntgu;
+  if (largescale)
+    h11pfun = @(z) sigx.*z - (1/fe)*At(A(z)) + 1/fe^2*(atr'*z)*atr;
+    [dx, cgres, cgiter] = cgsolve(h11pfun, w1p, cgtol, cgmaxiter, 0);
+    if (cgres > 1/2)
+      disp('Cannot solve system.  Returning previous iterate.  (See Section 4 of notes for more information.)');
+      xp = x;  up = u;
+      return
+    end
+    Adx = A(dx);
+  else
+    H11p = diag(sigx) - (1/fe)*AtA + (1/fe)^2*atr*atr';
+    opts.POSDEF = true; opts.SYM = true;
+    [dx,hcond] = linsolve(H11p, w1p, opts);
+    if (hcond < 1e-14)
+      disp('Matrix ill-conditioned.  Returning previous iterate.  (See Section 4 of notes for more information.)');
+      xp = x;  up = u;
+      return
+    end
+    Adx = A*dx;
+  end
+  du = (1./sig11).*ntgu - (sig12./sig11).*dx;  
+ 
+  % minimum step size that stays in the interior
+  ifu1 = find((dx-du) > 0); ifu2 = find((-dx-du) > 0);
+  aqe = Adx'*Adx;   bqe = 2*r'*Adx;   cqe = r'*r - epsilon^2;
+  smax = min(1,min([...
+    -fu1(ifu1)./(dx(ifu1)-du(ifu1)); -fu2(ifu2)./(-dx(ifu2)-du(ifu2)); ...
+    (-bqe+sqrt(bqe^2-4*aqe*cqe))/(2*aqe)
+    ]));
+  s = (0.99)*smax;
+  
+  % backtracking line search
+  suffdec = 0;
+  backiter = 0;
+  while (~suffdec)
+    xp = x + s*dx;  up = u + s*du;  rp = r + s*Adx;
+    fu1p = xp - up;  fu2p = -xp - up;  fep = 1/2*(rp'*rp - epsilon^2);
+    fp = sum(up) - (1/tau)*(sum(log(-fu1p)) + sum(log(-fu2p)) + log(-fep));
+    flin = f + alpha*s*(gradf'*[dx; du]);
+    suffdec = (fp <= flin);
+    s = beta*s;
+    backiter = backiter + 1;
+    if (backiter > 32)
+      disp('Stuck on backtracking line search, returning previous iterate.  (See Section 4 of notes for more information.)');
+      xp = x;  up = u;
+      return
+    end
+  end
+  
+  % set up for next iteration
+  x = xp; u = up;  r = rp;
+  fu1 = fu1p;  fu2 = fu2p;  fe = fep;  f = fp;
+  
+  lambda2 = -(gradf'*[dx; du]);
+  stepsize = s*norm([dx; du]);
+  niter = niter + 1;
+  done = (lambda2/2 < newtontol) | (niter >= newtonmaxiter);
+  
+  disp(sprintf('Newton iter = %d, Functional = %8.3f, Newton decrement = %8.3f, Stepsize = %8.3e', ...
+    niter, f, lambda2/2, stepsize));
+  if (largescale)
+    disp(sprintf('                CG Res = %8.3e, CG Iter = %d', cgres, cgiter));
+  else
+    disp(sprintf('                  H11p condition number = %8.3e', hcond));
+  end
+      
+end
+
+
+
+
diff --git a/Optimization/tvdantzig_logbarrier.m b/Optimization/tvdantzig_logbarrier.m
new file mode 100644
index 0000000..39c9463
--- /dev/null
+++ b/Optimization/tvdantzig_logbarrier.m
@@ -0,0 +1,120 @@
+% tvdantzig_logbarrier.m
+%
+% Solve the total variation Dantzig program
+%
+% min_x TV(x)  subject to  ||A'(Ax-b)||_\infty <= epsilon
+%
+% Recast as the SOCP
+% min sum(t) s.t.  ||D_{ij}x||_2 <= t,  i,j=1,...,n
+%                  <a_{ij},Ax - b> <= epsilon  i,j=1,...,n
+% and use a log barrier algorithm.
+%
+% Usage:  xp = tvdantzig_logbarrier(x0, A, At, b, epsilon, lbtol, mu, cgtol, cgmaxiter)
+%
+% x0 - Nx1 vector, initial point.
+%
+% A - Either a handle to a function that takes a N vector and returns a K 
+%     vector , or a KxN matrix.  If A is a function handle, the algorithm
+%     operates in "largescale" mode, solving the Newton systems via the
+%     Conjugate Gradients algorithm.
+%
+% At - Handle to a function that takes a K vector and returns an N vector.
+%      If A is a KxN matrix, At is ignored.
+%
+% b - Kx1 vector of observations.
+%
+% epsilon - scalar, constraint relaxation parameter
+%
+% lbtol - The log barrier algorithm terminates when the duality gap <= lbtol.
+%         Also, the number of log barrier iterations is completely
+%         determined by lbtol.
+%         Default = 1e-3.
+%
+% mu - Factor by which to increase the barrier constant at each iteration.
+%      Default = 10.
+%
+% cgtol - Tolerance for Conjugate Gradients; ignored if A is a matrix.
+%     Default = 1e-8.
+%
+% cgmaxiter - Maximum number of iterations for Conjugate Gradients; ignored
+%     if A is a matrix.
+%     Default = 200.
+%
+% Written by: Justin Romberg, Caltech
+% Email: jrom@acm.caltech.edu
+% Created: October 2005
+%
+
+function xp = tvdantzig_logbarrier(x0, A, At, b, epsilon, lbtol, mu, cgtol, cgmaxiter)  
+
+largescale = isa(A,'function_handle');
+
+if (nargin < 6), lbtol = 1e-3; end
+if (nargin < 7), mu = 10; end
+if (nargin < 8), cgtol = 1e-8; end
+if (nargin < 9), cgmaxiter = 200; end
+
+newtontol = lbtol;
+newtonmaxiter = 50;
+
+N = length(x0);
+n = round(sqrt(N));
+
+% create (sparse) differencing matrices for TV
+Dv = spdiags([reshape([-ones(n-1,n); zeros(1,n)],N,1) ...
+  reshape([zeros(1,n); ones(n-1,n)],N,1)], [0 1], N, N);
+Dh = spdiags([reshape([-ones(n,n-1) zeros(n,1)],N,1) ...
+  reshape([zeros(n,1) ones(n,n-1)],N,1)], [0 n], N, N);
+
+if (largescale)
+  if (norm(A(x0)-b) > epsilon)
+    disp('Starting point infeasible; using x0 = At*inv(AAt)*y.');
+    AAt = @(z) A(At(z));
+    [w, cgres] = cgsolve(AAt, b, cgtol, cgmaxiter, 0);
+    if (cgres > 1/2)
+      disp('A*At is ill-conditioned: cannot find starting point');
+      xp = x0;
+      return;
+    end
+    x0 = At(w);
+  end
+else
+  if (norm(A*x0-b) > epsilon)
+    disp('Starting point infeasible; using x0 = At*inv(AAt)*y.');
+    opts.POSDEF = true; opts.SYM = true;
+    [w, hcond] = linsolve(A*A', b, opts);
+    if (hcond < 1e-14)
+      disp('A*At is ill-conditioned: cannot find starting point');
+      xp = x0;
+      return;
+    end
+    x0 = A'*w;
+  end  
+end
+x = x0;
+Dhx = Dh*x;  Dvx = Dv*x;
+t = 1.05*sqrt(Dhx.^2 + Dvx.^2) + .01*max(sqrt(Dhx.^2 + Dvx.^2));
+
+% choose initial value of tau so that the duality gap after the first
+% step will be about the origial TV
+tau = 3*N/sum(sqrt(Dhx.^2+Dvx.^2));
+
+lbiter = ceil((log(3*N)-log(lbtol)-log(tau))/log(mu));
+disp(sprintf('Number of log barrier iterations = %d\n', lbiter));
+totaliter = 0;
+for ii = 1:lbiter
+  
+  [xp, tp, ntiter] = tvdantzig_newton(x, t, A, At, b, epsilon, tau, newtontol, newtonmaxiter, cgtol, cgmaxiter);
+  totaliter = totaliter + ntiter;
+  tvxp = sum(sqrt((Dh*xp).^2 + (Dv*xp).^2));
+  
+  disp(sprintf('\nLog barrier iter = %d, TV = %.3f, functional = %8.3f, tau = %8.3e, total newton iter = %d\n', ...
+    ii, tvxp, sum(tp), tau, totaliter));
+  
+  x = xp;
+  t = tp;
+  
+  tau = mu*tau;
+  
+end
+                   
\ No newline at end of file
diff --git a/Optimization/tvdantzig_newton.m b/Optimization/tvdantzig_newton.m
new file mode 100644
index 0000000..68d7148
--- /dev/null
+++ b/Optimization/tvdantzig_newton.m
@@ -0,0 +1,185 @@
+% tvdantzig_newton.m
+%
+% Newton iterations for TV Dantzig log-barrier subproblem.
+%
+% Usage : [xp, tp, niter] = tvdantzig_newton(x0, t0, A, At, b, epsilon, tau, 
+%                                          newtontol, newtonmaxiter, cgtol, cgmaxiter)
+%
+% x0,t0 - Nx1 vectors, initial points.
+%
+% A - Either a handle to a function that takes a N vector and returns a K 
+%     vector , or a KxN matrix.  If A is a function handle, the algorithm
+%     operates in "largescale" mode, solving the Newton systems via the
+%     Conjugate Gradients algorithm.
+%
+% At - Handle to a function that takes a K vector and returns an N vector.
+%      If A is a KxN matrix, At is ignored.
+%
+% b - Kx1 vector of observations.
+%
+% epsilon - scalar, constraint relaxation parameter
+%
+% tau - Log barrier parameter.
+%
+% newtontol - Terminate when the Newton decrement is <= newtontol.
+%
+% newtonmaxiter - Maximum number of iterations.
+%
+% cgtol - Tolerance for Conjugate Gradients; ignored if A is a matrix.
+%
+% cgmaxiter - Maximum number of iterations for Conjugate Gradients; ignored
+%     if A is a matrix.
+%
+% Written by: Justin Romberg, Caltech
+% Email: jrom@acm.caltech.edu
+% Created: October 2005
+%
+
+
+function [xp, tp, niter] = tvdantzig_newton(x0, t0, A, At, b, epsilon, tau, newtontol, newtonmaxiter, cgtol, cgmaxiter) 
+
+largescale = isa(A,'function_handle'); 
+
+alpha = 0.01;
+beta = 0.5;  
+
+N = length(x0);
+n = round(sqrt(N));
+
+% create (sparse) differencing matrices for TV
+Dv = spdiags([reshape([-ones(n-1,n); zeros(1,n)],N,1) ...
+  reshape([zeros(1,n); ones(n-1,n)],N,1)], [0 1], N, N);
+Dh = spdiags([reshape([-ones(n,n-1) zeros(n,1)],N,1) ...
+  reshape([zeros(n,1) ones(n,n-1)],N,1)], [0 n], N, N);
+
+% initial point
+x = x0;
+t = t0;
+if (largescale)
+  r = A(x) - b;
+  Atr = At(r);
+else  
+  AtA = A'*A;
+  r = A*x - b;
+  Atr = A'*r;
+end  
+Dhx = Dh*x;  Dvx = Dv*x;
+ft = 1/2*(Dhx.^2 + Dvx.^2 - t.^2);
+fe1 = Atr - epsilon;
+fe2 = -Atr - epsilon;
+f = sum(t) - (1/tau)*(sum(log(-ft)) + sum(log(-fe1)) + sum(log(-fe2)));
+
+niter = 0;
+done = 0;
+while (~done)
+  
+  if (largescale)
+    ntgx = Dh'*((1./ft).*Dhx) + Dv'*((1./ft).*Dvx) + At(A(1./fe1-1./fe2));
+  else
+    ntgx = Dh'*((1./ft).*Dhx) + Dv'*((1./ft).*Dvx) + AtA*(1./fe1-1./fe2);
+  end
+  ntgt = -tau - t./ft;
+  gradf = -(1/tau)*[ntgx; ntgt];
+  
+  sig22 = 1./ft + (t.^2)./(ft.^2);
+  sig12 = -t./ft.^2;
+  sigb = 1./ft.^2 - (sig12.^2)./sig22;
+  siga = 1./fe1.^2 + 1./fe2.^2;
+  
+  w11 = ntgx - Dh'*(Dhx.*(sig12./sig22).*ntgt) - Dv'*(Dvx.*(sig12./sig22).*ntgt);
+  if (largescale)
+    h11pfun = @(w) H11p(w, A, At, Dh, Dv, Dhx, Dvx, sigb, ft, siga);
+    [dx, cgres, cgiter] = cgsolve(h11pfun, w11, cgtol, cgmaxiter, 0);
+    if (cgres > 1/2)
+      disp('Cannot solve system.  Returning previous iterate. (See Section 4 of notes for more information.)');
+      xp = x;  tp = t;
+      return
+    end
+    Adx = A(dx);
+    AtAdx = At(Adx);
+  else
+    H11p =  Dh'*sparse(diag(-1./ft + sigb.*Dhx.^2))*Dh + ...
+      Dv'*sparse(diag(-1./ft + sigb.*Dvx.^2))*Dv + ...
+      Dh'*sparse(diag(sigb.*Dhx.*Dvx))*Dv + ...
+      Dv'*sparse(diag(sigb.*Dhx.*Dvx))*Dh + ...
+      AtA*sparse(diag(siga))*AtA;
+    opts.POSDEF = true; opts.SYM = true;
+    [dx,hcond] = linsolve(H11p, w11, opts);
+    if (hcond < 1e-14)
+      disp('Matrix ill-conditioned.  Returning previous iterate. (See Section 4 of notes for more information.)');
+      xp = x;  tp = t;
+      return
+    end
+    Adx = A*dx;
+    AtAdx = A'*Adx;
+  end
+  Dhdx = Dh*dx;  Dvdx = Dv*dx;
+  dt = (1./sig22).*(ntgt - sig12.*(Dhx.*Dhdx + Dvx.*Dvdx));
+  
+  % minimum step size that stays in the interior
+  ife1 = find(AtAdx > 0); ife2 = find(-AtAdx > 0); 
+  aqt = Dhdx.^2 + Dvdx.^2 - dt.^2;   
+  bqt = 2*(Dhdx.*Dhx + Dvdx.*Dvx - t.*dt);  
+  cqt = Dhx.^2 + Dvx.^2 - t.^2;
+  tsols = [(-bqt+sqrt(bqt.^2-4*aqt.*cqt))./(2*aqt); ...
+    (-bqt-sqrt(bqt.^2-4*aqt.*cqt))./(2*aqt) ];
+  indt = find([(bqt.^2 > 4*aqt.*cqt); (bqt.^2 > 4*aqt.*cqt)] & (tsols > 0));
+  smax = min(1, min([-fe1(ife1)./AtAdx(ife1); -fe2(ife2)./(-AtAdx(ife2)); tsols(indt)])); 
+  s = (0.99)*smax;
+  
+  % backtracking line search
+  suffdec = 0;
+  backiter = 0;
+  while (~suffdec)
+    xp = x + s*dx;  tp = t + s*dt;
+    rp = r + s*Adx;  Atrp = Atr + s*AtAdx;
+    Dhxp = Dhx + s*Dhdx;  Dvxp = Dvx + s*Dvdx;
+    ftp = 1/2*(Dhxp.^2 + Dvxp.^2 - tp.^2);
+    fe1p = Atrp - epsilon;
+    fe2p = -Atrp - epsilon;
+    fp = sum(tp) - (1/tau)*(sum(log(-ftp)) + sum(log(-fe1p)) + sum(log(-fe2p)));
+    flin = f + alpha*s*(gradf'*[dx; dt]);
+    suffdec = (fp <= flin);
+    s = beta*s;
+    backiter = backiter + 1;
+    if (backiter > 32)
+      disp('Stuck backtracking, returning previous iterate. (See Section 4 of notes for more information.)');
+      xp = x;  tp = t;
+      return
+    end
+  end
+  
+  % set up for next iteration
+  x = xp; t = tp;
+  r = rp;  Atr = Atrp;
+  Dvx = Dvxp;  Dhx = Dhxp; 
+  ft = ftp; fe1 = fe1p; fe2 = fe2p;  f = fp;
+  
+  lambda2 = -(gradf'*[dx; dt]);
+  stepsize = s*norm([dx; dt]);
+  niter = niter + 1;
+  done = (lambda2/2 < newtontol) | (niter >= newtonmaxiter);
+  
+  disp(sprintf('Newton iter = %d, Functional = %8.3f, Newton decrement = %8.3f, Stepsize = %8.3e', ...
+    niter, f, lambda2/2, stepsize));
+  if (largescale)
+    disp(sprintf('                  CG Res = %8.3e, CG Iter = %d', cgres, cgiter));
+  else
+    disp(sprintf('                  H11p condition number = %8.3e', hcond));
+  end
+      
+end
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% H11p auxiliary function
+function y = H11p(v, A, At, Dh, Dv, Dhx, Dvx, sigb, ft, siga)
+
+Dhv = Dh*v;
+Dvv = Dv*v;
+
+y = Dh'*((-1./ft + sigb.*Dhx.^2).*Dhv + sigb.*Dhx.*Dvx.*Dvv) + ...
+  Dv'*((-1./ft + sigb.*Dvx.^2).*Dvv + sigb.*Dhx.*Dvx.*Dhv) + ...
+  At(A(siga.*At(A(v)))); 
+                                                                                                                           
+
diff --git a/Optimization/tveq_logbarrier.m b/Optimization/tveq_logbarrier.m
new file mode 100644
index 0000000..617bf2e
--- /dev/null
+++ b/Optimization/tveq_logbarrier.m
@@ -0,0 +1,118 @@
+% tveq_logbarrier.m
+%
+% Solve equality constrained TV minimization
+% min TV(x)  s.t.  Ax=b.
+%
+% Recast as the SOCP
+% min sum(t) s.t.  ||D_{ij}x||_2 <= t,  i,j=1,...,n
+%                  Ax=b
+% and use a log barrier algorithm.
+%
+% Usage:  xp = tveq_logbarrier(x0, A, At, b, lbtol, mu, slqtol, slqmaxiter)
+%
+% x0 - Nx1 vector, initial point.
+%
+% A - Either a handle to a function that takes a N vector and returns a K 
+%     vector , or a KxN matrix.  If A is a function handle, the algorithm
+%     operates in "largescale" mode, solving the Newton systems via the
+%     Conjugate Gradients algorithm.
+%
+% At - Handle to a function that takes a K vector and returns an N vector.
+%      If A is a KxN matrix, At is ignored.
+%
+% b - Kx1 vector of observations.
+%
+% lbtol - The log barrier algorithm terminates when the duality gap <= lbtol.
+%         Also, the number of log barrier iterations is completely
+%         determined by lbtol.
+%         Default = 1e-3.
+%
+% mu - Factor by which to increase the barrier constant at each iteration.
+%      Default = 10.
+%
+% slqtol - Tolerance for SYMMLQ; ignored if A is a matrix.
+%     Default = 1e-8.
+%
+% slqmaxiter - Maximum number of iterations for SYMMLQ; ignored
+%     if A is a matrix.
+%     Default = 200.
+%
+% Written by: Justin Romberg, Caltech
+% Email: jrom@acm.caltech.edu
+% Created: October 2005
+%
+
+function xp = tveq_logbarrier(x0, A, At, b, lbtol, mu, slqtol, slqmaxiter)  
+
+largescale = isa(A,'function_handle'); 
+
+if (nargin < 5), lbtol = 1e-3; end
+if (nargin < 6), mu = 10; end
+if (nargin < 7), slqtol = 1e-8; end
+if (nargin < 8), slqmaxiter = 200; end
+
+newtontol = lbtol;
+newtonmaxiter = 50;
+
+N = length(x0);
+n = round(sqrt(N));
+
+% create (sparse) differencing matrices for TV
+Dv = spdiags([reshape([-ones(n-1,n); zeros(1,n)],N,1) ...
+  reshape([zeros(1,n); ones(n-1,n)],N,1)], [0 1], N, N);
+Dh = spdiags([reshape([-ones(n,n-1) zeros(n,1)],N,1) ...
+  reshape([zeros(n,1) ones(n,n-1)],N,1)], [0 n], N, N);
+
+% starting point --- make sure that it is feasible
+if (largescale)
+  if (norm(A(x0)-b)/norm(b) > slqtol)
+    disp('Starting point infeasible; using x0 = At*inv(AAt)*y.');
+    AAt = @(z) A(At(z));
+    [w,cgres] = cgsolve(AAt, b, slqtol, slqmaxiter, 0);
+    if (cgres > 1/2)
+      disp('A*At is ill-conditioned: cannot find starting point');
+      xp = x0;
+      return;
+    end
+    x0 = At(w);
+  end
+else
+  if (norm(A*x0-b)/norm(b) > slqtol)
+    disp('Starting point infeasible; using x0 = At*inv(AAt)*y.');
+    opts.POSDEF = true; opts.SYM = true;
+    [w, hcond] = linsolve(A*A', b, opts);
+    if (hcond < 1e-14)
+      disp('A*At is ill-conditioned: cannot find starting point');
+      xp = x0;
+      return;
+    end
+    x0 = A'*w;
+  end  
+end
+x = x0;
+Dhx = Dh*x;  Dvx = Dv*x;
+t = (0.95)*sqrt(Dhx.^2 + Dvx.^2) + (0.1)*max(sqrt(Dhx.^2 + Dvx.^2));
+
+% choose initial value of tau so that the duality gap after the first
+% step will be about the origial TV
+tau = N/sum(sqrt(Dhx.^2+Dvx.^2));
+
+lbiter = ceil((log(N)-log(lbtol)-log(tau))/log(mu));
+disp(sprintf('Number of log barrier iterations = %d\n', lbiter));
+totaliter = 0;
+for ii = 1:lbiter
+  
+  [xp, tp, ntiter] = tveq_newton(x, t, A, At, b, tau, newtontol, newtonmaxiter, slqtol, slqmaxiter);
+  totaliter = totaliter + ntiter;
+  
+  tvxp = sum(sqrt((Dh*xp).^2 + (Dv*xp).^2));
+  disp(sprintf('\nLog barrier iter = %d, TV = %.3f, functional = %8.3f, tau = %8.3e, total newton iter = %d\n', ...
+    ii, tvxp, sum(tp), tau, totaliter));
+  
+  x = xp;
+  t = tp;
+  
+  tau = mu*tau;
+  
+end
+                   
\ No newline at end of file
diff --git a/Optimization/tveq_newton.m b/Optimization/tveq_newton.m
new file mode 100644
index 0000000..9e71b73
--- /dev/null
+++ b/Optimization/tveq_newton.m
@@ -0,0 +1,180 @@
+% tveq_newton.m
+%
+% Newton algorithm for log-barrier subproblems for TV minimization
+% with equality constraints.
+%
+% Usage: 
+% [xp,tp,niter] = tveq_newton(x0, t0, A, At, b, tau, 
+%                             newtontol, newtonmaxiter, slqtol, slqmaxiter)
+%
+% x0,t0 - starting points
+%
+% A - Either a handle to a function that takes a N vector and returns a K 
+%     vector , or a KxN matrix.  If A is a function handle, the algorithm
+%     operates in "largescale" mode, solving the Newton systems via the
+%     Conjugate Gradients algorithm.
+%
+% At - Handle to a function that takes a K vector and returns an N vector.
+%      If A is a KxN matrix, At is ignored.
+%
+% b - Kx1 vector of observations.
+%
+% tau - Log barrier parameter.
+%
+% newtontol - Terminate when the Newton decrement is <= newtontol.
+%
+% newtonmaxiter - Maximum number of iterations.
+%
+% slqtol - Tolerance for SYMMLQ; ignored if A is a matrix.
+%
+% slqmaxiter - Maximum number of iterations for SYMMLQ; ignored
+%     if A is a matrix.
+%
+% Written by: Justin Romberg, Caltech
+% Email: jrom@acm.caltech.edu
+% Created: October 2005
+%
+
+function [xp, tp, niter] = tveq_newton(x0, t0, A, At, b, tau, newtontol, newtonmaxiter, slqtol, slqmaxiter) 
+
+largescale = isa(A,'function_handle'); 
+
+alpha = 0.01;
+beta = 0.5;  
+
+N = length(x0);
+n = round(sqrt(N));
+K = length(b);
+
+% create (sparse) differencing matrices for TV
+Dv = spdiags([reshape([-ones(n-1,n); zeros(1,n)],N,1) ...
+  reshape([zeros(1,n); ones(n-1,n)],N,1)], [0 1], N, N);
+Dh = spdiags([reshape([-ones(n,n-1) zeros(n,1)],N,1) ...
+  reshape([zeros(n,1) ones(n,n-1)],N,1)], [0 n], N, N);
+
+% auxillary matrices for preconditioning
+Mdv = spdiags([reshape([ones(n-1,n); zeros(1,n)],N,1) ...
+  reshape([zeros(1,n); ones(n-1,n)],N,1)], [0 1], N, N);
+Mdh = spdiags([reshape([ones(n,n-1) zeros(n,1)],N,1) ...
+  reshape([zeros(n,1) ones(n,n-1)],N,1)], [0 n], N, N);
+Mmd = reshape([ones(n-1,n-1) zeros(n-1,1); zeros(1,n)],N,1);
+
+
+% initial point
+x = x0;
+t = t0;
+Dhx = Dh*x;  Dvx = Dv*x;
+ft = 1/2*(Dhx.^2 + Dvx.^2 - t.^2);
+f = sum(t) - (1/tau)*(sum(log(-ft)));
+
+niter = 0;
+done = 0;
+while (~done)
+  
+  ntgx = Dh'*((1./ft).*Dhx) + Dv'*((1./ft).*Dvx);
+  ntgt = -tau - t./ft;
+  gradf = -(1/tau)*[ntgx; ntgt];
+   
+  sig22 = 1./ft + (t.^2)./(ft.^2);
+  sig12 = -t./ft.^2;
+  sigb = 1./ft.^2 - (sig12.^2)./sig22;
+  
+  w1p = ntgx - Dh'*(Dhx.*(sig12./sig22).*ntgt) - Dv'*(Dvx.*(sig12./sig22).*ntgt);
+  wp = [w1p; zeros(K,1)];
+  if (largescale)
+    % diagonal of H11p
+    dg11p = Mdh'*(-1./ft + sigb.*Dhx.^2) + Mdv'*(-1./ft + sigb.*Dvx.^2) + 2*Mmd.*sigb.*Dhx.*Dvx;
+    afac = max(dg11p);
+    hpfun = @(z) Hpeval(z, A, At, Dh, Dv, Dhx, Dvx, sigb, ft, afac);
+    [dxv,slqflag,slqres,slqiter] = symmlq(hpfun, wp, slqtol, slqmaxiter);
+    if (slqres > 1/2)
+      disp('Cannot solve system.  Returning previous iterate.  (See Section 4 of notes for more information.)');
+      xp = x;
+      return
+    end
+  else
+    H11p =  Dh'*sparse(diag(-1./ft + sigb.*Dhx.^2))*Dh + ...
+      Dv'*sparse(diag(-1./ft + sigb.*Dvx.^2))*Dv + ...
+      Dh'*sparse(diag(sigb.*Dhx.*Dvx))*Dv + ...
+      Dv'*sparse(diag(sigb.*Dhx.*Dvx))*Dh;
+    afac = max(diag(H11p));
+    Hp = full([H11p afac*A'; afac*A zeros(K)]);
+    %keyboard
+    opts.SYM = true;
+    [dxv, hcond] = linsolve(Hp, wp, opts);
+    if (hcond < 1e-14)
+      disp('Matrix ill-conditioned.  Returning previous iterate.  (See Section 4 of notes for more information.)');
+      xp = x;  tp = t;
+      return
+    end
+  end
+  dx = dxv(1:N);  
+  Dhdx = Dh*dx;  Dvdx = Dv*dx;
+  dt = (1./sig22).*(ntgt - sig12.*(Dhx.*Dhdx + Dvx.*Dvdx));
+  
+  % minimum step size that stays in the interior
+  aqt = Dhdx.^2 + Dvdx.^2 - dt.^2;   
+  bqt = 2*(Dhdx.*Dhx + Dvdx.*Dvx - t.*dt);  
+  cqt = Dhx.^2 + Dvx.^2 - t.^2;
+  tsols = [(-bqt+sqrt(bqt.^2-4*aqt.*cqt))./(2*aqt); ...
+    (-bqt-sqrt(bqt.^2-4*aqt.*cqt))./(2*aqt) ];
+  indt = find([(bqt.^2 > 4*aqt.*cqt); (bqt.^2 > 4*aqt.*cqt)] & (tsols > 0));
+  smax = min(1, min(tsols(indt))); 
+  s = (0.99)*smax;
+  
+  % line search
+  suffdec = 0;
+  backiter = 0;
+  while (~suffdec)
+    xp = x + s*dx;  tp = t + s*dt;
+    Dhxp = Dhx + s*Dhdx;  Dvxp = Dvx + s*Dvdx;
+    ftp = 1/2*(Dhxp.^2 + Dvxp.^2 - tp.^2);
+    fp = sum(tp) - (1/tau)*(sum(log(-ftp)));
+    flin = f + alpha*s*(gradf'*[dx; dt]);
+    suffdec = (fp <= flin);
+    s = beta*s;
+    backiter = backiter + 1;
+    if (backiter > 32)
+      disp('Stuck backtracking, returning last iterate. (See Section 4 of notes for more information.)');
+      xp = x;  tp = t;
+      return
+    end
+  end
+  
+  % set up for next iteration
+  x = xp; t = tp;
+  Dvx = Dvxp;  Dhx = Dhxp; 
+  ft = ftp; f = fp;
+  
+  lambda2 = -(gradf'*[dx; dt]);
+  stepsize = s*norm([dx; dt]);
+  niter = niter + 1;
+  done = (lambda2/2 < newtontol) | (niter >= newtonmaxiter);
+
+  disp(sprintf('Newton iter = %d, Functional = %8.3f, Newton decrement = %8.3f, Stepsize = %8.3e', ...
+    niter, f, lambda2/2, stepsize));
+  if (largescale)
+    disp(sprintf('                  SYMMLQ Res = %8.3e, SYMMLQ Iter = %d', slqres, slqiter));
+  else
+    disp(sprintf('                  H11p condition number = %8.3e', hcond));
+  end
+    
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Implicit application of Hessian
+function y = Hpeval(z, A, At, Dh, Dv, Dhx, Dvx, sigb, ft, afac)
+
+N = length(ft);
+K = length(z)-N;
+w = z(1:N);
+v = z(N+1:N+K);
+
+Dhw = Dh*w;
+Dvw = Dv*w;
+
+y1 = Dh'*((-1./ft + sigb.*Dhx.^2).*Dhw + sigb.*Dhx.*Dvx.*Dvw) + ...
+  Dv'*((-1./ft + sigb.*Dvx.^2).*Dvw + sigb.*Dhx.*Dvx.*Dhw) + afac*At(v);
+y2 = afac*A(w);
+
+y = [y1; y2];
diff --git a/Optimization/tvqc_logbarrier.m b/Optimization/tvqc_logbarrier.m
new file mode 100644
index 0000000..52cc6b1
--- /dev/null
+++ b/Optimization/tvqc_logbarrier.m
@@ -0,0 +1,121 @@
+% tvqc_logbarrier.m
+%
+% Solve quadractically constrained TV minimization
+% min TV(x)  s.t.  ||Ax-b||_2 <= epsilon.
+%
+% Recast as the SOCP
+% min sum(t) s.t.  ||D_{ij}x||_2 <= t,  i,j=1,...,n
+%                  ||Ax - b||_2 <= epsilon
+% and use a log barrier algorithm.
+%
+% Usage:  xp = tvqc_logbarrier(x0, A, At, b, epsilon, lbtol, mu, cgtol, cgmaxiter)
+%
+% x0 - Nx1 vector, initial point.
+%
+% A - Either a handle to a function that takes a N vector and returns a K 
+%     vector , or a KxN matrix.  If A is a function handle, the algorithm
+%     operates in "largescale" mode, solving the Newton systems via the
+%     Conjugate Gradients algorithm.
+%
+% At - Handle to a function that takes a K vector and returns an N vector.
+%      If A is a KxN matrix, At is ignored.
+%
+% b - Kx1 vector of observations.
+%
+% epsilon - scalar, constraint relaxation parameter
+%
+% lbtol - The log barrier algorithm terminates when the duality gap <= lbtol.
+%         Also, the number of log barrier iterations is completely
+%         determined by lbtol.
+%         Default = 1e-3.
+%
+% mu - Factor by which to increase the barrier constant at each iteration.
+%      Default = 10.
+%
+% cgtol - Tolerance for Conjugate Gradients; ignored if A is a matrix.
+%     Default = 1e-8.
+%
+% cgmaxiter - Maximum number of iterations for Conjugate Gradients; ignored
+%     if A is a matrix.
+%     Default = 200.
+%
+% Written by: Justin Romberg, Caltech
+% Email: jrom@acm.caltech.edu
+% Created: October 2005
+%
+
+
+function [xp, tp] = tvqc_logbarrier(x0, A, At, b, epsilon, lbtol, mu, cgtol, cgmaxiter)  
+
+largescale = isa(A,'function_handle'); 
+
+if (nargin < 6), lbtol = 1e-3; end
+if (nargin < 7), mu = 10; end
+if (nargin < 8), cgtol = 1e-8; end
+if (nargin < 9), cgmaxiter = 200; end
+
+newtontol = lbtol;
+newtonmaxiter = 50;
+
+N = length(x0);
+n = round(sqrt(N));
+
+% create (sparse) differencing matrices for TV
+Dv = spdiags([reshape([-ones(n-1,n); zeros(1,n)],N,1) ...
+  reshape([zeros(1,n); ones(n-1,n)],N,1)], [0 1], N, N);
+Dh = spdiags([reshape([-ones(n,n-1) zeros(n,1)],N,1) ...
+  reshape([zeros(n,1) ones(n,n-1)],N,1)], [0 n], N, N);
+
+% starting point --- make sure that it is feasible
+if (largescale)
+  if (norm(A(x0)-b) > epsilon)
+    disp('Starting point infeasible; using x0 = At*inv(AAt)*y.');
+    AAt = @(z) A(At(z));
+    [w, cgres] = cgsolve(AAt, b, cgtol, cgmaxiter, 0);
+    if (cgres > 1/2)
+      disp('A*At is ill-conditioned: cannot find starting point');
+      xp = x0;
+      return;
+    end
+    x0 = At(w);
+  end
+else
+  if (norm(A*x0-b) > epsilon)
+    disp('Starting point infeasible; using x0 = At*inv(AAt)*y.');
+    opts.POSDEF = true; opts.SYM = true;
+    [w, hcond] = linsolve(A*A', b, opts);
+    if (hcond < 1e-14)
+      disp('A*At is ill-conditioned: cannot find starting point');
+      xp = x0;
+      return;
+    end
+    x0 = A'*w;
+  end  
+end
+x = x0;
+Dhx = Dh*x;  Dvx = Dv*x;
+t = (0.95)*sqrt(Dhx.^2 + Dvx.^2) + (0.1)*max(sqrt(Dhx.^2 + Dvx.^2));
+
+% choose initial value of tau so that the duality gap after the first
+% step will be about the origial TV
+tau = (N+1)/sum(sqrt(Dhx.^2+Dvx.^2));
+
+lbiter = ceil((log((N+1))-log(lbtol)-log(tau))/log(mu));
+disp(sprintf('Number of log barrier iterations = %d\n', lbiter));
+totaliter = 0;
+for ii = 1:lbiter
+  
+  [xp, tp, ntiter] = tvqc_newton(x, t, A, At, b, epsilon, tau, newtontol, newtonmaxiter, cgtol, cgmaxiter);
+  totaliter = totaliter + ntiter;
+  
+  tvxp = sum(sqrt((Dh*xp).^2 + (Dv*xp).^2));
+  disp(sprintf('\nLog barrier iter = %d, TV = %.3f, functional = %8.3f, tau = %8.3e, total newton iter = %d\n', ...
+    ii, tvxp, sum(tp), tau, totaliter));
+  
+  x = xp;
+  t = tp;
+  
+  tau = mu*tau;
+  
+end
+                   
diff --git a/Optimization/tvqc_newton.m b/Optimization/tvqc_newton.m
new file mode 100644
index 0000000..febe8ff
--- /dev/null
+++ b/Optimization/tvqc_newton.m
@@ -0,0 +1,176 @@
+% tvqc_newton.m
+%
+% Newton algorithm for log-barrier subproblems for TV minimization
+% with quadratic constraints.
+%
+% Usage: 
+% [xp,tp,niter] = tvqc_newton(x0, t0, A, At, b, epsilon, tau, 
+%                             newtontol, newtonmaxiter, cgtol, cgmaxiter)
+%
+% x0,t0 - starting points
+%
+% A - Either a handle to a function that takes a N vector and returns a K 
+%     vector , or a KxN matrix.  If A is a function handle, the algorithm
+%     operates in "largescale" mode, solving the Newton systems via the
+%     Conjugate Gradients algorithm.
+%
+% At - Handle to a function that takes a K vector and returns an N vector.
+%      If A is a KxN matrix, At is ignored.
+%
+% b - Kx1 vector of observations.
+%
+% epsilon - scalar, constraint relaxation parameter
+%
+% tau - Log barrier parameter.
+%
+% newtontol - Terminate when the Newton decrement is <= newtontol.
+%
+% newtonmaxiter - Maximum number of iterations.
+%
+% cgtol - Tolerance for Conjugate Gradients; ignored if A is a matrix.
+%
+% cgmaxiter - Maximum number of iterations for Conjugate Gradients; ignored
+%     if A is a matrix.
+%
+% Written by: Justin Romberg, Caltech
+% Email: jrom@acm.caltech.edu
+% Created: October 2005
+%
+
+function [xp, tp, niter] = tvqc_newton(x0, t0, A, At, b, epsilon, tau, newtontol, newtonmaxiter, cgtol, cgmaxiter) 
+
+largescale = isa(A,'function_handle'); 
+
+alpha = 0.01;
+beta = 0.5;  
+
+N = length(x0);
+n = round(sqrt(N));
+
+% create (sparse) differencing matrices for TV
+Dv = spdiags([reshape([-ones(n-1,n); zeros(1,n)],N,1) ...
+  reshape([zeros(1,n); ones(n-1,n)],N,1)], [0 1], N, N);
+Dh = spdiags([reshape([-ones(n,n-1) zeros(n,1)],N,1) ...
+  reshape([zeros(n,1) ones(n,n-1)],N,1)], [0 n], N, N);
+
+if (~largescale),  AtA = A'*A;  end;
+
+% initial point
+x = x0;
+t = t0;
+if (largescale), r = A(x) - b;  else  r = A*x - b; end  
+Dhx = Dh*x;  Dvx = Dv*x;
+ft = 1/2*(Dhx.^2 + Dvx.^2 - t.^2);
+fe = 1/2*(r'*r - epsilon^2);
+f = sum(t) - (1/tau)*(sum(log(-ft)) + log(-fe));
+
+niter = 0;
+done = 0;
+while (~done)
+  
+  if (largescale),  Atr = At(r);  else  Atr = A'*r;  end
+  ntgx = Dh'*((1./ft).*Dhx) + Dv'*((1./ft).*Dvx) + 1/fe*Atr;
+  ntgt = -tau - t./ft;
+  gradf = -(1/tau)*[ntgx; ntgt];
+  
+  sig22 = 1./ft + (t.^2)./(ft.^2);
+  sig12 = -t./ft.^2;
+  sigb = 1./ft.^2 - (sig12.^2)./sig22;
+  
+  w1p = ntgx - Dh'*(Dhx.*(sig12./sig22).*ntgt) - Dv'*(Dvx.*(sig12./sig22).*ntgt);
+  if (largescale)
+    h11pfun = @(z) H11p(z, A, At, Dh, Dv, Dhx, Dvx, sigb, ft, fe, Atr);
+    [dx, cgres, cgiter] = cgsolve(h11pfun, w1p, cgtol, cgmaxiter, 0);
+    if (cgres > 1/2)
+      disp('Cannot solve system.  Returning previous iterate.  (See Section 4 of notes for more information.)');
+      xp = x;  tp = t;
+      return
+    end
+    Adx = A(dx);
+  else
+    H11p =  Dh'*sparse(diag(-1./ft + sigb.*Dhx.^2))*Dh + ...
+      Dv'*sparse(diag(-1./ft + sigb.*Dvx.^2))*Dv + ...
+      Dh'*sparse(diag(sigb.*Dhx.*Dvx))*Dv + ...
+      Dv'*sparse(diag(sigb.*Dhx.*Dvx))*Dh - ...
+      (1/fe)*AtA + (1/fe^2)*Atr*Atr';
+    opts.POSDEF = true; opts.SYM = true;
+    [dx,hcond] = linsolve(H11p, w1p, opts);
+    if (hcond < 1e-14)
+      disp('Matrix ill-conditioned.  Returning previous iterate.  (See Section 4 of notes for more information.)');
+      xp = x;  tp = t;
+      return
+    end
+    Adx = A*dx;
+  end
+  Dhdx = Dh*dx;  Dvdx = Dv*dx;
+  dt = (1./sig22).*(ntgt - sig12.*(Dhx.*Dhdx + Dvx.*Dvdx));
+  
+  % minimum step size that stays in the interior
+  aqt = Dhdx.^2 + Dvdx.^2 - dt.^2;   
+  bqt = 2*(Dhdx.*Dhx + Dvdx.*Dvx - t.*dt);  
+  cqt = Dhx.^2 + Dvx.^2 - t.^2;
+  tsols = [(-bqt+sqrt(bqt.^2-4*aqt.*cqt))./(2*aqt); ...
+    (-bqt-sqrt(bqt.^2-4*aqt.*cqt))./(2*aqt) ];
+  indt = find([(bqt.^2 > 4*aqt.*cqt); (bqt.^2 > 4*aqt.*cqt)] & (tsols > 0));
+  aqe = Adx'*Adx;   bqe = 2*r'*Adx;   cqe = r'*r - epsilon^2;
+  smax = min(1,min([...
+    tsols(indt); ...
+    (-bqe+sqrt(bqe^2-4*aqe*cqe))/(2*aqe)
+    ]));
+  s = (0.99)*smax;
+  
+  % backtracking line search
+  suffdec = 0;
+  backiter = 0;
+  while (~suffdec)
+    xp = x + s*dx;  tp = t + s*dt;
+    rp = r + s*Adx;  Dhxp = Dhx + s*Dhdx;  Dvxp = Dvx + s*Dvdx;
+    ftp = 1/2*(Dhxp.^2 + Dvxp.^2 - tp.^2);
+    fep = 1/2*(rp'*rp - epsilon^2);
+    fp = sum(tp) - (1/tau)*(sum(log(-ftp)) + log(-fep));
+    flin = f + alpha*s*(gradf'*[dx; dt]);
+    suffdec = (fp <= flin);
+    s = beta*s;
+    backiter = backiter + 1;
+    if (backiter > 32)
+      disp('Stuck on backtracking line search, returning previous iterate.  (See Section 4 of notes for more information.)');
+      xp = x;  tp = t;
+      return
+    end
+  end
+  
+  % set up for next iteration
+  x = xp; t = tp;
+  r = rp;  Dvx = Dvxp;  Dhx = Dhxp; 
+  ft = ftp; fe = fep; f = fp;
+  
+  lambda2 = -(gradf'*[dx; dt]);
+  stepsize = s*norm([dx; dt]);
+  niter = niter + 1;
+  done = (lambda2/2 < newtontol) | (niter >= newtonmaxiter);
+  
+  disp(sprintf('Newton iter = %d, Functional = %8.3f, Newton decrement = %8.3f, Stepsize = %8.3e', ...
+    niter, f, lambda2/2, stepsize));
+  if (largescale)
+    disp(sprintf('                  CG Res = %8.3e, CG Iter = %d', cgres, cgiter));
+  else
+    disp(sprintf('                  H11p condition number = %8.3e', hcond));
+  end
+ 
+end
+
+
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% H11p auxiliary function
+function y = H11p(v, A, At, Dh, Dv, Dhx, Dvx, sigb, ft, fe, atr)
+
+Dhv = Dh*v;
+Dvv = Dv*v;
+
+y = Dh'*((-1./ft + sigb.*Dhx.^2).*Dhv + sigb.*Dhx.*Dvx.*Dvv) + ...
+  Dv'*((-1./ft + sigb.*Dvx.^2).*Dvv + sigb.*Dhx.*Dvx.*Dhv) - ...
+  1/fe*At(A(v)) + 1/fe^2*(atr'*v)*atr;  
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%