% 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)  if (nargin < 6), lbtol = 1e-3; endif (nargin < 7), mu = 10; endif (nargin < 8), cgtol = 1e-8; endif (nargin < 9), cgmaxiter = 200; endnewtontol = lbtol;newtonmaxiter = 50;N = length(x0);n = round(sqrt(N));% create (sparse) differencing matrices for TVDv = 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);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 TVtau = (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