nikcleju@62: 
nikcleju@62: import numpy
nikcleju@62: import scipy.linalg
nikcleju@62: import math
nikcleju@62: 
nikcleju@63: class l1eqNotImplementedError(Exception):
nikcleju@62:   pass
nikcleju@62: 
nikcleju@62: #function xp = l1eq_pd(x0, A, At, b, pdtol, pdmaxiter, cgtol, cgmaxiter)
nikcleju@62: def l1eq_pd(x0, A, At, b, pdtol=1e-3, pdmaxiter=50, cgtol=1e-8, cgmaxiter=200, verbose=False):
nikcleju@62: 
nikcleju@62:     # Solve
nikcleju@62:     # min_x ||x||_1  s.t.  Ax = b
nikcleju@62:     #
nikcleju@62:     # Recast as linear program
nikcleju@62:     # min_{x,u} sum(u)  s.t.  -u <= x <= u,  Ax=b
nikcleju@62:     # and use primal-dual interior point method
nikcleju@62:     #
nikcleju@62:     # Usage: xp = l1eq_pd(x0, A, At, b, pdtol, pdmaxiter, cgtol, cgmaxiter)
nikcleju@62:     #
nikcleju@62:     # x0 - Nx1 vector, initial point.
nikcleju@62:     #
nikcleju@62:     # A - Either a handle to a function that takes a N vector and returns a K 
nikcleju@62:     #     vector , or a KxN matrix.  If A is a function handle, the algorithm
nikcleju@62:     #     operates in "largescale" mode, solving the Newton systems via the
nikcleju@62:     #     Conjugate Gradients algorithm.
nikcleju@62:     #
nikcleju@62:     # At - Handle to a function that takes a K vector and returns an N vector.
nikcleju@62:     #      If A is a KxN matrix, At is ignored.
nikcleju@62:     #
nikcleju@62:     # b - Kx1 vector of observations.
nikcleju@62:     #
nikcleju@62:     # pdtol - Tolerance for primal-dual algorithm (algorithm terminates if
nikcleju@62:     #     the duality gap is less than pdtol).  
nikcleju@62:     #     Default = 1e-3.
nikcleju@62:     #
nikcleju@62:     # pdmaxiter - Maximum number of primal-dual iterations.  
nikcleju@62:     #     Default = 50.
nikcleju@62:     #
nikcleju@62:     # cgtol - Tolerance for Conjugate Gradients; ignored if A is a matrix.
nikcleju@62:     #     Default = 1e-8.
nikcleju@62:     #
nikcleju@62:     # cgmaxiter - Maximum number of iterations for Conjugate Gradients; ignored
nikcleju@62:     #     if A is a matrix.
nikcleju@62:     #     Default = 200.
nikcleju@62:     #
nikcleju@62:     # Written by: Justin Romberg, Caltech
nikcleju@62:     # Email: jrom@acm.caltech.edu
nikcleju@62:     # Created: October 2005
nikcleju@62:     
nikcleju@62:     
nikcleju@62:     #---------------------
nikcleju@62:     # Original Matab code:    
nikcleju@62:         
nikcleju@62:     #largescale = isa(A,'function_handle');
nikcleju@62:     #
nikcleju@62:     #if (nargin < 5), pdtol = 1e-3;  end
nikcleju@62:     #if (nargin < 6), pdmaxiter = 50;  end
nikcleju@62:     #if (nargin < 7), cgtol = 1e-8;  end
nikcleju@62:     #if (nargin < 8), cgmaxiter = 200;  end
nikcleju@62:     #
nikcleju@62:     #N = length(x0);
nikcleju@62:     #
nikcleju@62:     #alpha = 0.01;
nikcleju@62:     #beta = 0.5;
nikcleju@62:     #mu = 10;
nikcleju@62:     #
nikcleju@62:     #gradf0 = [zeros(N,1); ones(N,1)];
nikcleju@62:     #
nikcleju@62:     ## starting point --- make sure that it is feasible
nikcleju@62:     #if (largescale)
nikcleju@62:     #  if (norm(A(x0)-b)/norm(b) > cgtol)
nikcleju@62:     #    disp('Starting point infeasible; using x0 = At*inv(AAt)*y.');
nikcleju@62:     #    AAt = @(z) A(At(z));
nikcleju@62:     #    [w, cgres, cgiter] = cgsolve(AAt, b, cgtol, cgmaxiter, 0);
nikcleju@62:     #    if (cgres > 1/2)
nikcleju@62:     #      disp('A*At is ill-conditioned: cannot find starting point');
nikcleju@62:     #      xp = x0;
nikcleju@62:     #      return;
nikcleju@62:     #    end
nikcleju@62:     #    x0 = At(w);
nikcleju@62:     #  end
nikcleju@62:     #else
nikcleju@62:     #  if (norm(A*x0-b)/norm(b) > cgtol)
nikcleju@62:     #    disp('Starting point infeasible; using x0 = At*inv(AAt)*y.');
nikcleju@62:     #    opts.POSDEF = true; opts.SYM = true;
nikcleju@62:     #    [w, hcond] = linsolve(A*A', b, opts);
nikcleju@62:     #    if (hcond < 1e-14)
nikcleju@62:     #      disp('A*At is ill-conditioned: cannot find starting point');
nikcleju@62:     #      xp = x0;
nikcleju@62:     #      return;
nikcleju@62:     #    end
nikcleju@62:     #    x0 = A'*w;
nikcleju@62:     #  end  
nikcleju@62:     #end
nikcleju@62:     #x = x0;
nikcleju@62:     #u = (0.95)*abs(x0) + (0.10)*max(abs(x0));
nikcleju@62:     #
nikcleju@62:     ## set up for the first iteration
nikcleju@62:     #fu1 = x - u;
nikcleju@62:     #fu2 = -x - u;
nikcleju@62:     #lamu1 = -1./fu1;
nikcleju@62:     #lamu2 = -1./fu2;
nikcleju@62:     #if (largescale)
nikcleju@62:     #  v = -A(lamu1-lamu2);
nikcleju@62:     #  Atv = At(v);
nikcleju@62:     #  rpri = A(x) - b;
nikcleju@62:     #else
nikcleju@62:     #  v = -A*(lamu1-lamu2);
nikcleju@62:     #  Atv = A'*v;
nikcleju@62:     #  rpri = A*x - b;
nikcleju@62:     #end
nikcleju@62:     #
nikcleju@62:     #sdg = -(fu1'*lamu1 + fu2'*lamu2);
nikcleju@62:     #tau = mu*2*N/sdg;
nikcleju@62:     #
nikcleju@62:     #rcent = [-lamu1.*fu1; -lamu2.*fu2] - (1/tau);
nikcleju@62:     #rdual = gradf0 + [lamu1-lamu2; -lamu1-lamu2] + [Atv; zeros(N,1)];
nikcleju@62:     #resnorm = norm([rdual; rcent; rpri]);
nikcleju@62:     #
nikcleju@62:     #pditer = 0;
nikcleju@62:     #done = (sdg < pdtol) | (pditer >= pdmaxiter);
nikcleju@62:     #while (~done)
nikcleju@62:     #  
nikcleju@62:     #  pditer = pditer + 1;
nikcleju@62:     #  
nikcleju@62:     #  w1 = -1/tau*(-1./fu1 + 1./fu2) - Atv;
nikcleju@62:     #  w2 = -1 - 1/tau*(1./fu1 + 1./fu2);
nikcleju@62:     #  w3 = -rpri;
nikcleju@62:     #  
nikcleju@62:     #  sig1 = -lamu1./fu1 - lamu2./fu2;
nikcleju@62:     #  sig2 = lamu1./fu1 - lamu2./fu2;
nikcleju@62:     #  sigx = sig1 - sig2.^2./sig1;
nikcleju@62:     #  
nikcleju@62:     #  if (largescale)
nikcleju@62:     #    w1p = w3 - A(w1./sigx - w2.*sig2./(sigx.*sig1));
nikcleju@62:     #    h11pfun = @(z) -A(1./sigx.*At(z));
nikcleju@62:     #    [dv, cgres, cgiter] = cgsolve(h11pfun, w1p, cgtol, cgmaxiter, 0);
nikcleju@62:     #    if (cgres > 1/2)
nikcleju@62:     #      disp('Cannot solve system.  Returning previous iterate.  (See Section 4 of notes for more information.)');
nikcleju@62:     #      xp = x;
nikcleju@62:     #      return
nikcleju@62:     #    end
nikcleju@62:     #    dx = (w1 - w2.*sig2./sig1 - At(dv))./sigx;
nikcleju@62:     #    Adx = A(dx);
nikcleju@62:     #    Atdv = At(dv);
nikcleju@62:     #  else
nikcleju@62:     #    w1p = -(w3 - A*(w1./sigx - w2.*sig2./(sigx.*sig1)));
nikcleju@62:     #    H11p = A*(sparse(diag(1./sigx))*A');
nikcleju@62:     #    opts.POSDEF = true; opts.SYM = true;
nikcleju@62:     #    [dv,hcond] = linsolve(H11p, w1p, opts);
nikcleju@62:     #    if (hcond < 1e-14)
nikcleju@62:     #      disp('Matrix ill-conditioned.  Returning previous iterate.  (See Section 4 of notes for more information.)');
nikcleju@62:     #      xp = x;
nikcleju@62:     #      return
nikcleju@62:     #    end
nikcleju@62:     #    dx = (w1 - w2.*sig2./sig1 - A'*dv)./sigx;
nikcleju@62:     #    Adx = A*dx;
nikcleju@62:     #    Atdv = A'*dv;
nikcleju@62:     #  end
nikcleju@62:     #  
nikcleju@62:     #  du = (w2 - sig2.*dx)./sig1;
nikcleju@62:     #  
nikcleju@62:     #  dlamu1 = (lamu1./fu1).*(-dx+du) - lamu1 - (1/tau)*1./fu1;
nikcleju@62:     #  dlamu2 = (lamu2./fu2).*(dx+du) - lamu2 - 1/tau*1./fu2;
nikcleju@62:     #  
nikcleju@62:     #  # make sure that the step is feasible: keeps lamu1,lamu2 > 0, fu1,fu2 < 0
nikcleju@62:     #  indp = find(dlamu1 < 0);  indn = find(dlamu2 < 0);
nikcleju@62:     #  s = min([1; -lamu1(indp)./dlamu1(indp); -lamu2(indn)./dlamu2(indn)]);
nikcleju@62:     #  indp = find((dx-du) > 0);  indn = find((-dx-du) > 0);
nikcleju@62:     #  s = (0.99)*min([s; -fu1(indp)./(dx(indp)-du(indp)); -fu2(indn)./(-dx(indn)-du(indn))]);
nikcleju@62:     #  
nikcleju@62:     #  # backtracking line search
nikcleju@62:     #  suffdec = 0;
nikcleju@62:     #  backiter = 0;
nikcleju@62:     #  while (~suffdec)
nikcleju@62:     #    xp = x + s*dx;  up = u + s*du; 
nikcleju@62:     #    vp = v + s*dv;  Atvp = Atv + s*Atdv; 
nikcleju@62:     #    lamu1p = lamu1 + s*dlamu1;  lamu2p = lamu2 + s*dlamu2;
nikcleju@62:     #    fu1p = xp - up;  fu2p = -xp - up;  
nikcleju@62:     #    rdp = gradf0 + [lamu1p-lamu2p; -lamu1p-lamu2p] + [Atvp; zeros(N,1)];
nikcleju@62:     #    rcp = [-lamu1p.*fu1p; -lamu2p.*fu2p] - (1/tau);
nikcleju@62:     #    rpp = rpri + s*Adx;
nikcleju@62:     #    suffdec = (norm([rdp; rcp; rpp]) <= (1-alpha*s)*resnorm);
nikcleju@62:     #    s = beta*s;
nikcleju@62:     #    backiter = backiter + 1;
nikcleju@62:     #    if (backiter > 32)
nikcleju@62:     #      disp('Stuck backtracking, returning last iterate.  (See Section 4 of notes for more information.)')
nikcleju@62:     #      xp = x;
nikcleju@62:     #      return
nikcleju@62:     #    end
nikcleju@62:     #  end
nikcleju@62:     #  
nikcleju@62:     #  
nikcleju@62:     #  # next iteration
nikcleju@62:     #  x = xp;  u = up;
nikcleju@62:     #  v = vp;  Atv = Atvp; 
nikcleju@62:     #  lamu1 = lamu1p;  lamu2 = lamu2p;
nikcleju@62:     #  fu1 = fu1p;  fu2 = fu2p;
nikcleju@62:     #  
nikcleju@62:     #  # surrogate duality gap
nikcleju@62:     #  sdg = -(fu1'*lamu1 + fu2'*lamu2);
nikcleju@62:     #  tau = mu*2*N/sdg;
nikcleju@62:     #  rpri = rpp;
nikcleju@62:     #  rcent = [-lamu1.*fu1; -lamu2.*fu2] - (1/tau);
nikcleju@62:     #  rdual = gradf0 + [lamu1-lamu2; -lamu1-lamu2] + [Atv; zeros(N,1)];
nikcleju@62:     #  resnorm = norm([rdual; rcent; rpri]);
nikcleju@62:     #  
nikcleju@62:     #  done = (sdg < pdtol) | (pditer >= pdmaxiter);
nikcleju@62:     #  
nikcleju@62:     #  disp(sprintf('Iteration = #d, tau = #8.3e, Primal = #8.3e, PDGap = #8.3e, Dual res = #8.3e, Primal res = #8.3e',...
nikcleju@62:     #    pditer, tau, sum(u), sdg, norm(rdual), norm(rpri)));
nikcleju@62:     #  if (largescale)
nikcleju@62:     #    disp(sprintf('                  CG Res = #8.3e, CG Iter = #d', cgres, cgiter));
nikcleju@62:     #  else
nikcleju@62:     #    disp(sprintf('                  H11p condition number = #8.3e', hcond));
nikcleju@62:     #  end
nikcleju@62:     #  
nikcleju@62:     #end
nikcleju@62:     
nikcleju@62:     # End of original Matab code
nikcleju@62:     #----------------------------    
nikcleju@62:     
nikcleju@65:     # Nic: check if b is 0; if so, return 0
nikcleju@65:     #    Otherwise it will break later
nikcleju@65:     if numpy.linalg.norm(b) < 1e-16:
nikcleju@65:         return numpy.zeros_like(x0)
nikcleju@65:     
nikcleju@62:     #largescale = isa(A,'function_handle');
nikcleju@62:     if hasattr(A, '__call__'):
nikcleju@62:         largescale = True
nikcleju@62:     else:
nikcleju@62:         largescale = False       
nikcleju@62:     
nikcleju@62:     #N = length(x0);
nikcleju@62:     N = x0.size
nikcleju@62:     
nikcleju@62:     alpha = 0.01
nikcleju@62:     beta = 0.5
nikcleju@62:     mu = 10
nikcleju@62:     
nikcleju@62:     #gradf0 = [zeros(N,1); ones(N,1)];
nikcleju@62:     gradf0 = numpy.hstack((numpy.zeros(N), numpy.ones(N)))
nikcleju@62:     
nikcleju@62:     # starting point --- make sure that it is feasible
nikcleju@62:     #if (largescale)
nikcleju@62:     if largescale:
nikcleju@63:       raise l1eqNotImplementedError('Largescale not implemented yet!')
nikcleju@62:     else:
nikcleju@62:       #if (norm(A*x0-b)/norm(b) > cgtol)
nikcleju@62:       if numpy.linalg.norm(numpy.dot(A,x0)-b) / numpy.linalg.norm(b) > cgtol:
nikcleju@62:         #disp('Starting point infeasible; using x0 = At*inv(AAt)*y.');
nikcleju@62:         if verbose:
nikcleju@62:           print 'Starting point infeasible; using x0 = At*inv(AAt)*y.'
nikcleju@62:         #opts.POSDEF = true; opts.SYM = true;
nikcleju@62:         #[w, hcond] = linsolve(A*A', b, opts);
nikcleju@62:         #if (hcond < 1e-14)
nikcleju@62:         #  disp('A*At is ill-conditioned: cannot find starting point');
nikcleju@62:         #  xp = x0;
nikcleju@62:         #  return;
nikcleju@62:         #end
nikcleju@62:         #x0 = A'*w;
nikcleju@62:         try:
nikcleju@62:             w = scipy.linalg.solve(numpy.dot(A,A.T), b, sym_pos=True)
nikcleju@65:             hcond = 1.0/numpy.linalg.cond(numpy.dot(A,A.T))
nikcleju@62:         except scipy.linalg.LinAlgError:
nikcleju@62:             if verbose:
nikcleju@62:               print 'A*At is ill-conditioned: cannot find starting point'
nikcleju@62:             xp = x0.copy()
nikcleju@62:             return xp
nikcleju@62:         if hcond < 1e-14:
nikcleju@62:             if verbose:
nikcleju@62:               print 'A*At is ill-conditioned: cannot find starting point'
nikcleju@62:             xp = x0.copy()
nikcleju@62:             return xp           
nikcleju@62:         x0 = numpy.dot(A.T, w)
nikcleju@62:       #end  
nikcleju@62:     #end
nikcleju@62:     x = x0.copy()
nikcleju@62:     #u = (0.95)*abs(x0) + (0.10)*max(abs(x0));
nikcleju@62:     u = (0.95)*numpy.abs(x0) + (0.10)*numpy.abs(x0).max()
nikcleju@62:     
nikcleju@62:     # set up for the first iteration
nikcleju@62:     fu1 = x - u
nikcleju@62:     fu2 = -x - u
nikcleju@62:     lamu1 = -1/fu1
nikcleju@62:     lamu2 = -1/fu2
nikcleju@62:     if (largescale):
nikcleju@62:       #v = -A(lamu1-lamu2);
nikcleju@62:       #Atv = At(v);
nikcleju@62:       #rpri = A(x) - b;
nikcleju@63:       raise l1eqNotImplementedError('Largescale not implemented yet!')
nikcleju@62:     else:
nikcleju@62:       #v = -A*(lamu1-lamu2);
nikcleju@62:       #Atv = A'*v;
nikcleju@62:       #rpri = A*x - b;
nikcleju@62:       v = numpy.dot(-A, lamu1-lamu2)
nikcleju@62:       Atv = numpy.dot(A.T, v)
nikcleju@62:       rpri = numpy.dot(A,x) - b
nikcleju@62:     #end
nikcleju@62:     
nikcleju@62:     #sdg = -(fu1'*lamu1 + fu2'*lamu2);
nikcleju@62:     sdg = -(numpy.dot(fu1,lamu1) + numpy.dot(fu2,lamu2))
nikcleju@62:     tau = mu*2*N/sdg
nikcleju@62:     
nikcleju@62:     #rcent = [-lamu1.*fu1; -lamu2.*fu2] - (1/tau);
nikcleju@62:     rcent = numpy.hstack((-numpy.dot(lamu1,fu1),  -numpy.dot(lamu2,fu2))) - (1/tau)
nikcleju@62:     #rdual = gradf0 + [lamu1-lamu2; -lamu1-lamu2] + [Atv; zeros(N,1)];
nikcleju@62:     rdual = gradf0 + numpy.hstack((lamu1-lamu2, -lamu1-lamu2)) + numpy.hstack((Atv, numpy.zeros(N)))
nikcleju@62:     #resnorm = norm([rdual; rcent; rpri]);
nikcleju@62:     resnorm = numpy.linalg.norm(numpy.hstack((rdual, rcent, rpri)))
nikcleju@62:     
nikcleju@62:     pditer = 0
nikcleju@62:     #done = (sdg < pdtol) | (pditer >= pdmaxiter);
nikcleju@62:     done = (sdg < pdtol) or (pditer >= pdmaxiter)
nikcleju@62:     #while (~done)
nikcleju@62:     while not done:
nikcleju@62:       
nikcleju@62:       pditer = pditer + 1
nikcleju@62:       
nikcleju@62:       #w1 = -1/tau*(-1./fu1 + 1./fu2) - Atv;
nikcleju@62:       w1 = -1/tau*(-1/fu1 + 1/fu2) - Atv
nikcleju@62:       w2 = -1 - 1/tau*(1/fu1 + 1/fu2)
nikcleju@62:       w3 = -rpri
nikcleju@62:       
nikcleju@62:       #sig1 = -lamu1./fu1 - lamu2./fu2;
nikcleju@62:       sig1 = -lamu1/fu1 - lamu2/fu2
nikcleju@62:       sig2 = lamu1/fu1 - lamu2/fu2
nikcleju@62:       #sigx = sig1 - sig2.^2./sig1;
nikcleju@62:       sigx = sig1 - sig2**2/sig1
nikcleju@62:       
nikcleju@62:       if largescale:
nikcleju@62:         #w1p = w3 - A(w1./sigx - w2.*sig2./(sigx.*sig1));
nikcleju@62:         #h11pfun = @(z) -A(1./sigx.*At(z));
nikcleju@62:         #[dv, cgres, cgiter] = cgsolve(h11pfun, w1p, cgtol, cgmaxiter, 0);
nikcleju@62:         #if (cgres > 1/2)
nikcleju@62:         #  disp('Cannot solve system.  Returning previous iterate.  (See Section 4 of notes for more information.)');
nikcleju@62:         #  xp = x;
nikcleju@62:         #  return
nikcleju@62:         #end
nikcleju@62:         #dx = (w1 - w2.*sig2./sig1 - At(dv))./sigx;
nikcleju@62:         #Adx = A(dx);
nikcleju@62:         #Atdv = At(dv);
nikcleju@63:         raise l1eqNotImplementedError('Largescale not implemented yet!')
nikcleju@62:       else:
nikcleju@62:         #w1p = -(w3 - A*(w1./sigx - w2.*sig2./(sigx.*sig1)));
nikcleju@62:         w1p = -(w3 - numpy.dot(A,(w1/sigx - w2*sig2/(sigx*sig1))))
nikcleju@62:         #H11p = A*(sparse(diag(1./sigx))*A');
nikcleju@62:         H11p = numpy.dot(A, numpy.dot(numpy.diag(1/sigx),A.T))
nikcleju@62:         #opts.POSDEF = true; opts.SYM = true;
nikcleju@62:         #[dv,hcond] = linsolve(H11p, w1p, opts);
nikcleju@62:         try:
nikcleju@62:           dv = scipy.linalg.solve(H11p, w1p, sym_pos=True)
nikcleju@62:           hcond = 1.0/numpy.linalg.cond(H11p)
nikcleju@62:         except scipy.linalg.LinAlgError:
nikcleju@62:             if verbose:
nikcleju@62:               print 'Matrix ill-conditioned.  Returning previous iterate.  (See Section 4 of notes for more information.)'
nikcleju@62:             xp = x.copy()
nikcleju@62:             return xp
nikcleju@62:         if hcond < 1e-14:
nikcleju@62:             if verbose:
nikcleju@62:               print 'Matrix ill-conditioned.  Returning previous iterate.  (See Section 4 of notes for more information.)'
nikcleju@62:             xp = x.copy()
nikcleju@62:             return xp            
nikcleju@62:         #if (hcond < 1e-14)
nikcleju@62:         #  disp('Matrix ill-conditioned.  Returning previous iterate.  (See Section 4 of notes for more information.)');
nikcleju@62:         #  xp = x;
nikcleju@62:         #  return
nikcleju@62:         #end
nikcleju@62:         
nikcleju@62:         #dx = (w1 - w2.*sig2./sig1 - A'*dv)./sigx;
nikcleju@62:         dx = (w1 - w2*sig2/sig1 - numpy.dot(A.T,dv))/sigx
nikcleju@62:         #Adx = A*dx;
nikcleju@62:         Adx = numpy.dot(A,dx)
nikcleju@62:         #Atdv = A'*dv;
nikcleju@62:         Atdv = numpy.dot(A.T,dv)
nikcleju@62:       #end
nikcleju@62:       
nikcleju@62:       #du = (w2 - sig2.*dx)./sig1;
nikcleju@62:       du = (w2 - sig2*dx)/sig1
nikcleju@62:       
nikcleju@62:       #dlamu1 = (lamu1./fu1).*(-dx+du) - lamu1 - (1/tau)*1./fu1;
nikcleju@62:       dlamu1 = (lamu1/fu1)*(-dx+du) - lamu1 - (1/tau)*1/fu1
nikcleju@62:       dlamu2 = (lamu2/fu2)*(dx+du) - lamu2 - 1/tau*1/fu2
nikcleju@62:       
nikcleju@62:       # make sure that the step is feasible: keeps lamu1,lamu2 > 0, fu1,fu2 < 0
nikcleju@62:       #indp = find(dlamu1 < 0);  indn = find(dlamu2 < 0);
nikcleju@62:       indp = numpy.nonzero(dlamu1 < 0)
nikcleju@62:       indn = numpy.nonzero(dlamu2 < 0)
nikcleju@62:       #s = min([1; -lamu1(indp)./dlamu1(indp); -lamu2(indn)./dlamu2(indn)]);
nikcleju@62:       s = numpy.min(numpy.hstack((numpy.array([1]), -lamu1[indp]/dlamu1[indp], -lamu2[indn]/dlamu2[indn])))
nikcleju@62:       #indp = find((dx-du) > 0);  indn = find((-dx-du) > 0);
nikcleju@62:       indp = numpy.nonzero((dx-du) > 0)
nikcleju@62:       indn = numpy.nonzero((-dx-du) > 0)
nikcleju@62:       #s = (0.99)*min([s; -fu1(indp)./(dx(indp)-du(indp)); -fu2(indn)./(-dx(indn)-du(indn))]);
nikcleju@62:       s = (0.99)*numpy.min(numpy.hstack((numpy.array([s]), -fu1[indp]/(dx[indp]-du[indp]), -fu2[indn]/(-dx[indn]-du[indn]))))
nikcleju@62:       
nikcleju@62:       # backtracking line search
nikcleju@62:       suffdec = 0
nikcleju@62:       backiter = 0
nikcleju@62:       #while (~suffdec)
nikcleju@62:       while not suffdec:
nikcleju@62:         #xp = x + s*dx;  up = u + s*du; 
nikcleju@62:         xp = x + s*dx
nikcleju@62:         up = u + s*du 
nikcleju@62:         #vp = v + s*dv;  Atvp = Atv + s*Atdv; 
nikcleju@62:         vp = v + s*dv
nikcleju@62:         Atvp = Atv + s*Atdv
nikcleju@62:         #lamu1p = lamu1 + s*dlamu1;  lamu2p = lamu2 + s*dlamu2;
nikcleju@62:         lamu1p = lamu1 + s*dlamu1
nikcleju@62:         lamu2p = lamu2 + s*dlamu2
nikcleju@62:         #fu1p = xp - up;  fu2p = -xp - up;  
nikcleju@62:         fu1p = xp - up
nikcleju@62:         fu2p = -xp - up
nikcleju@62:         #rdp = gradf0 + [lamu1p-lamu2p; -lamu1p-lamu2p] + [Atvp; zeros(N,1)];
nikcleju@62:         rdp = gradf0 + numpy.hstack((lamu1p-lamu2p, -lamu1p-lamu2p)) + numpy.hstack((Atvp, numpy.zeros(N)))
nikcleju@62:         #rcp = [-lamu1p.*fu1p; -lamu2p.*fu2p] - (1/tau);
nikcleju@62:         rcp = numpy.hstack((-lamu1p*fu1p, -lamu2p*fu2p)) - (1/tau)
nikcleju@62:         #rpp = rpri + s*Adx;
nikcleju@62:         rpp = rpri + s*Adx
nikcleju@62:         #suffdec = (norm([rdp; rcp; rpp]) <= (1-alpha*s)*resnorm);
nikcleju@62:         suffdec = (numpy.linalg.norm(numpy.hstack((rdp, rcp, rpp))) <= (1-alpha*s)*resnorm)
nikcleju@62:         s = beta*s
nikcleju@62:         backiter = backiter + 1
nikcleju@62:         if (backiter > 32):
nikcleju@62:           if verbose:
nikcleju@62:             print 'Stuck backtracking, returning last iterate.  (See Section 4 of notes for more information.)'
nikcleju@62:           xp = x.copy()
nikcleju@62:           return xp
nikcleju@62:         #end
nikcleju@62:       #end
nikcleju@62:       
nikcleju@62:       
nikcleju@62:       # next iteration
nikcleju@62:       #x = xp;  u = up;
nikcleju@62:       x = xp.copy()
nikcleju@62:       u = up.copy()
nikcleju@62:       #v = vp;  Atv = Atvp; 
nikcleju@62:       v = vp.copy()
nikcleju@62:       Atv = Atvp.copy()
nikcleju@62:       #lamu1 = lamu1p;  lamu2 = lamu2p;
nikcleju@62:       lamu1 = lamu1p.copy()
nikcleju@62:       lamu2 = lamu2p.copy()
nikcleju@62:       #fu1 = fu1p;  fu2 = fu2p;
nikcleju@62:       fu1 = fu1p.copy()
nikcleju@62:       fu2 = fu2p.copy()
nikcleju@62:       
nikcleju@62:       # surrogate duality gap
nikcleju@62:       #sdg = -(fu1'*lamu1 + fu2'*lamu2);
nikcleju@62:       sdg = -(numpy.dot(fu1,lamu1) + numpy.dot(fu2,lamu2))
nikcleju@62:       tau = mu*2*N/sdg
nikcleju@62:       rpri = rpp.copy()
nikcleju@62:       #rcent = [-lamu1.*fu1; -lamu2.*fu2] - (1/tau);
nikcleju@62:       rcent = numpy.hstack((-lamu1*fu1, -lamu2*fu2)) - (1/tau)
nikcleju@62:       #rdual = gradf0 + [lamu1-lamu2; -lamu1-lamu2] + [Atv; zeros(N,1)];
nikcleju@62:       rdual = gradf0 + numpy.hstack((lamu1-lamu2, -lamu1-lamu2)) + numpy.hstack((Atv, numpy.zeros(N)))
nikcleju@62:       #resnorm = norm([rdual; rcent; rpri]);
nikcleju@62:       resnorm = numpy.linalg.norm(numpy.hstack((rdual, rcent, rpri)))
nikcleju@62:       
nikcleju@62:       #done = (sdg < pdtol) | (pditer >= pdmaxiter);
nikcleju@62:       done = (sdg < pdtol) or (pditer >= pdmaxiter)
nikcleju@62:       
nikcleju@62:       if verbose:
nikcleju@62:         print 'Iteration =',pditer,', tau =',tau,', Primal =',numpy.sum(u),', PDGap =',sdg,', Dual res =',numpy.linalg.norm(rdual),', Primal res =',numpy.linalg.norm(rpri)
nikcleju@62:       if largescale:
nikcleju@62:         #disp(sprintf('                  CG Res = #8.3e, CG Iter = #d', cgres, cgiter));
nikcleju@63:         raise l1eqNotImplementedError('Largescale not implemented yet!')
nikcleju@62:       else:
nikcleju@62:         #disp(sprintf('                  H11p condition number = #8.3e', hcond));
nikcleju@62:         if verbose:
nikcleju@62:           print '                  H11p condition number =',hcond
nikcleju@62:       #end
nikcleju@62:       
nikcleju@62:     #end
nikcleju@62:     
nikcleju@62:     return xp