#!/usr/bin/python # Figure 4.11, page 185. # Regularized least-squares. import pylab from pickle import load from cvxopt import blas, lapack, matrix from cvxopt import solvers solvers.options['show_progress'] = 0 data = load(open("rls.bin",'r')) A, b = data['A'], data['b'] m, n = A.size # LS solution xls = +b lapack.gels(+A, xls) xls = xls[:n] # We compute the optimal values of # # minimize/maximize || A*x - b ||_2^2 # subject to x'*x = alpha # # via the duals. # # Lower bound: # # maximize -t - u*alpha # subject to [u*I, 0; 0, t] + [A, b]'*[A, b] >= 0 # # Upper bound: # # minimize t + u*alpha # subject to [u*I, 0; 0, t] - [A, b]'*[A, b] >= 0. # # Two variables (t, u). G = matrix(0.0, ((n+1)**2, 2)) G[-1, 0] = -1.0 # coefficient of t G[: (n+1)**2-1 : n+2, 1] = -1.0 # coefficient of u h = matrix( [ [ A.T * A, b.T * A ], [ A.T * b, b.T * b ] ] ) c = matrix(1.0, (2,1)) nopts = 40 alpha1 = [ 2.0/(nopts/2-1) * alpha for alpha in xrange(nopts/2) ] + \ [ 2.0 + (15.0 - 2.0)/(nopts/2) * alpha for alpha in xrange(1,nopts/2+1) ] lbnds = [ blas.nrm2(b)**2 ] for alpha in alpha1[1:]: c[1:] = alpha lbnds += [ -blas.dot(c, solvers.sdp(c, Gs=[G], hs=[h])['x']) ] nopts = 10 alpha2 = [ 1.0/(nopts-1) * alpha for alpha in xrange(nopts) ] ubnds = [ blas.nrm2(b)**2 ] for alpha in alpha2[1:]: c[1:] = alpha ubnds += [ blas.dot(c, solvers.sdp(c, Gs=[G], hs=[-h])['x']) ] pylab.figure(1, facecolor='w') pylab.plot(lbnds, alpha1, 'b-', ubnds, alpha2, 'b-') kmax = max([ k for k in xrange(len(alpha1)) if alpha1[k] < blas.nrm2(xls)**2 ]) pylab.plot( [ blas.nrm2(b)**2 ] + lbnds[:kmax] + [ blas.nrm2(A*xls-b)**2 ], [0.0] + alpha1[:kmax] + [ blas.nrm2(xls)**2 ], '-', linewidth=2) pylab.plot([ blas.nrm2(b)**2, blas.nrm2(A*xls-b)**2 ], [0.0, blas.nrm2(xls)**2], 'bo') pylab.fill(lbnds[-1::-1] + ubnds + [ubnds[-1]], alpha1[-1::-1] + alpha2+ [alpha1[-1]], facecolor = '#D0D0D0') pylab.axis([0, 15, -1.0, 15]) pylab.xlabel('||A*x-b||_2^2') pylab.ylabel('||x||_2^2') pylab.grid() pylab.title('Regularized least-squares (fig. 4.11)') pylab.show()