Code covered by the BSD License  

Highlights from
Sugeno-type FIS output tuning

Sugeno-type FIS output tuning

by

 

12 Aug 2010 (Updated )

Tunes linear parameters of Sugeno-type FIS output using various methods of solving linear equations

gradient(A, B, Xinit, options)
function X = gradient(A, B, Xinit, options)
%GRADIENT Uses gradient descent method to estimate output parameters
%   
%   It is to be called by SUGENOTUNE as method
%   for solving to linear problem A * X = B
%   with given initial Xinit.
% 
%   It also accepts additional argument variable 'options':
%       options(1) - maximum number of epochs (default: 10);
%       options(2) - learning paradigm: 0 for batch learning, 
%                    1 for pattern learning (default: 0);      
%       options(3) - initial step size (default: 0.01);
%       options(4) - step increasing rate (default: 1.1);
%       options(5) - step decreasing rate (default: 0.9);
%       options(6) - termination tolerance on step (default: 1E-5);
%       options(7) - info display during iteration (default: 1).
% 
%   See also SUGENOTUNE
% 
%   Example:
%       A = rand(5, 4);
%       B = rand(5, 2);
%       Xinit = zeros(4, 2);
%       X = gradient(A, B, Xinit);

%   Reference: Jang, J.-S. R., ANFIS: Adaptive-Network-based Fuzzy
%   Inference Systems, IEEE Transactions on Systems, Man, and Cybernetics,
%   Vol. 23, No. 3, pp. 665-685, 1993. 

%   Per Konstantin A. Sidelnikov, 2010.

L = size(B, 2);
M = size(A, 1);

% Set options
default_opts = [10, 0, 0.01, 1.1, 0.9, 1E-5, 1];
if nargin < 4
    options = default_opts;
else
    if length(options) < length(default_opts)
        tmp = default_opts;
        tmp(1 : length(options)) = options;
        options = tmp;
    end
    k = find(isnan(options));
    options(k) = default_opts(k);
end
zmax = options(1);    % maximum epoch
method = options(2);  % batch (0) or pattern (1) learning
alpha = options(3);   % step size
incrate = options(4); % step increasing rate
decrate = options(5); % step decreasing rate
steptol = options(6); % step tolerance
display = options(7); % display option

X = Xinit;
for l = 1 : L
    % Minimum error
    Emin = Inf;    
    % Best unknowns w.r.t Emin
    x_best = Xinit(:, l);
    % Current unknowns
    x = Xinit(:, l);    
    % Number of consequtive reductions and oscillations
    reduc = 0; oscil = 0;
    % Previous error values
    Eprev1 = 0; Eprev2 = 0;
        
    if method == 0
        b_calc = A * x; % needed to start batch learning
    end    
    
    % Main loop
    b = B(:, l);
    for z = 1 : zmax
        x_prev = x;
        
        % Estimate unknowns
        if method == 0 % off-line update
            e = b - b_calc;
            dE = -2 * A' * e;
            x = paramupdate(x, dE, alpha);
            b_calc = A * x;
        elseif method == 1 % on-line update
            for m = 1 : M
                a = A(m, :);
                bm = a * x;
                em = b(m) - bm;
                dem = -2 * a' * em;                
                x = paramupdate(x, dem, alpha);
            end
            b_calc = A * x;
        else
            error('Unknown learning paradigm.');
        end        
      
        % Calculate RMSE for both methods
        E = norm(b - b_calc) / sqrt(M);
                
        % Check if current unknowns are better 
        % than all previous unknowns then save them
        if E < Emin
            Emin = E;
            x_best = x;
        end
        
        % Check convergence based on
        % current and previous unknown values
        if norm(x - x_prev) / ...
                (1 + max(norm(x), norm(x_prev))) < steptol
            fprintf(['Convergence for output %d ', ...
                'at epoch %d is designated\n'], l, z);
            break;
        end
        
        if display
            fprintf(['%d. Epoch: %d of %d, Error: %g, ', ...
                'Step: %g\n'], l, z, zmax, E, alpha);
        end        

        % Update step size if necessary
        % on the basis of previous error dynamics
        sign1 = sign(Eprev1 - Eprev2);
        sign2 = sign(E - Eprev1);
        
        % If E < Eprev1
        if sign2 == -1
            if reduc == 3
                alpha = alpha * incrate;
                reduc = 0;
            else
                reduc = reduc + 1;
            end
        else
            reduc = 0;
        end
        
        % If (E < Eprev1) and (Eprev1 > Eprev2)
        % or (E > Eprev1) and (Eprev1 < Eprev2)
        if sign1 ~= sign2
            if oscil == 3
                alpha = alpha * decrate;
                oscil = 0;
            else
                oscil = oscil + 1;
            end
        else
            oscil = 0;
        end
        
        % Save previous error values
        Eprev2 = Eprev1;
        Eprev1 = E;
    end
    
    X(:, l) = x_best;
end

function x_new = paramupdate(x_old, dE, alpha)
% Update unknowns if length of dE is nonzero

len = norm(dE(:));
if len ~= 0
    etha = alpha / len;
    x_new = x_old - etha * dE;
else
    x_new = x_old;
end

Contact us