ScalarEpsSVR_SMO.m

Go to the documentation of this file.

namespace general{
namespace regression{


/* (Autoinserted by mtoc++)
 * This source code has been filtered by the mtoc++ executable,
 * which generates code that can be processed by the doxygen documentation tool.
 *
 * On the other hand, it can neither be interpreted by MATLAB, nor can it be compiled with a C++ compiler.
 * Except for the comments, the function bodies of your M-file functions are untouched.
 * Consequently, the FILTER_SOURCE_FILES doxygen switch (default in our Doxyfile.template) will produce
 * attached source files that are highly readable by humans.
 *
 * Additionally, links in the doxygen generated documentation to the source code of functions and class members refer to
 * the correct locations in the source code browser.
 * However, the line numbers most likely do not correspond to the line numbers in the original MATLAB source files.
 */

class ScalarEpsSVR_SMO
  :public general.regression.BaseScalarSVR {
    public: /* ( setObservable ) */

        Eps = .1;
        StopEps = 1e-4;
        Vis = 0;
        Version = 2;
        NNk = 10;
        MaxCount = 60000;
    public:

        integer LastIterations = "[]";
    public:

        function [ai , sf ] = regress(fxi,initialai) {
            if any(abs(fxi) > 1)
                error(" This SVR implementation requires fxi values in [-1,1] ");
            end
/*              initialai = reshape(initialai,1,[]); */
            if this.Version == 1
                /*  Call 1D-Optimizer */
                [ai, sf] = this.regress1D(fxi, initialai);
            else
                /*  Call 2D-Optimizer */
                [ai, sf] = this.regress2D(fxi, initialai);
            end
        }


        function copy = clone() {
            copy = general.regression.ScalarEpsSVR_SMO;
            copy = clone@general.regression.BaseScalarSVR(this, copy);
            copy.Eps= this.Eps;
            copy.StopEps= this.StopEps;
            copy.Vis= this.Vis;
            copy.Version= this.Version;
            copy.NNk= this.NNk;
            copy.MaxCount= this.MaxCount;
            copy.LastIterations= this.LastIterations;
        }

    
    private:


        function [ai , sf ] = regress1D(fxi,initialai) {
            
            if ~isempty(initialai)
                /*  Init - warm start */
                [a, dW, T] = I2_WarmStart(this, fxi, initialai);
            else
                /*  Init - cold start */
                [a, dW, T] = I0_ColdStart(this, fxi);

                /*  Init - kernel rule
 *                 [a, dW, T] = I1_ColdStartKernelRule(this, fxi); */
            end

            /* % Preps */
            cnt = 1;

            n = length(fxi);

            sgn = ones(1,2*n);
            sgn(n+1:end) = -1;

            S = n*this.C;

            if this.Vis > 1
                h = figure(1);
            end

            stop = this.StopEps/(2*this.Lambda);
            while S > stop && cnt < this.MaxCount+1

                if this.Vis > 1
                    afxi = (a(1:n)-a(n+1:end))*this.K;
                    figure(h);
                    subplot(1,2,1);
                    plot(1:n,fxi," r ",1:n,[fxi-this.Eps; fxi+this.Eps]," r-- ",1:n,afxi," b ");
                    legend(" f(x_i) training values "," +\epsilon "," -\epsilon "," approximation values ");
                    axis tight;
                end

                /* % Select working set */
                [idx, M] = this.WSS0_1D_GetMaxGainIndex(a, dW);

                /* [idx, M] = this.WSS1_1D_GetMaxGradientIndex(a, dW); */

                r = max(-a(idx),min(this.C-a(idx),dW(idx)));

                /* % Misc vars
                 * For index: Subtract n if index is for an \alpha^- (linear indexing) */
                idx1 = idx - n*(idx > n);

                if this.Vis > 1
                    subplot(1,2,1);
                    hold on;
                    plot(idx1,afxi(idx1)," . "," MarkerSize ",5);
                    hold off;
                end

                /* % update ai at max gain index */
                a(idx) = a(idx) + r;                  

                /* % Update T with old dW */
                T = T + r*(r - 2*dW(idx) + sgn(idx)*fxi(idx1) - this.Eps);

                /* % Update dW - only on side with change */
                Ki = this.K(idx1, :);
                dW = dW + sign(idx-n-.5) * r * [Ki -Ki];

                /* % Update E */
                E = this.C * sum(max(0,min(2-this.Eps,dW)));

                S = abs(T + E);
                cnt = cnt+1;
            end
            if this.Vis > 0
                disp(sprintf(" Finished after %d/%d iterations.\n ",cnt,this.MaxCount));/* #ok */

            end

            /*  \alpha = \alpha^+ - \alpha^- */
            ai = (a(1:n)-a(n+1:end))^t;

            sf = StopFlag.TOL_OK;
            if cnt == this.MaxCount
                sf = StopFlag.MAX_ITER;
            end
            this.LastIterations= cnt;
        }


        function [ai , sf ] = regress2D(fxi,initialai) {
            if ~isempty(initialai)
                /*  Init - warm start */
                [a, dW, T] = I2_WarmStart(this, fxi, initialai);
            else
                /*  Init - cold start */
                [a, dW, T] = I0_ColdStart(this, fxi);

                /*  Init - kernel rule
 *                 [a, dW, T] = I1_ColdStartKernelRule(this, fxi); */
            end

            if this.Vis > 1
                h = figure(1);
            end

            /* % Preps */
            n = length(fxi);
            sgn = ones(1,2*n);
            sgn(n+1:end) = -1;

            S = n*this.C;

            /*  Initialize i to -1 if MaxGainPlusLast is used. */
            i = -1;

            cnt = 1;
            stop = this.StopEps/(2*this.Lambda);
            while abs(S) > stop && cnt < this.MaxCount+1
                /* % Max gain computation - O(n^2) strategy */
                ij = [];/* #ok<*AGROW> */


/*                  ij = [ij; this.WSS0_BestGainIndices(a, dW, sgn); */

/*                  ij = [ij; this.WSS1024_RandomAdmissibleIndices(a); */

/*                  ij = [ij; this.WSS1_1DMaxGainPlusLast(a, dW, i)]; */

/*                  ij = [ij; this.WSS2_TwoSetMax1DGain(a, dW)]; */

/*                  ij = [ij; this.WSS4_1DMaxGainPluskNN(a, dW)]; */

                ij = [ij; this.WSS7_Combi(a, dW, i)];

                [r, s, idx] = this.getMax2DGainUpdatesForIndices(a, dW, ij);

                /*  Not usable with current scheme
                 *[i, j] = WSS128_ApproxBestGainIndices(this, a, dW, sgn); */

                i = ij(idx,1);
                j = ij(idx,2);
                i1 = i - (i > n)*n;
                j1 = j - (j > n)*n;

                if this.Vis > 1
                    afxi = (a(1:n)-a(n+1:end))*this.K;
                    figure(h);
                    subplot(1,2,1);
                    plot(1:n,fxi," r ",1:n,[fxi-this.Eps; fxi+this.Eps]" , "r--" ,1:n,afxi, "b^t);
                    hold on;
                    if (i1 == i)
                        plot(i,afxi(i)," o "," MarkerEdgeColor "," k ",...
                            " MarkerFaceColor ",[.1 .8 .1],...
                            " MarkerSize ",6);
                    else
                        plot(i1,afxi(i1)," o "," MarkerEdgeColor "," k ",...
                            " MarkerFaceColor ",[.8 .1 .1],...
                            " MarkerSize ",6);
                    end
                    if (j1 == j)
                        plot(j,afxi(j)," o "," MarkerEdgeColor "," k ",...
                            " MarkerFaceColor ",[.1 .8 .1],...
                            " MarkerSize ",6);
                    else
                        plot(j1,afxi(j1)," o "," MarkerEdgeColor "," k ",...
                            " MarkerFaceColor ",[.8 .1 .1],...
                            " MarkerSize ",6);
                    end

                    hold off;
                    axis tight;
                    if this.Vis > 2
                        fprintf(" alpha_{%d} change: %e, alpha_{%d} change: %e\n ",i,r,j,s);
                    end
                    pause;
                end

                a(i) = a(i) + r;
                a(j) = a(j) + s;

                /*  If difference is smaller than n, they are ++ or -- updates, meaning addition.
                 * otherwise, is |i-j| \geq n, we have +- or -+ cases, meaning subtraction.
                 * subtraction of .5 avoids having sign(0) = 0. */
                pm = sign(i - n - .5)*sign(j - n - .5);
                /* % Update T using old dW */
                T = T + r*(r - 2*dW(i) - this.Eps + sgn(i)*fxi(i1)) ...
                      + s*(s - 2*dW(j) - this.Eps + sgn(j)*fxi(j1)) ...
                      + 2*pm*r*s*this.K(i1,j1);

                /* % Update dW - only on side with change */
                Ki = this.K(i1, :);
                Kj = this.K(j1, :);
                dW = dW + sign(i-n-.5) * r * [Ki -Ki] ...
                        + sign(j-n-.5) * s * [Kj -Kj];

                /* % Get new E term */
                E = this.C * sum(max(0,min(2-this.Eps,dW)));

                S = T + E;
                cnt = cnt+1;
            end
            cnt = cnt-1;
            if this.Vis > 0 || KerMor.App.Verbose > 0
                fprintf(" ScalarEpsSVR_SMO: Finished after %d/%d iterations.\n ",cnt,this.MaxCount);
            end

            /*  \alpha = \alpha^+ - \alpha^- */
            ai = (a(1:n)-a(n+1:end))^t;  

            sf = StopFlag.TOL_OK;
            if cnt == this.MaxCount
                sf = StopFlag.MAX_ITER;
            end
            this.LastIterations= cnt;
        }


        function [i , M ] = WSS0_1D_GetMaxGainIndex(a,dW) {
            n = numel(a)/2;
            ind = 1:n;
            /*  Check which alpha values might be changed (only the ones with their partner
             * variable equal to zero) */
            ch = [ind(a(n+1:end) == 0) ind(a(1:n) == 0)+n];

            /*  Find optima for changeable alphas */
            r = max(-a(ch),min(this.C-a(ch),dW(ch)));

            /*  Get gain */
            g = r .* (dW(ch) - .5*r);

            [M, idxch] = max(g);
            i = ch(idxch);

            if this.Vis > 1
                subplot(1,2,2);
                r2 = zeros(size(a));
                r2(ch) = r;
                plot(1:2*n,r2," b ");
                hold on;
                g2 = zeros(size(a));
                g2(ch) = g;
                plot(1:2*n,g2," g ");
                plot(i,g(idxch)," r. "," MarkerSize ",6);
                plot([n n+eps],[0, max([r g])]," black ");
                hold off;
                title(sprintf(" Best gain index: %d, gain: %e, \\alpha_{%d} change: %e ",i, M, i, r(idxch)));
                legend(" \alpha difference "," gain ");
                axis tight;
            end
        }


        function [i , M ] = WSS1_1D_GetMaxGradientIndex(a,dW) {
             [M, i] = max(max([dW .* (a < this.C); -dW .* (a > 0)],[],1));
        }
        function ij = WSS0_AllIndices(a) {

            n = numel(a)/2;
            ind = 1:n;
            apch = ind(a(n+1:end) == 0);
            amch = ind(a(1:n) == 0);
            ch = [apch amch+n];
            /*  Create unique combinations of indices */
            A = repmat(ch,length(ch),1);
            I = triu(A,1);
            J = triu(A^t,1);
            ij = [I(I~=0) J(J~=0)];
            /*  Remove indices that would look at moving the same alpha up and down at the same time */
            ij(ij(:,1)-ij(:,2)-n == 0,:) = [];
        }
        function ij = WSS1_1DMaxGainPlusLast(a,dW,i) {
            in = this.WSS0_1D_GetMaxGainIndex(a, dW);
            if i ~= -1
                j = i;
                i = in;
            else
                [dWs, sortidx] = sort(abs(dW));
                pick = 1;
                while sortidx(pick) == in
                    pick = pick + 1;
                end
                i = in;
                j = sortidx(pick);
            end
            ij = [i j];
        }
        function ij = WSS2_TwoSetMax1DGain(a,dW) {
            n = numel(a)/2;
            len = round(n/2);
            set1 = 1:len;
            set2 = set1(end)+1:n;
            i = this.WSS0_1D_GetMaxGainIndex(a([set1 set1+n]), dW([set1 set1+n]));
            if i > len
                i = i+(n-len); /*  Move to \am index (add number of set2 elements), i.e. i = i+numel(set2) */

            end
            j = this.WSS0_1D_GetMaxGainIndex(a([set2 set2+n]), dW([set2 set2+n]));
            if j > numel(set2)
                j = j+len; /*  If larger than set2 size add number of set1 elements */

            end
            /*  Add +len as this is the offset for the second set in each case */
            ij = [i j+len];
        }


        function ij = WSS4_1DMaxGainPluskNN(a,dW) {
            i = this.WSS0_1D_GetMaxGainIndex(a, dW);

            n = numel(a)/2;
            ki = i;
            if ki > n
                ki = ki - n;
            end
            [dist, idx] = sort(sqrt(1-this.K(ki,:)));
            /*  Extract indices of alpha^+ that are "close" wrt to the kernel metric \sqrt{1-\K(x_i,x_j)}
             * Start at 2 as i itself is not an option for second index */
            nidx = idx(2: min(this.NNk+1,numel(idx)));
            /*  Of course also consider the \alpha^- */
            ij = [i*ones(size(nidx,2)*2,1) [nidx nidx+n]^t];
        }


        function ij = WSS7_Combi(a,dW,i) {
            ij = this.WSS2_TwoSetMax1DGain(a, dW);
            ij = [ij; this.WSS4_1DMaxGainPluskNN(a, dW)];

            /*  Implement Max1DGainPlusLast directly to save a call to the 1D method. */
            in = ij(end,1);
            if i ~= -1
                j = i;
                i = in;
            else
                [dWs, sortidx] = sort(abs(dW));
                pick = 1;
                while sortidx(pick) == in
                    pick = pick + 1;
                end
                i = in;
                j = sortidx(pick);
            end
            ij = [ij; i j];
        }


        function ij = WSS128_ApproxBestGainIndices(a,dW) {

            /*  Only consider changeable \alpha^{+,-}, i.e. the ones whos partner \alpha equals zero. */
            n = numel(a)/2;
            ind = 1:n;
            apch = ind(a(n+1:end) == 0);
            amch = ind(a(1:n) == 0);
            ch = [apch amch+n];
            /*  Create unique combinations of indices */
            A = repmat(ch,length(ch),1);
            I = triu(A,1);
            J = triu(A^t,1);
            ij = [I(I~=0) J(J~=0)];
            /*  Remove indices that would look at moving the same alpha up and down at the same time */
            ij(ij(:,1)-ij(:,2)-n == 0,:) = [];

            kijidx = ij;
            kijidxup = kijidx > n;
            kijidx(kijidxup) = kijidx(kijidxup)-n;
            kijidx = sub2ind(size(this.K),kijidx(:,1),kijidx(:,2));
            Kij = this.K(kijidx);

            sgn = 1-mod(sum(kijidxup,2),2)*2;
            dwi = dW(ij);
            aij = a(ij);

            div = 1 ./ (1-Kij.^2);

            /*  Compute update r */
            rs = [(dwi(:,1) - sgn .* Kij .* dwi(:,2)) .* div ...
                  (dwi(:,2) - sgn .* Kij .* dwi(:,1)) .* div];
            ub = rs > this.C - aij;
            lb = rs < - aij;
            rs(ub) = this.C-aij(ub);
            rs(lb) = -aij(lb);            

            /* % Compute gain for r_i,s_j combinations */
            g = sum(rs .* (dwi - .5*rs),2) - sgn.*rs(:,1).*rs(:,2).*Kij;

            /*  Extract i,j indices */
            [maxg,idx] = max(g);

            if this.Vis > 2
                fprintf(" Updates found inside WSS128: r=%f, s=%f\n ",rs(idx,1),rs(idx,2));
            end
            ij = ij(idx,:); 
            if this.Vis > 1    
                subplot(1,2,2);
                plot(1:length(g),max(0,g));
                hold on;
                plot(idx,g(idx)," black. "," MarkerSize ",8);
                hold off;
                title(sprintf(" Best gain index: (%d,%d), gain: %e ",ij(1),ij(2), maxg));
                axis tight;
            end
        }
        function ij = WSS1024_RandomAdmissibleIndices(a) {
            n = numel(a)/2;
            ind = 1:n;
            ch = [ind(a(n+1:end) == 0) ind(a(1:n) == 0)+n];
            i = ch(randi(length(ch)));
            j = ch(randi(length(ch)));
            while j==i || abs(j-i) == n
                j = ch(randi(length(ch)));
            end
            ij = [i j];
        }


        function [a , dW , T ] = I0_ColdStart(fxi) {
            T = 0;
            a = zeros(1,2*length(fxi)); /* 1..n = a^+, n+1..2n=a^- */

            dW = [fxi -fxi] - this.Eps;
        }


        function [ai , dW , T ] = I1_ColdStartKernelRule(fxi) {
            n = length(fxi);
            /* ai(1:n) = 1; ai(n+1:2*n) = 0; */
            ai = max(0,min(this.C,[fxi.*(fxi>=0) -fxi.*(fxi<0)]));
            ai = ai * 10 / n;
            a = ai(1:n)-ai(n+1:end);
            hlp = -a*this.K + fxi;
            dW = [hlp -hlp] - this.Eps;
            /* T = a*this.K*a' + this.Eps*sum(ai) - fxi*a'; */
            T = -sum(ai .* dW);
        }



        function [ai , dW , T ] = I2_WarmStart(fxi,initialai) {
            n = length(fxi);
            ap = zeros(1,n);
            am = ap;
            ap(initialai > 0) = initialai(initialai>0);
            am(initialai < 0) = initialai(initialai<0);
            ai = [ap am];
            hlp = -(ap-am)*this.K + fxi;
            dW = [hlp -hlp] - this.Eps;
            T = -sum(ai .* dW);
        }


        function [r , s , idx ] = getMax2DGainUpdatesForIndices(a,dW,ij) {
            n = numel(a)/2;
            kijidx = ij;
            kijidxup = kijidx > n;
            kijidx(kijidxup) = kijidx(kijidxup)-n;
            kijidx = sub2ind([n n],kijidx(:,1),kijidx(:,2));
            Kij = this.K(kijidx);

            sgn = 1-mod(sum(kijidxup,2),2)*2;
            dwi = dW(ij);
            aij = a(ij);

            div = 1 ./ (1-Kij.^2);

            /*  Compute update r */
            rs = [(dwi(:,1) - sgn .* Kij .* dwi(:,2)) .* div ...
                  (dwi(:,2) - sgn .* Kij .* dwi(:,1)) .* div];

            /*  Handle constraints */
            b = rs > this.C - aij | rs < - aij;

            ronly = find(b(:,1) & ~b(:,2));
            sonly = find(~b(:,1) & b(:,2));
            both = find(b(:,1) & b(:,2));

            /*  r is constrained only */
            if ~isempty(ronly)
                rs(ronly,1) = min(this.C - aij(ronly,1),max(- aij(ronly,1), rs(ronly,1)));
                rs(ronly,2) = min(this.C-aij(ronly,2), max(-aij(ronly,2), dwi(ronly,2) - rs(ronly,1) .* Kij(ronly)));
            end

            /*  s is constrained only */
            if ~isempty(sonly)
                rs(sonly,2) = min(this.C - aij(sonly,2),max(-aij(sonly,2), rs(sonly,2)));
                rs(sonly,1) = min(this.C-aij(sonly,1), max(-aij(sonly,1), dwi(sonly,1) - rs(sonly,2) .* Kij(sonly)));
            end

            /*  Case of both vars constrained
             * Update s2 and r2 for the cases of the other one is constrained */
            if ~isempty(both)
                Kijb = Kij(both);
                r1 = min(this.C - aij(both,1),max(- aij(both,1), rs(both,1)));
                s1 = min(this.C-aij(both,2), max(-aij(both,2), dwi(both,1) - r1 .* Kijb));
                s2 = min(this.C - aij(both,2),max(- aij(both,2), rs(both,2)));
                r2 = min(this.C-aij(both,1), max(-aij(both,1), dwi(both,2) - s2 .* Kijb));

                /*  Determine update via get higher gain
                 * Gain if r is constrained and s accordingly updated */
                g1 = r1 .* (dwi(both,1) -.5*r1) + s1 .*(dwi(both,2) -.5 * s1) + sgn(both) .* r1 .* s1.* Kijb;
                /*  Gain if s is constrained and r accordingly updated */
                g2 = r2 .* (dwi(both,1) -.5*r2) + s2 .*(dwi(both,2) -.5 * s2) + sgn(both) .* r2 .* s2 .* Kijb;

                /* % See which gain is higher. */
                cmp = g1 > g2;
                rs(both(cmp),1) = r1(cmp);
                rs(both(cmp),2) = s1(cmp);
                rs(both(~cmp),1) = r2(~cmp);
                rs(both(~cmp),2) = s2(~cmp);
            end

            /* % Compute gain */
            g = sum(rs .* (dwi - .5*rs),2) + sgn.*rs(:,1).*rs(:,2).*Kij;

            /*  Extract r,s updates */
            [dummy, idx] = max(g);
            r = rs(idx,1);
            s = rs(idx,2);
        }

    
    public: /* ( Static ) */

        static function res = test_ScalarEpsSVR_SMO() {
            x = -5:.1:5;
            fx = sinc(x)+.2*x;
            fx = fx ./ max(abs(fx));

            svr = general.regression.ScalarEpsSVR_SMO;
            svr.Eps= .1;
            svr.Lambda= 1/20;
            svr.Vis= 0;

            kernel = kernels.GaussKernel(.8);
            svr.K= kernel.evaluate(x,x);

            res = compute(1);
            res = res & compute(2);

            function res = compute(version)
                svr.Version= version;
                [ai, svidx] = svr.computeKernelCoefficients(fx,[]);
                sv = x(:,svidx);
                svfun = @(x)ai" *(kernel.evaluate(x,sv) ");
                fsvr = svfun(x);
                fdiff = abs(fsvr(svidx)-fx(svidx));
                res = isempty(find(fdiff > 1.01*svr.Eps,1));
            end
        }
};
}
}