# Classification using NAG Second-order Conic Programming via CVXPY

## Correct Rendering of this notebook

This notebook makes use of the `latex_envs` Jupyter extension for equations and references.  If the LaTeX is not rendering properly in your local installation of Jupyter , it may be because you have not installed this extension.  Details at https://jupyter-contrib-nbextensions.readthedocs.io/en/latest/nbextensions/latex_envs/README.html

## Installing the NAG library and running this notebook

This notebook depends on the NAG library for Python to run. Please read the instructions in the [Readme.md](https://github.com/numericalalgorithmsgroup/NAGPythonExamples/blob/master/local_optimization/Readme.md#install) file to download, install and obtain a licence for the library.

Instruction on how to run the notebook can be found [here](https://github.com/numericalalgorithmsgroup/NAGPythonExamples/blob/master/local_optimization/Readme.md#jupyter).

## Introduction

In this notebook, we demonstrate how NAG second-order conic programming (SOCP) solver can be used to find a quartic polynomial
\vspace{0.1cm}
\begin{equation}\label{poly}
p(x,y) = \sum_{i,j} a_{ij} x^i y^j ~\mbox{with}~ i+j\leq4
\end{equation}
\vspace{0.1cm}
that separates two sets of points in a plane robustly. The optimization model for this problem can be cast as
\vspace{0.1cm}
\begin{equation}\label{prob}
\begin{array}{ll}
\underset{a\in\Re^{15}}{\mbox{minimize}} & \|a\|\\[0.6ex]
\mbox{subject to} & p(x,y)\leq -1, \mbox{for x, y in set 1},\\[0.6ex]
     & p(x,y)\geq 1, \mbox{for x, y in set 2}.
\end{array}
\end{equation}
\vspace{0.1cm}


```python
# Import necessary packages
import numpy as np
import cvxpy as cvx
import matplotlib.pyplot as plt
```

## Data preparation
In this section, we generate two sets of random points used as training data for our model.


```python
# Fix random seed
np.random.seed(3)
# Generate two primary sets of points
# Number of points in each set
n = 80
set_1_base = np.random.uniform(-1.0, 1.0, (n,2))
set_2_base = np.random.uniform(-1.0, 1.0, (n,2))
# Scale the primary points to make set 1 surrounded by set 2
for i in range(n):
    set_1_base[i,0:] = 0.9 * set_1_base[i,0:] * np.random.rand() / \
                       np.linalg.norm(set_1_base[i,0:])
    set_2_base[i,0:] = set_2_base[i,0:] * (1.1 + np.random.rand() / \
                                np.linalg.norm(set_2_base[i,0:]))
# Further process the data to make abnormal shape of sets of points
maxnorm_set_1 = max(np.linalg.norm(set_1_base,axis=1))
set_2_pick = set_2_base[np.linalg.norm(set_2_base,axis=1)>maxnorm_set_1,0:]
set_1_pick = np.concatenate((set_1_base,
                             set_2_base[np.linalg.norm(set_2_base,axis=1)<=
                                        maxnorm_set_1,0:]))
# The shape of set 1 is primarily round, I punch set 1 from the left to make it abnormal
# Feel free to modify punch power to see the resulting graph
punch_power = 1.0
set_1 = set_1_pick[np.linalg.norm(set_1_pick-[-1.0,-0.5],axis=1)>punch_power,0:]
set_2 = np.concatenate((set_2_pick,
                        set_1_pick[np.linalg.norm(set_1_pick-[-1.0,-0.5],axis=1)<=
                                   punch_power,0:]))
```

Now the training data is ready to use. Visualize it.


```python
data = (set_1, set_2)
colors = ("red", "blue")
groups = ("set1", "set2")

fig = plt.figure()
ax = fig.add_subplot(1, 1, 1, alpha=1.0)
for data, color, group in zip(data, colors, groups):
    x = data[0:,0]
    y = data[0:,1]
    ax.scatter(x, y, alpha=0.8, c=color, edgecolors='none', s=10, label=group)
plt.xlim(-2.0, 2.0)
plt.ylim(-2.0, 2.0)
plt.title('Training Data')
plt.legend(loc=2)
plt.show()
```

The number of coefficient for a quartic polynomial is $15$, which is therefore our number of variables.


```python
# Number of variables
nvar = 15
```

Both constraints in (\ref{prob}) are linear with respect to variable $a$. Use set_1 and set_2 to set linear coefficients for those constraints.


```python
A_1 = np.zeros((set_1.shape[0], nvar))
A_2 = np.zeros((set_2.shape[0], nvar))

for i in range(set_1.shape[0]):
    counter = 0
    for j in range(5):
        for k in range(5-j):
            A_1[i,counter] = set_1[i,0]**j * set_1[i,1]**k
            counter = counter + 1
for i in range(set_2.shape[0]):
    counter = 0
    for j in range(5):
        for k in range(5-j):
            A_2[i,counter] = set_2[i,0]**j * set_2[i,1]**k
            counter = counter + 1
```

## Model our problem via CVXPY and solve


```python
# Define decision variables
x = cvx.Variable(nvar)
# Define objective function
objective = cvx.Minimize(cvx.norm(x))
# Define constraints
constraint = [A_1@x <= -1.0, A_2@x >= 1.0]
# Define the entire problem
problem = cvx.Problem(objective, constraint)
# Solve, Bing!
problem.solve(solver='NAG', verbose=True)
# Save result
coef = x.value
```

    
     ------------------------------------------------
      E04PT, Interior point method for SOCP problems
     ------------------------------------------------
    
     Begin of Options
         Print File                    =                  -1     * U
         Print Level                   =                   2     * d
         Print Options                 =                 Yes     * d
         Print Solution                =                  No     * d
         Monitoring File               =                   6     * U
         Monitoring Level              =                   2     * U
         Socp Monitor Frequency        =                   0     * d
    
         Infinite Bound Size           =         1.00000E+20     * d
         Task                          =            Minimize     * d
         Stats Time                    =                  No     * d
         Time Limit                    =         1.00000E+06     * d
    
         Socp Iteration Limit          =                 100     * d
         Socp Max Iterative Refinement =                   9     * d
         Socp Presolve                 =                 Yes     * d
         Socp Scaling                  =                None     * d
         Socp Stop Tolerance           =         1.05367E-08     * d
         Socp Stop Tolerance 2         =         1.05367E-08     * d
         Socp System Formulation       =                Auto     * d
     End of Options
    
     Problem Statistics
       No of variables                 32
         bounds               not defined
       No of lin. constraints         176
         nonzeroes                   2432
       No of quad.constraints           0
       No of cones                      1
         biggest cone size             16
       Objective function          Linear
    
     Presolved Problem Measures
       No of variables                191
       No of lin. constraints         175
         nonzeroes                   2590
       No of cones                      1
    
    
     ------------------------------------------------------------------------
      it|    pobj    |    dobj    |  p.inf  |  d.inf  |  d.gap  |   tau  | I
     ------------------------------------------------------------------------
       0  1.00000E+00  0.00000E+00  2.70E-02  5.35E-03  2.00E+00  1.0E+00
       1  5.13493E+01  8.32818E+01  1.29E-02  2.55E-03  9.55E-01  5.0E-02
       2  2.80351E+01  5.34418E+01  9.73E-03  1.93E-03  7.22E-01  3.6E-01
       3  1.11360E+02  2.17980E+02  4.46E-03  8.86E-04  3.31E-01  1.2E-01
       4  2.07571E+02  4.06965E+02  2.08E-03  4.13E-04  1.55E-01  5.9E-02
       5  3.02366E+02  5.91847E+02  1.38E-03  2.74E-04  1.03E-01  3.5E-02
       6  2.42309E+02  4.51372E+02  2.93E-04  5.82E-05  2.18E-02  2.0E-02
       7  1.29884E+02  1.56324E+02  8.30E-05  1.65E-05  6.16E-03  2.2E-02
       8  7.71822E+01  8.68099E+01  3.34E-05  6.63E-06  2.48E-03  3.4E-02
       9  6.67725E+01  7.13985E+01  1.75E-05  3.48E-06  1.30E-03  3.6E-02
      10  6.49722E+01  6.63180E+01  7.68E-06  1.52E-06  5.70E-04  3.2E-02
      11  6.27050E+01  6.28895E+01  2.58E-06  5.11E-07  1.91E-04  3.2E-02
      12  6.17692E+01  6.18968E+01  1.62E-06  3.21E-07  1.20E-04  2.9E-02
      13  6.14587E+01  6.14705E+01  4.56E-07  9.04E-08  3.38E-05  2.9E-02
      14  6.11965E+01  6.11981E+01  5.12E-08  1.02E-08  3.80E-06  2.8E-02
      15  6.11753E+01  6.11754E+01  2.23E-09  4.43E-10  1.66E-07  2.8E-02
      16  6.11747E+01  6.11747E+01  8.75E-11  1.72E-11  6.42E-09  2.8E-02
     ------------------------------------------------------------------------------
     Status: converged, an optimal solution found
     ------------------------------------------------------------------------------
     Final primal objective value         6.117468E+01
     Final dual objective value           6.117469E+01
     Absolute primal infeasibility        6.493819E-09
     Relative primal infeasibility        8.753236E-11
     Absolute dual infeasibility          3.212062E-09
     Relative dual infeasibility          1.718111E-11
     Absolute complementarity gap         6.424124E-09
     Relative complementarity gap         6.424124E-09
     Iterations                                     16


## Visualize the classifier


```python
# Generate a mesh
x = np.arange(-5.0,5.0,0.008)
y = np.arange(-5.0,5.0,0.008)
xx, yy = np.meshgrid(x,y)

# Ploynomial value on the mesh
polyval = np.zeros(xx.shape)
counter = 0
for i in range(5):
    for j in range(5-i):
        polyval = polyval + coef[counter]*np.power(xx,i)*np.power(yy,j)
        counter = counter + 1

# Plot the trained polynomial
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1, alpha=1.0)

ax.scatter(xx[polyval<=-1], yy[polyval<=-1], alpha=0.002)

data = (set_1, set_2)
colors = ("red", "blue")
groups = ("set1", "set2")
for data, color, group in zip(data, colors, groups):
    x = data[0:,0]
    y = data[0:,1]
    ax.scatter(x, y, alpha=0.8, c=color, edgecolors='none', s=10, label=group)
plt.xlim(-2.0, 2.0)
plt.ylim(-2.0, 2.0)
plt.title('Polynomial')
plt.legend(loc=2)
plt.show()
```

0ae1273a1c42309281bd95f3e445416e2a3c6766
ipynb
Jupyter Notebook
local_optimization/SOCP/cvxpy_classification.ipynb
Brunochris13/NAGPythonExamples
e57fc05ab9b27db66d06a52f9b9412205e984544
BSD-3-Clause
2018-12-06T20:20:01.000Z
2022-03-05T23:09:31.000Z
local_optimization/SOCP/cvxpy_classification.ipynb
kelly1208/NAGPythonExamples
bd20f719c176bbbbc878fea7d0962e5fa10d9d3e
BSD-3-Clause
2019-03-25T11:52:51.000Z
2021-04-12T14:08:31.000Z
local_optimization/SOCP/cvxpy_classification.ipynb
kelly1208/NAGPythonExamples
bd20f719c176bbbbc878fea7d0962e5fa10d9d3e
BSD-3-Clause
2019-01-22T13:30:57.000Z
2021-12-15T13:05:14.000Z
Qwen/Qwen-72B
1. YES
2. YES