# Support vector machine (SVM)¶

Support vector machine seeks an affine function that approximately classifies the two sets of points. The problem can be stated as

$\begin{array}{ll} \mbox{minimize} & \frac{1}{2} x^T x + \gamma \sum_{i=1}^{m} \max(0, b_i a_i^T x + 1), \end{array}$

where $$b_i \in \{ -1, +1 \}$$ is a set label, and $$a_i$$ is a vector of features for the $$i$$-th point. The problem has the following equivalent form

$\begin{split}\begin{array}{ll} \mbox{minimize} & \frac{1}{2} x^T x + \gamma \boldsymbol{1}^T t \\ \mbox{subject to} & t \ge {\rm diag}(b) Ax + 1 \\ & t \ge 0, \end{array}\end{split}$

where $${\rm diag}(b)$$ denotes the diagonal matrix with elements of $$b$$ on its diagonal.

## Python¶

import osqp
import numpy as np
import scipy as sp
from scipy import sparse

# Generate problem data
sp.random.seed(1)
n = 10
m = 1000
N = int(m / 2)
gamma = 1.0
b = np.hstack([np.ones(N), -np.ones(N)])
A_upp = sparse.random(N, n, density=0.5)
A_low = sparse.random(N, n, density=0.5)
A_upp / np.sqrt(n) + (A_upp != 0.).astype(float) / n,
A_low / np.sqrt(n) - (A_low != 0.).astype(float) / n
], format='csc')

# OSQP data
Im = sparse.eye(m)
P = sparse.block_diag([sparse.eye(n), sparse.csc_matrix((m, m))], format='csc')
q = np.hstack([np.zeros(n), gamma*np.ones(m)])
[None,                Im]], format='csc')
l = np.hstack([-np.inf*np.ones(m), np.zeros(m)])
u = np.hstack([-np.ones(m), np.inf*np.ones(m)])

# Create an OSQP object
prob = osqp.OSQP()

# Setup workspace
prob.setup(P, q, A, l, u)

# Solve problem
res = prob.solve()


## Matlab¶

% Generate problem data
rng(1)
n = 10;
m = 1000;
N = ceil(m/2);
gamma = 1;
A_upp = sprandn(N, n, 0.5);
A_low = sprandn(N, n, 0.5);
Ad = [A_upp / sqrt(n) + (A_upp ~= 0) / n;
A_low / sqrt(n) - (A_low ~= 0) / n];
b = [ones(N, 1); -ones(N,1)];

% OSQP data
P = blkdiag(speye(n), sparse(m, m));
q = [zeros(n,1); gamma*ones(m,1)];
sparse(m, n), speye(m)];
l = [-inf*ones(m, 1); zeros(m, 1)];
u = [-ones(m, 1); inf*ones(m, 1)];

% Create an OSQP object
prob = osqp;

% Setup workspace
prob.setup(P, q, A, l, u);

% Solve problem
res = prob.solve();


## CVXPY¶

from cvxpy import *
import numpy as np
import scipy as sp
from scipy import sparse

# Generate problem data
sp.random.seed(1)
n = 10
m = 1000
N = int(m / 2)
gamma = 1.0
b = np.hstack([np.ones(N), -np.ones(N)])
A_upp = sparse.random(N, n, density=0.5)
A_low = sparse.random(N, n, density=0.5)
A = sparse.vstack([
A_upp / np.sqrt(n) + (A_upp != 0.).astype(float) / n,
A_low / np.sqrt(n) - (A_low != 0.).astype(float) / n
], format='csc')

# Define problem
x = Variable(n)
objective = 0.5*sum_squares(x) + gamma*sum(pos(diag(b)@A@x + 1))

# Solve with OSQP
Problem(Minimize(objective)).solve(solver=OSQP)


## YALMIP¶

% Generate problem data
rng(1)
n = 10;
m = 1000;
N = ceil(m/2);
gamma = 1;
A_upp = sprandn(N, n, 0.5);
A_low = sprandn(N, n, 0.5);
A = [A_upp / sqrt(n) + (A_upp ~= 0) / n;
A_low / sqrt(n) - (A_low ~= 0) / n];
b = [ones(N, 1); -ones(N,1)];

% Define problem
x = sdpvar(n, 1);
objective = 0.5*norm(x)^2 + gamma*sum(max(diag(b)*A*x + 1, 0));

% Solve with OSQP
options = sdpsettings('solver','osqp');
optimize([],objective,options);