svcm_train.m

svm_increasement Matlab Program

                warning off % suppress div. by 0
                zlime = z0-g(inde)./gammae;
                warning on
                zlime(void) = polc*Inf;          % ... which don't enter the constraints
                [zmine, ie] = min(zlime*polc,[],1);
                imin = find(zlime==zmine);
                if length(imin)>1
                    [gmax, imax] = max(abs(gammae(imin)),[],1);
                    ie = imin(imax);
                end
                zlime = zmine*polc;              % pick tightest constraint
            end
        end
    end

    % ordinary vector constraints: g(indo)>=0 (only for those that already are)
    zlimo = Inf*polc;                            % by default, immaterial
    if lo>0
        gammao = gamma(indo);
        void = (indo==indc)|(g(indo)<0)|(gammao*polc>0)|(gammao==0);
                                                 % void c, g negative, g moving up, or zero gamma...
        if online
            void = void|~processed(indo);        % ... or, if online, points not seen previously,...
        end
        if any(any(~void))
            warning off % suppress div. by 0
            zlimo = z0-g(indo)./gammao;
            warning on 
            zlimo(void) = polc*Inf;              % ... which don't enter the constraints
            [zmino, io] = min(zlimo*polc,[],1);
            imin = find(zlimo==zmino);
            if length(imin)>1
                [gmax, imax] = max(abs(gammao(imin)),[],1);
                io = imin(imax);
            end
            zlimo = zmino*polc;                  % pick tightest constraint
        end
    end

    % find constraint-satisfying z
    [z,flag] = min([zlim;zlimc;zlims;zlime;zlimo]*polc);
    z = z*polc;
    if (z-z0)*polc<0
        fprintf('\nsvcm_train error: z-z0 of wrong polarity (Q not positive definite)\n\n')
        return
    end

    if verbose & ~leaveoneout & abs(z-z0)<eps
        fprintf('%g*', flag)            % procrastinating iteration!  no progress made
    end

    % update a, b, g and W from z-z0
    if ls>0                             % z = a(indc)
        a(indc) = z;
        b = b+(z-z0)*beta(1);
        a(inds) = a(inds)+(z-z0)*beta(2:ls+1);
        W = W+(z-z0)*(g(indc)'+0.5*(z-z0)*gammac);                % energy
    else                                % z = y(indc)*b
        b = y(indc)*z;
    end
    g = g+(z-z0)*gamma;                 % update g
    iter = iter+1;
    if visualize & ~leaveoneout
        atraj(1:L,iter)=a;              % record trajectory of a(:) over time
        gtraj(1:L,iter)=g;
        ctraj(iter)=indc;
    end 

    % bookkeeping: move elements across indc, inds, inde and indo, and update R and Qs
    converged = (flag<3);               % done with indc; no other changes in inds/inde
    incl_inds = 0;
    if flag==1                          % a(indc) reaches the limits 0 or C, stop moving
        if upc                                  % a(indc)=C, add to inde
            inde = [inde; indc];
            le = le+1;
            if keepe
                Qe = [Qe; Qc];
            end
            a(indc) = C;                        % should be OK, just to avoid round-off
        else % ~upc                                % a(indc)=0, indc stays in (or moves to) indo
            a(indc) = 0;                        % should be OK, just to avoid round-off
        end
    elseif flag==2                      % add indc to support vectors ...
        incl_inds = 1;
        indb = indc;                            % ... store in buffer indb for now
        Qb = Qc;
    elseif flag==3                      % one of support vectors becomes error or other vector
        indso = inds(is);                       % outgoing inds
        Qso = Qs(is+1,:);                       % could be reused later
        free_indc = (indc==indso);              % leave-indc-out: indc is part of inds
        if beta(is+1)*polc<0 | free_indc        % a(indso)=0 or indso=indc, move to indo
            if ~free_indc
                if keepr
                    indr = [indr; indso];       % also recycle into indr for later use
                    Qr = [Qr; Qs(is+1,:)];
                    lr = lr+1;
                end
                if a(indso)>C/2
                    fprintf('svcm_train error: a(indso)=%g; 0 anticipated\n', a(indso));
                end
                a(indso) = 0;                   % should be OK, just to avoid round-off
            end
            g(indso) = 0;                       % same
            indo = [indo; indso];
            lo = lo+1;
        else % beta(is+1)*polc>0 & ~free_indc   % a(indso)=C, move to inde
            if a(indso)<C/2
                fprintf('svcm_train error: a(indso)=%g; C anticipated\n', a(indso));
            end
            if keepe
                Qe = [Qe; Qs(is+1,:)];          % save to memory cache
            end
            a(indso) = C;                       % should be OK, just to avoid round-off
            g(indso) = 0;                       % same
            inde = [inde; indso];
            le = le+1;
        end
        inds = inds([1:is-1,is+1:ls]);          % remove from inds
        stripped = [1:is,is+2:ls+1];            % also from Qs and R ...
        Qs = Qs(stripped,:);
        ls = ls-1;
        if ls > 0
            if R(is+1,is+1)==0
                fprintf('\nsvcm_train error: divide by zero in R contraction\n')
                R(is+1,is+1)=1e-8;
            end
            R = R(stripped,stripped)-R(stripped,is+1)*R(is+1,stripped)/R(is+1,is+1);
        else % no support vectors left
            R = Inf;
        end
    elseif flag==4                       % one of error vectors becomes support/other vector
        indeo = inde(ie);                        % outgoing inde
        if indc==indeo                           % leave-indc-out
            indo = [indo; indeo];                % add inde(ie) to other vectors
            lo = lo+1;
        else
            incl_inds = 1;                       % add inde(ie) to support vectors ...
            indb = indeo;                        % ... store in buffer indb for now
            if keepe
                Qb = Qe(ie,:);                   % recover from Qe cache
            elseif indb==indso
                Qb = Qso;                        % recover from previous outgoing support vector
            else                                 % not in memory either way--- recompute
                Qb = (y(indb)*y').*kernel(x(indb,:),x);
                Qb(indb) = Qb(indb)+eps2;
                kernelcount = kernelcount+1;
            end
        end
        inde = inde([1:ie-1,ie+1:le]);           % remove from inde
        if keepe
            Qe = Qe([1:ie-1,ie+1:le],:);         % remove from Qe
        end
        le = le-1;
    elseif flag==5                       % one of other vectors becomes support vector
        indoo = indo(io);                        % outgoing indo
        incl_inds = 1;                           % add indo(io) to support vectors ...
        indb = indoo;                            % ... store in buffer indb for now
        if keepr & lr>0 & any(find(indr==indb))  % check for match among recycled vectors
            ir = find(indr==indb);                       % found, reuse
            Qb = Qr(ir,:);
            indr = indr([1:ir-1,ir+1:lr]);               % ... remove from indr
            Qr = Qr([1:ir-1,ir+1:lr],:);                 % ... and Qr
            lr = lr-1;
        elseif indb==indso
            Qb = Qso;                            % recover from previous outgoing support vector
        else                                     % not in memory either way--- recompute
            Qb = (y(indb)*y').*kernel(x(indb,:),x);
            Qb(indb) = Qb(indb)+eps2;
            kernelcount = kernelcount+1;
        end
        indo = indo([1:io-1,io+1:lo]);           % remove from indo
        lo = lo-1;
    end

    if incl_inds                        % move buffer indb into support vectors inds
        inds = [inds; indb];                     % move to inds ...
        ls = ls+1;
        Qs = [Qs; Qb];                           % and also Qs and R ...
        if ls==1                                 % compute R directly
            R = [-Qb(indb), y(indb); y(indb), 0];
        else                                     % compute R recursively
            if flag==2                           % from indc; use beta and gamma
                    pivot = gamma(indb);
            else % flag==4                       % from inde; compute beta and pivot
                beta=-R*Qs(1:ls,indb);
                pivot = [beta',1]*Qs(:,indb);
            end
            if pivot<eps2        % should be eps2 when kernel is singular (e.g., linear)
                fprintf('\nsvcm_train error: pivot = %g < %g in R expansion\n\n', pivot, eps2)
                pivot = eps2;
            end
            R = [R,zeros(ls,1);zeros(1,ls+1)]+[beta;1]*[beta',1]/pivot;
        end
    end

    % minor correction in R to avoid numerical instability in recursion when data is near-singular
    Qss = [[0;y(inds)],Qs(:,inds)];
    R = R+R'-R*Qss*R';

    % indc index adjustments (including leave-one-out)
    if converged & (upc|flag==2)                % indc is now part of inds or inde
        i = find(indo==indc);
        indo = indo([1:i-1,i+1:lo]);  % remove indc from indo
        lo = lo-1;
    elseif keepr
        indr = [indr; indc];                    % recycle again into indr for later use
        Qr = [Qr; Qc];
        lr = lr+1;
    end
    if leaveoneout
        indoc = indo(find(indo~=indc));         % indo other than indc
        free = a(indoc)>0|g(indoc)<0;           % candidate support/error vectors in indoc
        if visualize & ~any(free)
            gctraj = [gctraj,[a(indc);g(indc)]];
        end
        converged = converged & ~any(free);
        if converged                            % leave-indc-out reached a(indc)=0, all others settled
            indl = indl(find(indl~=indc));      % remove indc from indl
            if visualize
                figure(1)
                hold on
                if any(gctraj)
                    plot(gctraj(1,:),gctraj(2,:),'ok',gctraj(1,:),gctraj(2,:),'-k')
                end
                gctraj = [];                     % cleanup for next curve
            end
            if g(indc)<-1                       % if leave-indc-out generates an error ...
                indw = [indw; indc];            % ... store its index
                lw = lw+1;                      % ... and increment leave-one-out error count
            end
        end
    end

    % debugging mode: check for consistency; a, R, g and W
    if debug
        f = find(isnan(a));
        if any(f)
           fprintf('svcm_train error: a(%g) = %g\n', [f, a(f)]')
        end