% This script will illustrate various classification tasks related to Support Vectors (and Margins).
% To begin our discussion, consider a binary classification task.  The labeling of the two classes is arbitrary, but for the sake of discussion let us label them by -1 and 1.   Each datapoint will have a d dimensional feature vector x_i associated with it; for the sake of illustration, we will use d=2 here, but all equations will be presented for general numbers of dimensions.
% Now, if we are very lucky, our dataset will be perfectly separable by a linear structure.  That is, it might be the case that the feature vectors for class -1 and the feature vectors for class 1 are sufficiently separated in feature space such that there is at least one plane that can be drawn such that all of the class -1 points are on one side of the plane and all of the class 1 points are on the other.  Below is some code that (should) do this in 2D:

nn=100; % number of negative and positive points
np=100;
n=np+nn;
sd=.2; % a value for the standard deviation of the two normals, controls whether or not the groups overlap in feature space.
xn=mvnrnd([-1,0],sd*eye(2),nn); % generate data, using two different means to make the groups somewhat distinct
xp=mvnrnd([2,2],sd*eye(2),np);
scatter(xn(:,1),xn(:,2)) % visualize the data, negative is blue and positive is red
hold on;
scatter(xp(:,1),xp(:,2))
hold off

yvals=[-ones(nn,1);ones(np,1)];
xvals=[xn;xp];

% Though the above code does not choose any particular line to draw to separate these two clusters of datapoints, it should be clear that at least one (and probably and infinite number) of such lines exist.
%
% Now, recall back to our discussion of logistic regression.  We noted then that in scenarios such as these, logistic regression is actually unable to identify a unique linear structure to maximize the log-likelihood.  However, recall also that the LDA method constructs linear decision boundaries between classes, and LDA will have no problem determining a single, optimal (under the assumptions of LDA, at least) line to separate the two classes perfectly.
% But, suppose I don't believe the assumptions of LDA hold for my dataset, but I still would like to determine a unique "optimal" linear structure to separate my perfectly separable data.  To do so, we must make a definition on what be believe "optimal" should mean.  One such approach is to pick a line that will maximize what is called the margin.  We define the margin by first noting that a general linear structure in d dimensional space can be written as the linear equation:
% \beta_0+\beta_1x_1+\beta_2x_2+\ldots+\beta_dx_d=\beta_0+x^T\beta=0
% for some constants \beta_j, j=0,\ldots,d. Further, we note that if we were to scale all of the \beta_j values above by any nonzero constant, it would not alter the linear structure (the equation would still hold if every term were divided by the same nonzero constant).  Hence, we can always choose to represent our linear structure by a set of constants such that the vector of values not including \beta_0 is a unit vector; that is
% \sum_{j=1}^d\beta_j^2=1.
% This is a convenient choice as, under this representation, one can show that for a given feature vector x_i, the quantity \beta_0+x_i^T\beta is the signed perpendicular distance from the linear structure to that feature vector.  We say the distance is signed since it could be positive or negative, indicating which side of the structure ("above" or "below", whatever those mean) the datapoint is on.
% Now, back to our dataset and margins.  Suppose we choose some linear structure and its corresponding \beta_j values that does indeed perfectly separate our data: all of the class 1 points are "above" the line (and have positive signed distance) and all of the class -1 points are "below" the line (with negative signed distance); note that since the labeling of the classes is arbitrary, we can alwasy make this true for any linear structure that perfectly separates the data.  We can then easily find the unsigned distance (so, just the magnitude) of each point to the linear structure by just doing
% y_i\left(\beta_0+x_i^T\beta\right).
% Then, for that linear structure, the margin Mis defined to be the smallest unsigned distance observed over all of the datapoints:
% M=\min_i y_i\left(\beta_0+x_i^T\beta\right).
% Of course, given its definition as the minimum of some set, the following inequality must hold:
% y_i\left(\beta_0+x_i^T\beta\right)\geq M~\forall i.
% Based on all of these definitions, we can finally introduce the maximal margin classifier.  This is the linear structure (defined by its \beta_j) that, for a perfectly separable dataset, maximizes the margin obtained.  That is, it separates the data perfectly while being as far as possible (in a perpendicular sense) from those datapoints that are closest to it.  The datapoints that are closest to the maximal margin classifier (and therefore determine what the margin actually is) are called support vectors.  To find the maximal margin classifier, one needs to solve the following optimization problem:
% \max_{M,\beta_0,\beta} M\\
% \textrm{subject to  } \beta^T\beta=1\\
% y_i\left(\beta_0+x_i^T\beta\right)\geq M~\forall i
% There are actually more efficient ways of representing this same problem, but we'll ignore those.
% The code below will fit a maximal margin classifier to our data.

