Simple SVM in MATLAB

Question

I'm studying SVMs and wrote a demo in MATLAB (because I couldn't get a quadratic programming package to work correctly in Python). Right now it's simple and can only do linearly-separable cases (nonlinear kernels will be implemented later). This is the first MATLAB program I've written outside of a MATLAB seminar I'm taking, so any thoughts would be appreciated.

% simple SVM demo
clc;clear all;

% input data; y's are labels for data
y = [ -1; -1; -1; -1; -1; 1; 1; 1; 1; 1; 1 ];
data = [
    0 3;
    1 2;
    2 1;
    3 3;
    5 1;
    5 5;
    6 5;
    6 7;
    7 3;
    8 6;
    8 1;
];

% separate positive and negative datapoints
posDps = data(find(y == 1), :);
negDps = data(find(y == -1), :);

% space is number of dimensions (n-space)
% right now, b/c being plotted in 2D, can only be 2-space
space = size(data, 2);
% n is number of input variables
n = length(data);
[ y1, y2 ] = meshgrid(y, y);
[ i, j ] = meshgrid(1:n, 1:n);

% generate matrices for minimization
P(i, j) = y1(i, j) .* y2(i, j) .* (data(i,:) * data(j,:)');
q = -1 * ones(n, 1);

% generate matrices for inequality constraint (alpha >= 0)
% also can use LB parameter
A = -1 * eye(n);
b = zeros(n, 1);

% generate matrices for equality constraint
Aeq = y';
beq = [ 0 ];

% use quadprog package to generate alphas
alpha = quadprog(P, q, A, b, Aeq, beq, [], [], [], optimoptions('quadprog', 'Display', 'off'));
% calculate w from alpha
w = (data' .* repmat(y', [space 1])) * alpha;

% finding b parameter
% (y_n)(x_n * w + b) = 1 for support vector, so b = y_n - w * x_n
threshold = 1e-5;
svIndices = find(alpha > threshold);
b = y(svIndices(1)) - data(svIndices(1),:) * w;

% display points
figure;
hold on;
scatter(posDps(:,1), posDps(:,2));
scatter(negDps(:,1), negDps(:,2));

% draw line (only 2D for now)
margin = 1;
domain = (min(data(:,1)) - margin):(max(data(:,1)) + margin);
plot(domain, (w(1) .* domain + b)/(-w(2)));

% plotting gutters
plot(domain, (-1 + w(1) .* domain + b)/(-w(2)), 'g:');
plot(domain, (1 + w(1) .* domain + b)/(-w(2)), 'g:');

I'd like any feedback, but here are some specific questions I had in mind:

1. For storing one-dimensional data (arrays), is there usually a preference between column and row matrices?

2. Is there an easy way to deal with nested matrices? I tried to combine the labels and attribute vectors like so:

   data = [
       -1 [ 0 3 ];
       -1 [ 1 2 ];
       % ...
   ]
   

   and use more indices (e.g., data(1,2,2)) but that didn't work.

3. There are a ton of different ways to create matrices, and I'm not sure if the ones I am using are the most appropriate (e.g., via literals, find(), meshgrid(), ones() and zeros(), repmat(), etc.).
4. What is there to do about arbitrary imprecision? For example, I had to put in an arbitrary value threshold to find non-zero values because using find(A > 0) wasn't working. (Since the values were on the magnitude of about 1e-12, this was larger than the value of eps() so that function didn't seem very helpful.)
5. Naming conventions -- is there one for MATLAB? I don't believe I've seen a consistent style.

Answer 1

First to answer your specific questions:

1. For storing one-dimensional data (arrays), is there usually a preference between column and row matrices?

It doesn't matter, use what is most convenient. Often horizontal vectors are used because they are easy to create: 1:10.

2. Is there an easy way to deal with nested matrices?

The best way to collect matrices is using a cell array. A cell array's elements are arbitrary values (matrices, other cell arrays, objects, function handles, whatever you can assign to a variable).

