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

Additional Clean up of Mac dirs and tex generated files

1 files changed, 118 insertions, 0 deletions

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