#!/usr/bin/python # Figure 6.2, page 297. # Penalty approximation. # # The problem data are not the same as in the book figure. import pylab, numpy from cvxopt import lapack, solvers, matrix, spdiag, log, div, normal from cvxopt.modeling import variable, op, max, sum solvers.options['show_progress'] = 0 m, n = 100, 30 A = normal(m,n) b = normal(m,1) b /= (1.1 * max(abs(b))) # Make x = 0 feasible for log barrier. # l1 approximation # # minimize || A*x + b ||_1 x = variable(n) op(sum(abs(A*x+b))).solve() x1 = x.value pylab.figure(1, facecolor='w', figsize=(10,10)) pylab.subplot(411) nbins = 100 bins = [-1.5 + 3.0/(nbins-1)*k for k in xrange(nbins)] pylab.hist( A*x1+b , numpy.array(bins)) nopts = 200 xs = -1.5 + 3.0/(nopts-1) * matrix(range(nopts)) pylab.plot(xs, (35.0/1.5) * abs(xs), 'g-') pylab.axis([-1.5, 1.5, 0, 40]) pylab.ylabel('l1') pylab.title('Penalty function approximation (fig. 6.2)') # l2 approximation # # minimize || A*x + b ||_2 x = matrix(0.0, (m,1)) lapack.gels(+A, x) x2 = x[:n] pylab.subplot(412) pylab.hist(A*x2+b, numpy.array(bins)) pylab.plot(xs, (8.0/1.5**2) * xs**2 , 'g-') pylab.ylabel('l2') pylab.axis([-1.5, 1.5, 0, 10]) # Deadzone approximation # # minimize sum(max(abs(A*x+b)-0.5, 0.0)) x = variable(n) op(sum(max(abs(A*x+b)-0.5, 0.0))).solve() xdz = x.value pylab.subplot(413) pylab.hist(A*xdz+b, numpy.array(bins)) pylab.plot(xs, 15.0/1.0 * matrix([ max(abs(xk)-0.5, 0.0) for xk in xs ]), 'g-') pylab.ylabel('Deadzone') pylab.axis([-1.5, 1.5, 0, 20]) # Log barrier # # minimize -sum (log ( 1.0 - A*x+b)**2) def F(x=None, z=None): if x is None: return 0, matrix(0.0, (n,1)) y = A*x+b if max(abs(y)) >= 1.0: return None f = -sum(log(1.0 - y**2)) gradf = 2.0 * A.T * div(y, 1-y**2) if z is None: return f, gradf.T H = A.T * spdiag(2.0 * z[0] * div( 1.0+y**2, (1.0 - y**2)**2 )) * A return f, gradf.T, H xlb = solvers.cp(F)['x'] pylab.subplot(414) pylab.hist(A*xlb+b, numpy.array(bins)) nopts = 200 pylab.plot(xs, (8.0/1.5**2) * xs**2, 'g--') xs2 = -0.99999 + (2*0.99999 /(nopts-1)) * matrix(range(nopts)) pylab.plot(xs2, -3.0 * log(1.0 - abs(xs2)**2), 'g-') pylab.ylabel('Log barrier') pylab.xlabel('residual') pylab.axis([-1.5, 1.5, 0, 10]) pylab.show()