nikcleju@66: # -*- coding: utf-8 -*-
nikcleju@66: """
nikcleju@66: 
nikcleju@66: Author: Nicolae Cleju
nikcleju@66: """
nikcleju@66: __author__ = "Nicolae Cleju"
nikcleju@66: __license__ = "GPL"
nikcleju@66: __email__ = "nikcleju@gmail.com"
nikcleju@66: 
nikcleju@66: 
nikcleju@66: import numpy
nikcleju@66: import scipy
nikcleju@66: import math
nikcleju@67: #import cvxpy
nikcleju@66: 
nikcleju@66: class EllipseProjDaiError(Exception):
nikcleju@66: #  def __init__(self, e):
nikcleju@66: #    self._e = e
nikcleju@66:   pass
nikcleju@66: 
nikcleju@66: def ellipse_proj_cvxpy(A,b,p,epsilon):
nikcleju@66:   b = b.reshape((b.size,1))
nikcleju@66:   p = p.reshape((p.size,1))
nikcleju@66:   A = cvxpy.matrix(A)
nikcleju@66:   b = cvxpy.matrix(b)
nikcleju@66:   p = cvxpy.matrix(p)
nikcleju@66: 
nikcleju@66:   x = cvxpy.variable(p.shape[0],1)                 # Create optimization variable
nikcleju@66:   prog = cvxpy.program(cvxpy.minimize(cvxpy.norm2(p - x)),              # Create problem instance
nikcleju@66:             [cvxpy.leq(cvxpy.norm2(b-A*x),epsilon)])
nikcleju@66:   #prog.solve(quiet=True)                              # Get optimal value
nikcleju@66:   prog.solve()                              # Get optimal value
nikcleju@66:   return numpy.array(x.value).flatten()
nikcleju@66: 
nikcleju@66: def ellipse_proj_dai(A,x,s,eps):
nikcleju@66:   
nikcleju@66:   A_pinv = numpy.linalg.pinv(A)
nikcleju@66:   
nikcleju@66:   AA = numpy.dot(A.T, A)
nikcleju@66:   bb = -numpy.dot(x.T,A)
nikcleju@66:   alpha = eps*eps - numpy.inner(x,x)
nikcleju@66:   
nikcleju@66:   direction = numpy.dot(A_pinv,(numpy.dot(A,s)-x))
nikcleju@66:   s0 = s - (1.0 - eps/numpy.linalg.norm(numpy.dot(A,direction))) * direction
nikcleju@66:   
nikcleju@66:   a = s
nikcleju@66:   xk = s0
nikcleju@66:   m1 = 1
nikcleju@66:   m2 = 1
nikcleju@66:   c1 = 0.1
nikcleju@66:   c2 = 0.8
nikcleju@66:   done = False
nikcleju@66:   k = 0
nikcleju@66:   while not done and k < 50:
nikcleju@66:     
nikcleju@66:     uk = numpy.dot(AA,xk) + bb
nikcleju@66:     vk = a - xk
nikcleju@66:     
nikcleju@66:     zk = uk - (numpy.inner(uk,vk)/numpy.inner(vk,vk))*vk
nikcleju@66:     vtav = numpy.inner(vk,numpy.dot(AA,vk))
nikcleju@66:     vtv = numpy.inner(vk,vk)
nikcleju@66:     ztaz = numpy.inner(zk,numpy.dot(AA,zk))
nikcleju@66:     ztz = numpy.inner(zk,zk)
nikcleju@66:     vtaz = numpy.inner(vk,numpy.dot(AA,zk))
nikcleju@66:     gammak1 = 1.0/(0.5 * ( vtav/vtv + ztaz/ztz + math.sqrt((vtaz/vtv - ztaz/ztz)**2 + 4*vtaz**2/vtv/ztz)))
nikcleju@66:     
nikcleju@66:     pk = vk / numpy.linalg.norm(vk)
nikcleju@66:     qk = zk / numpy.linalg.norm(zk)
nikcleju@66:     Hk = numpy.zeros((pk.size,2))
nikcleju@66:     Hk[:,0] = pk
nikcleju@66:     Hk[:,1] = qk
nikcleju@66:     Ak = numpy.dot(Hk.T, numpy.dot(AA,Hk))
nikcleju@66:     bk = numpy.dot(Hk.T,uk)
nikcleju@66:     
nikcleju@66:     #al = numpy.dot(Hk, numpy.dot(Hk.T, a))
nikcleju@66:     al = numpy.array([numpy.linalg.norm(vk), 0])
nikcleju@66:     D,Q = numpy.linalg.eig(Ak)
nikcleju@66:     Q = Q.T
nikcleju@66:     ah = numpy.dot(Q,al) + (1.0/D) * numpy.dot(Q,bk)
nikcleju@66:     
nikcleju@66:     l1 = D[0]
nikcleju@66:     l2 = D[1]
nikcleju@66:     Qbk = numpy.dot(Q,bk)
nikcleju@66:     beta = numpy.dot(Qbk.T, (1.0/D) * Qbk)
nikcleju@66:     hstar1s = numpy.roots(numpy.array([ (l1-l2)**2*l1, 
nikcleju@66:                                   2*(l1-l2)*l2*ah[0]*l1,
nikcleju@66:                                   -(l1-l2)**2*beta + l2**2*ah[0]**2*l1 + l1**2*l2*ah[1]**2,
nikcleju@66:                                   -2*beta*(l1-l2)*l2*ah[0],
nikcleju@66:                                   -beta*l2**2*ah[0]**2]))
nikcleju@66:     hstar2s = numpy.zeros_like(hstar1s)
nikcleju@66:     i_s = []
nikcleju@66:     illcond = False  # flag if ill conditioned problem (numerical errros)
nikcleju@66:     for i in range(hstar1s.size):
nikcleju@66:       
nikcleju@66:       # Protect against bad conditioning (ratio of two very small values)
nikcleju@66:       if numpy.abs(l1*ah[1]) > 1e-6 and numpy.abs((l1-l2) + l2*ah[0]/hstar1s[i]) > 1e-6:
nikcleju@66:         
nikcleju@66:         # Ignore small imaginary parts
nikcleju@66:         if numpy.abs(numpy.imag(hstar1s[i])) < 1e-10:
nikcleju@66:           hstar1s[i] = numpy.real(hstar1s[i])
nikcleju@66:         
nikcleju@66:         hstar2s[i] = l1*ah[1] / ((l1-l2)*hstar1s[i] + l2*ah[0]) * hstar1s[i]
nikcleju@66:         
nikcleju@66:         # Ignore small imaginary parts
nikcleju@66:         if numpy.abs(numpy.imag(hstar2s[i])) < 1e-10:
nikcleju@66:           hstar2s[i] = numpy.real(hstar2s[i])
nikcleju@66:           
nikcleju@66:         if (ah[0] - hstar1s[i]) / (l1*hstar1s[i]) > 0 and (ah[1] - hstar2s[i]) / (l2*hstar2s[i]) > 0:
nikcleju@66:           i_s.append(i)
nikcleju@66:         
nikcleju@66:       else:
nikcleju@66:         # Cannot rely on hstar2s[i] calculation since it is the ratio of two small numbers
nikcleju@66:         # Do a vertical projection instead
nikcleju@66:           hstar1 = ah[0]
nikcleju@66:           hstar2 = numpy.sign(ah[1]) * math.sqrt((beta - l1*ah[0]**2)/l2)
nikcleju@66:           illcond = True  # Flag, so we don't take i_s[] into account anymore
nikcleju@66:     
nikcleju@66:     if illcond:
nikcleju@66:       print "Ill conditioned problem, do vertical projection instead"
nikcleju@66:       # hstar1 and hstar2 are already set above, nothing to do here
nikcleju@66:     else:
nikcleju@66:       if len(i_s) > 1:
nikcleju@66:         hstar1 = hstar1s[i_s[0]].real
nikcleju@66:         hstar2 = hstar2s[i_s[0]].real
nikcleju@66:       elif len(i_s) == 0:
nikcleju@66:         # Again do vertical projection
nikcleju@66:         hstar1 = ah[0]
nikcleju@66:         hstar2 = numpy.sign(ah[1]) * math.sqrt((beta - l1*ah[0]**2)/l2)
nikcleju@66:       else:
nikcleju@66:         # Everything is ok
nikcleju@66:         hstar1 = hstar1s[i_s].real
nikcleju@66:         hstar2 = hstar2s[i_s].real
nikcleju@66:       
nikcleju@66:     ahstar = numpy.array([hstar1, hstar2]).flatten()
nikcleju@66:     alstar = numpy.dot(Q.T, ahstar - (1.0/D) * numpy.dot(Q,bk))
nikcleju@66:     alstar1 = alstar[0]
nikcleju@66:     alstar2 = alstar[1]
nikcleju@66:     etak = 1 - alstar1/numpy.linalg.norm(vk) + numpy.inner(uk,vk)/vtv*alstar2/numpy.linalg.norm(zk)
nikcleju@66:     gammak2 = -alstar2/numpy.linalg.norm(zk)/etak
nikcleju@66:     if k % (m1+m2) < m1:
nikcleju@66:       gammak = gammak2
nikcleju@66:     else:
nikcleju@66:       gammak = c1*gammak1 + c2*gammak2
nikcleju@66:       
nikcleju@66:     # Safeguard:
nikcleju@66:     if gammak < 0:
nikcleju@66:       gammak = gammak1
nikcleju@66: 
nikcleju@66:     ck = xk - gammak * uk
nikcleju@66:     wk = ck - a
nikcleju@66:     
nikcleju@66:     ga = numpy.dot(AA,a) + bb
nikcleju@66:     qa = numpy.inner(a,numpy.dot(AA,a)) + 2*numpy.inner(bb,a)
nikcleju@66:     # Check if delta < 0 but very small (possibly numerical errors)
nikcleju@66:     if (numpy.inner(ga,wk)/numpy.inner(wk, numpy.dot(AA,wk)))**2 - (qa-alpha)/numpy.inner(wk, numpy.dot(AA,wk)) < 0 and abs((numpy.inner(ga,wk)/numpy.inner(wk, numpy.dot(AA,wk)))**2 - (qa-alpha)/numpy.inner(wk, numpy.dot(AA,wk))) < 1e-10:
nikcleju@66:       etak = -numpy.inner(ga,wk)/numpy.inner(wk, numpy.dot(AA,wk))
nikcleju@66:     else:
nikcleju@66:       assert ((numpy.inner(ga,wk)/numpy.inner(wk, numpy.dot(AA,wk)))**2 - (qa-alpha)/numpy.inner(wk, numpy.dot(AA,wk)) >= 0)
nikcleju@66:       etak = -numpy.inner(ga,wk)/numpy.inner(wk, numpy.dot(AA,wk)) - math.sqrt((numpy.inner(ga,wk)/numpy.inner(wk, numpy.dot(AA,wk)))**2 - (qa-alpha)/numpy.inner(wk, numpy.dot(AA,wk)))
nikcleju@66:     xk = a + etak * wk
nikcleju@66:     
nikcleju@66:     if (1 - numpy.inner(uk,vk) / numpy.linalg.norm(uk) / numpy.linalg.norm(vk)) <= 0.01:
nikcleju@66:       done = True
nikcleju@66:     k = k+1
nikcleju@66:   
nikcleju@66:   return xk
nikcleju@66: 
nikcleju@66:   
nikcleju@66: def ellipse_proj_proj(A,x,s,eps,L2):
nikcleju@66:   A_pinv = numpy.linalg.pinv(A)
nikcleju@66:   u,singvals,v  = numpy.linalg.svd(A, full_matrices=0)
nikcleju@66:   singvals = numpy.flipud(singvals)
nikcleju@66:   s_orig = s
nikcleju@66:   
nikcleju@66:   for j in range(L2):
nikcleju@66:     direction = numpy.dot(A_pinv,(numpy.dot(A,s)-x))
nikcleju@66: 
nikcleju@66:     if (numpy.linalg.norm(numpy.dot(A,direction)) >= eps):
nikcleju@66:       
nikcleju@66:       P = numpy.dot(numpy.dot(u,v), s)
nikcleju@66:       sproj = s - (1.0 - eps/numpy.linalg.norm(numpy.dot(A,direction))) * direction
nikcleju@66:       P0 = numpy.dot(numpy.dot(u,v), sproj)
nikcleju@66:       
nikcleju@66:       tangent = (P0 * (singvals**2)).reshape((1,P.size))
nikcleju@66:       uu,ss,vv = numpy.linalg.svd(tangent)
nikcleju@66:       svd = vv[1:,:]
nikcleju@66:       P1 = numpy.linalg.solve(numpy.vstack((tangent,svd)), numpy.vstack((numpy.array([[eps]]), numpy.dot(svd, P).reshape((svd.shape[0],1)))))
nikcleju@66:       
nikcleju@66:       # Take only a smaller step
nikcleju@66:       #P1 = P0.reshape((P0.size,1)) + 0.1*(P1-P0.reshape((P0.size,1)))
nikcleju@66:       
nikcleju@66:       #s = numpy.dot(A_pinv,P1).flatten() + numpy.dot(A_pinv,numpy.dot(A,s)).flatten()
nikcleju@66:       s = numpy.dot(numpy.linalg.pinv(numpy.dot(u,v)),P1).flatten() + (s-numpy.dot(A_pinv,numpy.dot(A,s)).flatten())
nikcleju@66:       
nikcleju@66:       #assert(numpy.linalg.norm(x - numpy.dot(A,s)) < eps + 1e-6)          
nikcleju@66:   
nikcleju@66:   direction = numpy.dot(A_pinv,(numpy.dot(A,s)-x))
nikcleju@66:   if (numpy.linalg.norm(numpy.dot(A,direction)) >= eps):
nikcleju@66:     s = s - (1.0 - eps/numpy.linalg.norm(numpy.dot(A,direction))) * direction  
nikcleju@66:     
nikcleju@66:   return s
nikcleju@66: 
nikcleju@66: def ellipse_proj_logbarrier(A,b,p,epsilon,verbose=False):
nikcleju@66:   return l1qc_logbarrier(numpy.zeros_like(p), p, A, A.T, b, epsilon, lbtol=1e-3, mu=10, cgtol=1e-8, cgmaxiter=200, verbose=verbose)
nikcleju@66: 
nikcleju@66: def l1qc_newton(x0, p, A, At, b, epsilon, tau, newtontol, newtonmaxiter, cgtol, cgmaxiter, verbose=False):
nikcleju@66:     # Newton algorithm for log-barrier subproblems for l1 minimization
nikcleju@66:     # with quadratic constraints.
nikcleju@66:     #
nikcleju@66:     # Usage: 
nikcleju@66:     # [xp,up,niter] = l1qc_newton(x0, u0, A, At, b, epsilon, tau, 
nikcleju@66:     #                             newtontol, newtonmaxiter, cgtol, cgmaxiter)
nikcleju@66:     #
nikcleju@66:     # x0,u0 - starting points
nikcleju@66:     #
nikcleju@66:     # A - Either a handle to a function that takes a N vector and returns a K 
nikcleju@66:     #     vector , or a KxN matrix.  If A is a function handle, the algorithm
nikcleju@66:     #     operates in "largescale" mode, solving the Newton systems via the
nikcleju@66:     #     Conjugate Gradients algorithm.
nikcleju@66:     #
nikcleju@66:     # At - Handle to a function that takes a K vector and returns an N vector.
nikcleju@66:     #      If A is a KxN matrix, At is ignored.
nikcleju@66:     #
nikcleju@66:     # b - Kx1 vector of observations.
nikcleju@66:     #
nikcleju@66:     # epsilon - scalar, constraint relaxation parameter
nikcleju@66:     #
nikcleju@66:     # tau - Log barrier parameter.
nikcleju@66:     #
nikcleju@66:     # newtontol - Terminate when the Newton decrement is <= newtontol.
nikcleju@66:     #         Default = 1e-3.
nikcleju@66:     #
nikcleju@66:     # newtonmaxiter - Maximum number of iterations.
nikcleju@66:     #         Default = 50.
nikcleju@66:     #
nikcleju@66:     # cgtol - Tolerance for Conjugate Gradients; ignored if A is a matrix.
nikcleju@66:     #     Default = 1e-8.
nikcleju@66:     #
nikcleju@66:     # cgmaxiter - Maximum number of iterations for Conjugate Gradients; ignored
nikcleju@66:     #     if A is a matrix.
nikcleju@66:     #     Default = 200.
nikcleju@66:     #
nikcleju@66:     # Written by: Justin Romberg, Caltech
nikcleju@66:     # Email: jrom@acm.caltech.edu
nikcleju@66:     # Created: October 2005
nikcleju@66:     #
nikcleju@66:     
nikcleju@66:     # check if the matrix A is implicit or explicit
nikcleju@66:     #largescale = isa(A,'function_handle');
nikcleju@66:     if hasattr(A, '__call__'):
nikcleju@66:         largescale = True
nikcleju@66:     else:
nikcleju@66:         largescale = False    
nikcleju@66:     
nikcleju@66:     # line search parameters
nikcleju@66:     alpha = 0.01
nikcleju@66:     beta = 0.5
nikcleju@66:     
nikcleju@66:     #if (~largescale), AtA = A'*A; end
nikcleju@66:     if not largescale:
nikcleju@66:         AtA = numpy.dot(A.T,A)
nikcleju@66:     
nikcleju@66:     # initial point
nikcleju@66:     x = x0.copy()
nikcleju@66:     #u = u0.copy()
nikcleju@66:     #if (largescale), r = A(x) - b; else  r = A*x - b; end
nikcleju@66:     if largescale:
nikcleju@66:         r = A(x) - b
nikcleju@66:     else:
nikcleju@66:         r = numpy.dot(A,x) - b
nikcleju@66:         
nikcleju@66:     #fu1 = x - u
nikcleju@66:     #fu2 = -x - u
nikcleju@66:     fe = 1.0/2*(numpy.vdot(r,r) - epsilon**2)
nikcleju@66:     #f = u.sum() - (1.0/tau)*(numpy.log(-fu1).sum() + numpy.log(-fu2).sum() + math.log(-fe))
nikcleju@66:     f = numpy.linalg.norm(p-x)**2 - (1.0/tau) * math.log(-fe)
nikcleju@66:     
nikcleju@66:     niter = 0
nikcleju@66:     done = 0
nikcleju@66:     while not done:
nikcleju@66:       
nikcleju@66:       #if (largescale), atr = At(r); else  atr = A'*r; end
nikcleju@66:       if largescale:
nikcleju@66:           atr = At(r)
nikcleju@66:       else:
nikcleju@66:           atr = numpy.dot(A.T,r)
nikcleju@66:       
nikcleju@66:       ##ntgz = 1./fu1 - 1./fu2 + 1/fe*atr;
nikcleju@66:       #ntgz = 1.0/fu1 - 1.0/fu2 + 1.0/fe*atr
nikcleju@66:       ntgz = 1.0/fe*atr
nikcleju@66:       ##ntgu = -tau - 1./fu1 - 1./fu2;
nikcleju@66:       #ntgu = -tau - 1.0/fu1 - 1.0/fu2
nikcleju@66:       ##gradf = -(1/tau)*[ntgz; ntgu];
nikcleju@66:       #gradf = -(1.0/tau)*numpy.concatenate((ntgz, ntgu),0)
nikcleju@66:       gradf = 2*x + 2*p -(1.0/tau)*ntgz
nikcleju@66:       
nikcleju@66:       ##sig11 = 1./fu1.^2 + 1./fu2.^2;
nikcleju@66:       #sig11 = 1.0/(fu1**2) + 1.0/(fu2**2)
nikcleju@66:       ##sig12 = -1./fu1.^2 + 1./fu2.^2;
nikcleju@66:       #sig12 = -1.0/(fu1**2) + 1.0/(fu2**2)
nikcleju@66:       ##sigx = sig11 - sig12.^2./sig11;
nikcleju@66:       #sigx = sig11 - (sig12**2)/sig11
nikcleju@66:         
nikcleju@66:       ##w1p = ntgz - sig12./sig11.*ntgu;
nikcleju@66:       #w1p = ntgz - sig12/sig11*ntgu
nikcleju@66:       #w1p = -tau * (2*x + 2*p) + ntgz
nikcleju@66:       w1p = -gradf
nikcleju@66:       if largescale:
nikcleju@66:         #h11pfun = @(z) sigx.*z - (1/fe)*At(A(z)) + 1/fe^2*(atr'*z)*atr;
nikcleju@66:         h11pfun = lambda z: sigx*z - (1.0/fe)*At(A(z)) + 1.0/(fe**2)*numpy.dot(numpy.dot(atr.T,z),atr)
nikcleju@66:         dx,cgres,cgiter = cgsolve(h11pfun, w1p, cgtol, cgmaxiter, 0)
nikcleju@66:         if (cgres > 1.0/2):
nikcleju@66:           if verbose:
nikcleju@66:             print 'Cannot solve system.  Returning previous iterate.  (See Section 4 of notes for more information.)'
nikcleju@66:           xp = x.copy()
nikcleju@66:           up = u.copy()
nikcleju@66:           return xp,up,niter
nikcleju@66:         #end
nikcleju@66:         Adx = A(dx)
nikcleju@66:       else:
nikcleju@66:         ##H11p = diag(sigx) - (1/fe)*AtA + (1/fe)^2*atr*atr';
nikcleju@66:         # Attention: atr is column vector, so atr*atr' means outer(atr,atr)
nikcleju@66:         #H11p = numpy.diag(sigx) - (1.0/fe)*AtA + (1.0/fe)**2*numpy.outer(atr,atr)
nikcleju@66:         H11p = 2 * numpy.eye(x.size) - (1.0/fe)*AtA + (1.0/fe)**2*numpy.outer(atr,atr)
nikcleju@66:         #opts.POSDEF = true; opts.SYM = true;
nikcleju@66:         #[dx,hcond] = linsolve(H11p, w1p, opts);
nikcleju@66:         #if (hcond < 1e-14)
nikcleju@66:         #  disp('Matrix ill-conditioned.  Returning previous iterate.  (See Section 4 of notes for more information.)');
nikcleju@66:         #  xp = x;  up = u;
nikcleju@66:         #  return
nikcleju@66:         #end
nikcleju@66:         try:
nikcleju@66:             dx = scipy.linalg.solve(H11p, w1p, sym_pos=True)
nikcleju@66:             #dx = numpy.linalg.solve(H11p, w1p)
nikcleju@66:             hcond = 1.0/numpy.linalg.cond(H11p)
nikcleju@66:         except scipy.linalg.LinAlgError:
nikcleju@66:             if verbose:
nikcleju@66:               print 'Matrix ill-conditioned.  Returning previous iterate.  (See Section 4 of notes for more information.)'
nikcleju@66:             xp = x.copy()
nikcleju@66:             #up = u.copy()
nikcleju@66:             #return xp,up,niter
nikcleju@66:             return xp,niter
nikcleju@66:         if hcond < 1e-14:
nikcleju@66:             if verbose:
nikcleju@66:               print 'Matrix ill-conditioned.  Returning previous iterate.  (See Section 4 of notes for more information.)'
nikcleju@66:             xp = x.copy()
nikcleju@66:             #up = u.copy()
nikcleju@66:             #return xp,up,niter
nikcleju@66:             return xp,niter
nikcleju@66:         
nikcleju@66:         #Adx = A*dx;
nikcleju@66:         Adx = numpy.dot(A,dx)
nikcleju@66:       #end
nikcleju@66:       ##du = (1./sig11).*ntgu - (sig12./sig11).*dx;  
nikcleju@66:       #du = (1.0/sig11)*ntgu - (sig12/sig11)*dx;
nikcleju@66:      
nikcleju@66:       # minimum step size that stays in the interior
nikcleju@66:       ##ifu1 = find((dx-du) > 0); ifu2 = find((-dx-du) > 0);
nikcleju@66:       #ifu1 = numpy.nonzero((dx-du)>0)
nikcleju@66:       #ifu2 = numpy.nonzero((-dx-du)>0)
nikcleju@66:       #aqe = Adx'*Adx;   bqe = 2*r'*Adx;   cqe = r'*r - epsilon^2;
nikcleju@66:       aqe = numpy.dot(Adx.T,Adx)
nikcleju@66:       bqe = 2*numpy.dot(r.T,Adx)
nikcleju@66:       cqe = numpy.vdot(r,r) - epsilon**2
nikcleju@66:       #smax = min(1,min([...
nikcleju@66:       #  -fu1(ifu1)./(dx(ifu1)-du(ifu1)); -fu2(ifu2)./(-dx(ifu2)-du(ifu2)); ...
nikcleju@66:       #  (-bqe+sqrt(bqe^2-4*aqe*cqe))/(2*aqe)
nikcleju@66:       #  ]));
nikcleju@66:       #smax = min(1,numpy.concatenate( (-fu1[ifu1]/(dx[ifu1]-du[ifu1]) , -fu2[ifu2]/(-dx[ifu2]-du[ifu2]) , numpy.array([ (-bqe + math.sqrt(bqe**2-4*aqe*cqe))/(2*aqe) ]) ) , 0).min())
nikcleju@66:       smax = min(1,numpy.array([ (-bqe + math.sqrt(bqe**2-4*aqe*cqe))/(2*aqe) ] ).min())
nikcleju@66:       
nikcleju@66:       s = 0.99 * smax
nikcleju@66:       
nikcleju@66:       # backtracking line search
nikcleju@66:       suffdec = 0
nikcleju@66:       backiter = 0
nikcleju@66:       while not suffdec:
nikcleju@66:         #xp = x + s*dx;  up = u + s*du;  rp = r + s*Adx;
nikcleju@66:         xp = x + s*dx
nikcleju@66:         #up = u + s*du
nikcleju@66:         rp = r + s*Adx
nikcleju@66:         #fu1p = xp - up;  fu2p = -xp - up;  fep = 1/2*(rp'*rp - epsilon^2);
nikcleju@66:         #fu1p = xp - up
nikcleju@66:         #fu2p = -xp - up
nikcleju@66:         fep = 0.5*(numpy.vdot(rp,rp) - epsilon**2)
nikcleju@66:         ##fp = sum(up) - (1/tau)*(sum(log(-fu1p)) + sum(log(-fu2p)) + log(-fep));
nikcleju@66:         #fp = up.sum() - (1.0/tau)*(numpy.log(-fu1p).sum() + numpy.log(-fu2p).sum() + math.log(-fep))
nikcleju@66:         fp = numpy.linalg.norm(p-xp)**2 - (1.0/tau) * math.log(-fep)
nikcleju@66:         #flin = f + alpha*s*(gradf'*[dx; du]);
nikcleju@66:         flin = f + alpha*s*numpy.dot(gradf.T , dx)
nikcleju@66:         #suffdec = (fp <= flin);
nikcleju@66:         if fp <= flin:
nikcleju@66:             suffdec = True
nikcleju@66:         else:
nikcleju@66:             suffdec = False
nikcleju@66:         
nikcleju@66:         s = beta*s
nikcleju@66:         backiter = backiter + 1
nikcleju@66:         if (backiter > 132):
nikcleju@66:           if verbose:
nikcleju@66:             print 'Stuck on backtracking line search, returning previous iterate.  (See Section 4 of notes for more information.)'
nikcleju@66:           xp = x.copy()
nikcleju@66:           #up = u.copy()
nikcleju@66:           #return xp,up,niter
nikcleju@66:           return xp,niter
nikcleju@66:         #end
nikcleju@66:       #end
nikcleju@66:       
nikcleju@66:       # set up for next iteration
nikcleju@66:       ##x = xp; u = up;  r = rp;
nikcleju@66:       x = xp.copy()
nikcleju@66:       #u = up.copy()
nikcleju@66:       r = rp.copy()
nikcleju@66:       ##fu1 = fu1p;  fu2 = fu2p;  fe = fep;  f = fp;
nikcleju@66:       #fu1 = fu1p.copy()
nikcleju@66:       #fu2 = fu2p.copy()
nikcleju@66:       fe = fep
nikcleju@66:       f = fp
nikcleju@66:       
nikcleju@66:       ##lambda2 = -(gradf'*[dx; du]);
nikcleju@66:       #lambda2 = -numpy.dot(gradf.T , numpy.concatenate((dx,du),0))
nikcleju@66:       lambda2 = -numpy.dot(gradf.T , dx)
nikcleju@66:       ##stepsize = s*norm([dx; du]);
nikcleju@66:       #stepsize = s * numpy.linalg.norm(numpy.concatenate((dx,du),0))
nikcleju@66:       stepsize = s * numpy.linalg.norm(dx)
nikcleju@66:       niter = niter + 1
nikcleju@66:       #done = (lambda2/2 < newtontol) | (niter >= newtonmaxiter);
nikcleju@66:       if lambda2/2.0 < newtontol or niter >= newtonmaxiter:
nikcleju@66:           done = 1
nikcleju@66:       else:
nikcleju@66:           done = 0
nikcleju@66:       
nikcleju@66:       #disp(sprintf('Newton iter = #d, Functional = #8.3f, Newton decrement = #8.3f, Stepsize = #8.3e', ...
nikcleju@66:       if verbose:
nikcleju@66:         print 'Newton iter = ',niter,', Functional = ',f,', Newton decrement = ',lambda2/2.0,', Stepsize = ',stepsize
nikcleju@66: 
nikcleju@66:       if verbose:
nikcleju@66:         if largescale:
nikcleju@66:           print '                CG Res = ',cgres,', CG Iter = ',cgiter
nikcleju@66:         else:
nikcleju@66:           print '                  H11p condition number = ',hcond
nikcleju@66:       #end
nikcleju@66:           
nikcleju@66:     #end
nikcleju@66:     #return xp,up,niter
nikcleju@66:     return xp,niter
nikcleju@66: 
nikcleju@66: #function xp = l1qc_logbarrier(x0, A, At, b, epsilon, lbtol, mu, cgtol, cgmaxiter)
nikcleju@66: def l1qc_logbarrier(x0, p, A, At, b, epsilon, lbtol=1e-3, mu=10, cgtol=1e-8, cgmaxiter=200, verbose=False):
nikcleju@66:     # Solve quadratically constrained l1 minimization:
nikcleju@66:     # min ||x||_1   s.t.  ||Ax - b||_2 <= \epsilon
nikcleju@66:     #
nikcleju@66:     # Reformulate as the second-order cone program
nikcleju@66:     # min_{x,u}  sum(u)   s.t.    x - u <= 0,
nikcleju@66:     #                            -x - u <= 0,
nikcleju@66:     #      1/2(||Ax-b||^2 - \epsilon^2) <= 0
nikcleju@66:     # and use a log barrier algorithm.
nikcleju@66:     #
nikcleju@66:     # Usage:  xp = l1qc_logbarrier(x0, A, At, b, epsilon, lbtol, mu, cgtol, cgmaxiter)
nikcleju@66:     #
nikcleju@66:     # x0 - Nx1 vector, initial point.
nikcleju@66:     #
nikcleju@66:     # A - Either a handle to a function that takes a N vector and returns a K 
nikcleju@66:     #     vector , or a KxN matrix.  If A is a function handle, the algorithm
nikcleju@66:     #     operates in "largescale" mode, solving the Newton systems via the
nikcleju@66:     #     Conjugate Gradients algorithm.
nikcleju@66:     #
nikcleju@66:     # At - Handle to a function that takes a K vector and returns an N vector.
nikcleju@66:     #      If A is a KxN matrix, At is ignored.
nikcleju@66:     #
nikcleju@66:     # b - Kx1 vector of observations.
nikcleju@66:     #
nikcleju@66:     # epsilon - scalar, constraint relaxation parameter
nikcleju@66:     #
nikcleju@66:     # lbtol - The log barrier algorithm terminates when the duality gap <= lbtol.
nikcleju@66:     #         Also, the number of log barrier iterations is completely
nikcleju@66:     #         determined by lbtol.
nikcleju@66:     #         Default = 1e-3.
nikcleju@66:     #
nikcleju@66:     # mu - Factor by which to increase the barrier constant at each iteration.
nikcleju@66:     #      Default = 10.
nikcleju@66:     #
nikcleju@66:     # cgtol - Tolerance for Conjugate Gradients; ignored if A is a matrix.
nikcleju@66:     #     Default = 1e-8.
nikcleju@66:     #
nikcleju@66:     # cgmaxiter - Maximum number of iterations for Conjugate Gradients; ignored
nikcleju@66:     #     if A is a matrix.
nikcleju@66:     #     Default = 200.
nikcleju@66:     #
nikcleju@66:     # Written by: Justin Romberg, Caltech
nikcleju@66:     # Email: jrom@acm.caltech.edu
nikcleju@66:     # Created: October 2005
nikcleju@66:     #
nikcleju@66: 
nikcleju@66:     
nikcleju@66:     # Check if epsilon > 0. If epsilon is 0, the algorithm fails. You should run the algo with equality constraint instead  
nikcleju@66:     if epsilon == 0:
nikcleju@66:       raise l1qcInputValueError('Epsilon should be > 0!')       
nikcleju@66:     
nikcleju@66:     #largescale = isa(A,'function_handle');
nikcleju@66:     if hasattr(A, '__call__'):
nikcleju@66:         largescale = True
nikcleju@66:     else:
nikcleju@66:         largescale = False
nikcleju@66:     
nikcleju@66:     #    if (nargin < 6), lbtol = 1e-3; end
nikcleju@66:     #    if (nargin < 7), mu = 10; end
nikcleju@66:     #    if (nargin < 8), cgtol = 1e-8; end
nikcleju@66:     #    if (nargin < 9), cgmaxiter = 200; end
nikcleju@66:     # Nic: added them as optional parameteres
nikcleju@66:     
nikcleju@66:     newtontol = lbtol
nikcleju@66:     newtonmaxiter = 150
nikcleju@66:     
nikcleju@66:     #N = length(x0);
nikcleju@66:     N = x0.size
nikcleju@66:     
nikcleju@66:     # starting point --- make sure that it is feasible
nikcleju@66:     if largescale:
nikcleju@66:       if numpy.linalg.norm(A(x0) - b) > epsilon:
nikcleju@66:         if verbose:
nikcleju@66:           print 'Starting point infeasible; using x0 = At*inv(AAt)*y.'
nikcleju@66:         #AAt = @(z) A(At(z));
nikcleju@66:         AAt = lambda z: A(At(z))
nikcleju@66:         # TODO: implement cgsolve
nikcleju@66:         w,cgres,cgiter = cgsolve(AAt, b, cgtol, cgmaxiter, 0)
nikcleju@66:         if (cgres > 1.0/2):
nikcleju@66:           if verbose:
nikcleju@66:             print 'A*At is ill-conditioned: cannot find starting point'
nikcleju@66:           xp = x0.copy()
nikcleju@66:           return xp
nikcleju@66:         #end
nikcleju@66:         x0 = At(w)
nikcleju@66:       #end
nikcleju@66:     else:
nikcleju@66:       if numpy.linalg.norm( numpy.dot(A,x0) - b ) > epsilon:
nikcleju@66:         if verbose:
nikcleju@66:           print 'Starting point infeasible; using x0 = At*inv(AAt)*y.'
nikcleju@66:         #opts.POSDEF = true; opts.SYM = true;
nikcleju@66:         #[w, hcond] = linsolve(A*A', b, opts);
nikcleju@66:         #if (hcond < 1e-14)
nikcleju@66:         #  disp('A*At is ill-conditioned: cannot find starting point');
nikcleju@66:         #  xp = x0;
nikcleju@66:         #  return;
nikcleju@66:         #end
nikcleju@66:         try:
nikcleju@66:             w = scipy.linalg.solve(numpy.dot(A,A.T), b, sym_pos=True)
nikcleju@66:             #w = numpy.linalg.solve(numpy.dot(A,A.T), b)
nikcleju@66:             hcond = 1.0/numpy.linalg.cond(numpy.dot(A,A.T))
nikcleju@66:         except scipy.linalg.LinAlgError:
nikcleju@66:             if verbose:
nikcleju@66:               print 'A*At is ill-conditioned: cannot find starting point'
nikcleju@66:             xp = x0.copy()
nikcleju@66:             return xp
nikcleju@66:         if hcond < 1e-14:
nikcleju@66:             if verbose:
nikcleju@66:               print 'A*At is ill-conditioned: cannot find starting point'
nikcleju@66:             xp = x0.copy()
nikcleju@66:             return xp           
nikcleju@66:         #x0 = A'*w;
nikcleju@66:         x0 = numpy.dot(A.T, w)
nikcleju@66:       #end  
nikcleju@66:     #end
nikcleju@66:     x = x0.copy()
nikcleju@66:     #u = (0.95)*numpy.abs(x0) + (0.10)*numpy.abs(x0).max()
nikcleju@66:     
nikcleju@66:     #disp(sprintf('Original l1 norm = #.3f, original functional = #.3f', sum(abs(x0)), sum(u)));
nikcleju@66:     if verbose:
nikcleju@66:       #print 'Original l1 norm = ',numpy.abs(x0).sum(),'original functional = ',u.sum()
nikcleju@66:       print 'Original distance ||p-x|| = ',numpy.linalg.norm(p-x)
nikcleju@66:     # choose initial value of tau so that the duality gap after the first
nikcleju@66:     # step will be about the origial norm
nikcleju@66:     #tau = max(((2*N+1.0)/numpy.abs(x0).sum()), 1)
nikcleju@66:     # Nic:
nikcleju@66:     tau = max(((2*N+1.0)/numpy.linalg.norm(p-x0)**2), 1)
nikcleju@66:                                                                                                                               
nikcleju@66:     lbiter = math.ceil((math.log(2*N+1)-math.log(lbtol)-math.log(tau))/math.log(mu))
nikcleju@66:     #disp(sprintf('Number of log barrier iterations = #d\n', lbiter));
nikcleju@66:     if verbose:
nikcleju@66:       print 'Number of log barrier iterations = ',lbiter
nikcleju@66:     
nikcleju@66:     totaliter = 0
nikcleju@66:     
nikcleju@66:     # Added by Nic, to fix some crashing
nikcleju@66:     if lbiter == 0:
nikcleju@66:         xp = numpy.zeros(x0.size)
nikcleju@66:     #end
nikcleju@66:     
nikcleju@66:     #for ii = 1:lbiter
nikcleju@66:     for ii in numpy.arange(lbiter):
nikcleju@66:     
nikcleju@66:       #xp,up,ntiter = l1qc_newton(x, u, A, At, b, epsilon, tau, newtontol, newtonmaxiter, cgtol, cgmaxiter)
nikcleju@66:       xp,ntiter = l1qc_newton(x, p, A, At, b, epsilon, tau, newtontol, newtonmaxiter, cgtol, cgmaxiter, verbose=verbose)
nikcleju@66:       totaliter = totaliter + ntiter
nikcleju@66:       
nikcleju@66:       #disp(sprintf('\nLog barrier iter = #d, l1 = #.3f, functional = #8.3f, tau = #8.3e, total newton iter = #d\n', ...
nikcleju@66:       #  ii, sum(abs(xp)), sum(up), tau, totaliter));
nikcleju@66:       if verbose:
nikcleju@66:         #print 'Log barrier iter = ',ii,', l1 = ',numpy.abs(xp).sum(),', functional = ',up.sum(),', tau = ',tau,', total newton iter = ',totaliter
nikcleju@66:         print 'Log barrier iter = ',ii,', l1 = ',numpy.abs(xp).sum(),', functional = ',numpy.linalg.norm(p-xp),', tau = ',tau,', total newton iter = ',totaliter
nikcleju@66:       x = xp.copy()
nikcleju@66:       #u = up.copy()
nikcleju@66:      
nikcleju@66:       tau = mu*tau
nikcleju@66:       
nikcleju@66:     #end
nikcleju@66:     return xp