trybeetle

take it slow!!

SVMについて


Posted on July 28, 2018, 4:31 a.m.



SVM(support vector machine)について記載します。
CourseraのMachine Learningコースを元にしています。

SVM(support vector machine)では、境界線とtraining datasetとの距離を最大化させる(マージン最大化)する重みθを求めます。図としては、logistic関数を少し大雑把にした形です。

上図のとおり、目的関数を最小化するために、下記の2つを条件とします。
①赤色で囲った部分が最小値となるように算出する(赤色の部分が第二の目的関数となる)
②青色で囲った部分を小さく(y=1のときcost1≥1、y=0のときcost0≦-1)

※logistic関数のときと異なり、regulationのパラメータlambdaが、重みθ^2に対して掛かるのではなく、SVMでは、前半部分に掛けます(C=1/lambda)。
Cはlambdaを分母としているので、Cの値が大きくなれば、相対的にθ^2の影響力は少なくなります。(regulation効果は小さくなります。)
一方で、Cの値が小さくなれば、相対的にθ^2の影響力は大きくなります。(regulation効果は大きくなります。)


データの読み込み

今回のdatasetは、feature数が2のデータに対し、label=0 or 1となっています。


    load('ex6data1.mat');
SVMのトレーニング用関数

SVMのトレーニングを実行し、重みθを求める関数です。今回は、svmtrainを使用しています。(なお、θはwとして記載してあります。)返り値は、modelとして返します。


    function [model] = svmTrain(X, Y, C, kernelFunction, ...
                                tol, max_passes)
    %SVMTRAIN Trains an SVM classifier using a simplified version of the SMO
    %algorithm.
    %   [model] = SVMTRAIN(X, Y, C, kernelFunction, tol, max_passes) trains an
    %   SVM classifier and returns trained model. X is the matrix of training
    %   examples.  Each row is a training example, and the jth column holds the
    %   jth feature.  Y is a column matrix containing 1 for positive examples
    %   and 0 for negative examples.  C is the standard SVM regularization
    %   parameter.  tol is a tolerance value used for determining equality of
    %   floating point numbers. max_passes controls the number of iterations
    %   over the dataset (without changes to alpha) before the algorithm quits.
    %
    % Note: This is a simplified version of the SMO algorithm for training
    %       SVMs. In practice, if you want to train an SVM classifier, we
    %       recommend using an optimized package such as:
    %
    %           LIBSVM   (http://www.csie.ntu.edu.tw/~cjlin/libsvm/)
    %           SVMLight (http://svmlight.joachims.org/)
    %
    %

    if ~exist('tol', 'var') || isempty(tol)
        tol = 1e-3;
    end

    if ~exist('max_passes', 'var') || isempty(max_passes)
        max_passes = 5;
    end

    % Data parameters
    m = size(X, 1);
    n = size(X, 2);

    % Map 0 to -1
    Y(Y==0) = -1;

    % Variables
    alphas = zeros(m, 1);
    b = 0;
    E = zeros(m, 1);
    passes = 0;
    eta = 0;
    L = 0;
    H = 0;

    % Pre-compute the Kernel Matrix since our dataset is small
    % (in practice, optimized SVM packages that handle large datasets
    %  gracefully will _not_ do this)
    %
    % We have implemented optimized vectorized version of the Kernels here so
    % that the svm training will run faster.
    if strcmp(func2str(kernelFunction), 'linearKernel')
        % Vectorized computation for the Linear Kernel
        % This is equivalent to computing the kernel on every pair of examples
        K = X*X';
    elseif strfind(func2str(kernelFunction), 'gaussianKernel')
        % Vectorized RBF Kernel
        % This is equivalent to computing the kernel on every pair of examples
        X2 = sum(X.^2, 2);
        K = bsxfun(@plus, X2, bsxfun(@plus, X2', - 2 * (X * X')));
        K = kernelFunction(1, 0) .^ K;
    else
        % Pre-compute the Kernel Matrix
        % The following can be slow due to the lack of vectorization
        K = zeros(m);
        for i = 1:m
            for j = i:m
                 K(i,j) = kernelFunction(X(i,:)', X(j,:)');
                 K(j,i) = K(i,j); %the matrix is symmetric
            end
        end
    end

    % Train
    fprintf('\nTraining ...');
    dots = 12;
    while passes < max_passes,

        num_changed_alphas = 0;
        for i = 1:m,

            % Calculate Ei = f(x(i)) - y(i) using (2).
            % E(i) = b + sum (X(i, :) * (repmat(alphas.*Y,1,n).*X)') - Y(i);
            E(i) = b + sum (alphas.*Y.*K(:,i)) - Y(i);

            if ((Y(i)*E(i) < -tol && alphas(i) < C) || (Y(i)*E(i) > tol && alphas(i) > 0)),

                % In practice, there are many heuristics one can use to select
                % the i and j. In this simplified code, we select them randomly.
                j = ceil(m * rand());
                while j == i,  % Make sure i \neq j
                    j = ceil(m * rand());
                end

                % Calculate Ej = f(x(j)) - y(j) using (2).
                E(j) = b + sum (alphas.*Y.*K(:,j)) - Y(j);

                % Save old alphas
                alpha_i_old = alphas(i);
                alpha_j_old = alphas(j);

                % Compute L and H by (10) or (11).
                if (Y(i) == Y(j)),
                    L = max(0, alphas(j) + alphas(i) - C);
                    H = min(C, alphas(j) + alphas(i));
                else
                    L = max(0, alphas(j) - alphas(i));
                    H = min(C, C + alphas(j) - alphas(i));
                end

                if (L == H),
                    % continue to next i.
                    continue;
                end

                % Compute eta by (14).
                eta = 2 * K(i,j) - K(i,i) - K(j,j);
                if (eta >= 0),
                    % continue to next i.
                    continue;
                end

                % Compute and clip new value for alpha j using (12) and (15).
                alphas(j) = alphas(j) - (Y(j) * (E(i) - E(j))) / eta;

                % Clip
                alphas(j) = min (H, alphas(j));
                alphas(j) = max (L, alphas(j));

                % Check if change in alpha is significant
                if (abs(alphas(j) - alpha_j_old) < tol),
                    % continue to next i.
                    % replace anyway
                    alphas(j) = alpha_j_old;
                    continue;
                end

                % Determine value for alpha i using (16).
                alphas(i) = alphas(i) + Y(i)*Y(j)*(alpha_j_old - alphas(j));

                % Compute b1 and b2 using (17) and (18) respectively.
                b1 = b - E(i) ...
                     - Y(i) * (alphas(i) - alpha_i_old) *  K(i,j)' ...
                     - Y(j) * (alphas(j) - alpha_j_old) *  K(i,j)';
                b2 = b - E(j) ...
                     - Y(i) * (alphas(i) - alpha_i_old) *  K(i,j)' ...
                     - Y(j) * (alphas(j) - alpha_j_old) *  K(j,j)';

                % Compute b by (19).
                if (0 < alphas(i) && alphas(i) < C),
                    b = b1;
                elseif (0 < alphas(j) && alphas(j) < C),
                    b = b2;
                else
                    b = (b1+b2)/2;
                end

                num_changed_alphas = num_changed_alphas + 1;

            end

        end

        if (num_changed_alphas == 0),
            passes = passes + 1;
        else
            passes = 0;
        end

        fprintf('.');
        dots = dots + 1;
        if dots > 78
            dots = 0;
            fprintf('\n');
        end
        if exist('OCTAVE_VERSION')
            fflush(stdout);
        end
    end
    fprintf(' Done! \n\n');

    % Save the model
    idx = alphas > 0;
    model.X= X(idx,:);
    model.y= Y(idx);
    model.kernelFunction = kernelFunction;
    model.b= b;
    model.alphas= alphas(idx);
    model.w = ((alphas.*Y)'*X)';

    end
トレーニングモデルの実行

上記のsvmtrainを実行します。regulation値C=1としています。(Cが小さいほど、regulation効果が大きい)


    C = 1;
    model = svmTrain(X, y, C, @linearKernel, 1e-3, 20);

これで、重みθの値が、model.wとしてもとまりました。

続いて、上記で求めた重みを元に境界線をグラフ化します。


    function visualizeBoundaryLinear(X, y, model)
    w = model.w;
    b = model.b;
    xp = linspace(min(X(:,1)), max(X(:,1)), 100);
    yp = - (w(1)*xp + b)/w(2);
    plotData(X, y);
    hold on;
    plot(xp, yp, '-b');
    hold off
    end

グラフ化する関数を実行します。


    visualizeBoundaryLinear(X, y, model);

Gaussian kernel

SVMにはカーネルを設定使用し、polynomialな計算式でnon linearな境界線を求めることもできます。

Gaussian kernelでは、feature(x) と landmark(l) の間の距離からsimilarityを求め(f)とします。求めたfをxの代わりに使用することで、polynomialな式にも対応することができます。

f(i) = sim(x,l) = exp(-(|| x-l ||^2) / (2*sigma^2))


Datasetの読み込み

今回のdatasetは、feature数が2のデータに対し、label=0 or 1となっています。


    load('ex6data2.mat');

SVMのトレーニング用関数

Gaussian kernelの関数を用意します。


    function sim = gaussianKernel(x1, x2, sigma)
    x1 = x1(:); x2 = x2(:);
    sim = 0;
    sim = exp(-(norm(x1 - x2).^2) / (2 * (sigma^2)));
    end

linear regressionと同様に、svmtrainを実行します。regulation値C=1としています。
(Cが小さいほど、regulation効果が大きい)
(sigmaが大きいほど、regulationの効果が大きい)


    C = 1;
    sigma = 0.1;
    model= svmTrain(X, y, C, @(x1, x2) gaussianKernel(x1, x2, sigma));

これで、重みθの値が、model.wとしてもとまりました。

続いて、予測を行うsvmPredictの関数です。


    function pred = svmPredict(model, X)
    %SVMPREDICT returns a vector of predictions using a trained SVM model
    %(svmTrain). 
    %   pred = SVMPREDICT(model, X) returns a vector of predictions using a 
    %   trained SVM model (svmTrain). X is a mxn matrix where there each 
    %   example is a row. model is a svm model returned from svmTrain.
    %   predictions pred is a m x 1 column of predictions of {0, 1} values.
    %

    % Check if we are getting a column vector, if so, then assume that we only
    % need to do prediction for a single example
    if (size(X, 2) == 1)
        % Examples should be in rows
        X = X';
    end

    % Dataset 
    m = size(X, 1);
    p = zeros(m, 1);
    pred = zeros(m, 1);

    if strcmp(func2str(model.kernelFunction), 'linearKernel')
        % We can use the weights and bias directly if working with the 
        % linear kernel
        p = X * model.w + model.b;
    elseif strfind(func2str(model.kernelFunction), 'gaussianKernel')
        % Vectorized RBF Kernel
        % This is equivalent to computing the kernel on every pair of examples
        X1 = sum(X.^2, 2);
        X2 = sum(model.X.^2, 2)';
        K = bsxfun(@plus, X1, bsxfun(@plus, X2, - 2 * X * model.X'));
        K = model.kernelFunction(1, 0) .^ K;
        K = bsxfun(@times, model.y', K);
        K = bsxfun(@times, model.alphas', K);
        p = sum(K, 2);
    else
        % Other Non-linear kernel
        for i = 1:m
            prediction = 0;
            for j = 1:size(model.X, 1)
                prediction = prediction + ...
                    model.alphas(j) * model.y(j) * ...
                    model.kernelFunction(X(i,:)', model.X(j,:)');
            end
            p(i) = prediction + model.b;
        end
    end

    % Convert predictions into 0 / 1
    pred(p >= 0) =  1;
    pred(p <  0) =  0;

    end

上記で求めた予測結果を元に境界線をグラフ化します。


    function visualizeBoundary(X, y, model, varargin)
    plotData(X, y)
    x1plot = linspace(min(X(:,1)), max(X(:,1)), 100)';
    x2plot = linspace(min(X(:,2)), max(X(:,2)), 100)';
    [X1, X2] = meshgrid(x1plot, x2plot);
    vals = zeros(size(X1));
    for i = 1:size(X1, 2)
       this_X = [X1(:, i), X2(:, i)];
       vals(:, i) = svmPredict(model, this_X);
    endfor
    hold on
    contour(X1, X2, vals, [0.5 0.5], 'b');
    hold off;

    end

グラフ化する関数を実行します。


    visualizeBoundary(X, y, model);

Cとδの最適値について

Cとδの最適値を求めるために、下記のように、候補となる値を列挙しておき、for loopで組み合わせを全て試し、コスト関数を最小にできた値が、Cとδの最適値です。


    paramlist = [0.01,0.03,0.1,0.3,1,3,10,30];

以上になります。


Category:ML
Tag: ML Octave
July 28, 2018, 4:31 a.m.

Comments