View Single Post
  #49  
Old 05-23-2012, 01:02 AM
sakumar sakumar is offline
Member
 
Join Date: Apr 2012
Posts: 40
Default Re: What Quadratic Programming Package?

Quote:
Originally Posted by kurts View Post
sakumar,
I am an Octave newbie - can you show me the commands you used to generate those plots?
Sure. Here it is. It is hacked up a bit because I cut and paste to get Figure 1 and Figure 2. It just runs one iteration and generates two plots.

Code:
iterations = 1;
n = 100;

test_n = 10000; % number of points used to determine Eout

svm_threshold = 10^(-5); % if alpha value is greater than this, count it as a support vector


function w = GetWFrom2Points(x, y)
    w(3) = - 1.0;
    w(2) = (x(2) - y(2)) / (x(1) - y(1));
    w(1) = x(2) - w(2) * x(1);
end

X = zeros(3,n);
testX = zeros(3, test_n);


min_val_vec = zeros(n, 1);
min_val_vec = svm_threshold;

svm_winner_count = 0;
sv_count = 0;

for m = 1:iterations

    svm_fraction_incorrect = 0.0;
    perceptron_fraction_incorrect = 0.0;

    do % we need to check to make sure the n points are not all on one side of the line
       % asked to discard datapoint if that happens and start again
        a = unifrnd(-1.0, 1.0, 1, 2);
        b = unifrnd(-1.0, 1.0, 1, 2);

        fw = GetWFrom2Points(a, b); % This is the target function
        fw = fw';

        % now generate n points -- these are the in-sample points
        X(2:3,1:n) = unifrnd(-1.0, 1.0, 2, n);
        X(1,1:n) = 1.0;

        % get true value for in-sample points
        ty = sign(fw'*X);
    until abs(sum(ty)) != n % if all the y's added up equal to n, it means all the points are on one side of the line


    % now we run the perceptron
    % initialize hypothesis ("all-zero vector")
    w = zeros(3,1);

    while (1)
        y = sign(w'*X);
        comp = (y == ty);

        if (comp) % if all comp values are 1, then we have converged
            break;
        end

        % get random mismatched index
        ilist = find(!comp);
        m_indx = ilist(unidrnd(sum(!comp)));

        w = w + ty(m_indx)*X(:, m_indx);
    end

    % now get probability that f(x) is not equal to g(x)
    testX(2:3,1:test_n) = unifrnd(-1.0, 1.0, 2, test_n);
    testX(1,1:test_n) = 1.0;

    test_ty = sign(fw'*testX);
    test_y = sign(w'*testX);

    comp = (test_y == test_ty);

    perceptron_fraction_incorrect += sum(!comp)/test_n;


    % Now on to the quadratic programming. Am going to try the "qp" function in Octave

    % initialize alpha0 vector to zeros (??)
    alpha0 = ones(n, 1);

    % the "q" matrix
    q = ones(n,1)*-1;

    % the A matrix (actually vector)
    A = ty;

    % the b value
    b = 0;

    % the lb vector (all zeros)
    lb = zeros(n, 1);

    X_svm = X(2:3, :);
    % build the H (called Q in lecture) matrix
    H = (ty'*ty).*(X_svm'*X_svm);

    tempH = H + eye(100).*10^-15;
    [alpha0, obj, info, lambda] = qp([], tempH, q, A, b, lb, []);

    options.MaxIter = 1000;
    [alpha, obj, info, lambda] = qp(alpha0, H, q, A, b, lb, [], options);

    'INFO.INFO'
    info.info

    % count number of support vectors
    comp = ge(alpha, min_val_vec);
    sv_index = find(comp);
    num_sv = sum(comp);

    sv_count += num_sv;


    svm_w = zeros(3, 1);

    svm_w(2) = sum(X_svm(1, :).*(alpha'.*ty));
    svm_w(3) = sum(X_svm(2, :).*(alpha'.*ty));

    % get index of the maximum alpha
    [maxval, maxindx] = max(alpha);

    svm_w(1) = 1/ty(maxindx) - (svm_w(2:3,1))'*X_svm(:, maxindx);


    negmark = lt(ty, min_val_vec);
    negindex = find(negmark);

    posmark = gt(ty, min_val_vec);
    posindex = find(posmark);

    figure(1);
    plot(X(2,negindex), X(3,negindex), "@11", X(2, posindex), X(3, posindex), "@22", X(2, sv_index), X(3, sv_index), "o");

    printf("NUMBER OF SUPPORT VECTORS %i\n", length(sv_index));

    options.MaxIter = 1000;
    [alpha, obj, info, lambda] = qp([], H, q, A, b, lb, [], options);

    'INFO.INFO'
    info.info

    % count number of support vectors
    comp = ge(alpha, min_val_vec);
    sv_index = find(comp);
    num_sv = sum(comp);

    sv_count += num_sv;


    svm_w = zeros(3, 1);

    svm_w(2) = sum(X_svm(1, :).*(alpha'.*ty));
    svm_w(3) = sum(X_svm(2, :).*(alpha'.*ty));

    % get index of the minimum alpha
    [maxval, maxindx] = max(alpha);

    svm_w(1) = 1/ty(maxindx) - (svm_w(2:3,1))'*X_svm(:, maxindx);


    negmark = lt(ty, min_val_vec);
    negindex = find(negmark);

    posmark = gt(ty, min_val_vec);
    posindex = find(posmark);

    figure(2);
    plot(X(2,negindex), X(3,negindex), "@11", X(2, posindex), X(3, posindex), "@22", X(2, sv_index), X(3, sv_index), "o");

    printf("NUMBER OF SUPPORT VECTORS %i\n", length(sv_index));


end
Reply With Quote