%find maximal margin classifier.  The variable z used below in the
%objective function and constraints is [M, \beta_0, \beta^T]
conmat=[ones(n,1),-yvals]; %begin constructing the matrix that the inequality constraint uses
for d=1:2
    conmat=[conmat,-yvals.*xvals(:,d)];
end

myobj =@(z) -z(1); %the objective function to minimize, which is just -M

mmc=fmincon(myobj,[0,0,1,0],conmat,zeros(n,1),[],[],[],[],@mycons) %matlab constrained optimization.  See function mycons below for the constraint that beta^Tbeta=1.
xplot=min(xvals(:,1)):.1:max(xvals(:,1));
slope=-mmc(3)/mmc(4);
int=-mmc(2)/mmc(4);
yplot=int+slope*xplot; %the y values of the obtain MMC linear structure
hold on
plot(xplot,yplot)

%this bit finds the support vectors and plots them in green
mm=mmc(1);
unsdist=mm-conmat*mmc';
tmp=1:n;
svloc=abs(unsdist-mm)<10^-3;
nsv=sum(svloc);
sv=tmp(svloc);
scatter(xvals(svloc,1),xvals(svloc,2),'go',"filled")

D=mm*sqrt(1+slope^2);
plot(xplot,yplot+D,'g-',xplot,yplot-D,'g-');
hold off
wrong=sum(unsdist<0)

