Steepest Ascent

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Ausarbeitung zu Steepest Ascent in Matlab 7.1 mit ausführlichen Kommentar
% Vorlesung: Optimierung
% Copyright 2007 Thomas Ludwig
% This file is distributed under the terms of the GNU General Public License, see http://www.gnu.org/copyleft/gpl.txt
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Steepest Ascent - Steilster Anstieg
function Steepest_Ascent(c);
 
    % Startwerte
    k= 0;           % Iterationen
    x0= [-1,1];     % Startpunkt
    xk= x0;     
    eps= 0.00001;   % Abbruchbedingung
    tol= 1e-6;      % Toleranzgrenze
    
    % Konstante 0 < tk < 0.5 in Armijo Goldstein Regel
    tk= 0.5;        % Schrittweite
    tk_halb= 0.5;   % Fuer die halbe Schrittweite (tk/2)
    
    % Beispiel laden
    fx= gradient_beispiele(c);
    
    % Ausgabe der Startwerte
    fprintf('\nFunktionsparameter \n');
    fprintf('-> Startpunkt: f(%g,%g)= %f \n', xk(1), xk(2), fnc(xk, fx));
    fprintf('-> Anfangsschrittweite: %g \n', tk);
    fprintf('-> Abbruchbedingung: eps= %g \n', eps);
    fprintf('\n-----------------------------------------------------------\n');
    
    % Solange die Abbruchbedingung eps nicht erfüllt ist 
    % und solange tk >= 0 ist
    while  norm(gradienten(xk, fx)) > tol 
    
            if (abbruch(eps, xk, fx) == 0);
                
                % Gradienten berechnen
                gradf= gradienten(xk, fx);
 
                % Hesse und/oder Eigenwerte verwenden
                fx= hesse(xk, fx);
                fx= eigenwerte(fx);
 
                % Richtungsvektor dk (in Richtung Optimum bestimmen)
                if (fx.Richtung == 'Maximiere'); dk= gradf; end  % Steepest Ascent
                if (fx.Richtung == 'Minimiere'); dk= -gradf; end % Steepest Descent
 
                % Ausgabe
                fprintf('Iteration %g: x1= %g  | x2= %g | tk= %g | f(x,y)= %g | g= (%g,%g) \n', k, xk(1), xk(2), tk, fnc(xk, fx), gradf(1), gradf(2));
 
                % Schrittweitenkontrolle nach der Armijo Goldstein Regel
                while ( ( fnc(xk+(tk*dk), fx) - fnc(xk, fx) ) > (tk*(gradf'*dk) ) );
                   tk = tk*tk_halb;
                end
                
                xk= xk + (tk*dk);   % Neuen Punkt setzen
                k= k+1;             % Naechste Iteration
            
            else
                break;
            end
            
    end
 
end

% Quadratische Funktion
function f= fnc(x, fx);
    f= fx.a*x(1).^2 + fx.b*x(1)*x(2) + fx.c*x(2).^2 + fx.d*x(1) + fx.e*x(2) + fx.f;         
end

% Gradienten
function gradf= gradienten(x, fx);
    gradf= [2*fx.a*x(1)+fx.b*x(2)+fx.d , 2*fx.c*x(2)+fx.b*x(1)+ fx.e];
end

% Um einen Sattelpunkt auszuschließen wird die Hesse-Matrix auf positive Definitheit überprüft.
function fx= hesse(x, fx);
 
    fx.Hesse= [2*fx.a, fx.b; fx.b, 2*fx.c];
    fx.det_Hesse= det(fx.Hesse);
    
%     % Positiv Semidefinit -  Relatives Minimum
%     if ( (fx.Hesse(1, 1) >= 0) & (fx.det_Hesse >= 0) ), fx.Richtung= 'Minimiere';
% 
%          % Streng Konvex - Streng Positiv Semidefinit
%         elseif ( (fx.Hesse(1, 1) >= 0) & (fx.det_Hesse > 0) ), fx.Richtung= 'Minimiere'; 
% 
%         % Konkav - Negativ Semidefinit
%         elseif ( (fx.Hesse(1, 1) <= 0) & (fx.det_Hesse >= 0) ), fx.Richtung= 'Maximiere';  
% 
%         % Streng Negativ Semidefinit - rechtsgekrümmt                              
%         % Sattelpunkt oder relatives Maximum
%         elseif ( (fx.Hesse(1, 1) < 0) & (fx.det_Hesse > 0) ), fx.Richtung= 'Maximiere';  
% 
%         % Sattelpunkt - Indefinit
%         elseif ( (fx.Hesse(1, 1) > 0) & (fx.det_Hesse < 0) ), fx.Richtung= 'Unbestimmt';  
%             
%     end
 
end

function fx= eigenwerte(fx)
 
    if min(eig(fx.Hesse)) > 0 
        fx.Richtung= 'Minimiere';   % steepest ascent
    else
        fx.Richtung= 'Maximiere';  % steepest descent
    end    
    
end

function s= abbruch(eps, xk, fx);
 
    gradf= gradienten(xk, fx);
    gradf_betrag= sqrt(gradf(1)^2 + gradf(2)^2);
 
    if (gradf_betrag < eps);
       s= 1;    % Abbrechen
    elseif (gradf_betrag >= eps);
       s= 0;    % Weiter
    end
 
end

%
% Koeffizienten fuer verschiedene Berechnungen
% Parameter c waehlt die gewuenschte Aufgabe aus
%
function fx= gradient_beispiele(c);
 
    switch c
        
        case 1
            % x² und y² sind hier positiv
            % Somit ist die Funktion nach oben geoeffnet und besitzt
            % ein Minimum. Deshalb suchen wir ein lokales Minimum     
            % Dies ist dann eine Quelle (das lokale Minimum)
            fx.a= 1.5;
            fx.b= 2;
            fx.c= 1;
            fx.d= -1;
            fx.e= -2;  
            fx.f= 0;
            
        case 2
            % x² und y² sind hier negativ
            % Somit ist die Funktion nach unten geoeffnet und besitzt
            % ein Maximum. Deshalb suchen wir ein lokales Maximum
            fx.a= -2.5;
            fx.b= -4;
            fx.c= -3;
            fx.d= -4;
            fx.e= -2.5;  
            fx.f= 0;
            
        case 3
            % x² und y² sind hier positiv
            % Somit ist die Funktion nach oben geoeffnet und besitzt
            % ein Minimum. Deshalb suchen wir ein lokales Minimum     
            % Dies ist dann eine Quelle (das lokale Minimum)
            fx.a= 1.5;
            fx.b= 2;
            fx.c= 1;
            fx.d= 0.5;
            fx.e= 1;  
            fx.f= 0;
            
        case 4
            % x² und y² sind hier negativ
            % Somit ist die Funktion nach unten geoeffnet und besitzt
            % ein Maximum. Deshalb suchen wir ein lokales Maximum
            fx.a= -2.5;
            fx.b= 4;
            fx.c= -3;
            fx.d= 2;
            fx.e= 1.25;  
            fx.f= 0;   
 
    end
end