Homec4science
Diffusion EMPoWER (master) Browse

R1252/simulation-code/cvx/examples/cvxbook/Ch08_geometric_probs/html/quad_discr.html
de096106547amaster

Last Change
Raw File

quad_discr.html
View Options

<!DOCTYPE HTML>
<html>
<head>
<meta charset="UTF-8">
<title>Quadratic discrimination (separating ellipsoid)</title>
<link rel="canonical" href="http://cvxr.com/cvx/examples/cvxbook/Ch08_geometric_probs/html/quad_discr.html">
<link rel="stylesheet" href="../../../examples.css" type="text/css">
</head>
<body>
<div id="header">
<h1>Quadratic discrimination (separating ellipsoid)</h1>
Jump to:&nbsp;&nbsp;&nbsp;&nbsp;
<a href="#source">Source code</a>&nbsp;&nbsp;&nbsp;&nbsp;
<a href="#output">Text output</a>
&nbsp;&nbsp;&nbsp;&nbsp;
<a href="#plots">Plots</a>
&nbsp;&nbsp;&nbsp;&nbsp;<a href="../../../index.html">Library index</a>
</div>
<div id="content">
<a id="source"></a>
<pre class="codeinput">
<span class="comment">% Section 8.6.2, Boyd &amp; Vandenberghe "Convex Optimization"</span>
<span class="comment">% Original by Lieven Vandenberghe</span>
<span class="comment">% Adapted for CVX by Joelle Skaf - 10/16/05</span>
<span class="comment">% (a figure is generated)</span>
<span class="comment">%</span>
<span class="comment">% The goal is to find an ellipsoid that contains all the points</span>
<span class="comment">% x_1,...,x_N but none of the points y_1,...,y_M. The equation of the</span>
<span class="comment">% ellipsoidal surface is: z'*P*z + q'*z + r =0</span>
<span class="comment">% P, q and r can be obtained by solving the SDP feasibility problem:</span>
<span class="comment">% minimize 0</span>
<span class="comment">% s.t. x_i'*P*x_i + q'*x_i + r &gt;= 1 for i = 1,...,N</span>
<span class="comment">% y_i'*P*y_i + q'*y_i + r &lt;= -1 for i = 1,...,M</span>
<span class="comment">% P &lt;= -I</span>
<span class="comment">% data generation</span>
n = 2;
rand(<span class="string">'state'</span>,0); randn(<span class="string">'state'</span>,0);
N=50;
X = randn(2,N); X = X*diag(0.99*rand(1,N)./sqrt(sum(X.^2)));
Y = randn(2,N); Y = Y*diag((1.02+rand(1,N))./sqrt(sum(Y.^2)));
T = [1 -1; 2 1]; X = T*X; Y = T*Y;
<span class="comment">% Solution via CVX</span>
fprintf(1,<span class="string">'Find the optimal ellipsoid that seperates the 2 classes...'</span>);
cvx_begin <span class="string">sdp</span>
variable <span class="string">P(n,n)</span> <span class="string">symmetric</span>
variables <span class="string">q(n)</span> <span class="string">r(1)</span>
P &lt;= -eye(n);
sum((X'*P).*X',2) + X'*q + r &gt;= +1;
sum((Y'*P).*Y',2) + Y'*q + r &lt;= -1;
cvx_end
fprintf(1,<span class="string">'Done! \n'</span>);
<span class="comment">% Displaying results</span>
r = -r; P = -P; q = -q;
c = 0.25*q'*inv(P)*q - r;
xc = -0.5*inv(P)*q;
nopts = 1000;
angles = linspace(0,2*pi,nopts);
ell = inv(sqrtm(P/c))*[cos(angles); sin(angles)] + repmat(xc,1,nopts);
graph=plot(X(1,:),X(2,:),<span class="string">'o'</span>, Y(1,:), Y(2,:),<span class="string">'o'</span>, ell(1,:), ell(2,:),<span class="string">'-'</span>);
set(graph(2),<span class="string">'MarkerFaceColor'</span>,[0 0.5 0]);
set(gca,<span class="string">'XTick'</span>,[]); set(gca,<span class="string">'YTick'</span>,[]);
title(<span class="string">'Quadratic discrimination'</span>);
<span class="comment">% print -deps ellips.eps</span>
</pre>
<a id="output"></a>
<pre class="codeoutput">
Find the optimal ellipsoid that seperates the 2 classes...
Calling SDPT3: 103 variables, 6 equality constraints
For improved efficiency, SDPT3 is solving the dual problem.
------------------------------------------------------------
num. of constraints = 6
dim. of sdp var = 2, num. of sdp blk = 1
dim. of linear var = 100
*******************************************************************
SDPT3: Infeasible path-following algorithms
*******************************************************************
version predcorr gam expon scale_data
HKM 1 0.000 1 0
it pstep dstep pinfeas dinfeas gap prim-obj dual-obj cputime
-------------------------------------------------------------------
0|0.000|0.000|2.9e+03|9.6e+00|6.4e+04| 2.626414e+03 0.000000e+00| 0:0:00| chol 1 1
1|0.636|0.827|1.0e+03|1.7e+00|3.0e+04| 7.649394e+02 0.000000e+00| 0:0:00| chol 1 1
2|0.797|1.000|2.1e+02|1.6e-03|6.7e+03| 1.293218e+02 0.000000e+00| 0:0:00| chol 1 1
3|0.987|1.000|2.7e+00|1.6e-04|8.4e+01| 1.627774e+00 0.000000e+00| 0:0:00| chol 1 1
4|0.989|1.000|2.9e-02|1.6e-05|9.3e-01| 1.791336e-02 0.000000e+00| 0:0:00| chol 1 1
5|0.989|1.000|3.2e-04|5.9e-03|1.0e-02| 1.968772e-04 0.000000e+00| 0:0:00| chol 1 1
6|0.989|1.000|3.6e-06|6.5e-05|1.1e-04| 2.176872e-06 0.000000e+00| 0:0:00| chol 1 1
7|0.989|1.000|3.9e-08|7.3e-07|1.2e-06| 2.392715e-08 0.000000e+00| 0:0:00| chol 1 1
8|0.983|1.000|6.9e-10|7.9e-09|2.2e-08| 4.410963e-10 0.000000e+00| 0:0:00| chol 1 1
9|0.985|1.000|1.1e-11|1.4e-10|3.4e-10| 6.901849e-12 0.000000e+00| 0:0:00|
stop: max(relative gap, infeasibilities) &lt; 1.49e-08
-------------------------------------------------------------------
number of iterations = 9
primal objective value = 6.90184871e-12
dual objective value = 0.00000000e+00
gap := trace(XZ) = 3.40e-10
relative gap = 3.40e-10
actual relative gap = 6.90e-12
rel. primal infeas = 1.06e-11
rel. dual infeas = 1.38e-10
norm(X), norm(y), norm(Z) = 9.9e-13, 1.2e+02, 1.1e+03
norm(A), norm(b), norm(C) = 7.8e+01, 1.0e+00, 6.4e+01
Total CPU time (secs) = 0.13
CPU time per iteration = 0.01
termination code = 0
DIMACS: 1.1e-11 0.0e+00 4.5e-10 0.0e+00 6.9e-12 3.4e-10
-------------------------------------------------------------------
------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): -6.90185e-12
Done!
</pre>
<a id="plots"></a>
<div id="plotoutput">
<img src="quad_discr__01.png" alt="">
</div>
</div>
</body>
</html>
c4science · Help