Polynomial discrimination
rand('state',0);
N = 100;
M = 120;
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);
p1 = [0,0,1,0,1,2,0,1,2,3,0,1,2,3,4]';
p2 = [0,1,1,2,2,2,3,3,3,3,4,4,4,4,4]'-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));
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;
norm(a) <= 1;
cvx_end
fprintf(1,'Done! \n');
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('Optimal order-4 polynomial that separates the 2 classes')
Finding the optimal polynomial of order 4 that separates the 2 classes...
Calling SDPT3: 227 variables, 16 equality constraints
For improved efficiency, SDPT3 is solving the dual problem.
------------------------------------------------------------
num. of constraints = 16
dim. of socp var = 16, num. of socp blk = 1
dim. of linear var = 211
*******************************************************************
SDPT3: Infeasible path-following algorithms
*******************************************************************
version predcorr gam expon scale_data
NT 1 0.000 1 0
it pstep dstep pinfeas dinfeas gap mean(obj) cputime
-------------------------------------------------------------------
0|0.000|0.000|4.3e+03|3.8e+02|1.6e+05| 4.000000e+00| 0:0:00| chol 1 1
1|0.989|1.000|4.8e+01|1.0e-01|1.9e+03|-2.388888e+01| 0:0:00| chol 1 1
2|0.881|1.000|5.7e+00|1.0e-02|2.6e+02|-2.457008e+01| 0:0:00| chol 1 1
3|1.000|0.577|3.4e-08|4.8e-03|5.2e+01|-2.322700e+01| 0:0:00| chol 1 1
4|0.938|0.930|2.3e-07|4.3e-04|3.8e+00| 2.108478e-01| 0:0:00| chol 1 1
5|0.481|0.362|1.2e-07|2.8e-04|3.0e+00|-2.693317e-01| 0:0:00| chol 1 1
6|1.000|0.358|2.7e-09|1.8e-04|2.3e+00|-4.524383e-01| 0:0:00| chol 1 1
7|0.602|1.000|1.6e-09|1.0e-07|1.1e+00|-2.193218e-02| 0:0:00| chol 1 1
8|0.752|0.531|3.9e-10|5.3e-08|8.0e-01|-2.139053e-02| 0:0:00| chol 1 1
9|0.492|0.924|2.0e-10|5.0e-09|6.2e-01| 5.195829e-03| 0:0:00| chol 1 1
10|0.757|0.981|4.9e-11|2.3e-10|3.0e-01| 9.452139e-02| 0:0:00| chol 1 1
11|0.962|1.000|1.8e-12|2.0e-11|1.3e-01| 6.607111e-02| 0:0:00| chol 1 1
12|1.000|1.000|6.3e-16|2.0e-12|7.1e-02| 5.511548e-02| 0:0:00| chol 1 1
13|1.000|1.000|6.9e-16|1.1e-12|2.5e-02| 4.750589e-02| 0:0:00| chol 1 1
14|0.969|0.951|1.2e-15|1.1e-12|8.2e-03| 4.376829e-02| 0:0:00| chol 1 1
15|1.000|0.868|4.0e-15|1.1e-12|4.7e-03| 4.294058e-02| 0:0:00| chol 1 1
16|0.936|0.912|5.7e-16|1.1e-12|7.1e-04| 4.127616e-02| 0:0:00| chol 1 1
17|1.000|0.998|7.3e-14|1.0e-12|1.6e-04| 4.106513e-02| 0:0:00| chol 1 1
18|1.000|1.000|3.3e-14|1.0e-12|2.6e-05| 4.100493e-02| 0:0:00| chol 1 1
19|1.000|1.000|2.2e-13|1.0e-12|2.0e-06| 4.099416e-02| 0:0:01| chol 1 1
20|1.000|1.000|4.0e-14|1.0e-12|3.9e-08| 4.099324e-02| 0:0:01| chol 1 1
21|0.999|0.998|2.7e-13|1.0e-12|4.3e-10| 4.099323e-02| 0:0:01|
stop: max(relative gap, infeasibilities) < 1.49e-08
-------------------------------------------------------------------
number of iterations = 21
primal objective value = 4.09932258e-02
dual objective value = 4.09932254e-02
gap := trace(XZ) = 4.28e-10
relative gap = 3.96e-10
actual relative gap = 3.95e-10
rel. primal infeas = 2.74e-13
rel. dual infeas = 1.00e-12
norm(X), norm(y), norm(Z) = 5.6e-01, 1.0e+00, 2.1e+01
norm(A), norm(b), norm(C) = 1.0e+02, 2.0e+00, 2.0e+00
Total CPU time (secs) = 0.6
CPU time per iteration = 0.0
termination code = 0
DIMACS: 2.7e-13 0.0e+00 1.0e-12 0.0e+00 4.0e-10 4.0e-10
-------------------------------------------------------------------
------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): -0.0409932
Done!