# Polynomial discrimination

Jump to: Source code Text output Plots Library index

    % Section 8.6.2, Boyd & Vandenberghe "Convex Optimization"
    % Original by Lieven Vandenberghe
    % Adapted for CVX by Joelle Skaf - 10/23/05
    % (a figure is generated)
    %
    % The goal is to find the polynomial of degree 3 on R^n that separates
    % two sets of points {x_1,...,x_N} and {y_1,...,y_N}. We are trying to find
    % the coefficients of an order-3-polynomial P(x) that would satisfy:
    %           minimize    t
    %               s.t.    P(x_i) <= t  for i = 1,...,N
    %                       P(y_i) >= t   for i = 1,...,M

    % Data generation
    rand('state',0);
    N = 100;
    M = 120;

    % The points X lie within a circle of radius 0.9, with a wedge of points
    % near [1.1,0] removed. The points Y lie outside a circle of radius 1.1,
    % with a wedge of points near [1.1,0] added. The wedges are precisely what
    % makes the separation difficult and interesting.
    X = 2 * rand(2,N) - 1;
    X = X * diag(0.9*rand(1,N)./sqrt(sum(X.^2)));
    Y = 2 * rand(2,M) - 1;
    Y = Y * diag((1.1+rand(1,M))./sqrt(sum(Y.^2)));
    d = sqrt(sum((X-[1.1;0]*ones(1,N)).^2));
    Y = [ Y, X(:,d<0.9) ];
    X = X(:,d>1);
    N = size(X,2);
    M = size(Y,2);

    % Construct Vandermonde-style monomial matrices
    p1   = [0,0,1,0,1,2,0,1,2,3]';
    p2   = [0,1,1,2,2,2,3,3,3,3]'-p1;
    np   = length(p1);
    op   = ones(np,1);
    monX = X(op,:) .^ p1(:,ones(1,N)) .* X(2*op,:) .^ p2(:,ones(1,N));
    monY = Y(op,:) .^ p1(:,ones(1,M)) .* Y(2*op,:) .^ p2(:,ones(1,M));

    % Solution via CVX
    fprintf(1,'Finding the optimal polynomial of order 4 that separates the 2 classes...');

    cvx_begin
        variables a(np) t(1)
        minimize ( t )
        a'*monX <= t;
        a'*monY >= -t;
    cvx_end

    fprintf(1,'Done! \n');

    % Displaying results
    nopts = 2000;
    angles = linspace(0,2*pi,nopts);
    cont = zeros(2,nopts);
    for i=1:nopts
       v = [cos(angles(i)); sin(angles(i))];
       l = 0;  u = 1;
       while ( u - l > 1e-3 )
          s = (u+l)/2;
          x = s * v;
          if a' * ( x(op,:) .^ p1 .* x(2*op) .^ p2 ) > 0,
              u = s;
          else
              l = s;
          end
       end;
       s = (u+l)/2;
       cont(:,i) = s*v;
    end;

    graph=plot(X(1,:),X(2,:),'o', Y(1,:), Y(2,:),'o', cont(1,:), cont(2,:), '-');
    set(graph(2),'MarkerFaceColor',[0 0.5 0]);
    title('No cubic polynomial can separate the 2 classes')

    % Results
    disp('-----------------------------------------------------------------');
    disp('As seen on the figure, the 2 sets of points are not separated.   ');
    disp('There exists no cubic polynomial that can separate these 2 sets.');

    Finding the optimal polynomial of order 4 that separates the 2 classes...
    Calling sedumi: 211 variables, 11 equality constraints
       For improved efficiency, sedumi is solving the dual problem.
    ------------------------------------------------------------
    SeDuMi 1.21 by AdvOL, 2005-2008 and Jos F. Sturm, 1998-2003.
    Alg = 2: xz-corrector, Adaptive Step-Differentiation, theta = 0.250, beta = 0.500
    eqs m = 11, order n = 212, dim = 212, blocks = 1
    nnz(A) = 2321 + 0, nnz(ADA) = 121, nnz(L) = 66
     it :     b*y       gap    delta  rate   t/tP*  t/tD*   feas cg cg  prec
      0 :            4.26E+02 0.000
      1 :  -3.39E-01 1.69E+02 0.000 0.3957 0.9000 0.9000   1.99  1  1  9.7E+01
      2 :  -5.79E-02 1.23E+02 0.000 0.7309 0.9000 0.9000   9.56  1  1  1.3E+01
      3 :  -3.90E-03 5.32E+01 0.000 0.4318 0.9000 0.9000   9.12  1  1  1.6E+00
      4 :  -2.19E-03 2.29E+01 0.000 0.4309 0.9000 0.9000   0.93  1  1  7.3E-01
      5 :  -1.10E-03 1.10E+01 0.000 0.4807 0.9000 0.9000   1.56  1  1  2.7E-01
      6 :  -3.97E-04 5.29E+00 0.000 0.4804 0.9000 0.9000   1.77  1  1  9.5E-02
      7 :  -1.54E-04 2.84E+00 0.000 0.5359 0.9000 0.9000   1.99  1  1  3.6E-02
      8 :  -3.21E-05 8.44E-01 0.000 0.2975 0.9000 0.9000   1.89  1  1  7.1E-03
      9 :  -6.12E-07 1.86E-02 0.000 0.0221 0.9900 0.9900   1.29  1  1  1.4E-04
     10 :  -3.35E-12 1.03E-07 0.088 0.0000 1.0000 1.0000   1.01  1  1
    iter seconds digits       c*x               b*y
     10      0.2  17.6  0.0000000000e+00 -3.3652985170e-26
    |Ax-b| =   5.2e-16, [Ay-c]_+ =   1.7E-22, |x|=  4.8e-01, |y|=  2.8e-23

    Detailed timing (sec)
       Pre          IPM          Post
    3.000E-02    1.700E-01    1.000E-02
    Max-norms: ||b||=1, ||c|| = 0,
    Cholesky |add|=0, |skip| = 0, ||L.L|| = 13.8426.
    ------------------------------------------------------------
    Status: Solved
    Optimal value (cvx_optval): +3.3653e-26
    Done!
    -----------------------------------------------------------------
    As seen on the figure, the 2 sets of points are not separated.
    There exists no cubic polynomial that can separate these 2 sets.