% The maximal margin classifier is interesting, in that in the end it is completely determined by the support vectors and nothing else.  That is, if we took our dataset and changed the location of any of the points that are not suppor vectors, the classifier would be unaffected (assuming we didn't move the point to within the current margin, in which case it would become a new support vector and replace some existing one(s)).  Contrast this with LDA, whose linear structure depends on the mean and covariance of the full dataset(s), and so is influenced by all data points in some way.
% Another interesting thing about the maximal margin classifier is that, if our data is not perfectly separable, it will tell us that by finding an M that is negative: that is, any given line in that case will result in some datapoints having a negative "unsigned" distance.  You can explore this by changing the dataset above to have a higher variance so that the two groups overlap, then running the MMC part again.
%
% The maximal margin classifier is quite "unstable": adding a single additional datapoint that would lie within the margin of the current classifier (or be misclassified) can have a huge impact on the decision boundary.  You can see this in action by adding an additional datapoint that would lie on the wrong side of the line shown above and rerunning the classifier.

n=n+1;
newx=[0,3]; %some new datapoint.  Try a positive (red) point on the wrong side of the classifier or within the margin, but still leaving a separable dataset
hold on
scatter(newx(1),newx(2),'ro');
hold off

xvals=[xvals;newx];
yvals=[yvals;1];
%go back and rerun the classifier

%
% How can we make a classifier that has the flavor of maximal margin, but is perhaps less rigid in its insistence of having no incorrectly classified datapoints, and therefore possibly less unstable to changes in data?  One way is the Support Vector Classifier.  The only difference between the maximal margin classifier and the support vector classifier is in the final constraint.  In the MMC, the constrain forces all points to be more distant than the margin (and implicitly wants the margin to be positive; i.e, no incorrectly classified points).  The SVC on the other hand allows for some "slack" for each datapoint, where this slack will take into account the ability of points to be inside the margin or incorrectly classified.  In one precise formulation, we introduce n additional slack parameters \epsilon_i\geq 0, and rewrite the final constraint to be
% y_i\left(\beta_0+x_i^T\beta\right)\geq M(1-\epsilon_i)~\forall i
% In this way, any given datapoint could be correctly classified but within the margin (0<\epsilon_i<1) or could be incorrectly classified (\epsilon_i>1), while still allowing for M to be positive.  Note that this problem is not well posed, though, as we could make M as large as we want by just making the slack variables all arbitrarily large.  So, we add another constraint, that limits the total amount of slack that can be used:
% \sum_{i=1}^n\epsilon_i\leq C.
% Here, the value C sets our overall "budget" for slack, and does not allows arbitrarily large M values.  Note that one immediate consequence of a slack budget C is that there can be at most C points misclassified (since those points have \epsilon_i>1).
% If turns out that there are more efficient ways of writing the SVC than that presented above.

%find SV classifier.  The variable z used below in the
%objective function and constraints is [\beta_0, \beta^T, epsilon^T]
Cbudget=10;

conmat=-yvals; %begin constructing the matrix that the inequality constraint uses
for d=1:2
    conmat=[conmat,-yvals.*xvals(:,d)];
end
for i=1:n
    tmp=zeros(n,1);
    tmp(i)=-1;
    conmat=[conmat,tmp];
end
for i=1:n %this forces all epsilon positive
    tmp=zeros(1,n+3);
    tmp(i+3)=-1;
    conmat=[conmat;tmp];
end
tmp=ones(1,n+3); %this forces sum of epsilon <=Cbudget
tmp(1:3)=[0,0,0];
conmat=[conmat;tmp];


myobj =@(z) z(2)^2+z(3)^2; %the objective function to minimize, which is just beta^T beta

options=optimset('MaxFunEvals',1000000,'MaxIter',10000);
svc=fmincon(myobj,[0,1,0,zeros(1,n)],conmat,[-ones(n,1);zeros(n,1);Cbudget],[],[],[],[],[],options) %matlab constrained optimization.  See function mycons below for the constraint that beta^Tbeta=1.

slope=-svc(2)/svc(3);
int=-svc(1)/svc(3);
yplot=int+slope*xplot; %the y values of the obtain MMC linear structure


% scatter(xn(:,1),xn(:,2)) % visualize the data, negative is blue and positive is red
% hold on;
% scatter(xp(:,1),xp(:,2))
hold on
plot(xplot,yplot,'c-')

msvc=1/sqrt(svc(2)^2+svc(3)^2);
D=msvc*sqrt(1+slope^2);
plot(xplot,yplot+D,'k-',xplot,yplot-D,'k-');

hold off

eps=svc(4:end);
inmargin=sum(eps<1&eps>.00001)
wrong=sum(eps>1)

% The SVC is much more resilient to adding extra points that are within the margin or even incorrectlly classified.  Try this out by adding such a point above and then rerunning this part.  Similarly to the MMC, the SVC result does not in any way depend on points that are correctly classified and outside of the margin.  You can also verify this by adding such a point and then checking that the result does not change.
%
% Support Vector Machines are just like Support Vector Classifiers, except that SVMs transform the regular feature vectors via some choice of basis functions, so that the decision boundaries can be nonlinear.  The mathematics behind this is quite interesting and is discussed in class, but beyond the scope of this livescript.  However, Matlab has a built in SVM algorithm, so we can still illustrate their power below.  This example creates a dataset where the classes are fully separable, but not by a linear boundary.  The SVM has no problem with this.  The resulting plot shows several contours that correspond to different expected values for the points there; the actual decision boundary in this case is the contour at level 0, which is one of the ones shown.


nn=100; % number of negative and positive points
np=100;
n=np+nn;
sd=.25; % a value for the standard deviation of the two normals, controls whether or not the groups overlap in feature space.
xn=mvnrnd([-1,0],sd*eye(2),nn/2); % generate data, using two different means to make the groups somewhat distinct
xn=[xn;mvnrnd([5,4],sd*eye(2),nn/2)];
xp=mvnrnd([2,2],sd*eye(2),np);
scatter(xn(:,1),xn(:,2)) % visualize the data, negative is blue and positive is red
hold on;
scatter(xp(:,1),xp(:,2))

yvals=[-ones(nn,1);ones(np,1)];
xvals=[xn;xp];

svm=fitcsvm(xvals,yvals,"KernelFunction","rbf","KernelScale","auto");

[X1,X2]=meshgrid(min(xvals(:,1)):.01:max(xvals(:,1)),min(xvals(:,2)):.01:max(xvals(:,2)));
[~,score]=predict(svm,[X1(:),X2(:)]);
scoregrid=reshape(score(:,2),size(X1,1),size(X2,2));
contour(X1,X2,scoregrid);
contour(X1,X2,scoregrid,[0,0])
hold off

%

function [c,ceq]=mycons(x) %function needed for MMC

c=0;
ceq=x(3:end)*x(3:end)'-1;

return
end