3. There are a ton of different ways to create matrices, and I'm not sure if the ones I am using are the most appropriate (e.g., via literals, find(), meshgrid(), ones() and zeros(), repmat(), etc.).

Yes, MATLAB has a huge library and it is impossible to know if you are using the best possible function in each case. Experience will help.

4. What is there to do about arbitrary imprecision? For example, I had to put in an arbitrary value threshold to find non-zero values because using find(A > 0) wasn't working. (Since the values were on the magnitude of about 1e-12, this was larger than the value of eps() so that function didn't seem very helpful.)

Yes, you need to pick some threshold there. It depends on the source of the data what an appropriate value is. eps is the distance between 1.0 and the next representable floating-point value. That is, there are no possible values in between 1.0 and 1.0+eps. eps is not necessarily an appropriate threshold to find "approximately zero" values.

5. Naming conventions -- is there one for MATLAB? I don't believe I've seen a consistent style.

Not that I know of. Early versions of MATLAB were not case sensitive, so all code was typically written in lower case. Currently some built-in functionality uses an upper case for the first letter, and some use camelCase. There is no PEP8 for MATLAB. :)

---

Now I'll go through the bits of code that can be improved:

clc;clear all;

Please use the MATLAB Editor, and heeds its warnings. clear all is not useful here, it not only clears variables, it also clears functions from memory. These will need to be read in and parsed again, slowing down your program. clear suffices to remove variables from your workspace.

posDps = data(find(y == 1), :);
negDps = data(find(y == -1), :);

Here you can remove find. The MATLAB Editor also warns about this.

posDps = data(y == 1, :);
negDps = data(y == -1, :);

find is usually not necessary when indexing, only when you need to store an index to an array element.

space = size(data, 2);
n = length(data);

It is always better to specify directly what size you need. length(data) is defined as max(size(data)), I consider this dangerous. Instead:

[n,space] = size(data);

or

space = size(data, 2);
n = size(data, 1);

You should also know numel, which returns the number of elements in a matrix. This is much more efficient than prod(size(data)).

[ y1, y2 ] = meshgrid(y, y);
[ i, j ] = meshgrid(1:n, 1:n);

These are not necessary at all. Since MATLAB R2016b we have implicit singleton expansion, but before that you could use bsxfun to operate on two arrays of different sizes. So none of the following indexing operations are really necessary.

P(i, j) = y1(i, j) .* y2(i, j) .* (data(i,:) * data(j,:)');

can be rewritten as

P = (y .* y.') .* (data * data.');

Note that I use .' here, where you had '. ' is the complex conjugate transpose. Unless you want to compute a dot product, you should avoid it. .' transposes the matrix without changing any of its values. For real-valued data they're the same, but it's good to get used to using the right form, it will prevent difficult to find errors down the road.

Aeq = y';

Same here: y.'.

beq = [ 0 ];

Again, the MATLAB Editor warns here about the useless [] brackets. Heed the warnings!

There is actually no difference in MATLAB between a matrix and a scalar. A scalar is a matrix with a single value. The [...] operator is for concatenation. It concatenates two or more matrices into a single matrix. Every time you have a single matrix inside [], the MATLAB Editor will warn you about it.

Similarly, I often see things like [1:10]. 1:10 is a matrix, the [] don't do anything there.

w = (data' .* repmat(y', [space 1])) * alpha;

Same here: data.', y.'.

svIndices = find(alpha > threshold);
b = y(svIndices(1)) - data(svIndices(1),:) * w;

Here you can use the second argument to find to tell it to give you only one index:

svIndex = find(alpha > threshold,1);
b = y(svIndex) - data(svIndex,:) * w;

domain = (min(data(:,1)) - margin):(max(data(:,1)) + margin);

You don't need this many parentheses here, but they don't hurt either. I actually think they help readability. But the : has a much lower precedence than + and -. This should also work:

domain = min(data(:,1))-margin : max(data(:,1))+margin;